import os import sys import numpy as np import pylab as pl import string import h5py dir_pyartdeco = os.environ['DIR_PYARTDECO'] sys.path.append(dir_pyartdeco+'tools') import pyartdeco_utils as ad dir_artdeco = os.environ['DIR_ARTDECO'] dir_opt = dir_artdeco+'lib/opt/' # # This routine allows for mixing of particle # optical properties # # This will access opt_[...].dat files in # ARTDECO libraries and mix the optical # properties given a particle number mixing ratio # # NOTE: prop_[...].dat files must exist # for all components # # NOTE II : opt MUST have been computed using # MIE code in ARTDECO # delta_sumf11 = 1e-3 def compute_mix_opt(wvl, fraction, ptcles, ptcle_name): nptcle = len(ptcles) nwvl = len(wvl) # normalization de fraction fraction = fraction / np.sum(fraction) # first compute all components optical properties keywords = 'verbose opt_only' saveroot = 'mix_opt' Qext = np.zeros((nptcle, nwvl)) ssa = np.zeros((nptcle, nwvl)) g = np.zeros((nptcle, nwvl)) Qext_eff = np.zeros(nwvl) Qabs_eff = np.zeros(nwvl) Qsca_eff = np.zeros(nwvl) ssa_eff = np.zeros(nwvl) g_eff = np.zeros(nwvl) for iptcle in xrange(nptcle): artdeco_in=[] artdeco_in.append('# \n') artdeco_in.append('# MAIN INPUT FILE FOR ARTDECO PROGRAM \n') artdeco_in.append('# \n') artdeco_in.append('######################## \n') artdeco_in.append('# keywords (Ex: none, verbose,...) \n') artdeco_in.append(keywords+' \n') artdeco_in.append('######################## \n') artdeco_in.append('# outfiles root name \n') artdeco_in.append(saveroot+' \n') artdeco_in.append('######################## \n') artdeco_in.append('# mode \n') artdeco_in.append(' mono \n') artdeco_in.append('#######################\n') artdeco_in.append('# filters\n') artdeco_in.append(' none \n') artdeco_in.append('######################## \n') artdeco_in.append('# Wavelengths (microns) \n') s= '' for iwl in xrange(nwvl): s = s+' %.5f'%wvl[iwl] artdeco_in.append(s+' \n') artdeco_in.append('###################### \n') artdeco_in.append('# Particles \n') artdeco_in.append('# type interp. H-G Tau_550 vertical distribution type vdist parameters (km) \n') artdeco_in.append(ptcles[iptcle]+' no \n') artdeco_in.append('# \n') # write artdeco_in.dat f = open(dir_artdeco+'input/artdeco_in_mix_opt.dat', 'w') for j in xrange(len(artdeco_in)): f.write(artdeco_in[j]) f.close() # run ARTDECO one time os.system(dir_artdeco.strip()+"src/artdeco artdeco_in_mix_opt.dat") # retrieve optical properties for iwvl in xrange(nwvl): print wvl[iwvl] ang_tmp, p11_tmp, p44_tmp, p21_tmp, p34_tmp, \ Qext[iptcle,iwvl], ssa[iptcle,iwvl], g[iptcle,iwvl] = ad.get_opt(ptcles[iptcle], wvl[iwvl]) print 'done' if (iptcle == 0) and (iwvl == 0) : nang = len(ang_tmp) ang = np.zeros(nang) p11 = np.zeros((nang, nptcle, nwvl)) p44 = np.zeros((nang, nptcle, nwvl)) p21 = np.zeros((nang, nptcle, nwvl)) p34 = np.zeros((nang, nptcle, nwvl)) p11_eff = np.zeros((nang, nwvl)) p44_eff = np.zeros((nang, nwvl)) p21_eff = np.zeros((nang, nwvl)) p34_eff = np.zeros((nang, nwvl)) ang[:] = ang_tmp p11[:,iptcle,iwvl] = p11_tmp p44[:,iptcle,iwvl] = p44_tmp p21[:,iptcle,iwvl] = p21_tmp p34[:,iptcle,iwvl] = p34_tmp # check if angle def is the same for all wavelength and all particle size dist. if (nang!=len(ang_tmp)): print '' print 'There is a problem with nang' print '' exit(0) for iang in xrange(nang): if ang_tmp[iang]!=ang[iang]: print '' print 'There is a problem with ang' print '' exit(0) Qsca = ssa * Qext Qabs = (1.0-ssa) * Qext for iwvl in xrange(nwvl): Qsca_eff[iwvl] = np.sum(Qsca[:,iwvl] * fraction) Qabs_eff[iwvl] = np.sum(Qabs[:,iwvl] * fraction) Qext_eff[iwvl] = np.sum(Qext[:,iwvl] * fraction) ssa_eff[iwvl] = Qsca_eff[iwvl] / Qext_eff[iwvl] frac_sca = fraction * Qsca[:,iwvl] frac_sca = frac_sca / np.sum(frac_sca) for iang in xrange(nang): p11_eff[iang,iwvl] = np.sum(frac_sca * p11[iang,:,iwvl]) p44_eff[iang,iwvl] = np.sum(frac_sca * p44[iang,:,iwvl]) p21_eff[iang,iwvl] = np.sum(frac_sca * p21[iang,:,iwvl]) p34_eff[iang,iwvl] = np.sum(frac_sca * p34[iang,:,iwvl]) sumf11 = np.trapz(p11_eff[:,iwvl], np.cos(ang*np.pi/180.0)) if abs(sumf11 - 2.00) > delta_sumf11 : print 'sumf 11 = ', sumf11, ' for wvl=', wvl[iwvl] #pl.plot(ang, p11_eff[:,0]) #pl.show() # Write opt_ file to be used for th betal expension fopt = open(dir_artdeco+'lib/opt/opt_'+ptcle_name+'.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') for iwvl in xrange(nwvl): s = ' %.7e'%wvl[iwvl]+' %i'%nang+' %.7e'%Qext_eff[iwvl]+' %.7e'%ssa_eff[iwvl]+' -32768 \n' fopt.write(s) fopt.write('# Phase matrix \n') for iwvl in xrange(nwvl): fopt.write('# lambda = %.4f'%wvl[iwvl]+' \n') fopt.write('# u F11 F44 F33 F21 F34 \n') for iang in xrange(nang): s = ' %16.8e'%(np.cos(ang[iang]*np.pi/180.))+' %15.7e'%p11_eff[iang, iwvl]+' %15.7e'%p44_eff[iang, iwvl]+' %15.7e'%p21_eff[iang, iwvl]+' %15.7e'%p34_eff[iang, iwvl]+' \n' fopt.write(s) fopt.close() return if __name__=='__main__': #wvl = [0.410, 0.443, 0.490, 0.550, 0.555, 0.650, 0.670, 0.752, 0.763, 0.765, 0.865, 0.910, 0.914, 1.24, 1.365, 1.370, 1.375, 1.63, 1.650, 2.13, 2.25, 3.74, 3.959, 4.04, 4.05, 6.725, 7.325, 8.54, 10.69, 12.02, 13.345] # wavelength that will be computed wvl_min = 0.40 wvl_max = 0.70 delta_wvl = 0.05 nwvl = int( (wvl_max - wvl_min) / delta_wvl) + 2 wvl = np.zeros(nwvl) for iwvl in xrange(nwvl): wvl[iwvl] = wvl_min + delta_wvl*iwvl #print iwvl, wvl[iwvl] print wvl # OPAC Cumulus Maritime ptcle_name = 'opac_cumulus_maritime' fraction = np.array([1.0]) # this must be a fraction IN NUMBER # this will be normalized to 1 in the # code ptcles = ['opac_cumuma'] compute_mix_opt(wvl, fraction, ptcles, ptcle_name) # ## See OPAC (Hess, 1998 paper) # ptcle_name = 'opac_maritime_clean' # fraction = np.array([1500, 20, 3.2e-3]) # this must be a fraction IN NUMBER # # this will be normalized to 1 in the # # code # ptcles = ['opac_waso95','opac_ssam95','opac_sscm95'] # compute_mix_opt(wvl, fraction, ptcles, ptcle_name) # ptcle_name = 'opac_cont_avg' # fraction = np.array([ 7000, 0.4, 8300]) # this must be a fraction IN NUMBER # # this will be normalized to 1 in the # # code # ptcles = ['opac_waso50', 'opac_inso', 'opac_soot'] # compute_mix_opt(wvl, fraction, ptcles, ptcle_name) # ptcle_name = 'opac_urban' # fraction = np.array([28000, 1.5, 130000]) # this must be a fraction IN NUMBER # # this will be normalized to 1 in the # # code # ptcles = ['opac_waso50', 'opac_inso', 'opac_soot'] # compute_mix_opt(wvl, fraction, ptcles, ptcle_name) # ptcle_name = 'opac_desert' # fraction = np.array([2000, 269.5, 30.5, 0.142]) # this must be a fraction IN NUMBER # # this will be normalized to 1 in the # # code # ptcles = ['opac_waso50', 'opac_minm', 'opac_miam', 'opac_micm'] # compute_mix_opt(wvl, fraction, ptcles, ptcle_name)