import numpy as np from ncutils import read_nc import pylab as plt import h5py cdffile = "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/waso.mie.cdf" read_nc(cdffile,'') wvl_orig = read_nc(cdffile,'wavelen') nwvl_orig = len(wvl_orig) reff_orig = read_nc(cdffile,'reff') nreff_orig = len(reff_orig ) nteta = read_nc(cdffile,'ntheta') teta = read_nc(cdffile,'theta') phase = read_nc(cdffile,'phase') ext = read_nc(cdffile,'ext') ssa = read_nc(cdffile,'ssa') nteta = nteta[0,0,0] ext = ext ssa = ssa wvl = wvl_orig nwvl = 1 teta = teta[0,0,0,0:nteta] P1 = phase[0,0,0,0:nteta] P2 = phase[0,0,1,0:nteta] P3 = phase[0,0,2,0:nteta] P4 = phase[0,0,3,0:nteta] # plt.figure() # plt.xlim([0,180.]) # plt.ylim([0.1,10]) # plt.semilogy(teta,P1) # plt.figure() # plt.xlim([0,180.]) # plt.plot(teta,P2/P1) # plt.figure() # plt.xlim([0,180.]) # plt.plot(teta,P3/P1) # plt.figure() # plt.xlim([0,180.]) # plt.plot(teta,P4/P1) # plt.show() dir_opt = "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/" # Write opt_ file to be used for th betal expension fopt = open(dir_opt+'opt_waso_iprt.dat', 'w') fopt.write('# Optical properties to be used by ARTDECO \n') fopt.write('# \n') fopt.write('# Used model to obtain that properties is: \n') fopt.write(' unknown \n') fopt.write('# Used material is: \n') fopt.write(' unknown \n') fopt.write('# Number of phase matrix elements : \n') fopt.write(' 4 \n') fopt.write('# number of wavelengths \n') fopt.write('%i'%nwvl+' \n') fopt.write('# lambda(microns) nteta Cext (microns^2) SSA g \n') s = ' %.7e'%wvl+' %i'%nteta+' %.7e'%ext+' %.7e'%ssa+' -32768 \n' fopt.write(s) fopt.write('# Phase matrix \n') fopt.write('# lambda = %.4f'%wvl+' \n') fopt.write('# u F11 F44 F21 F34 \n') for iang in range(nteta): s = ' %17.9e'%(np.cos(teta[iang]*np.pi/180.))+' %17.9e'%P1[iang]+' %17.9e'%P3[iang]+' %17.9e'%P2[iang]+' %17.9e'%P4[iang]+' \n' fopt.write(s) fopt.close() dir_opt = "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/" # Write opt_ file to be used for th betal expension fopt = h5py.File(dir_opt+'opt_waso_iprt.h5', 'w') grp = fopt.require_group("waso_iprt") subgrp = grp.require_group("axis") dset = subgrp.require_dataset("wavelengths", (len(wvl),), dtype='float64') dset[...] = wvl[:] dset = subgrp.require_dataset("phase_angles", (len(teta),), dtype='float64') dset[...] = teta[:] dset = subgrp.require_dataset("mu", (len(teta),), dtype='float64') dset[...] = np.cos(teta[:]*np.pi/180.) subgrp = grp.require_group("data") dset = subgrp.require_dataset("single_scattering_albedo", (len(wvl),), dtype='float64') dset[...] = ssa[:] dset.attrs.create("dimensions", np.string_("wavelengths")) dset = subgrp.require_dataset("Cext", (len(wvl),), dtype='float64') dset[...] = ext[:] dset.attrs.create("dimensions", np.string_("wavelengths")) dset = subgrp.require_dataset("p11_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0] = P1 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) dset = subgrp.require_dataset("p21_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0] = P2 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) dset = subgrp.require_dataset("p34_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0] = P4 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) dset = subgrp.require_dataset("p44_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0]= P3 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) fopt.close() cdffile = "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/watercloud.mie.cdf" read_nc(cdffile,'') wvl_orig = read_nc(cdffile,'wavelen') nwvl_orig = len(wvl_orig) reff_orig = read_nc(cdffile,'reff') nreff_orig = len(reff_orig ) nteta = read_nc(cdffile,'ntheta') teta = read_nc(cdffile,'theta') phase = read_nc(cdffile,'phase') ext = read_nc(cdffile,'ext') ssa = read_nc(cdffile,'ssa') nteta = nteta[0,0,0] ext = ext ssa = ssa wvl = wvl_orig nwvl = 1 teta = teta[0,0,0,0:nteta] P1 = phase[0,0,0,0:nteta] P2 = phase[0,0,1,0:nteta] P3 = phase[0,0,2,0:nteta] P4 = phase[0,0,3,0:nteta] # plt.figure() # plt.xlim([0,180.]) # plt.ylim([0.01,10000]) # plt.semilogy(teta,P1) # plt.figure() # plt.xlim([0,180.]) # plt.plot(teta,P2/P1) # plt.figure() # plt.xlim([0,180.]) # plt.plot(teta,P3/P1) # plt.figure() # plt.xlim([0,180.]) # plt.plot(teta,P4/P1) # plt.show() dir_opt = "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/" # Write opt_ file to be used for th betal expension fopt = open(dir_opt+'opt_watercloud_iprt.dat', 'w') fopt.write('# Optical properties to be used by ARTDECO \n') fopt.write('# \n') fopt.write('# Used model to obtain that properties is: \n') fopt.write(' unknown \n') fopt.write('# Used material is: \n') fopt.write(' unknown \n') fopt.write('# Number of phase matrix elements : \n') fopt.write(' 4 \n') fopt.write('# number of wavelengths \n') fopt.write('%i'%nwvl+' \n') fopt.write('# lambda(microns) nteta Cext (microns^2) SSA g \n') s = ' %.7e'%wvl+' %i'%nteta+' %.7e'%ext+' %.7e'%ssa+' -32768 \n' fopt.write(s) fopt.write('# Phase matrix \n') fopt.write('# lambda = %.4f'%wvl+' \n') fopt.write('# u F11 F44 F21 F34 \n') for iang in range(nteta): s = ' %17.9e'%(np.cos(teta[iang]*np.pi/180.))+' %17.9e'%P1[iang]+' %17.9e'%P3[iang]+' %17.9e'%P2[iang]+' %17.9e'%P4[iang]+' \n' fopt.write(s) fopt.close() dir_opt = "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/" # Write opt_ file to be used for th betal expension fopt = h5py.File(dir_opt+'opt_watercloud_iprt.h5', 'w') grp = fopt.require_group("watercloud_iprt") subgrp = grp.require_group("axis") dset = subgrp.require_dataset("wavelengths", (len(wvl),), dtype='float64') dset[...] = wvl[:] dset = subgrp.require_dataset("phase_angles", (len(teta),), dtype='float64') dset[...] = teta[:] dset = subgrp.require_dataset("mu", (len(teta),), dtype='float64') dset[...] = np.cos(teta[:]*np.pi/180.) subgrp = grp.require_group("data") dset = subgrp.require_dataset("single_scattering_albedo", (len(wvl),), dtype='float64') dset[...] = ssa[:] dset.attrs.create("dimensions", np.string_("wavelengths")) dset = subgrp.require_dataset("Cext", (len(wvl),), dtype='float64') dset[...] = ext[:] dset.attrs.create("dimensions", np.string_("wavelengths")) dset = subgrp.require_dataset("p11_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0] = P1 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) dset = subgrp.require_dataset("p21_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0] = P2 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) dset = subgrp.require_dataset("p34_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0] = P4 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) dset = subgrp.require_dataset("p44_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0]= P3 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) fopt.close() cdffile = "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/sizedistr_spheroid.cdf" read_nc(cdffile,'') wvl_orig = read_nc(cdffile,'wavelen') nwvl_orig = len(wvl_orig) reff_orig = read_nc(cdffile,'reff') nreff_orig = len(reff_orig ) nteta = read_nc(cdffile,'ntheta') teta = read_nc(cdffile,'theta') phase = read_nc(cdffile,'phase') ext = read_nc(cdffile,'ext') ssa = read_nc(cdffile,'ssa') nteta = nteta[0,0,0] ext = ext ssa = ssa wvl = wvl_orig nwvl = 1 teta = teta[0,0,0,0:nteta] P1 = phase[0,0,0,0:nteta] P2 = phase[0,0,1,0:nteta] P3 = phase[0,0,2,0:nteta] P4 = phase[0,0,3,0:nteta] P5 = phase[0,0,4,0:nteta] P6 = phase[0,0,5,0:nteta] # plt.figure() # plt.xlim([0,180.]) # plt.ylim([0.01,10000]) # plt.semilogy(teta,P1) # plt.figure() # plt.xlim([0,180.]) # plt.plot(teta,P2/P1) # plt.figure() # plt.xlim([0,180.]) # plt.plot(teta,P3/P1) # plt.figure() # plt.xlim([0,180.]) # plt.plot(teta,P4/P1) # plt.figure() # plt.xlim([0,180.]) # plt.plot(teta,P5/P1) # plt.figure() # plt.xlim([0,180.]) # plt.plot(teta,P6/P1) # plt.show() dir_opt = "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/" # Write opt_ file to be used for th betal expension fopt = open(dir_opt+'opt_sizedistr_spheroid_iprt.dat', 'w') fopt.write('# Optical properties to be used by ARTDECO \n') fopt.write('# \n') fopt.write('# Used model to obtain that properties is: \n') fopt.write(' unknown \n') fopt.write('# Used material is: \n') fopt.write(' unknown \n') fopt.write('# Number of phase matrix elements : \n') fopt.write(' 6 \n') fopt.write('# number of wavelengths \n') fopt.write('%i'%nwvl+' \n') fopt.write('# lambda(microns) nteta Cext (microns^2) SSA g \n') s = ' %.7e'%wvl+' %i'%nteta+' %.7e'%ext+' %.7e'%ssa+' -32768 \n' fopt.write(s) fopt.write('# Phase matrix \n') fopt.write('# lambda = %.4f'%wvl+' \n') fopt.write('# u F11 F22 F33 F44 F21 F34 \n') for iang in range(nteta): s = ' %17.9e'%(np.cos(teta[iang]*np.pi/180.))+' %17.9e'%P1[iang]+' %17.9e'%P5[iang]+' %17.9e'%P3[iang]+' %17.9e'%P6[iang]+' %17.9e'%P2[iang]+' %17.9e'%P4[iang]+' \n' fopt.write(s) fopt.close() dir_opt = "/home/mathieu/work/RTM/artdeco/pyartdeco/f2py_utils/benchmark/iprt/" # Write opt_ file to be used for th betal expension fopt = h5py.File(dir_opt+'opt_sizedistr_spheroid_iprt.h5', 'w') grp = fopt.require_group("sizedistr_spheroid_iprt") subgrp = grp.require_group("axis") dset = subgrp.require_dataset("wavelengths", (len(wvl),), dtype='float64') dset[...] = wvl[:] dset = subgrp.require_dataset("phase_angles", (len(teta),), dtype='float64') dset[...] = teta[:] dset = subgrp.require_dataset("mu", (len(teta),), dtype='float64') dset[...] = np.cos(teta[:]*np.pi/180.) subgrp = grp.require_group("data") dset = subgrp.require_dataset("single_scattering_albedo", (len(wvl),), dtype='float64') dset[...] = ssa[:] dset.attrs.create("dimensions", np.string_("wavelengths")) dset = subgrp.require_dataset("Cext", (len(wvl),), dtype='float64') dset[...] = ext[:] dset.attrs.create("dimensions", np.string_("wavelengths")) dset = subgrp.require_dataset("p11_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0] = P1 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) dset = subgrp.require_dataset("p22_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0] = P5 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) dset = subgrp.require_dataset("p33_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0] = P3 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) dset = subgrp.require_dataset("p21_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0] = P2 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) dset = subgrp.require_dataset("p34_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0] = P4 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) dset = subgrp.require_dataset("p44_phase_function", (len(teta),len(wvl)), dtype='float64') dset[:,0]= P6 dset.attrs.create("dimensions", np.string_("mu,wavelengths")) fopt.close()