#!/usr/bin/env python from __future__ import print_function, division, absolute_import from builtins import range import sys import os import string import numpy as np import h5py from ncutils import read_nc from scipy.interpolate import interp1d from itertools import product import trunc dir_opt = '/rfs/proj/artdeco_lib/opt/' dir_data = dir_opt delta_sumf11 = 1e-3 trunc_model ='dm' ######################################################################################################### wl_interp = 'linear' flag_interp_ang = False # where slinear, quadratic and cubic refer to a spline interpolation # of first, second or third order) or as an integer specifying the order of # the spline interpolator to use. Default is linear. ##################################################### baum_dir = '/rfs/data/baum_opt/' cdffile = 'GeneralHabitMixture_SeverelyRough_AllWavelengths_FullPhaseMatrix.nc' ang_s = read_nc(baum_dir+cdffile,'phase_angles') ang_s = ang_s.astype(np.float64) ang = np.append(ang_s[ang_s <= 80.0], np.linspace(80.5, 175.5, num=191, endpoint=True)) ang = np.append(ang, ang_s[ang_s >= 176.0]) ang = np.sort(ang)[::-1] mu = np.cos(ang*np.pi/180.0) nang = len(ang) ##################################################### def doit(kind, wvl_bound, nbetal, res_suffixe=""): n_nbetal = len(nbetal) nbetalmax = np.amax(nbetal) if kind == "ghm": cdffile = '/rfs/data/baum_opt/GeneralHabitMixture_SeverelyRough_AllWavelengths_FullPhaseMatrix.nc' elif kind == "sc": cdffile = '/rfs/data/baum_opt/SolidColumns_SeverelyRough_AllWavelengths_FullPhaseMatrix.nc' elif kind == "asc": cdffile = '/rfs/data/baum_opt/AggregateSolidColumns_SeverelyRough_AllWavelengths_FullPhaseMatrix.nc' res_file = dir_opt+'opt_baum_ice_'+kind if len(res_suffixe)>0: res_file = res_file+"_"+res_suffixe res_file = res_file+'.h5' if os.path.isfile(res_file): f = h5py.File(res_file, 'r') if kind in f.keys(): print( kind+" already present in res file") f.close return else: f.close read_nc(cdffile,'') # -> wavelengths - (445,) float32 # -> effective_diameter - (23,) float32 # -> ice_water_content - (23, 445) float32 # -> total_area - (23, 445) float32 # -> asymmetry_parameter - (23, 445) float32 # -> single_scattering_albedo - (23, 445) float32 # -> extinction_efficiency - (23, 445) float32 # -> extinction_coefficient_over_iwc - (23, 445) float32 # -> p11_phase_function - (498, 23, 445) float32 # -> p21_phase_function - (498, 23, 445) float32 # -> p22_phase_function - (498, 23, 445) float32 # -> p33_phase_function - (498, 23, 445) float32 # -> p43_phase_function - (498, 23, 445) float32 # -> p44_phase_function - (498, 23, 445) float32 # -> phase_angles - (498,) float32 wvl_orig = read_nc(cdffile,'wavelengths') deff_orig = read_nc(cdffile,'effective_diameter') ndeff_orig = len(deff_orig ) nreff_orig = ndeff_orig reff_orig = deff_orig / 2.0 iwc_orig = read_nc(cdffile,'ice_water_content') tot_aera_orig = read_nc(cdffile,'total_area') g_orig = read_nc(cdffile,'asymmetry_parameter') ssa_orig = read_nc(cdffile,'single_scattering_albedo') Qext_orig = read_nc(cdffile,'extinction_efficiency') Sext_o_IWC_orig = read_nc(cdffile,'extinction_coefficient_over_iwc') p11_orig = read_nc(cdffile,'p11_phase_function') p21_orig = read_nc(cdffile,'p21_phase_function') p21_orig = -p21_orig p22_orig = read_nc(cdffile,'p22_phase_function') p33_orig = read_nc(cdffile,'p33_phase_function') p43_orig = read_nc(cdffile,'p43_phase_function') p44_orig = read_nc(cdffile,'p44_phase_function') ang_orig = read_nc(cdffile,'phase_angles') ang_orig = ang_orig.astype(np.float64) #print ang_orig #print np.cos(ang_orig*np.pi/180.0) #print type(ang_orig) nang_orig = len(ang_orig ) good = np.where((wvl_orig > wvl_bound[0]) & (wvl_orig < wvl_bound[1])) if wvl_bound[0] > 0.55: iwvl0p55 = [np.squeeze(np.argwhere(wvl_orig==0.55))] print(iwvl0p55) wvl_0p55 = wvl_orig[iwvl0p55] iwc_0p55 = iwc_orig[:,iwvl0p55] tot_aera_0p55 = tot_aera_orig[:,iwvl0p55] g_0p55 = g_orig[:,iwvl0p55] ssa_0p55 = ssa_orig[:,iwvl0p55] Qext_0p55 = Qext_orig[:,iwvl0p55] Sext_o_IWC_0p55 = Sext_o_IWC_orig[:,iwvl0p55] p11_0p55 = p11_orig[:,:,iwvl0p55] p21_0p55 = p21_orig[:,:,iwvl0p55] p22_0p55 = p22_orig[:,:,iwvl0p55] p33_0p55 = p33_orig[:,:,iwvl0p55] p43_0p55 = p43_orig[:,:,iwvl0p55] p44_0p55 = p44_orig[:,:,iwvl0p55] iwvl1 = good[0][0] iwvl2 = good[0][-1]+1 wvl_orig = wvl_orig[iwvl1:iwvl2] iwc_orig = iwc_orig[:,iwvl1:iwvl2] tot_aera_orig = tot_aera_orig[:,iwvl1:iwvl2] g_orig = g_orig[:,iwvl1:iwvl2] ssa_orig = ssa_orig[:,iwvl1:iwvl2] Qext_orig = Qext_orig[:,iwvl1:iwvl2] Sext_o_IWC_orig = Sext_o_IWC_orig[:,iwvl1:iwvl2] p11_orig = p11_orig[:,:,iwvl1:iwvl2] p21_orig = p21_orig[:,:,iwvl1:iwvl2] p22_orig = p22_orig[:,:,iwvl1:iwvl2] p33_orig = p33_orig[:,:,iwvl1:iwvl2] p43_orig = p43_orig[:,:,iwvl1:iwvl2] p44_orig = p44_orig[:,:,iwvl1:iwvl2] if wvl_bound[0] > 0.55: print(iwc_0p55.shape) print(iwc_orig.shape) wvl_orig = np.append(0.55, wvl_orig) iwc_orig = np.append(iwc_0p55, iwc_orig, axis=1) tot_aera_orig = np.append(tot_aera_0p55,tot_aera_orig, axis=1) g_orig = np.append(g_0p55,g_orig, axis=1) ssa_orig = np.append(ssa_0p55,ssa_orig, axis=1) Qext_orig = np.append(Qext_0p55,Qext_orig,axis=1) Sext_o_IWC_orig = np.append(Sext_o_IWC_0p55,Sext_o_IWC_orig, axis=1) p11_orig = np.append(p11_0p55,p11_orig, axis=2) p21_orig = np.append(p21_0p55,p21_orig, axis=2) p22_orig = np.append(p22_0p55,p22_orig, axis=2) p33_orig = np.append(p33_0p55,p33_orig, axis=2) p43_orig = np.append(p43_0p55,p43_orig, axis=2) p44_orig = np.append(p44_0p55,p44_orig, axis=2) nwvl_orig = len(wvl_orig ) p22_orig = p22_orig * p11_orig p33_orig = p33_orig * p11_orig p44_orig = p44_orig * p11_orig p43_orig = p43_orig * p11_orig p21_orig = p21_orig * p11_orig p11_orig = p11_orig[::-1,:,:] p22_orig = p22_orig[::-1,:,:] p33_orig = p33_orig[::-1,:,:] p44_orig = p44_orig[::-1,:,:] p43_orig = p43_orig[::-1,:,:] p21_orig = p21_orig[::-1,:,:] ang_orig = ang_orig[::-1] # interpolate to a new angular grid p11_orig = interp1d(ang_orig, p11_orig, axis=0)(ang) p22_orig = interp1d(ang_orig, p22_orig, axis=0)(ang) p33_orig = interp1d(ang_orig, p33_orig, axis=0)(ang) p44_orig = interp1d(ang_orig, p44_orig, axis=0)(ang) p43_orig = interp1d(ang_orig, p43_orig, axis=0)(ang) p21_orig = interp1d(ang_orig, p21_orig, axis=0)(ang) ang_orig = ang nang_orig = nang Cext_orig = np.zeros((ndeff_orig, nwvl_orig)) trunc_coeff = np.zeros((n_nbetal, ndeff_orig, nwvl_orig)) p11_recomp_integrate = np.zeros((n_nbetal, ndeff_orig, nwvl_orig)) p11_integrate = np.zeros((ndeff_orig, nwvl_orig)) alpha1 = np.zeros((nbetalmax+1, n_nbetal, ndeff_orig, nwvl_orig)) alpha2 = np.zeros((nbetalmax+1, n_nbetal, ndeff_orig, nwvl_orig)) alpha3 = np.zeros((nbetalmax+1, n_nbetal, ndeff_orig, nwvl_orig)) alpha4 = np.zeros((nbetalmax+1, n_nbetal, ndeff_orig, nwvl_orig)) beta1 = np.zeros((nbetalmax+1, n_nbetal, ndeff_orig, nwvl_orig)) beta2 = np.zeros((nbetalmax+1, n_nbetal, ndeff_orig, nwvl_orig)) p11_recomp = np.zeros((nang_orig, n_nbetal, ndeff_orig, nwvl_orig)) p22_recomp = np.zeros((nang_orig, n_nbetal, ndeff_orig, nwvl_orig)) p33_recomp = np.zeros((nang_orig, n_nbetal, ndeff_orig, nwvl_orig)) p44_recomp = np.zeros((nang_orig, n_nbetal, ndeff_orig, nwvl_orig)) p21_recomp = np.zeros((nang_orig, n_nbetal, ndeff_orig, nwvl_orig)) p34_recomp = np.zeros((nang_orig, n_nbetal, ndeff_orig, nwvl_orig)) #print("Write opt files...") p11_integrate[:,:] = np.trapz(p11_orig[:,:,:], x=np.cos(ang_orig*np.pi/180.0), axis=0) for ideff in range(ndeff_orig): Cext_orig[ideff,:] = Qext_orig[ideff,:] * np.pi * reff_orig[ideff] * reff_orig[ideff] # for ideff in range(ndeff_orig): # print("Reff = %.2f"%reff_orig[ideff]) # # Write opt_ file to be used for th betal expension # fopt = open(dir_opt+'opt_baum_'+kind+'_ice_reff_%.2f'%reff_orig[ideff]+'.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_orig+' \n') # fopt.write('# lambda(microns) nteta Cext (microns^2) SSA g \n') # for iwvl in range(nwvl_orig): # s = ' %.7e'%wvl_orig[iwvl]+' %i'%nang_orig+' %.7e'%Cext_orig[ideff, iwvl]+' %.7e'%ssa_orig[ideff, iwvl]+' %.7e'%g_orig[ideff, iwvl]+' \n' # fopt.write(s) # fopt.write('# Phase matrix \n') # for iwvl in range(nwvl_orig): # fopt.write('# lambda = %.4f'%wvl_orig[iwvl]+' \n') # fopt.write('# u F11 F22 F33 F44 F21 F34 \n') # for iang in range(nang_orig): # s = ' %17.9e'%(np.cos(ang_orig[iang]*np.pi/180.))+' %17.9e'%p11_orig[iang, ideff, iwvl]+' %17.9e'%p22_orig[iang, ideff, iwvl]+' %17.9e'%p33_orig[iang, ideff, iwvl]+' %17.9e'%p44_orig[iang, ideff, iwvl]+' %17.9e'%p21_orig[iang, ideff, iwvl]+' %17.9e'%(-p43_orig[iang, ideff, iwvl])+' \n' # fopt.write(s) # fopt.close() print('Compute Betal ...') F = np.zeros((6,len(ang)), order="F") F_recomp = np.zeros((6,len(ang)), order="F") coeftr = np.zeros((1), order="F") # F(1) = F11 # F(2) = F22 # F(3) = F33 # F(4) = F44 # F(5) = F12 # F(6) = F34 gauss_mu, gauss_wght = np.polynomial.legendre.leggauss(len(ang)) p11_gauss = interp1d( mu, p11_orig[:,:,:], axis=0)(gauss_mu) p22_gauss = interp1d( mu, p22_orig[:,:,:], axis=0)(gauss_mu) p33_gauss = interp1d( mu, p33_orig[:,:,:], axis=0)(gauss_mu) p44_gauss = interp1d( mu, p44_orig[:,:,:], axis=0)(gauss_mu) p21_gauss = interp1d( mu, p21_orig[:,:,:], axis=0)(gauss_mu) p43_gauss = interp1d( mu, p43_orig[:,:,:], axis=0)(gauss_mu) for ib,nb in enumerate(nbetal): print(nb) betal = np.zeros((6,nb+1), order="F") for ireff in range(len(reff_orig)): print( reff_orig[ireff] ) for iwvl in range(len(wvl_orig)): F[0,:]= p11_gauss[:,ireff,iwvl] F[1,:]= p22_gauss[:,ireff,iwvl] F[2,:]= p33_gauss[:,ireff,iwvl] F[3,:]= p44_gauss[:,ireff,iwvl] F[4,:]= p21_gauss[:,ireff,iwvl] F[5,:]=-p43_gauss[:,ireff,iwvl] trunc.mtrunc.trunc_dm(4, len(ang), gauss_mu, gauss_wght, F, nb, betal, coeftr) F_recomp[:,:] = 0.0 trunc.mtrunc.get_f_recomp(4, nb, len(ang), mu, betal, F_recomp) alpha1[:nb+1,ib,ireff,iwvl] = betal[0,:] alpha2[:nb+1,ib,ireff,iwvl] = betal[1,:] alpha3[:nb+1,ib,ireff,iwvl] = betal[2,:] alpha4[:nb+1,ib,ireff,iwvl] = betal[3,:] beta1[:nb+1,ib,ireff,iwvl] = betal[4,:] beta2[:nb+1,ib,ireff,iwvl] = betal[5,:] p11_recomp[:,ib,ireff,iwvl] = F_recomp[0,:] p22_recomp[:,ib,ireff,iwvl] = F_recomp[1,:] p33_recomp[:,ib,ireff,iwvl] = F_recomp[2,:] p44_recomp[:,ib,ireff,iwvl] = F_recomp[3,:] p21_recomp[:,ib,ireff,iwvl] = F_recomp[4,:] p34_recomp[:,ib,ireff,iwvl] = F_recomp[5,:] trunc_coeff[ib,ireff,iwvl] = coeftr[0] p11_recomp_integrate[:,:,:] = np.trapz( p11_recomp[:,:,:,:] , x=mu, axis=0) ############## # WRITE A RESULT FILE (HDF5) f = h5py.File(res_file, "a") if kind == "ghm": grp = f.require_group("ghm") elif kind == "sc": grp = f.require_group("sc") elif kind == "asc": grp = f.require_group("asc") subgrp = grp.require_group("axis") dset = subgrp.require_dataset("wavelengths", (nwvl_orig,), dtype='float64') dset[...] = wvl_orig dset = subgrp.require_dataset("reff", (nreff_orig,), dtype='float64') dset[...] = reff_orig dset = subgrp.require_dataset("phase_angles", (nang_orig,), dtype='float64') dset[...] = ang_orig dset = subgrp.require_dataset("mu", (nang_orig,), dtype='float64') dset[...] = np.cos(ang_orig*np.pi/180.0) dset = subgrp.require_dataset("nbetal", (n_nbetal,), dtype='i') dset[...] = nbetal dset = subgrp.require_dataset("ibetal", ( nbetalmax+1 ,), dtype='i') dset[...] = np.arange(nbetalmax+1) subgrp = grp.require_group("data") dset = subgrp.require_dataset("single_scattering_albedo", (nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("reff,wavelengths")) dset[...] = ssa_orig dset = subgrp.require_dataset("Cext", (nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("reff,wavelengths")) dset.attrs.create("unit", np.string_("micron^2")) dset[...] = Cext_orig dset = subgrp.require_dataset("extinction_coefficient_over_wc", (nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("reff,wavelengths")) dset.attrs.create("unit", np.string_("m^2g-1")) dset[...] = Sext_o_IWC_orig dset = subgrp.require_dataset("integrate_p11", (nreff_orig,nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("reff,wavelengths")) dset[...] = p11_integrate dset = subgrp.require_dataset("integrate_p11_recomp", (n_nbetal,nreff_orig,nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("nbetal,reff,wavelengths")) dset[...] = p11_recomp_integrate dset = subgrp.require_dataset("p11_phase_function", (nang_orig,nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,reff,wavelengths")) dset[...] = p11_orig dset = subgrp.require_dataset("p21_phase_function", (nang_orig,nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,reff,wavelengths")) dset[...] = p21_orig dset = subgrp.require_dataset("p34_phase_function", (nang_orig,nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,reff,wavelengths")) dset[...] = -p43_orig dset = subgrp.require_dataset("p22_phase_function", (nang_orig,nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,reff,wavelengths")) dset[...] = p22_orig dset = subgrp.require_dataset("p33_phase_function", (nang_orig,nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,reff,wavelengths")) dset[...] = p33_orig dset = subgrp.require_dataset("p44_phase_function", (nang_orig,nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,reff,wavelengths")) dset[...] = p44_orig dset = subgrp.require_dataset("truncation_coefficient", (n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = trunc_coeff dset = subgrp.require_dataset("alpha1_betal", (nbetalmax+1, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = alpha1 dset = subgrp.require_dataset("alpha2_betal", (nbetalmax+1, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = alpha2 dset = subgrp.require_dataset("alpha3_betal", (nbetalmax+1, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = alpha3 dset = subgrp.require_dataset("alpha4_betal", (nbetalmax+1, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = alpha4 dset = subgrp.require_dataset("beta1_betal", (nbetalmax+1, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = beta1 dset = subgrp.require_dataset("beta2_betal", (nbetalmax+1, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("ibetal,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = beta2 dset = subgrp.require_dataset("recomposed_p11_phase_function", (nang_orig, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = p11_recomp dset = subgrp.require_dataset("recomposed_p22_phase_function", (nang_orig, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = p22_recomp dset = subgrp.require_dataset("recomposed_p33_phase_function", (nang_orig, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = p33_recomp dset = subgrp.require_dataset("recomposed_p21_phase_function", (nang_orig, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = p21_recomp dset = subgrp.require_dataset("recomposed_p34_phase_function", (nang_orig, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = p34_recomp dset = subgrp.require_dataset("recomposed_p44_phase_function", (nang_orig, n_nbetal, nreff_orig, nwvl_orig), dtype='float32') dset.attrs.create("dimensions", np.string_("mu,nbetal,reff,wavelengths")) dset.attrs.create("truncation", np.string_("dm")) dset[...] = p44_recomp f.close if __name__=='__main__': nbetal = np.array([4, 8, 16, 32], dtype = int) wvl_bound = [0.1, 100.0] doit("asc", wvl_bound, nbetal, res_suffixe="") doit("ghm", wvl_bound, nbetal, res_suffixe="") doit("sc", wvl_bound, nbetal, res_suffixe="")