diff --git a/sice.py b/sice.py index 9cd8844..40ecf25 100644 --- a/sice.py +++ b/sice.py @@ -48,7 +48,7 @@ # area specific surface area (kg/m/m) # al effective absorption length(mm) # r0 reflectance of a semi-infinite non-absorbing snow layer - +# # plane BroadBand Albedo (BBA) # rp1 visible(0.3-0.7micron) # rp2 near-infrared (0.7-2.4micron) @@ -73,7 +73,7 @@ # quad_func calculation of quadratic parameters # trapzd trapezoidal rule for integral calculation # funp snow spectral planar and spherical albedo function -# ==================================== +# ==================================== import numpy as np from numpy import genfromtxt @@ -423,3 +423,4 @@ def func_solv(albedo): print(OutputFolder) print('---') pySICE(InputPath,OutputFolder) + diff --git a/sice_lib.py b/sice_lib.py index eb6a8a3..e60e5c2 100644 --- a/sice_lib.py +++ b/sice_lib.py @@ -83,6 +83,7 @@ from constants import w, bai, xa, ya, f0, f1, f2, bet, gam, coef1, coef2, coef3, coef4 import numpy as np # ================================================ + # tozon [i_channel] spectral ozone vertical optical depth at the fixed ozone concentration 404.59DU ( wavelength, VOD) # voda[i_channel] spectral water vapour vertical optical depth at the fixed concentration 3.847e+22 molecules per square sm @@ -92,31 +93,37 @@ # totadu ECMWF total column ozone in Dobson Unit # toa_cor_03 ozone-corrected OLCI toa relfectances -def molecular_absorption(ozone,tozon,sza,vza,toa): - scale = np.arccos(-1.)/180. # rad per degree + +def ozone_scattering(ozone, tozon, sza, vza, toa): + scale = np.arccos(-1.) / 180. # rad per degree eps = 1.55 # ecmwf ozone from OLCI file (in Kg.m-2) to DOBSON UNITS # 1 kg O3 / m2 = 46696.24 DOBSON Unit (DU) - totadu = 46696.24 * ozone + totadu = 46729. * ozone - inv_cos_za = 1./np.cos(sza*scale)+1./np.cos(vza*scale) + amf = 1. / np.cos(sza * scale) + 1. / np.cos(vza * scale) - BX=(toa[20]**(1.-eps)) * (toa[16]**eps) / toa[6] - BXXX=np.log(BX)/1.11e-4/inv_cos_za - BXXX[BXXX>500] = 999 - BXXX[BXXX<0] = 999 - - # Correcting TOA reflectance for ozone absorption - # bav 09-02-2020: now water scattering not accounted for - # kg/m**2. transfer to mol/cm**2 - # roznov = 2.99236e-22 # 1 moles Ozone = 47.9982 grams - # water vapor optical depth - # vap = water/roznov - # AKOWAT = vap/3.847e+22# tvoda = np.exp(inv_cos_za*voda*AKOWAT) - tvoda=tozon*0+1 - toa_cor_o3=toa*np.nan; + BX = (toa[20,:,:]**(1. - eps)) * (toa[16,:,:]**eps) / toa[6,:,:] + BXXX = np.log(BX) / 1.11e-4 / amf + BXXX[BXXX > 500] = 999 + BXXX[BXXX < 0] = 999 + + # Correcting TOA reflectance for ozone and water scattering + + # bav 09-02-2020: now water scattering not accounted for + # kg/m**2. transfer to mol/cm**2 + # roznov = 2.99236e-22 # 1 moles Ozone = 47.9982 grams + # water vapor optical depth + # vap = water/roznov + # AKOWAT = vap/3.847e+22# tvoda = np.exp(amf*voda*AKOWAT) + + tvoda = tozon * 0 + 1 + toa_cor_o3 = toa * np.nan + for i in range(21): - toa_cor_o3[i,:,:] = toa[i,:,:]*tvoda[i]*np.exp(inv_cos_za*tozon[i]*totadu/404.59) + + toa_cor_o3[i, :, :] = toa[i, :, :] * tvoda[i] \ + * np.exp(amf * tozon[i] * totadu / 404.59) return BXXX, toa_cor_o3 # viewing characteristics and aerosol properties @@ -128,29 +135,23 @@ def molecular_absorption(ozone,tozon,sza,vza,toa): # aot threshold value on aerosol optical thickness (aot) at 500nm # height height of underlying surface(meters) + def view_geometry(vaa, saa, sza, vza, aot, height): # transfer of OLCI relative azimuthal angle to the definition used in # radiative transfer code - # raa relative azimuth angle - # sza solar zenith angle - # vza viewing zenith angle - # cos_sa cosine of the scattering angle - # ak1 - # ak2 - - raa=180.-(vaa-saa) - sin_sza=np.sin(sza*np.pi/180.) - sin_vza=np.sin(vza*np.pi/180.) + raa = 180. - (vaa - saa) + as1 = np.sin(sza * np.pi / 180.) + as2 = np.sin(vza * np.pi / 180.) - cos_sza=np.cos(sza*np.pi/180.) - cos_vza=np.cos(vza*np.pi/180.) + am1 = np.cos(sza * np.pi / 180.) + am2 = np.cos(vza * np.pi / 180.) - ak1=3.*(1.+2.*cos_sza)/7. - ak2=3.*(1.+2.*cos_vza)/7. + ak1 = 3. * (1. + 2. * am1) / 7. + ak2 = 3. * (1. + 2. * am2) / 7. - cos_raa =np.cos(raa*np.pi/180.) - inv_cos_za=1./cos_sza+1./cos_vza - cos_sa=-cos_sza*cos_vza + sin_sza*sin_vza*cos_raa + cofi = np.cos(raa * np.pi / 180.) + amf = 1. / am1 + 1. / am2 + co = -am1 * am2 + as1 * as2 * cofi return raa, cos_sza, cos_vza, ak1, ak2, inv_cos_za, cos_sa # @@ -179,20 +180,22 @@ def aerosol_properties(aot, height, cos_sa): # molecular optical thickness taumol[i,:,:] = ak*0.00877*w[i]**(-4.05) tau[i,:,:] = tauaer[i] + taumol[i,:,:] + # aerosol asymmetry parameter - g[i,:,:]=tauaer[i]*gaer[i]/tau[i,:,:] + g[i, :, :] = tauaer[i] * gaer[i] / tau[i, :, :] # HG phase function for aerosol - pa[i,:,:]=(1-g[i,:,:]**2)/(1.-2.*g[i,:,:]*cos_sa+g[i,:,:]**2)**1.5 + pa[i, :, :] = (1 - g[i, :, :] ** 2) \ + / (1. - 2. * g[i, :, :] * co + g[i, :, :] ** 2) ** 1.5 - p[i,:,:]=(taumol[i,:,:]*pmol + tauaer[i]*pa[i,:,:])/tau[i,:,:] - - return tau, p, g, gaer,taumol,tauaer + p[i,:,:]=(taumol[i,:,:]*pmol + tauaer[i]*pa[i,:,:])/tau[i,:,:] + return tau, p, g, gaer, taumol, tauaer + +# %% snow properties -# snow properties def snow_properties(toa, ak1, ak2): - # retrieval of snow properties ( R_0, size of grains from OLCI channels 865[17] and 1020nm[21] + # retrieval of snow properties ( R_0, size of grains from OLCI channels 865[17] and 1020nm[21] # assumed not influenced by atmospheric scattering and absorption processes) # imaginary part of the ice refractive index at 1020nm @@ -209,9 +212,9 @@ def snow_properties(toa, ak1, ak2): # consequently: 1-eps = 1/(1-np.sqrt(alpha2/alpha1)) # reflectivity of nonabsorbing snow layer - rr1=toa[16,:,:] - rr2=toa[20,:,:] - r0 = (rr1**eps)*(rr2**(1.-eps)) + rr1 = toa[16, :, :] + rr2 = toa[20, :, :] + r0 = (rr1 ** eps) * (rr2 ** (1. - eps)) # effective absorption length(mm) bal = (np.log(rr2/r0)/(ak1*ak2/r0))**2 /alpha2 @@ -220,9 +223,14 @@ def snow_properties(toa, ak1, ak2): # effective grain size(mm):diameter # xi/(1-g) = 9.2 D=al/(9.2*16/9) + # snow specific area ( dimension: m*m/kg) - area= 6./D/0.917 - return D, area, al, r0, bal + area = 6. / D / 0.917 + + return D, area, al, r0, bal + +# %% ================================================= + # ================================================= def prepare_coef(tau, g, p, cos_sza, cos_vza, inv_cos_za, gaer, taumol, tauaer): @@ -232,16 +240,16 @@ def prepare_coef(tau, g, p, cos_sza, cos_vza, inv_cos_za, gaer, taumol, tauaer): t2=tau*np.nan # SOBOLEV - oskar=4.+3.*(1.-g)*tau - b1=1.+1.5*cos_sza+(1.-1.5*cos_sza)*np.exp(-tau/cos_sza) - b2=1.+1.5*cos_vza+(1.-1.5*cos_vza)*np.exp(-tau/cos_vza) + oskar = 4. + 3. * (1. - g) * tau + b1 = 1. + 1.5 * am1 + (1. - 1.5 * am1) * np.exp(-tau / am1) + b2 = 1. + 1.5 * am2 + (1. - 1.5 * am2) * np.exp(-tau / am2) - wa1=1.10363 - wa2=-6.70122 - wx0=2.19777 - wdx=0.51656 - bex=np.exp ( (g-wx0)/wdx ) - sssss= (wa1-wa2)/(1.+bex)+wa2 + wa1 = 1.10363 + wa2 = -6.70122 + wx0 = 2.19777 + wdx = 0.51656 + bex = np.exp((g - wx0) / wdx) + sssss = (wa1 - wa2) / (1. + bex) + wa2 for i in range(21): M[i,:,:]=(1.-np.exp(-tau[i,:,:]*inv_cos_za))/(cos_sza+cos_vza)/4. @@ -284,59 +292,62 @@ def atm_spher_alb(tau, g): # snow_imputirities def snow_impurities(alb_sph, bal): - # analysis of snow impurities + # analysis of snow impurities # ( the concentrations below 0.0001 are not reliable ) # bf normalized absorption coefficient of pollutants ay 1000nm ( in inverse mm) # bm Angstroem absorption coefficient of pollutants ( around 1 - for soot, 3-7 for dust) - bm=np.nan*bal - bf=bm + bm = np.nan * bal + bf = bm p1 = bm p2 = bm - ind_nonan = np.logical_and(np.logical_not(np.isnan(alb_sph[0,:,:])), - np.logical_not(np.isnan(alb_sph[1,:,:]))) - p1[ind_nonan]=np.log(alb_sph[0,ind_nonan])*np.log(alb_sph[0,ind_nonan]) - p2[ind_nonan]=np.log(alb_sph[1,ind_nonan])*np.log(alb_sph[1,ind_nonan]) - bm[ind_nonan]=np.log( p1[ind_nonan]/p2[ind_nonan])/np.log(w[1]/w[0]) + ind_nonan = np.logical_and(np.logical_not(np.isnan(alb_sph[0, :, :])), + np.logical_not(np.isnan(alb_sph[1, :, :]))) + + p1[ind_nonan] = np.log(alb_sph[0, ind_nonan]) * np.log(alb_sph[0, ind_nonan]) + p2[ind_nonan] = np.log(alb_sph[1, ind_nonan]) * np.log(alb_sph[1, ind_nonan]) + bm[ind_nonan] = np.log(p1[ind_nonan] / p2[ind_nonan]) / np.log(w[1] / w[0]) # type of pollutants - ntype=np.nan*bal - ntype[bm <= 1.2]=1 # soot - ntype[bm > 1.2]=2 # dust + ntype = np.nan * bal + ntype[bm <= 1.2] = 1 # soot + ntype[bm > 1.2] = 2 # dust - soda = bm*np.nan - soda[bm>=0.1]=(w[0])**bm[bm>=0.1] - bf=soda*p1/bal + soda = bm * np.nan + soda[bm >= 0.1] = (w[0]) ** bm[bm >= 0.1] + bf = soda * p1 / bal # normalized absorption coefficient of pollutants at the wavelength 1000nm - k_abs_0=p1/bal + bff = p1 / bal # bal -effective absorption length in microns - B_soot=1.6 # enhancement factors for soot - B_ice= 0.9 # enhancement factors for ice grains - alfa_soot=4.*np.pi*0.47/w[0] # bulk soot absorption coefficient at 1000nm - k_dust=0.01 # volumetric absorption coefficient of dust - - conc = bal*np.nan - conc[ntype == 1] = B_soot*k_abs_0[ntype == 1]/B_ice/alfa_soot - conc[ntype == 2] = B_soot*k_abs_0[ntype == 2]/k_dust - ntype[bm <= 0.5] = 3 # type is other or mixture - ntype[bm >= 10.] = 4 # type is other or mixture + BBBB = 1.6 # enhancement factors for soot + FFFF = 0.9 # enhancement factors for ice grains + alfa = 4. * np.pi * 0.47 / w[0] # bulk soot absorption coefficient at 1000nm + DUST = 0.01 # volumetric absorption coefficient of dust + + conc = bal * np.nan + conc[ntype == 1] = BBBB * bff[ntype == 1] / FFFF / alfa + conc[ntype == 2] = BBBB * bff[ntype == 2] / DUST + ntype[bm <= 0.5] = 3 # type is other or mixture + ntype[bm >= 10.] = 4 # type is other or mixture + return ntype, bf, conc # =========================================================================== def alb2rtoa(a, t1, t2, r0, ak1, ak2, ratm, r): -# Function that calculates the theoretical reflectance from a snow spherical albedo a -# This function can then be solved to find optimal snow albedo -# Inputs: -# a Surface albedo -# r0 reflectance of a semi-infinite non-absorbing snow layer -# -# Outputs: -# rs surface reflectance at specific channel - surf = t1*t2*r0*a**(ak1*ak2/r0)/(1-a*ratm) - rs=r + surf + # Function that calculates the theoretical reflectance from a snow spherical albedo a + # This function can then be solved to find optimal snow albedo + # Inputs: + # a Surface albedo + # r0 reflectance of a semi-infinite non-absorbing snow layer + # + # Outputs: + # rs surface reflectance at specific channel + surf = t1 * t2 * r0 * a ** (ak1 * ak2 / r0) / (1 - a * ratm) + rs = r + surf + return rs # =========================================================================== @@ -349,21 +360,24 @@ def salbed(tau, g): # g asymetry coefficient # outputs: # salbed spherical albedo - a_s = (.18016, -0.18229, 0.15535, -0.14223) - bs = (.58331, -0.50662, -0.09012, 0.0207) - cs = (0.21475, -0.1, 0.13639, -0.21948) - als = (0.16775, -0.06969, 0.08093, -0.08903) - bets = (1.09188, 0.08994, 0.49647, -0.75218) - - a = a_s[0]*g**0 + a_s[1]*g**1 + a_s[2]*g**2 + a_s[3]*g**3 - b = bs[0]*g**0 + bs[1]*g**1 + bs[2]*g**2 + bs[3]*g**3 - c = cs[0]*g**0 + cs[1]*g**1 + cs[2]*g**2 + cs[3]*g**3 - al= als[0]*g**0 + als[1]*g**1 + als[2]*g**2 + als[3]*g**3 - bet= bets[0]*g**0 + bets[1]*g**1 + bets[2]*g**2 + bets[3]*g**3 + a_s = (.18016, -0.18229, 0.15535, -0.14223) + bs = (.58331, -0.50662, -0.09012, 0.0207) + cs = (0.21475, -0.1, 0.13639, -0.21948) + als = (0.16775, -0.06969, 0.08093, -0.08903) + bets = (1.09188, 0.08994, 0.49647, -0.75218) + + a = a_s[0] * g ** 0 + a_s[1] * g ** 1 + a_s[2] * g ** 2 + a_s[3] * g ** 3 + b = bs[0] * g ** 0 + bs[1] * g ** 1 + bs[2] * g ** 2 + bs[3] * g ** 3 + c = cs[0] * g ** 0 + cs[1] * g ** 1 + cs[2] * g ** 2 + cs[3] * g ** 3 + al = als[0] * g ** 0 + als[1] * g ** 1 + als[2] * g ** 2 + als[3] * g ** 3 + bet = bets[0] * g ** 0 + bets[1] * g ** 1 + bets[2] * g ** 2 + bets[3] * g ** 3 - salbed = tau*(a*np.exp(-tau/al)+b*np.exp(-tau/bet)+c) + salbed = tau * (a * np.exp(-tau / al) + b * np.exp(-tau / bet) + c) + return salbed +# %% ===================================================================== + # ===================================================================== def zbrent(f, x0, x1, max_iter=100, tolerance=1e-6): @@ -371,7 +385,7 @@ def zbrent(f, x0, x1, max_iter=100, tolerance=1e-6): # https://en.wikipedia.org/wiki/Brent%27s_method # Brent’s is essentially the Bisection method augmented with Inverse # Quadratic Interpolation whenever such a step is safe. At it’s worst case - #it converges linearly and equal to Bisection, but in general it performs + # it converges linearly and equal to Bisection, but in general it performs # superlinearly; it combines the robustness of Bisection with the speedy # convergence and inexpensive computation of Quasi-Newtonian methods. # Because of this, you’re likely to find Brent’s as a default root-finding @@ -382,10 +396,10 @@ def zbrent(f, x0, x1, max_iter=100, tolerance=1e-6): fx0 = f(x0) fx1 = f(x1) -# print(str(fx0) + ", " + str(fx1)) + # print(str(fx0) + ", " + str(fx1)) if ((fx0 * fx1) > 0): -# print("Root not bracketed "+str(fx0)+", "+str(fx1)) -# assert ((fx0 * fx1) <= 0), ("-----Root not bracketed"+str(fx0)+", "+str(fx1)) + # print("Root not bracketed "+str(fx0)+", "+str(fx1)) + # assert ((fx0 * fx1) <= 0), ("-----Root not bracketed"+str(fx0)+", "+str(fx1)) return -999 if abs(fx0) < abs(fx1): @@ -397,7 +411,7 @@ def zbrent(f, x0, x1, max_iter=100, tolerance=1e-6): mflag = True steps_taken = 0 - while steps_taken < max_iter and abs(x1-x0) > tolerance: + while steps_taken < max_iter and abs(x1 - x0) > tolerance: fx0 = f(x0) fx1 = f(x1) fx2 = f(x2) @@ -409,13 +423,13 @@ def zbrent(f, x0, x1, max_iter=100, tolerance=1e-6): new = L0 + L1 + L2 else: - new = x1 - ( (fx1 * (x1 - x0)) / (fx1 - fx0) ) + new = x1 - ((fx1 * (x1 - x0)) / (fx1 - fx0)) - if ((new < ((3 * x0 + x1) / 4) or new > x1) or - (mflag == True and (abs(new - x1)) >= (abs(x1 - x2) / 2)) or - (mflag == False and (abs(new - x1)) >= (abs(x2 - d) / 2)) or - (mflag == True and (abs(x1 - x2)) < tolerance) or - (mflag == False and (abs(x2 - d)) < tolerance)): + if ((new < ((3 * x0 + x1) / 4) or new > x1) + or (mflag and (abs(new - x1)) >= (abs(x1 - x2) / 2)) + or (mflag == False and (abs(new - x1)) >= (abs(x2 - d) / 2)) + or (mflag and (abs(x1 - x2)) < tolerance) + or (mflag == False and (abs(x2 - d)) < tolerance)): new = (x0 + x1) / 2 mflag = True @@ -455,22 +469,38 @@ def funp(x, al, sph_calc, ak1): # bav 2020 # using numpy interpolation - y = np.interp(x,xa,ya) + y = np.interp(x, xa, ya) - dega = 1000.* al * 4.*np.pi*y/x + dega = 1000. * al * 4. * np.pi * y / x pow = np.sqrt(dega) - if (pow >= 1.e-6): rsd = np.exp(-pow) - else: rsd=1. + if (pow >= 1.e-6): + rsd = np.exp(-pow) + else: + rsd = 1. - if (sph_calc == 0): rs = rsd**ak1 - elif (sph_calc == 1): rs = rsd + if (sph_calc == 0): + rs = rsd**ak1 + elif (sph_calc == 1): + rs = rsd if (x < 0.4): x = 0.4 - funcs = f0+ f1*np.exp(-x*bet)+ f2*np.exp(-x*gam) + funcs = f0 + f1 * np.exp(-x * bet) + f2 * np.exp(-x * gam) - return rs*funcs + return rs * funcs + +# %% Approximation functions for BBA integration + + +def plane_albedo_sw_approx(D, am1): + anka = 0.7389 - 0.1783 * am1 + 0.0484 * am1 ** 2. + banka = 0.0853 + 0.0414 * am1 - 0.0127 * am1 ** 2. + canka = 0.1384 + 0.0762 * am1 - 0.0268 * am1 ** 2. + diam1 = 187.89 - 69.2636 * am1 + 40.4821 * am1 ** 2. + diam2 = 2687.25 - 405.09 * am1 + 94.5 * am1 ** 2. + return anka + banka * np.exp(-1000 * D / diam1) + canka \ + * np.exp(-1000 * D / diam2) # Approximation functions for BBA integration def plane_albedo_sw_approx(D,cos_sza): @@ -493,67 +523,76 @@ def spher_albedo_sw_approx(D): def BBA_calc_clean(al, ak1): # for clean snow # plane albedo - sph_calc = 0 # planar + sph_calc = 0 # planar # visible(0.3-0.7micron) + def func_integ(x): return funp(x, al, sph_calc, ak1) - p1 = qsimp(func_integ,0.3,0.7) + p1 = qsimp(func_integ, 0.3, 0.7) # near-infrared (0.7-2.4micron) -# p2 = trapzd(func_integ,0.7,2.4, 20) - p2 = qsimp(func_integ,0.7,2.4) + # p2 = trapzd(func_integ,0.7,2.4, 20) + p2 = qsimp(func_integ, 0.7, 2.4) # spherical albedo - sph_calc = 1 # spherical calculation + sph_calc = 1 # spherical calculation def func_integ(x): return funp(x, al, sph_calc, ak1) # visible(0.3-0.7micron) -# s1 = trapzd(func_integ,0.3,0.7, 20) - s1 = qsimp(func_integ,0.3,0.7) + # s1 = trapzd(func_integ,0.3,0.7, 20) + s1 = qsimp(func_integ, 0.3, 0.7) # near-infrared (0.7-2.4micron) -# s2 = trapzd(func_integ,0.7,2.4, 20) - s2 = qsimp(func_integ,0.7,2.4) + # s2 = trapzd(func_integ,0.7,2.4, 20) + s2 = qsimp(func_integ, 0.7, 2.4) # shortwave(0.3-2.4 micron) # END of clean snow bba calculation - return p1,p2,s1,s2 + + return p1, p2, s1, s2 + +# %% =============================== + # =============================== def qsimp(func,a,b): # integrate function between a and b using simpson's method. # works as fast as scipy.integrate quad - eps=1.e-3 - jmax=20 - ost=-1.e30 - os= -1.e30 + eps = 1.e-3 + jmax = 20 + ost = -1.e30 + os = -1.e30 + for j in range(jmax): - if (j==0): - st=0.5*(b-a)*(func(a)+func(b)) + + if (j == 0): + st = 0.5 * (b - a) * (func(a) + func(b)) else: - it=2**(j-1) - tnm=it - delta=(b-a)/tnm - x=a+0.5*delta - sum=0. + it = 2 ** (j - 1) + tnm = it + delta = (b - a) / tnm + x = a + 0.5 * delta + sum = 0. + for jj in range(it): - sum=sum+func(x) - x=x+delta - st=0.5*(st+(b-a)*sum/tnm) - s=(4.*st-ost)/3. - if (j>4): - if (abs(s-os) 4): + if (abs(s - os) < eps * abs(os)): return s - if (s==0) and (os==0.): + if (s == 0) and (os == 0.): return s - os=s - ost=st + os = s + ost = st print("Max iteration reached") + return s # Calculation f BBA for polluted snow - def BBA_calc_pol(alb, asol, sol1_pol, sol2, sol3_pol): # polluted snow # NEW CODE FOR BBA OF BARE ICE @@ -565,42 +604,42 @@ def BBA_calc_pol(alb, asol, sol1_pol, sol2, sol3_pol): # segment 1 # QUADRATIC POLYNOMIal for the range 400-709nm # input wavelength -# alam2=w[0] -# alam3=w[5] -# alam5=w[10] -# alam6=w[11] -# alam7=w[16] -# alam8=w[20] - - alam2=0.4 - alam3=0.56 - alam5=0.709 - alam6=0.753 - alam7=0.865 - alam8=1.02 + # alam2=w[0] + # alam3=w[5] + # alam5=w[10] + # alam6=w[11] + # alam7=w[16] + # alam8=w[20] + + alam2 = 0.4 + alam3 = 0.56 + alam5 = 0.709 + alam6 = 0.753 + alam7 = 0.865 + alam8 = 1.02 - #input reflectances - r2 =alb[0,:] - r3 =alb[5,:] - r5 =alb[10,:] - r6=alb[11,:] - r7=alb[16,:] - r8=alb[20,:] - - sa1, a1, b1, c1 = quad_func(alam2,alam3,alam5, r2 ,r3,r5) - ajx1 = a1*sol1_pol - ajx2 = b1*coef1 - ajx3 = c1*coef2 + # input reflectances + r2 = alb[0, :] + r3 = alb[5, :] + r5 = alb[10, :] + r6 = alb[11, :] + r7 = alb[16, :] + r8 = alb[20, :] + + sa1, a1, b1, c1 = quad_func(alam2, alam3, alam5, r2, r3, r5) + ajx1 = a1 * sol1_pol + ajx2 = b1 * coef1 + ajx3 = c1 * coef2 aj1 = ajx1 + ajx2 + ajx3 # segment 2.1 # QUADRATIC POLYNOMIal for the range 709-865nm - sa1, a2, b2, c2 = quad_func(alam5,alam6,alam7,r5,r6,r7) - ajx1 = a2*asol - ajx2 = b2*coef3 - ajx3 = c2*coef4 + sa1, a2, b2, c2 = quad_func(alam5, alam6, alam7, r5, r6, r7) + ajx1 = a2 * asol + ajx2 = b2 * coef3 + ajx3 = c2 * coef4 - aj2 = ajx1 + ajx2 + ajx3 # segment 2.2 + aj2 = ajx1 + ajx2 + ajx3 # segment 2.2 # exponential approximation for the range 865- 2400 nm z1=0.865 z2=2.4 @@ -625,13 +664,14 @@ def quad_func(x0,x1,x2,y0,y1,y2): # quadratic function used for the polluted snow BBA calculation # see BBA_calc_pol # compatible with arrays - d1=(x0-x1)*(x0-x2) - d2=(x1-x0)*(x1-x2) - d3=(x2-x0)*(x2-x1) - - a1=x1*x2*y0/d1+x0*x2*y1/d2+x0*x1*y2/d3 - b1=-(x1+x2)*y0/d1-(x0+x2)*y1/d2-(x0+x1)*y2/d3 - c1=y0/d1+y1/d2+y2/d3 - x=x1 - sa=a1+b1*x+c1*x*x + d1 = (x0 - x1) * (x0 - x2) + d2 = (x1 - x0) * (x1 - x2) + d3 = (x2 - x0) * (x2 - x1) + + a1 = x1 * x2 * y0 / d1 + x0 * x2 * y1 / d2 + x0 * x1 * y2 / d3 + b1 = -(x1 + x2) * y0 / d1 - (x0 + x2) * y1 / d2 - (x0 + x1) * y2 / d3 + c1 = y0 / d1 + y1 / d2 + y2 / d3 + x = x1 + sa = a1 + b1 * x + c1 * x * x + return sa, a1, b1, c1