#!/usr/bin/env python # This was developped by Mathieu Compiegne at HYGEOS from __future__ import print_function, division, absolute_import from builtins import range import os import sys import numpy as np import pylab as plt import time ######################################## dir_artdeco = "/home/mathieu/work/RTM/artdeco/fortran/" dir_pyartdeco = "/home/mathieu/work/RTM/artdeco/pyartdeco/" sys.path.append(dir_pyartdeco+'f2py_utils/') import pyartdeco_runlib as pyartdeco from f2py_utils import run_artdeco ######################################## def find_nearest(array, value): array = np.asarray(array) idx = (np.abs(array - value)).argmin() return idx def get_atm_iprt_b1b4(file_path, res_path): res = np.genfromtxt(file_path, names=['alt','tau']) cumul_tau = np.cumsum(res["tau"]) # print(cumul_tau) # print(cumul_tau.shape) # print(res["tau"].shape) p = cumul_tau / np.max(cumul_tau) * 1013.0 t = np.full_like(p,0.0) nair = np.full_like(p,0.0) print(np.max(cumul_tau)) f = open(res_path,"w") f.write("# number of sampled altitude : \n") f.write(" %i \n"%len(res['alt'])) f.write("# number of sampled gas : \n") f.write(" 0 \n") f.write("# sampled gas :\n") f.write(" none \n") f.write("# z(km) p(mb) T(K) air(cm-3) \n") for ialt in range(len(res["tau"])): s = " %12.6f"%res["alt"][ialt]+" %18.8f"%p[ialt]+" %12.6f"%t[ialt]+" %12.6f"%nair[ialt] f.write(s+ "\n") f.close() return def get_atm_iprt_b2(file_path, file_path_abs, res_path): res = np.genfromtxt(file_path, names=['alt','dtau']) cumul_tau = np.cumsum(res["dtau"]) p = cumul_tau / np.max(cumul_tau) * 1013.0 t = np.full_like(p,0.0) nair = np.full_like(p,0.0) print(np.max(cumul_tau)) resabs = np.genfromtxt(file_path_abs, names=['alt','dtau']) f = open(res_path,"w") f.write("# number of sampled altitude : \n") f.write(" %i \n"%len(res['alt'])) f.write("# number of wavelengths : \n") f.write(" 1 \n") f.write("# wavelengths : \n") f.write(" 0.325 \n") f.write("# z(km) p(mb) T(K) dtauabs \n") for ialt in range(len(res["dtau"])): s = " %12.6f"%res["alt"][ialt]+" %18.8f"%p[ialt]+" %12.6f"%t[ialt]+" %12.6e"%resabs["dtau"][ialt] f.write(s+ "\n") f.close() return def get_atm_iprt_b3(file_path, file_path_abs, res_path): res = np.genfromtxt(file_path, names=['alt','dtau']) cumul_tau = np.cumsum(res["dtau"]) p = cumul_tau / np.max(cumul_tau) * 1013.0 t = np.full_like(p,0.0) nair = np.full_like(p,0.0) print(np.max(cumul_tau)) resabs = np.genfromtxt(file_path_abs, names=['alt','dtau']) f = open(res_path,"w") f.write("# number of sampled altitude : \n") f.write(" %i \n"%len(res['alt'])) f.write("# number of wavelengths : \n") f.write(" 1 \n") f.write("# wavelengths : \n") f.write(" 0.350 \n") f.write("# z(km) p(mb) T(K) dtauabs \n") for ialt in range(len(res["dtau"])): s = " %12.6f"%res["alt"][ialt]+" %18.8f"%p[ialt]+" %12.6f"%t[ialt]+" %12.6e"%resabs["dtau"][ialt] f.write(s+ "\n") f.close() return def get_vprof_aer_iprt_b3(file_path): res = np.genfromtxt(file_path, names=['alt','dtau']) cumul_tau = np.sum(res["dtau"]) # alt = np.zeros( len(res["alt"])-1 ) # vprof = np.zeros( len(res["alt"])-1 ) # for ialt in range( len(alt) ): # alt[ialt] = (res["alt"][ialt] + res["alt"][ialt+1]) / 2.0 # vprof[ialt] = res["dtau"][ialt+1] # alt = np.append(alt, np.max(alt)+0.001 ) # vprof = np.append(vprof, 0.0 ) # alt = np.append(alt, 100.0 ) # vprof = np.append(vprof, 0.0 ) # alt = np.append(alt, np.min(alt)-0.001 ) # vprof = np.append(vprof, 0.0 ) # alt = np.append(alt, -10.0 ) # vprof = np.append(vprof, 0.0 ) # ind = np.argsort(alt) dz = np.abs(np.diff(res["alt"])) tau = np.zeros_like(dz) for itau in range(len(dz)): tau[itau] = res["dtau"][itau+1] / dz[itau] alt = np.linspace(-0.01, 31.0, num=1000) vprof = np.zeros_like( alt ) for ialt in range( len(alt) ): if (alt[ialt] > np.max(res["alt"])) or (alt[ialt] < np.min(res["alt"])): vprof[ialt] = 0.0 else: ii = find_nearest(res["alt"], alt[ialt]) if alt[ialt] >= res["alt"][ii]: vprof[ialt] = tau[ii-1] else: vprof[ialt] = tau[ii] ind = np.argsort(alt) # import matplotlib.pyplot as plt # plt.figure() # plt.plot(vprof[ind], alt[ind]) # plt.show() return cumul_tau, alt[ind], vprof[ind] def plot_polar_contour(values, azimuths, zeniths, barname='Pixel reflectance', cut=False): """Plot a polar contour plot, with 0 degrees at the North. Arguments: * `values` -- A list (or other iterable - eg. a NumPy array) of the values to plot on the contour plot (the `z` values) * `azimuths` -- A list of azimuths (in degrees) * `zeniths` -- A list of zeniths (that is, radii) The shapes of these lists are important, and are designed for a particular use case (but should be more generally useful). The values list should be `len(azimuths) * len(zeniths)` long with data for the first azimuth for all the zeniths, then the second azimuth for all the zeniths etc. This is designed to work nicely with data that is produced using a loop as follows: values = [] for azimuth in azimuths: for zenith in zeniths: # Do something and get a result values.append(result) After that code the azimuths, zeniths and values lists will be ready to be passed into this function. """ theta = np.radians(azimuths) zeniths = np.array(zeniths) values = np.array(values) val = values.reshape(len(azimuths), len(zeniths)) cut_min = 0.1 cut_max = 99.9 ncontour = 50 if cut: minval = np.percentile(val[~np.isnan(val)], cut_min) maxval = np.percentile(val[~np.isnan(val)], cut_max) print(minval) print(maxval) extend = "both" # if abs(minval - np.min(val[~np.isnan(val)])) < diff_lim : # extend = "max" # minval = np.min(val[~np.isnan(val)]) # if abs(maxval - np.max(val[~np.isnan(val)])) < diff_lim : # maxval = np.max(val[~np.isnan(val)]) # if extend == "max": # extend = "neither" # else: # extend = "min" levels = np.linspace(minval, maxval, num=ncontour, endpoint = True) r, theta = np.meshgrid(zeniths, np.radians(azimuths)) fig, ax = plt.subplots(subplot_kw=dict(projection='polar')) ax.set_theta_zero_location("E") ax.set_theta_direction(-1) #autumn() if cut: cax = ax.contourf(theta, r, val, levels, extend=extend) else: cax = ax.contourf(theta, r, val, ncontour) #autumn() cb = fig.colorbar(cax) cb.set_label(barname) return fig, ax, cax def get_kokha_ref(dir_kokha_data, case, nmat, vza, vaa): nvaa = np.size(vaa) nvza = np.size(vza) if nvaa == 1 : vaa = np.array([vaa]) if nvza == 1 : vza = np.array([vza]) for i in range (nvaa): if ((vaa[i] != 0.0) & (vaa[i] != 90.0) & (vaa[i] != 180.0)) : print('View azimuthal angle must be 0, 90 or 180 degree') exit() nvza_kokha = 90 vza_kokha = np.zeros(nvza_kokha) SvR_kokha = np.zeros((nvza_kokha, 3, 4)) if case == 'aer': f = open( dir_kokha_data+'aerosol_refl_N_240.dat', 'r') if case == 'cl': f = open( dir_kokha_data+'cloud_refl_N_360.dat', 'r') if case == 'ray': f = open( dir_kokha_data+'Rayleigh_refl_N_60.dat', 'r') for i in range(nvza_kokha): tmp = f.readline() vza_kokha[i] = float(tmp.split()[0]) SvR_kokha[i,0,0] = float(tmp.split()[1]) # I SvR_kokha[i,0,1] = -float(tmp.split()[2]) # Q SvR_kokha[i,0,2] = float(tmp.split()[3]) # U SvR_kokha[i,1,0] = float(tmp.split()[5]) # I SvR_kokha[i,1,1] = -float(tmp.split()[6]) # Q SvR_kokha[i,1,2] = float(tmp.split()[7]) # U SvR_kokha[i,1,3] = -float(tmp.split()[8]) # V SvR_kokha[i,2,0] = float(tmp.split()[9]) # I SvR_kokha[i,2,1] = -float(tmp.split()[10]) # Q SvR_kokha[i,2,2] = float(tmp.split()[11]) # U SvR_kokha[i,2,3] = -float(tmp.split()[12]) # V f.close() if nmat==1: SvR = np.zeros((nvza, nvaa)) dSvR = np.zeros((nvza, nvaa)) for k in range(nvaa): if (vaa[k] == 0.0): SvR[:,k] = np.interp(vza, vza_kokha, SvR_kokha[:, 0, 0]) if (vaa[k] == 90.0): SvR[:,k] = np.interp(vza, vza_kokha, SvR_kokha[:, 1, 0]) if (vaa[k] == 180.0): SvR[:,k] = np.interp(vza, vza_kokha, SvR_kokha[:, 2, 0]) # if (nvaa == 1) & (nvza == 1) : # SvR = SvR[0,0] # else: # if nvaa == 1 : # SvR = SvR[:,0] # if nvza == 1 : # SvR = SvR[0,:] else: SvR = np.zeros((nvza, nvaa, nmat)) dSvR = np.zeros((nvza, nvaa, nmat)) for k in range(nvaa): for u in range(nmat): if (vaa[k] == 0.0): SvR[:,k,u] = np.interp(vza, vza_kokha, SvR_kokha[:, 0, u]) if (vaa[k] == 90.0): SvR[:,k,u] = np.interp(vza, vza_kokha, SvR_kokha[:, 1, u]) if (vaa[k] == 180.0): SvR[:,k,u] = np.interp(vza, vza_kokha, SvR_kokha[:, 2, u]) # if (nvaa == 1) & (nvza == 1) : # SvR = SvR[0,0,:] # else: # if nvaa == 1 : # SvR = SvR[:,0,:] # if nvza == 1 : # SvR = SvR[0,:,:] return SvR #################################################### def kokha_rayleigh(nstreams_init): trans = False # Kokhanovsky Rayleigh ########################################### # artdeco_in technical parameters structure keywords = ['nowarn', 'od_no_check'] mode = 'mono' filters = ['none'] wavel = [0.40961879958331204] rt_model = 'doad' thermal = False corint = False nmat = 3 trunc_method = 'none' artdeco_in = pyartdeco.artdeco_in(keywords, mode, filters, wavel, \ trunc_method, nstreams_init, rt_model, \ corint, thermal, nmat) #################### # read kdis coeff kdis = pyartdeco.kdis_coeff(artdeco_in, dir_artdeco+'lib/kdis/', 'ascii') ###################### # load solar TOA flux solrad = pyartdeco.solar_irradiance(artdeco_in, dir_artdeco+'lib/solrad/', 'ascii') #################### # load atmosphere gas = ['none'] ppmv = [-1] atm = 'atm_kokha.dat' dir_data = dir_artdeco+'lib/atm/' wldepol = np.array([0.1, 200.0]) # microns depol = np.array([0.0, 0.0]) atmos = pyartdeco.atmosphere(artdeco_in, dir_data, atm, gas, ppmv, 'ascii', wldepol, depol, interp=True, kdis=kdis) ############################### # set ptcle structure ptcle_def = [] opt_interp = False wlref = None wlptcle = None dir_data = None ptcle = pyartdeco.particle_old(artdeco_in, dir_data, ptcle_def, wlref, wlptcle, opt_interp) #print ptcle.__dict__.keys() #print ptcle.__dict__ ######################## # set surface structure name = "black" family = "lambert" kind = "land" # wavelength definition for the surface wl = np.array([0.1, 200.0]) # microns alb = np.array([0.0, 0.0]) surface = pyartdeco.surface(name, family, kind, wl=wl, alb=alb, interp=True) ######################## # Geometry sza = np.array([60.0]) nvza = 90 if trans : min_vza = 91. max_vza = 180. else: min_vza = 0 max_vza = 89. dvza = (max_vza - min_vza) / (nvza-1) vza = np.zeros(nvza) for j in range(nvza): vza[j] = max_vza - (j * dvza) print(vza[j]) vaa = np.array([0.0, 90.0, 180.0]) geom = pyartdeco.geometry(sza, vza, vaa) ############### # run ARTDECO lamb, rad, rad_levels, flux, alt = run_artdeco(artdeco_in, atmos, surface, solrad, ptcle, geom, kdis=kdis, verbose=True) ##(artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda) ##(artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nlambda) if trans: rad = rad_levels[:,:,:,:,-1,:] * np.pi / np.cos(geom.sza[0] * np.pi / 180.) else: # noramalisation as Kokhanovsky rad = rad * np.pi / np.cos(geom.sza[0] * np.pi / 180.) ##################### # Plot the result # get Kokha ref SvR_kokha = get_kokha_ref('/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/kokha/' ,'ray', artdeco_in.nmat, geom.vza, geom.vaa) color = ['r','g','b','c','m','y','k','w'] line = ['-','--'] sign = ['d', 'd'] fsize=15 plt.rcParams.update({'font.size': fsize}) lw = 2 if trans: xlim=[90,180] else: xlim=[0,90] if artdeco_in.nmat == 1 : plt.figure(figsize=(6,6)) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$I_R$') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, rad[0,:,k,0,0],color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$ ARTDECO') if not trans: plt.plot(vza, SvR_kokha[:,k],color[k]+line[1],label='$\phi_v ='+str(vaa[k])+'^o$ Kokha') plt.legend(loc="best", prop={'size':fsize}) # plt.savefig('/home/mathieu/tmp/kokha_test_notms.png') elif artdeco_in.nmat == 3 : plt.figure(figsize=(6,12)) plt.subplots_adjust(bottom=0.05, top=0.97, left=0.17, right=0.95) plt.subplot(3,1,1) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$I_R$') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, rad[0,:,k,0,0],color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$ ARTDECO', lw=lw) if not trans: plt.plot(vza, SvR_kokha[:,k,0],color[k]+line[1],label='$\phi_v ='+str(vaa[k])+'^o$ Kokha', lw=lw) plt.legend(loc="best", prop={'size':12}) plt.subplot(3,1,2) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$-Q_R$') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, -rad[0,:,k,1,0],color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$',lw=lw) if not trans: plt.plot(vza, -SvR_kokha[:,k,1],color[k]+line[1],label='$\phi_v ='+str(vaa[k])+'^o$ Kokha',lw=lw) plt.subplot(3,1,3) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$ U_R$') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, -rad[0,:,k,2,0],color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$',lw=lw) if not trans: plt.plot(vza, SvR_kokha[:,k,2],color[k]+line[1],label='$\phi_v ='+str(vaa[k])+'^o$ Kokha',lw=lw) #plt.tight_layout #plt.savefig('/home/mathieu/tmp/kokha_test.png') if not trans: # diff if artdeco_in.nmat == 1 : plt.figure(figsize=(6,6)) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$\delta I_R / $ (%)') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, (rad[0,:,k,0,0] - SvR_kokha[:,k]) / SvR_kokha[:,k] * 100., color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$') plt.legend(loc="best", prop={'size':fsize}) # plt.savefig('/home/mathieu/tmp/kokha_test_delta_notms.png') elif artdeco_in.nmat == 3 : plt.figure(figsize=(6,12)) plt.subplots_adjust(bottom=0.05, top=0.97, left=0.17, right=0.95) plt.subplot(3,1,1) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$\delta I_R / $ (%)') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, (rad[0,:,k,0,0] - SvR_kokha[:,k,0]) / SvR_kokha[:,k, 0] * 100., color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$',lw=lw) plt.legend(loc="best", prop={'size':10}) plt.subplot(3,1,2) plt.grid(True) plt.ylim([-15.0,15.0]) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$\delta Q_R / $ (%)') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, (rad[0,:,k,1,0] - SvR_kokha[:,k,1]) / SvR_kokha[:,k, 1] * 100., color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$',lw=lw) plt.subplot(3,1,3) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$\delta U_R / $ (%)') plt.xlim(xlim) k = 1 plt.plot(vza, (rad[0,:,k,2,0] - -SvR_kokha[:,k,2]) / SvR_kokha[:,k, 2] * 100., color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$',lw=lw) #plt.savefig('/home/mathieu/tmp/kokha_test_delta.png') plt.show() # if artdeco_in.nmat == 1 : # plt.grid(True) # plt.xlabel('$\Theta_v$ (deg.)') # plt.ylabel('$I_R$') # plt.xlim([0,90]) # for k in range(len(geom.vaa)): # plt.plot(geom.vza, rad[0,:,k,0,0],color[k]+line[0],label='$\phi_v ='+str(geom.vaa[k])+'^o$ ARTDECO') # plt.plot(geom.vza, SvR_kokha[:,k],color[k]+line[1],label='$\phi_v ='+str(geom.vaa[k])+'^o$ Kokha') # plt.legend(loc="best", prop={'size':10}) # plt.show() # elif artdeco_in.nmat == 3 : # plt.figure(figsize=(6,12)) # plt.subplots_adjust(bottom=0.05, top=0.97, left=0.17, right=0.95) # plt.subplot(3,1,1) # plt.grid(True) # plt.xlabel('$\Theta_v$ (deg.)') # plt.ylabel('$I_R$') # plt.xlim([0,90]) # for k in range(len(geom.vaa)): # plt.plot(geom.vza, rad[0,:,k,0,0],color[k]+line[0],label='$\phi_v ='+str(geom.vaa[k])+'^o$ ARTDECO') # plt.plot(geom.vza, SvR_kokha[:,k,0],color[k]+line[1],label='$\phi_v ='+str(geom.vaa[k])+'^o$ Kokha') # plt.legend(loc="best", prop={'size':10}) # plt.subplot(3,1,2) # plt.grid(True) # plt.xlabel('$\Theta_v$ (deg.)') # plt.ylabel('$-Q_R$') # plt.xlim([0,90]) # for k in range(len(geom.vaa)): # plt.plot(geom.vza, -rad[0,:,k,1,0],color[k]+line[0],label='$\phi_v ='+str(geom.vaa[k])+'^o$') # plt.plot(geom.vza, -SvR_kokha[:,k,1],color[k]+line[1],label='$\phi_v ='+str(geom.vaa[k])+'^o$ Kokha') # plt.subplot(3,1,3) # plt.grid(True) # plt.xlabel('$\Theta_v$ (deg.)') # plt.ylabel('$ U_R$') # plt.xlim([0,90]) # for k in range(len(geom.vaa)): # plt.plot(geom.vza, -rad[0,:,k,2,0],color[k]+line[0],label='$\phi_v ='+str(geom.vaa[k])+'^o$') # plt.plot(geom.vza, SvR_kokha[:,k,2],color[k]+line[1],label='$\phi_v ='+str(geom.vaa[k])+'^o$ Kokha') # # diff # if artdeco_in.nmat == 1 : # plt.figure() # plt.grid(True) # plt.xlabel('$\Theta_v$ (deg.)') # plt.ylabel('$\delta I_R / $ (%)') # plt.xlim([0,90]) # for k in range(len(geom.vaa)): # plt.plot(geom.vza, (rad[0,:,k,0,0] - SvR_kokha[:,k]) / SvR_kokha[:,k] * 100., color[k]+line[0],label='$\phi_v ='+str(geom.vaa[k])+'^o$') # plt.legend(loc="best", prop={'size':10}) # plt.show() # elif artdeco_in.nmat == 3 : # plt.figure(figsize=(6,12)) # plt.subplots_adjust(bottom=0.05, top=0.97, left=0.17, right=0.95) # plt.subplot(3,1,1) # plt.grid(True) # plt.xlabel('$\Theta_v$ (deg.)') # plt.ylabel('$\delta I_R / $ (%)') # plt.xlim([0,90]) # for k in range(len(geom.vaa)): # plt.plot(geom.vza, (rad[0,:,k,0,0] - SvR_kokha[:,k,0]) / SvR_kokha[:,k, 0] * 100., color[k]+line[0],label='$\phi_v ='+str(geom.vaa[k])+'^o$') # plt.legend(loc="best", prop={'size':10}) # plt.subplot(3,1,2) # plt.grid(True) # plt.xlabel('$\Theta_v$ (deg.)') # plt.ylabel('$\delta Q_R / $ (%)') # plt.xlim([0,90]) # for k in range(len(geom.vaa)): # plt.plot(geom.vza, (rad[0,:,k,1,0] - SvR_kokha[:,k,1]) / SvR_kokha[:,k, 1] * 100., color[k]+line[0],label='$\phi_v ='+str(geom.vaa[k])+'^o$') # plt.subplot(3,1,3) # plt.grid(True) # plt.xlabel('$\Theta_v$ (deg.)') # plt.ylabel('$\delta U_R / $ (%)') # plt.xlim([0,90]) # k = 1 # plt.plot(geom.vza, (rad[0,:,k,2,0] - -SvR_kokha[:,k,2]) / SvR_kokha[:,k, 2] * 100., color[k]+line[0],label='$\phi_v ='+str(geom.vaa[k])+'^o$') # plt.show() def kokha_cl_aer(nstreams_init, case): trans = False ########################################### # artdeco_in technical parameters structure keywords = ['od_no_check', "no_rayleigh"] mode = 'mono' filters = ['none'] wavel = [0.412] rt_model = 'doad' thermal = False corint = True nmat = 3 trunc_method = 'dm' artdeco_in = pyartdeco.artdeco_in(keywords, mode, filters, wavel, \ trunc_method, nstreams_init, rt_model, \ corint, thermal, nmat) #################### # read kdis coeff kdis = pyartdeco.kdis_coeff(artdeco_in, dir_artdeco+'lib/kdis/', 'ascii') ###################### # load solar TOA flux solrad = pyartdeco.solar_irradiance(artdeco_in, dir_artdeco+'lib/solrad/', 'ascii') #################### # load atmosphere gas = ['none'] ppmv = [-1] atm = 'atm_noatm.dat' dir_data = "/rfs/proj/artdeco_lib/atm/" wldepol = np.array([0.1, 200.0]) # microns depol = np.array([0.0, 0.0]) atmos = pyartdeco.atmosphere(artdeco_in, dir_data, atm, gas, ppmv, 'ascii', wldepol, depol, interp=True, kdis=kdis) ############################### # set ptcle structure if case == 'cl': ptcle_def = [ {"file_path":"/rfs/proj/artdeco_lib/opt/opt_kokha.h5", "name":"kokha_cl", "tau":5.0, "alt_distrib":"layer", "z":1.0 } ] elif case == 'aer': ptcle_def = [ {"file_path":"/rfs/proj/artdeco_lib/opt/opt_kokha.h5", "name":"kokha_aer", "tau":0.3262, "alt_distrib":"layer", "z":1.0} ] # if case == 'cl': # ptcle_def = [ {"file":"opt_kokha.h5", "name":"kokha_cl", "tau":5.0, "alt_distrib":"layer", "z":1.0 } ] # elif case == 'aer': # ptcle_def = [ {"file":"opt_kokha.h5", "name":"kokha_aer", "tau":0.3262, "alt_distrib":"layer", "z":1.0} ] opt_interp = True wlref = 0.412 # wavelengths to read in files wlptcle = [0.1, 100.0] dir_data = dir_artdeco+'lib/opt/' # ptcle = pyartdeco.particle_old(artdeco_in, dir_data, ptcle_def, wlref, wlptcle, opt_interp) ptcle_opt = pyartdeco.ptcle_optical_properties(ptcle_def, artdeco_in.nstreams, wlptcle, wlref, read_betal=False, opt_interp=True) ptcle = pyartdeco.particle(artdeco_in, ptcle_def, ptcle_opt) #print ptcle.__dict__.keys() #print ptcle.__dict__ ######################## # set surface structure name = "black" family = "lambert" kind = "land" # wavelength definition for the surface wl = np.array([0.1, 200.0]) # microns alb = np.array([0.0, 0.0]) surface = pyartdeco.surface(name, family, kind, wl=wl, alb=alb, interp=True) ######################## # Geometry sza = np.array([60.0]) nvza = 90 if trans : min_vza = 91. max_vza = 180. else: min_vza = 0 max_vza = 89. dvza = (max_vza - min_vza) / (nvza-1) vza = np.zeros(nvza) for j in range(nvza): vza[j] = max_vza - (j * dvza) #print(vza[j]) vaa = np.array([0.0, 90.0, 180.0]) ########################################## # a single call for all geom geom = pyartdeco.geometry(sza, vza, vaa) lamb, rad, rad_levels, flux, alt = run_artdeco(artdeco_in, atmos, surface, solrad, ptcle, geom, kdis=kdis, verbose=True) ##(artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda) ##(artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nvaa, artdeco.mcommon.nmat, artdeco.mcommon.nlambda) if trans: rad = rad_levels[:,:,:,:,-1,:] * np.pi / np.cos(geom.sza[0] * np.pi / 180.) else: rad = rad * np.pi / np.cos(geom.sza[0] * np.pi / 180.) # ################################ # # separate call for vza ans sza # flag_final_array = False # verbose = True # for isza in range(len(sza)): # for ivza in range(len(vza)): # geom = pyartdeco.geometry(np.array([sza[isza]]), np.array([vza[ivza]]), vaa) # lamb_tmp, rad_tmp, flux_tmp, alt_tmp = run_artdeco(artdeco_in, atmos, surface, solrad, ptcle, geom, kdis=kdis, verbose=verbose) # verbose = False # if not flag_final_array: # nlamb = rad_tmp.shape[4] # rad = np.zeros((len(sza), len(vza), len(vaa), artdeco_in.nmat, nlamb)) # flag_final_array = True # rad[isza,ivza,:,:,:] = rad_tmp[0,0,:,:,:] * np.pi / np.cos(sza[isza] * np.pi / 180.) ##################### # Plot the result # get Kokha ref SvR_kokha = get_kokha_ref('benchmark/kokha/', case, artdeco_in.nmat, vza, vaa) color = ['r','g','b','c','m','y','k','w'] line = ['-','--'] sign = ['d', 'd'] fsize=15 plt.rcParams.update({'font.size': fsize}) lw = 2 if trans: xlim=[90,180] else: xlim=[0,90] if artdeco_in.nmat == 1 : plt.figure(figsize=(6,6)) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$I_R$') plt.xlim(xlim) for k in range(len(vaa)): if not trans: plt.plot(vza, rad[0,:,k,0,0],color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$ ARTDECO') plt.plot(vza, SvR_kokha[:,k],color[k]+line[1],label='$\phi_v ='+str(vaa[k])+'^o$ Kokha') else: plt.semilogy(vza, rad[0,:,k,0,0],color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$ ARTDECO') plt.legend(loc="best", prop={'size':fsize}) figname = 'kokha_'+case if corint: figname = figname+"_corint" figname=figname+".png" plt.savefig(figname) elif artdeco_in.nmat == 3 : plt.figure(figsize=(6,12)) plt.subplots_adjust(bottom=0.05, top=0.97, left=0.17, right=0.95) plt.subplot(3,1,1) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$I_R$') plt.xlim(xlim) for k in range(len(vaa)): if not trans: plt.plot(vza, rad[0,:,k,0,0],color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$ ARTDECO', lw=lw) plt.plot(vza, SvR_kokha[:,k,0],color[k]+line[1],label='$\phi_v ='+str(vaa[k])+'^o$ Kokha', lw=lw) else: plt.semilogy(vza, rad[0,:,k,0,0],color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$ ARTDECO', lw=lw) plt.legend(loc="best", prop={'size':12}) plt.subplot(3,1,2) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$-Q_R$') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, -rad[0,:,k,1,0],color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$',lw=lw) if not trans: plt.plot(vza, -SvR_kokha[:,k,1],color[k]+line[1],label='$\phi_v ='+str(vaa[k])+'^o$ Kokha',lw=lw) plt.subplot(3,1,3) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$ U_R$') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, -rad[0,:,k,2,0],color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$',lw=lw) if not trans: plt.plot(vza, SvR_kokha[:,k,2],color[k]+line[1],label='$\phi_v ='+str(vaa[k])+'^o$ Kokha',lw=lw) plt.tight_layout figname = 'kokha_'+case if corint: figname = figname+"_corint" figname=figname+".png" plt.savefig(figname) if not trans: # diff if artdeco_in.nmat == 1 : plt.figure(figsize=(6,6)) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$\delta I_R / $ (%)') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, (rad[0,:,k,0,0] - SvR_kokha[:,k]) / SvR_kokha[:,k] * 100., color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$') plt.legend(loc="best", prop={'size':fsize}) figname = 'kokha_delta_'+case if corint: figname = figname+"_corint" figname=figname+".png" plt.savefig(figname) elif artdeco_in.nmat == 3 : plt.figure(figsize=(6,12)) plt.subplots_adjust(bottom=0.05, top=0.97, left=0.17, right=0.95) plt.subplot(3,1,1) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$\delta I_R / $ (%)') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, (rad[0,:,k,0,0] - SvR_kokha[:,k,0]) / SvR_kokha[:,k, 0] * 100., color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$',lw=lw) plt.legend(loc="best", prop={'size':10}) plt.subplot(3,1,2) plt.grid(True) plt.ylim([-15.0,15.0]) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$\delta Q_R / $ (%)') plt.xlim(xlim) for k in range(len(vaa)): plt.plot(vza, (rad[0,:,k,1,0] - SvR_kokha[:,k,1]) / SvR_kokha[:,k, 1] * 100., color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$',lw=lw) plt.subplot(3,1,3) plt.grid(True) plt.xlabel('$\Theta_v$ (deg.)') plt.ylabel('$\delta U_R / $ (%)') plt.xlim(xlim) k = 1 plt.plot(vza, (rad[0,:,k,2,0] - -SvR_kokha[:,k,2]) / SvR_kokha[:,k, 2] * 100., color[k]+line[0],label='$\phi_v ='+str(vaa[k])+'^o$',lw=lw) figname = 'kokha_delta_'+case if corint: figname = figname+"_corint" figname=figname+".png" plt.savefig(figname) plt.show() def iprt_a(nstreams_init, trans, case, nmat, rt_model, corint, trunc_method, test_od_ims=False): if "a5"in case: sub_case = case.split('_')[1] case = case.split('_')[0] elif "a1" in case: sub_case = int(case.split('_')[1]) case = case.split('_')[0] ########################################### # artdeco_in technical parameters structure if case in ["a1","a2","a6"]: keywords = ['nowarn', 'od_no_check'] elif case in ["a3","a4","a5"]: keywords = ["no_rayleigh", 'nowarn', 'od_no_check'] mode = 'mono' filters = ['none'] if case=="a1": wavel = [0.3] elif case in ["a2","a6"]: wavel = [0.5] elif case in ["a3","a4"]: wavel = [0.35] elif case=="a5": wavel = [0.8] thermal = False potter_theta_min = 15.0 potter_theta_max = 20.0 if test_od_ims: artdeco_in = pyartdeco.artdeco_in(keywords, mode, filters, wavel, \ trunc_method, nstreams_init, rt_model, \ corint, thermal, nmat, \ od_nstr=8, od_deltam=True, od_corint=True, od_do_secsca=True) else: artdeco_in = pyartdeco.artdeco_in(keywords, mode, filters, wavel, \ trunc_method, nstreams_init, rt_model, \ corint, thermal, nmat, potter_theta_min=potter_theta_min, \ potter_theta_max=potter_theta_max) #################### # read kdis coeff kdis = pyartdeco.kdis_coeff(artdeco_in, dir_artdeco+'lib/kdis/', 'ascii') ###################### # load solar TOA flux solrad = pyartdeco.solar_irradiance(artdeco_in, dir_artdeco+'lib/solrad/', 'ascii') #################### # load atmosphere gas = ['none'] ppmv = [-1] if case in ["a1","a2","a3","a5","a4","a6"]: atm = 'atm_kokha.dat' dir_data = dir_artdeco+'lib/atm/' wldepol = np.array([0.1, 200.0]) # microns if case == "a1": taur = 0.5 if sub_case == 1: depol = 0.0 elif sub_case == 2: depol = 0.03 elif sub_case == 3: depol = 0.1 elif case in ["a2","a6"]: taur = 0.1 depol = 0.03 elif case in ["a3","a5","a4"]: depol = 0.0 taur = 0.0 atmos = pyartdeco.atmosphere(artdeco_in, dir_data, atm, gas, ppmv, 'ascii', wldepol, np.array([depol, depol]), interp=True, kdis=kdis, tau_ray = { "wvl":np.array([0.1,200.]), "tau":np.array([taur,taur]) }) ############################### # set ptcle structure if case in ["a1","a2","a6"]: ptcle_def = [] wlref = None wlptcle = None elif case in ["a3"]: ptcle_def = [ {"file_path":"/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/opt_waso_iprt.h5", "name":"waso_iprt", "tau":0.2, "alt_distrib":"layer", "z":1.0} ] wlref = 0.35 wlptcle = [0.1, 100.0] elif case in ["a4"]: ptcle_def = [ {"file_path":"/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/opt_sizedistr_spheroid_iprt.h5", "name":"sizedistr_spheroid_iprt", "tau":0.2, "alt_distrib":"layer", "z":1.0} ] wlref = 0.35 wlptcle = [0.1, 100.0] elif case in ["a5"]: ptcle_def = [ {"file_path":"/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/opt_watercloud_iprt.h5", "name":"watercloud_iprt", "tau":5.0, "alt_distrib":"layer", "z":1.0} ] wlref = 0.8 wlptcle = [0.1, 100.0] ptcle_opt = pyartdeco.ptcle_optical_properties(ptcle_def, artdeco_in.nstreams, wlptcle, wlref, opt_interp=False, read_betal=False) ptcle = pyartdeco.particle(artdeco_in, ptcle_def, ptcle_opt) #print ptcle.__dict__.keys() #print ptcle.__dict__ ######################## # set surface structure if case in ["a3","a1","a5","a4","a2"]: name = "toto" family = "lambert" kind = "land" # wavelength definition for the surface wl = np.array([0.1, 200.0]) # microns if case in ["a3","a1","a5","a4"]: alb = np.array([0.0, 0.0]) elif case == "a2": alb = np.array([0.3, 0.3]) surface = pyartdeco.surface(name, family, kind, wl=wl, alb=alb, interp=True) elif case=="a6": surface_name = "surface" surface_interp = False surface_family = "brdf" surface_kind = "ocean" surface_wl = np.array([0.1, 100.0]) # microns surface_nwvl = len(surface_wl) surface = pyartdeco.surface(surface_name, surface_family, surface_kind, interp=surface_interp, \ wdspd = 2.0, shadow=True, test_glitter=True) ######################## # Geometry if case == "a1": if sub_case == 1: sza = 0.0 saa = 65.0 elif sub_case == 2: sza = 30.0 saa = 0 elif sub_case == 3: sza = 30.0 saa = 65 nvza = 17 nraa = 73 min_raa = 0.0 max_raa = 360.0 elif case == "a2": sza = 50.0 saa = 0.0 nvza = 17 nraa = 37 min_raa = 0.0 max_raa = 180.0 elif case in ["a3","a4"]: sza = 40.0 saa = 0.0 nvza = 17 nraa = 37 min_raa = 0.0 max_raa = 180.0 elif case == "a5": sza = 50.0 saa = 0.0 nvza = 81 if sub_case=="pp": nraa = 2 min_raa = 0.0 max_raa = 180.0 elif case == "a6": sza = 45.0 saa = 0.0 nvza = 17 nraa = 37 min_raa = 0.0 max_raa = 180.0 if trans : min_vza = 100. max_vza = 180. else: min_vza = 0 max_vza = 80. dvza = (max_vza - min_vza) / (nvza-1) vza = np.zeros(nvza) for j in range(nvza): vza[j] = max_vza - (j * dvza) #print(vza[j]) draa = (max_raa - min_raa) / (nraa-1) raa = np.zeros(nraa) vaa = np.zeros(nraa) for j in range(nraa): raa[j] = max_raa - (j * draa) vaa[j] = saa - (360.0 - raa[j]) if vaa[j] < 0: vaa[j] = vaa[j] + 360.0 #print(vaa[j], raa[j]) vaa = vaa[np.argsort(raa)] raa = raa[np.argsort(raa)] geom = pyartdeco.geometry(np.array([sza]), vza, raa) ############### # run ARTDECO lamb, rad, rad_levels, flux, alt = run_artdeco(artdeco_in, atmos, surface, solrad, ptcle, geom, kdis=kdis, verbose=True) ##(artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nraa, artdeco.mcommon.nmat, artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda) ##(artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nraa, artdeco.mcommon.nmat, artdeco.mcommon.nlambda) if trans: rad = rad_levels[0,:,:,:,-1,0] else: # normalisation as Kokhanovsky rad = rad[0,:,:,:,0] ##################### # Plot the result print("Plotting...") indvaa = np.argsort(vaa) rad = rad[:,indvaa,:] vaa = vaa[indvaa] code_ref = "ipol" if case =="a5": path_file = "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/iprt_case_"+case+"_"+sub_case+"_"+code_ref+".dat" else: path_file = "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/iprt_case_"+case+"_"+code_ref+".dat" print(path_file) iprt_val = np.genfromtxt(path_file, names=["depol","alt","sza","saa","vza","vaa","I","Q","U","V"]) # names=["depol","alt","sza","saa","vza","vaa","I","Q","U","V"]) I_ref = np.zeros((nvza,nraa)) Q_ref = np.zeros((nvza,nraa)) U_ref = np.zeros((nvza,nraa)) V_ref = np.zeros((nvza,nraa)) if trans: sza_iprt = iprt_val["sza"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] vza_iprt = iprt_val["vza"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] saa_iprt = iprt_val["saa"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] vaa_iprt = iprt_val["vaa"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] I_iprt = iprt_val["I"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] Q_iprt = iprt_val["Q"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] U_iprt = iprt_val["U"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] V_iprt = iprt_val["V"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] else: sza_iprt = iprt_val["sza"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] vza_iprt = iprt_val["vza"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] saa_iprt = iprt_val["saa"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] vaa_iprt = iprt_val["vaa"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] I_iprt = iprt_val["I"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] Q_iprt = iprt_val["Q"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] U_iprt = iprt_val["U"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] V_iprt = iprt_val["V"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] for ivaa in range(nraa): for ivza in range(nvza): #print(vza[ivza],vaa[ivaa]) #print([ (vza_iprt == vza[ivza]) & (vaa_iprt == vaa[ivaa]) ]) I_ref[ivza,ivaa] = I_iprt[ (vza_iprt == 180.0-vza[ivza]) & (vaa_iprt == vaa[ivaa]) ] Q_ref[ivza,ivaa] = Q_iprt[ (vza_iprt == 180.0-vza[ivza]) & (vaa_iprt == vaa[ivaa]) ] U_ref[ivza,ivaa] = U_iprt[ (vza_iprt == 180.0-vza[ivza]) & (vaa_iprt == vaa[ivaa]) ] V_ref[ivza,ivaa] = V_iprt[ (vza_iprt == 180.0-vza[ivza]) & (vaa_iprt == vaa[ivaa]) ] if trans: vza = 180-vza I_ref = np.transpose(I_ref) Q_ref = np.transpose(Q_ref) U_ref = np.transpose(U_ref) V_ref = np.transpose(V_ref) I = rad[:,:,0] I = np.transpose(I) if nmat>1: Q = rad[:,:,1] Q = np.transpose(Q) U = rad[:,:,2] U = np.transpose(U) if nmat>3: V = rad[:,:,3] V = np.transpose(V) # if case in ["a1","a2","a3","a4"]: # figI, axI, caxI = plot_polar_contour(I, vaa, vza, barname='$I$', cut=False) # #figI_ref, axI_ref, caxI_ref = plot_polar_contour(I_ref, vaa, vza, barname='$I_{%s}$'%code_ref, cut=False) # figdI, axdI, caxdI = plot_polar_contour(I-I_ref, vaa, vza, barname='$I-I_{%s}$'%code_ref, cut=False) # figQ, axQ, caxQ = plot_polar_contour(Q, vaa, vza, barname='$Q$', cut=False) # # figQ_ref, axQ_ref, caxQ_ref = plot_polar_contour(Q_ref, vaa, vza, barname='$Q_{%s}$'%code_ref, cut=False) # figdQ, axdQ, caxdQ = plot_polar_contour(Q-Q_ref, vaa, vza, barname='$Q-Q_{%s}$'%code_ref, cut=False) # figU, axU, caxU = plot_polar_contour(U, vaa, vza, barname='$U$', cut=False) # # figU_ref, axU_ref, caxU_ref = plot_polar_contour(U_ref, vaa, vza, barname='$U_{%s}$'%code_ref, cut=False) # figdU, axdU, caxdU = plot_polar_contour(U-U_ref, vaa, vza, barname='$U-U_{%s}$'%code_ref, cut=False) # plt.show() fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) ax1.set_ylabel("I") ax1.set_xlim([np.min(-vza),np.max(vza)]) for ivaa in range(nraa): if raa[ivaa] == 0.0: ax1.plot(vza,I_ref[ivaa,:],'b--', label = code_ref) ax1.plot(vza,I[ivaa,:],'b.-', label="ARTDECO") if raa[ivaa] == 180.0: ax1.plot(-vza,I_ref[ivaa,:],'b--') ax1.plot(-vza,I[ivaa,:],'b.-') if raa[ivaa] == 90.0: ax1.plot(vza,I_ref[ivaa,:],'r--') ax1.plot(vza,I[ivaa,:],'r.-') if raa[ivaa] == 270.0: ax1.plot(-vza,I_ref[ivaa,:],'r--') ax1.plot(-vza,I[ivaa,:],'r.-') if case == "a6": if raa[ivaa] == 10.0: ax1.plot(vza,I_ref[ivaa,:],'g--') ax1.plot(vza,I[ivaa,:],'g.-') if raa[ivaa] == 190.0: ax1.plot(-vza,I_ref[ivaa,:],'g--') ax1.plot(-vza,I[ivaa,:],'g.-') ax1.legend() # ax2.set_ylabel('$I-I_{%s}$'%code_ref) # ax2.set_xlabel("$VZA$") # ax2.plot([np.min(-vza),np.max(vza)],[0.0,0.0],'k') # for ivaa in range(nraa): # if raa[ivaa] == 0.0: # ax2.plot(vza,I[ivaa,:]-I_ref[ivaa,:],'b', label = "principal plane") # if raa[ivaa] == 180.0: # ax2.plot(-vza,I[ivaa,:]-I_ref[ivaa,:],'b') # if raa[ivaa] == 90.0: # ax2.plot(vza,I[ivaa,:]-I_ref[ivaa,:],'r', label = "Phi=90") # if raa[ivaa] == 270.0: # ax2.plot(-vza,I[ivaa,:]-I_ref[ivaa,:],'r') # if case=="a6": # if raa[ivaa] == 10.0: # ax2.plot(vza,I[ivaa,:]-I_ref[ivaa,:],'g', label = "Phi=10") # if raa[ivaa] == 190.0: # ax2.plot(-vza,I[ivaa,:]-I_ref[ivaa,:],'g') ax2.set_ylabel('($I-I_{%s}$)'%code_ref+'$/I_{%s}$'%code_ref) ax2.set_xlabel("$VZA$") ax2.plot([np.min(-vza),np.max(vza)],[0.0,0.0],'k') for ivaa in range(nraa): if raa[ivaa] == 0.0: ax2.plot(vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'b', label = "principal plane") if raa[ivaa] == 180.0: ax2.plot(-vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'b') if raa[ivaa] == 90.0: ax2.plot(vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'r', label = "Phi=90") if raa[ivaa] == 270.0: ax2.plot(-vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'r') if case=="a6": if raa[ivaa] == 10.0: ax2.plot(vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'g', label = "Phi=10") if raa[ivaa] == 190.0: ax2.plot(-vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'g') ax2.legend() plt.tight_layout() if nmat>1: fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) ax1.set_ylabel("Q") for ivaa in range(nraa): if raa[ivaa] == 0.0: ax1.plot(vza,Q_ref[ivaa,:],'b--', label = code_ref) ax1.plot(vza,Q[ivaa,:],'b', label="ARTDECO") if raa[ivaa] == 180.0: ax1.plot(-vza,Q_ref[ivaa,:],'b--') ax1.plot(-vza,Q[ivaa,:],'b') if raa[ivaa] == 90.0: ax1.plot(vza,Q_ref[ivaa,:],'r--') ax1.plot(vza,Q[ivaa,:],'r') if raa[ivaa] == 270.0: ax1.plot(-vza,Q_ref[ivaa,:],'r--') ax1.plot(-vza,Q[ivaa,:],'r') if case=="a6": if raa[ivaa] == 10.0: ax1.plot(vza,Q_ref[ivaa,:],'g--') ax1.plot(vza,Q[ivaa,:],'g') if raa[ivaa] == 190.0: ax1.plot(-vza,Q_ref[ivaa,:],'g--') ax1.plot(-vza,Q[ivaa,:],'g') ax1.legend() ax2.set_ylabel('$Q-Q_{%s}$'%code_ref) ax2.set_xlabel("$VZA$") for ivaa in range(nraa): if raa[ivaa] == 0.0: ax2.plot(vza,Q[ivaa,:]-Q_ref[ivaa,:],'b', label = "principal plane") if raa[ivaa] == 180.0: ax2.plot(-vza,Q[ivaa,:]-Q_ref[ivaa,:],'b') if raa[ivaa] == 90.0: ax2.plot(vza,Q[ivaa,:]-Q_ref[ivaa,:],'r', label = "Phi=90") if raa[ivaa] == 270.0: ax2.plot(-vza,Q[ivaa,:]-Q_ref[ivaa,:],'r') if case=="a6": if raa[ivaa] == 10.0: ax2.plot(vza,Q[ivaa,:]-Q_ref[ivaa,:],'g', label = "Phi=10") if raa[ivaa] == 190.0: ax2.plot(-vza,Q[ivaa,:]-Q_ref[ivaa,:],'g') ax2.legend() plt.tight_layout() fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) ax1.set_ylabel("U") for ivaa in range(nraa): if raa[ivaa] == 0.0: ax1.plot(vza,U_ref[ivaa,:],'b--', label = code_ref) ax1.plot(vza,U[ivaa,:],'b', label="ARTDECO") if raa[ivaa] == 180.0: ax1.plot(-vza,U_ref[ivaa,:],'b--') ax1.plot(-vza,U[ivaa,:],'b') if raa[ivaa] == 90.0: ax1.plot(vza,U_ref[ivaa,:],'r--') ax1.plot(vza,U[ivaa,:],'r') if raa[ivaa] == 270.0: ax1.plot(-vza,U_ref[ivaa,:],'r--') ax1.plot(-vza,U[ivaa,:],'r') if case=="a6": if raa[ivaa] == 10.0: ax1.plot(vza,U_ref[ivaa,:],'g--') ax1.plot(vza,U[ivaa,:],'g') if raa[ivaa] == 190.0: ax1.plot(-vza,U_ref[ivaa,:],'g--') ax1.plot(-vza,U[ivaa,:],'g') ax1.legend() ax2.set_ylabel('$U-U_{%s}$'%code_ref) ax2.set_xlabel("$VZA$") for ivaa in range(nraa): if raa[ivaa] == 0.0: ax2.plot(vza,U[ivaa,:]-U_ref[ivaa,:],'b', label = "principal plane") if raa[ivaa] == 180.0: ax2.plot(-vza,U[ivaa,:]-U_ref[ivaa,:],'b') if raa[ivaa] == 90.0: ax2.plot(vza,U[ivaa,:]-U_ref[ivaa,:],'r', label = "Phi=90") if raa[ivaa] == 270.0: ax2.plot(-vza,U[ivaa,:]-U_ref[ivaa,:],'r') if case=="a6": if raa[ivaa] == 10.0: ax2.plot(vza,U[ivaa,:]-U_ref[ivaa,:],'g', label = "Phi=10") if raa[ivaa] == 190.0: ax2.plot(-vza,U[ivaa,:]-U_ref[ivaa,:],'g') ax2.legend() plt.tight_layout() plt.show() def iprt_b(nstreams_init, trans, case_in, nmat, rt_model, corint, trunc_method, test_od_ims=False, hard=False): ########################################### # artdeco_in technical parameters structure keywords=[] wavel=[] mode = 'mono' if "b1" in case_in: keywords = ['nowarn', 'od_no_check'] sub_case = int(case_in.split('_')[1]) case = case_in.split('_')[0] elif "b2" in case_in: keywords = ['nowarn', 'od_no_check'] mode = 'lbl' sub_case = int(case_in.split('_')[1]) case = case_in.split('_')[0] elif "b3" in case_in: keywords = ['nowarn', 'od_no_check'] mode = 'lbl' sub_case = int(case_in.split('_')[1]) case = case_in.split('_')[0] elif "b4" in case_in: keywords = ['nowarn', 'od_no_check'] sub_case = int(case_in.split('_')[1]) case = case_in.split('_')[0] filters = ['none'] if case=="b1": wavel = [0.450] elif case=="b2": wavel = [-1,-1] elif case=="b3": wavel = [-1,-1] elif case=="b4": wavel = [0.800] thermal = False if test_od_ims: artdeco_in = pyartdeco.artdeco_in(keywords, mode, filters, wavel, \ trunc_method, nstreams_init, rt_model, \ corint, thermal, nmat, \ od_nstr=8, od_deltam=True, od_corint=True, od_do_secsca=True) else: artdeco_in = pyartdeco.artdeco_in(keywords, mode, filters, wavel, \ trunc_method, nstreams_init, rt_model, \ corint, thermal, nmat) #################### # read kdis coeff kdis = pyartdeco.kdis_coeff(artdeco_in, dir_artdeco+'lib/kdis/', 'ascii') ###################### # load solar TOA flux solrad = pyartdeco.solar_irradiance(artdeco_in, dir_artdeco+'lib/solrad/', 'ascii') #################### # load atmosphere gas = ['none'] ppmv = [-1] atm = '' dir_data = '' fmt_atm = "ascii" if case in ["b1"]: atm = 'atm_iprt_b1.dat' dir_data = '/home/mathieu/work/RTM/artdeco/pyartdeco/misc/benchmark/iprt/' fmt_atm = "ascii" elif case in ["b2"]: atm = 'atm_iprt_b2.dat' dir_data = '/home/mathieu/work/RTM/artdeco/pyartdeco/misc/benchmark/iprt/' fmt_atm = "lbl_iprt" elif case in ["b3"]: atm = 'atm_iprt_b3.dat' dir_data = '/home/mathieu/work/RTM/artdeco/pyartdeco/misc/benchmark/iprt/' fmt_atm = "lbl_iprt" elif case in ["b4"]: atm = 'atm_iprt_b4.dat' dir_data = '/home/mathieu/work/RTM/artdeco/pyartdeco/misc/benchmark/iprt/' fmt_atm = "ascii" wldepol = np.array([0.1, 200.0]) # microns taur = np.nan depol = np.nan if case == "b1": taur = 0.2188074939 depol = 0.03 if case in ["b2"]: taur = 0.8540404312 depol = 0.03 if case in ["b3"]: taur = 0.6230 depol = 0.03 if case in ["b4"]: taur = 0.02098185032 depol = 0.03 atmos = pyartdeco.atmosphere(artdeco_in, dir_data, atm, gas, ppmv, fmt_atm, wldepol, np.array([depol, depol]), interp=True, kdis=kdis, tau_ray = { "wvl":np.array([0.1,200.]), "tau":np.array([taur,taur]) }) ############################### # set ptcle structure if case in ["b1", "b2"]: ptcle_def = [] wlref = None wlptcle = None elif case in ["b3"]: tauaer, alt_aer, vprof_aer = get_vprof_aer_iprt_b3("/home/mathieu/work/RTM/artdeco/pyartdeco/misc/benchmark/iprt/tau_aerosol.dat") #print(alt_aer) #print(vprof_aer) ptcle_def = [ {"file_path":"/home/mathieu/work/RTM/artdeco/pyartdeco/misc/benchmark/iprt/opt_sizedistr_spheroid_iprt.h5", "name":"sizedistr_spheroid_iprt", "tau":tauaer, "alt_distrib":"user", "user_vdist":vprof_aer, "user_alt":alt_aer } ] wlref = 0.35 wlptcle = [0.1, 100.0] elif case in ["b4"]: ptcle_def = [ {"file_path":"/home/mathieu/work/RTM/artdeco/pyartdeco/misc/benchmark/iprt/opt_watercloud_iprt.h5", "name":"watercloud_iprt", "tau":5.0, "alt_distrib":"homogeneous", "z":3.0, "dz":1.0 } ] wlref = 0.800 wlptcle = [0.1, 100.0] ptcle_opt = pyartdeco.ptcle_optical_properties(ptcle_def, artdeco_in.nstreams, wlptcle, wlref, opt_interp=False, read_betal=False) ptcle = pyartdeco.particle(artdeco_in, ptcle_def, ptcle_opt) #print ptcle.__dict__.keys() #print ptcle.__dict__ ######################## # set surface structure if case in ["b1", "b2", "b3"]: name = "toto" family = "lambert" kind = "land" # wavelength definition for the surface wl = np.array([0.1, 200.0]) # microns alb = np.array([0.0, 0.0]) surface = pyartdeco.surface(name, family, kind, wl=wl, alb=alb, interp=True) if case in ["b4"]: surface_name = "surface" surface_interp = False surface_family = "brdf" surface_kind = "ocean" surface_wl = np.array([0.1, 100.0]) # microns surface_nwvl = len(surface_wl) surface = pyartdeco.surface(surface_name, surface_family, surface_kind, interp=surface_interp, \ wdspd = 2.0, shadow=True, test_glitter=True) ######################## # Geometry if case in ["b1","b2","b4"]: sza = 60.0 saa = 0.0 if case in ["b3"]: sza = 30.0 saa = 0.0 if case in ["b1","b2","b3","b4"]: nraa = 37 min_raa = 0.0 max_raa = 180.0 if sub_case == 1: nvza = 18 min_vza = 95. max_vza = 180. if not trans : print("not possible") return elif sub_case == 2: nvza = 18 if trans : min_vza = 95. max_vza = 180. else: min_vza = 0. max_vza = 85. elif sub_case == 3: nvza = 18 min_vza = 0. max_vza = 85. if trans : print("not possible") return dvza = (max_vza - min_vza) / (nvza-1) vza = np.zeros(nvza) for j in range(nvza): vza[j] = max_vza - (j * dvza) #print(vza[j]) draa = (max_raa - min_raa) / (nraa-1) raa = np.zeros(nraa) vaa = np.zeros(nraa) for j in range(nraa): raa[j] = max_raa - (j * draa) vaa[j] = saa - (360.0 - raa[j]) if vaa[j] < 0: vaa[j] = vaa[j] + 360.0 #print(vaa[j], raa[j]) vaa = vaa[np.argsort(raa)] raa = raa[np.argsort(raa)] geom = pyartdeco.geometry(np.array([sza]), vza, raa) ############### # run ARTDECO lamb, rad, rad_levels, flux, alt = run_artdeco(artdeco_in, atmos, surface, solrad, ptcle, geom, kdis=kdis, verbose=True) print("number of NaN values in rad = ",np.count_nonzero(np.isnan(rad))) ##(artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nraa, artdeco.mcommon.nmat, artdeco.mcommon.nalt_atm, artdeco.mcommon.nlambda) ##(artdeco.mcommon.nsza, artdeco.mcommon.nvza, artdeco.mcommon.nraa, artdeco.mcommon.nmat, artdeco.mcommon.nlambda) if trans: rad = rad_levels[0,:,:,:,-1,0] else: # normalisation as Kokhanovsky rad = rad[0,:,:,:,0] ##################### # Plot the result print("Plotting...") indvaa = np.argsort(vaa) rad = rad[:,indvaa,:] vaa = vaa[indvaa] code_ref = "ipol" #code_ref = "3dmcpol" #code_ref = "mystic" path_file = "/home/mathieu/work/RTM/artdeco/pyartdeco/misc/benchmark/iprt/iprt_case_"+case+"_"+code_ref+".dat" print(path_file) iprt_val = np.genfromtxt(path_file, names=["depol","alt","sza","saa","vza","vaa","I","Q","U","V"]) # names=["depol","alt","sza","saa","vza","vaa","I","Q","U","V"]) I_ref = np.zeros((nvza,nraa)) Q_ref = np.zeros((nvza,nraa)) U_ref = np.zeros((nvza,nraa)) V_ref = np.zeros((nvza,nraa)) if trans: if code_ref == "mystic": print(" check file format for trans and mystic...") return else: sza_iprt = iprt_val["sza"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] vza_iprt = iprt_val["vza"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] saa_iprt = iprt_val["saa"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] vaa_iprt = iprt_val["vaa"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] I_iprt = iprt_val["I"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] Q_iprt = iprt_val["Q"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] U_iprt = iprt_val["U"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] V_iprt = iprt_val["V"][(iprt_val["alt"]==0) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] else: if code_ref == "mystic": sza_iprt = iprt_val["sza"][(iprt_val["alt"]==30) & (iprt_val["sza"]==sza) ] vza_iprt = iprt_val["vza"][(iprt_val["alt"]==30) & (iprt_val["sza"]==sza) ] saa_iprt = iprt_val["saa"][(iprt_val["alt"]==30) & (iprt_val["sza"]==sza) ] vaa_iprt = iprt_val["vaa"][(iprt_val["alt"]==30) & (iprt_val["sza"]==sza)] I_iprt = iprt_val["I"][(iprt_val["alt"]==30) & (iprt_val["sza"]==sza)] Q_iprt = iprt_val["Q"][(iprt_val["alt"]==30) & (iprt_val["sza"]==sza)] U_iprt = iprt_val["U"][(iprt_val["alt"]==30) & (iprt_val["sza"]==sza)] V_iprt = iprt_val["V"][(iprt_val["alt"]==30) & (iprt_val["sza"]==sza)] else: sza_iprt = iprt_val["sza"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] vza_iprt = iprt_val["vza"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] saa_iprt = iprt_val["saa"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] vaa_iprt = iprt_val["vaa"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] I_iprt = iprt_val["I"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] Q_iprt = iprt_val["Q"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] U_iprt = iprt_val["U"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] V_iprt = iprt_val["V"][(iprt_val["alt"]==1) & (iprt_val["sza"]==sza) & (iprt_val["depol"]==depol)] if code_ref == "mystic" or code_ref == "3dmcpol": U_iprt = -U_iprt V_iprt = -V_iprt for ivaa in range(nraa): for ivza in range(nvza): #print(vza[ivza],vaa[ivaa]) #print([ (vza_iprt == vza[ivza]) & (vaa_iprt == vaa[ivaa]) ]) I_ref[ivza,ivaa] = I_iprt[ (vza_iprt == 180.0-vza[ivza]) & (vaa_iprt == vaa[ivaa]) ] Q_ref[ivza,ivaa] = Q_iprt[ (vza_iprt == 180.0-vza[ivza]) & (vaa_iprt == vaa[ivaa]) ] U_ref[ivza,ivaa] = U_iprt[ (vza_iprt == 180.0-vza[ivza]) & (vaa_iprt == vaa[ivaa]) ] V_ref[ivza,ivaa] = V_iprt[ (vza_iprt == 180.0-vza[ivza]) & (vaa_iprt == vaa[ivaa]) ] if trans: vza = 180-vza I_ref = np.transpose(I_ref) Q_ref = np.transpose(Q_ref) U_ref = np.transpose(U_ref) V_ref = np.transpose(V_ref) I = rad[:,:,0] I = np.transpose(I) if nmat>1: Q = rad[:,:,1] Q = np.transpose(Q) U = rad[:,:,2] U = np.transpose(U) if nmat>3: V = rad[:,:,3] V = np.transpose(V) fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) ax1.set_ylabel("I") ax1.set_xlim([np.min(-vza),np.max(vza)]) for ivaa in range(nraa): if raa[ivaa] == 0.0: ax1.plot(vza,I_ref[ivaa,:],'b--', label = code_ref) ax1.plot(vza,I[ivaa,:],'b.-', label="ARTDECO") if raa[ivaa] == 180.0: ax1.plot(-vza,I_ref[ivaa,:],'b--') ax1.plot(-vza,I[ivaa,:],'b.-') if raa[ivaa] == 90.0: ax1.plot(vza,I_ref[ivaa,:],'r--') ax1.plot(vza,I[ivaa,:],'r.-') if raa[ivaa] == 270.0: ax1.plot(-vza,I_ref[ivaa,:],'r--') ax1.plot(-vza,I[ivaa,:],'r.-') if case == "a6": if raa[ivaa] == 10.0: ax1.plot(vza,I_ref[ivaa,:],'g--') ax1.plot(vza,I[ivaa,:],'g.-') if raa[ivaa] == 190.0: ax1.plot(-vza,I_ref[ivaa,:],'g--') ax1.plot(-vza,I[ivaa,:],'g.-') ax1.legend() # ax2.set_ylabel('$I-I_{%s}$'%code_ref) # ax2.set_xlabel("$VZA$") # ax2.plot([np.min(-vza),np.max(vza)],[0.0,0.0],'k') # for ivaa in range(nraa): # if raa[ivaa] == 0.0: # ax2.plot(vza,I[ivaa,:]-I_ref[ivaa,:],'b', label = "principal plane") # if raa[ivaa] == 180.0: # ax2.plot(-vza,I[ivaa,:]-I_ref[ivaa,:],'b') # if raa[ivaa] == 90.0: # ax2.plot(vza,I[ivaa,:]-I_ref[ivaa,:],'r', label = "Phi=90") # if raa[ivaa] == 270.0: # ax2.plot(-vza,I[ivaa,:]-I_ref[ivaa,:],'r') # if case=="a6": # if raa[ivaa] == 10.0: # ax2.plot(vza,I[ivaa,:]-I_ref[ivaa,:],'g', label = "Phi=10") # if raa[ivaa] == 190.0: # ax2.plot(-vza,I[ivaa,:]-I_ref[ivaa,:],'g') ax2.set_ylabel('($I-I_{%s}$)'%code_ref+'$/I_{%s}$'%code_ref+"($\%$)") ax2.set_xlabel("$VZA$") ax2.plot([np.min(-vza),np.max(vza)],[0.0,0.0],'k') for ivaa in range(nraa): if raa[ivaa] == 0.0: ax2.plot(vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'b', label = "principal plane") if raa[ivaa] == 180.0: ax2.plot(-vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'b') if raa[ivaa] == 90.0: ax2.plot(vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'r', label = "Phi=90") if raa[ivaa] == 270.0: ax2.plot(-vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'r') if case=="a6": if raa[ivaa] == 10.0: ax2.plot(vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'g', label = "Phi=10") if raa[ivaa] == 190.0: ax2.plot(-vza,(I[ivaa,:]-I_ref[ivaa,:])/I_ref[ivaa,:]*100.,'g') ax2.legend() plt.tight_layout() if hard: plt.savefig("/tmp/fig_I_"+case_in+".png",dpi=150) if nmat>1: fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) ax1.set_ylabel("Q") for ivaa in range(nraa): if raa[ivaa] == 0.0: ax1.plot(vza,Q_ref[ivaa,:],'b--', label = code_ref) ax1.plot(vza,Q[ivaa,:],'b', label="ARTDECO") if raa[ivaa] == 180.0: ax1.plot(-vza,Q_ref[ivaa,:],'b--') ax1.plot(-vza,Q[ivaa,:],'b') if raa[ivaa] == 90.0: ax1.plot(vza,Q_ref[ivaa,:],'r--') ax1.plot(vza,Q[ivaa,:],'r') if raa[ivaa] == 270.0: ax1.plot(-vza,Q_ref[ivaa,:],'r--') ax1.plot(-vza,Q[ivaa,:],'r') if case=="a6": if raa[ivaa] == 10.0: ax1.plot(vza,Q_ref[ivaa,:],'g--') ax1.plot(vza,Q[ivaa,:],'g') if raa[ivaa] == 190.0: ax1.plot(-vza,Q_ref[ivaa,:],'g--') ax1.plot(-vza,Q[ivaa,:],'g') ax1.legend() ax2.set_ylabel('$Q-Q_{%s}$'%code_ref) ax2.set_xlabel("$VZA$") for ivaa in range(nraa): if raa[ivaa] == 0.0: ax2.plot(vza,Q[ivaa,:]-Q_ref[ivaa,:],'b', label = "principal plane") if raa[ivaa] == 180.0: ax2.plot(-vza,Q[ivaa,:]-Q_ref[ivaa,:],'b') if raa[ivaa] == 90.0: ax2.plot(vza,Q[ivaa,:]-Q_ref[ivaa,:],'r', label = "Phi=90") if raa[ivaa] == 270.0: ax2.plot(-vza,Q[ivaa,:]-Q_ref[ivaa,:],'r') if case=="a6": if raa[ivaa] == 10.0: ax2.plot(vza,Q[ivaa,:]-Q_ref[ivaa,:],'g', label = "Phi=10") if raa[ivaa] == 190.0: ax2.plot(-vza,Q[ivaa,:]-Q_ref[ivaa,:],'g') ax2.legend() plt.tight_layout() if hard: plt.savefig("/tmp/fig_Q_"+case_in+".png",dpi=150) fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) ax1.set_ylabel("U") for ivaa in range(nraa): if raa[ivaa] == 0.0: ax1.plot(vza,U_ref[ivaa,:],'b--', label = code_ref) ax1.plot(vza,U[ivaa,:],'b', label="ARTDECO") if raa[ivaa] == 180.0: ax1.plot(-vza,U_ref[ivaa,:],'b--') ax1.plot(-vza,U[ivaa,:],'b') if raa[ivaa] == 90.0: ax1.plot(vza,U_ref[ivaa,:],'r--') ax1.plot(vza,U[ivaa,:],'r') if raa[ivaa] == 270.0: ax1.plot(-vza,U_ref[ivaa,:],'r--') ax1.plot(-vza,U[ivaa,:],'r') if case=="a6": if raa[ivaa] == 10.0: ax1.plot(vza,U_ref[ivaa,:],'g--') ax1.plot(vza,U[ivaa,:],'g') if raa[ivaa] == 190.0: ax1.plot(-vza,U_ref[ivaa,:],'g--') ax1.plot(-vza,U[ivaa,:],'g') ax1.legend() ax2.set_ylabel('$U-U_{%s}$'%code_ref) ax2.set_xlabel("$VZA$") for ivaa in range(nraa): if raa[ivaa] == 0.0: ax2.plot(vza,U[ivaa,:]-U_ref[ivaa,:],'b', label = "principal plane") if raa[ivaa] == 180.0: ax2.plot(-vza,U[ivaa,:]-U_ref[ivaa,:],'b') if raa[ivaa] == 90.0: ax2.plot(vza,U[ivaa,:]-U_ref[ivaa,:],'r', label = "Phi=90") if raa[ivaa] == 270.0: ax2.plot(-vza,U[ivaa,:]-U_ref[ivaa,:],'r') if case=="a6": if raa[ivaa] == 10.0: ax2.plot(vza,U[ivaa,:]-U_ref[ivaa,:],'g', label = "Phi=10") if raa[ivaa] == 190.0: ax2.plot(-vza,U[ivaa,:]-U_ref[ivaa,:],'g') ax2.legend() plt.tight_layout() if hard: plt.savefig("/tmp/fig_U_"+case_in+".png",dpi=150) if not hard: plt.show() if __name__=='__main__': # nstr = 8 # #kokha_rayleigh(nstr) # kokha_cl_aer(nstr, "cl") # #kokha_cl_aer(nstr, "aer") get_atm_iprt_b1b4("benchmark/iprt/tau_rayleigh_450.dat", "benchmark/iprt/atm_iprt_b1.dat") nstr = 8 case = "b3_3" trans = False rt_model = "doad" nmat = 3 trunc_method = "dm" corint = True iprt_b(nstr, trans, case, nmat, rt_model, corint, trunc_method, hard=True) get_atm_iprt_b1b4("benchmark/iprt/tau_rayleigh_450.dat", "benchmark/iprt/atm_iprt_b1.dat") nstr = 8 case = "b4_3" trans = False rt_model = "doad" nmat = 3 trunc_method = "dm" corint = True iprt_b(nstr, trans, case, nmat, rt_model, corint, trunc_method, hard=True) exit() # get_atm_iprt_b2("/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/tau_rayleigh_325.dat", "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/tau_molabs_325.dat", "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/atm_iprt_b2.dat") # nstr = 8 # case = "b2_2" # trans = False # rt_model = "doad" # nmat = 3 # trunc_method="dm" # corint = False # iprt_b(nstr, trans, case, nmat, rt_model, corint, trunc_method) #get_atm_iprt_b3("/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/tau_rayleigh_350.dat", "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/tau_molabs_350.dat", "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/atm_iprt_b3.dat") nstr = 32 case = "b3_2" trans = False rt_model = "doad" nmat = 3 trunc_method="dm" corint = True iprt_b(nstr, trans, case, nmat, rt_model, corint, trunc_method, hard=True) #get_atm_iprt_b1b4("/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/tau_rayleigh_800.dat", "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/atm_iprt_b4.dat") nstr = 32 case = "b4_2" trans = False rt_model = "doad" nmat = 3 trunc_method="dm" corint = True iprt_b(nstr, trans, case, nmat, rt_model, corint, trunc_method, hard=True) # nstr = 32 # case = "a1_2" # trans = False # rt_model = "doad" # nmat = 3 # trunc_method="dm" # corint = False # iprt_a(nstr, trans, case, nmat, rt_model, corint, trunc_method) # nstr = 32 # case = "a6" # trans = False # rt_model = "doad" # nmat = 3 # trunc_method="dm" # corint = True # iprt_a(nstr, trans, case, nmat, rt_model, corint, trunc_method) # nstr = 16 # case = "a4" # trans = True # rt_model = "disort" # nmat = 1 # trunc_method="dm" # corint = True # iprt_a(nstr, trans, case, nmat, rt_model, corint, trunc_method) # nstr = 8 # case = "a4" # trans = False # rt_model = "doad" # nmat = 3 # trunc_method="dm" # corint = True # iprt_a(nstr, trans, case, nmat, rt_model, corint, trunc_method) # nstr = 32 # case = "a4" # trans = False # rt_model = "doad" # nmat = 3 # trunc_method="dm" # corint = True # iprt_a(nstr, trans, case, nmat, rt_model, corint, trunc_method) # nstr = 8 # case = "a5_pp" # trans = False # rt_model = "disort" # nmat = 1 # trunc_method="dm" # corint = True # iprt_a(nstr, trans, case, nmat, rt_model, corint, trunc_method) # nstr = 8 # case = "a5_pp" # trans = True # rt_model = "disort" # nmat = 1 # trunc_method="dm" # corint = True # iprt_a(nstr, trans, case, nmat, rt_model, corint, trunc_method) # nstr = 361 # case = "a5_pp" # trans = True # rt_model = "disort" # nmat = 1 # trunc_method="none" # corint = False # iprt_a(nstr, trans, case, nmat, rt_model, corint, trunc_method, test_od_ims=True) # nstr = 16 # case = "a5_pp" # trans = True # rt_model = "disort" # nmat = 1 # trunc_method ="dm" # corint =True # iprt_a(nstr, trans, case, nmat, rt_model, corint, trunc_method)