diff --git a/sice.py b/sice.py
index 8a570d4..46fcc87 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
@@ -81,313 +81,304 @@
 import rasterio as rio
 import time
 import sys
-import os
-import pandas as pd
 from constants import w, bai, sol1_clean, sol2, sol3_clean, sol1_pol, sol3_pol, asol
 np.seterr(invalid='ignore')
 
 start_time = time.process_time()
-print(sys.argv[1])
-
-#%% input text file
-if os.path.isfile(sys.argv[1]):  
-    InputFolder = os.path.dirname(os.path.dirname(sys.argv[1]))+'/'
-    print('\nText file input')
-    # data_in = pd.read_csv(sys.argv[1])
-    data_in = pd.read_csv(sys.argv[1])
-    toa = np.expand_dims(data_in[[c for c in data_in.columns if c.find('reflec')>=0]].to_numpy().transpose(), axis=2)
-    
-        
-    ozone = np.expand_dims(data_in['total_ozone'], axis=1)
-    water = np.expand_dims(data_in['total_columnar_water_vapour'], axis=1)
-    sza = np.expand_dims(data_in['sza'], axis=1)
-    saa = np.expand_dims(data_in['saa'], axis=1)
-    vza = np.expand_dims(data_in['vza'], axis=1)
-    vaa = np.expand_dims(data_in['vaa'], axis=1)
-    height = np.expand_dims(data_in['altitude'], axis=1)
-    
-    sza[np.isnan(toa[0,:,:])] = np.nan
-    saa[np.isnan(toa[0,:,:])] = np.nan
-    vza[np.isnan(toa[0,:,:])] = np.nan
-    vaa[np.isnan(toa[0,:,:])] = np.nan   
-
-#%% ========= input tif ===============
-elif os.path.isdir(sys.argv[1]):
-    InputFolder =  sys.argv[1] + '/'
-    print("\nTiff files input")  
-    Oa01 = rio.open(InputFolder+'r_TOA_01.tif')
-    meta = Oa01.meta
-    with rio.Env():    
-        meta.update(compress='DEFLATE')
-    
-    def WriteOutput(var,var_name,in_folder):
-        # this functions write tif files based on a model file, here "Oa01"
-        # opens a file for writing
-    
-        with rio.open(in_folder+var_name+'.tif', 'w+', **meta) as dst:
-            dst.write(var.astype('float32'),1)
-        
-    toa = np.tile(Oa01.read(1)*np.nan, (21,1,1))
-    
-    for i in range(21):
-        dat = rio.open((InputFolder+'r_TOA_'+str(i+1).zfill(2)+'.tif'))
-        toa[i,:,:] = dat.read(1)
-        
-    ozone = rio.open(InputFolder+'O3.tif').read(1)
-    water = rio.open(InputFolder+'WV.tif').read(1)
-    sza = rio.open(InputFolder+'SZA.tif').read(1)
-    saa = rio.open(InputFolder+'SAA.tif').read(1)
-    vza = rio.open(InputFolder+'OZA.tif').read(1)
-    vaa = rio.open(InputFolder+'OAA.tif').read(1)
-    height = rio.open(InputFolder+'height.tif').read(1)
+
+InputFolder = sys.argv[1] + '/'
+
+# %% ========= input tif ================
+
+Oa01 = rio.open(InputFolder + 'r_TOA_01.tif')
+meta = Oa01.meta
+
+with rio.Env():    
+    meta.update(compress='DEFLATE')
+
+
+def WriteOutput(var, var_name, in_folder):
+    # this functions write tif files based on a model file, here "Oa01"
+    # opens a file for writing
+
+    with rio.open(in_folder + var_name + '.tif', 'w+', **meta) as dst:
+        dst.write(var.astype('float32'), 1)
+  
+toa = np.tile(Oa01.read(1).astype('float32') * np.nan, (21, 1, 1))
+
+for i in range(21):
     
-    sza[np.isnan(toa[0,:,:])] = np.nan
-    saa[np.isnan(toa[0,:,:])] = np.nan
-    vza[np.isnan(toa[0,:,:])] = np.nan
-    vaa[np.isnan(toa[0,:,:])] = np.nan
+    try:
+        dat = rio.open((InputFolder + 'r_TOA_' + str(i + 1).zfill(2) + '.tif'))
+        toa[i, :, :] = dat.read(1).astype('float32')
+    except:
+        toa[i, :, :] = np.nan
     
-else:  
-    print("\n Input path neither file or directory" )
+ozone = rio.open(InputFolder + 'O3.tif').read(1).astype('float32')
+water = rio.open(InputFolder + 'WV.tif').read(1).astype('float32')
+sza = rio.open(InputFolder + 'SZA.tif').read(1).astype('float32')
+saa = rio.open(InputFolder + 'SAA.tif').read(1).astype('float32')
+vza = rio.open(InputFolder + 'OZA.tif').read(1).astype('float32')
+vaa = rio.open(InputFolder + 'OAA.tif').read(1).astype('float32')
+height = rio.open(InputFolder + 'height.tif').read(1).astype('float32')
+
+sza[np.isnan(toa[0, :, :])] = np.nan
+saa[np.isnan(toa[0, :, :])] = np.nan
+vza[np.isnan(toa[0, :, :])] = np.nan
+vaa[np.isnan(toa[0, :, :])] = np.nan
 
-# %% water and ozone spectral optical density    
 water_vod = genfromtxt('./tg_water_vod.dat', delimiter='   ')
-voda = water_vod[range(21),1]
+voda = water_vod[range(21), 1]
 
 ozone_vod = genfromtxt('./tg_vod.dat', delimiter='   ')
-tozon = ozone_vod[range(21),1]
+tozon = ozone_vod[range(21), 1]
 aot = 0.1
 
-#%%   declaring variables
-BXXX, isnow, D, area, al, r0, isnow, conc, ntype, rp1, rp2, rp3, rs1, rs2, rs3 =   vaa*np.nan, vaa*np.nan, vaa*np.nan, vaa*np.nan, vaa*np.nan, vaa*np.nan,  vaa*np.nan, vaa*np.nan, vaa*np.nan, vaa*np.nan, vaa*np.nan, vaa*np.nan, vaa*np.nan, vaa*np.nan, vaa*np.nan
+# %%   declaring variables
 
-alb_sph, rp, refl =  toa*np.nan, toa*np.nan, toa*np.nan
+BXXX, isnow, D, area, al, r0, isnow, conc, ntype, rp1, rp2, rp3, rs1, rs2, rs3 =  \
+    vaa * np.nan, vaa * np.nan, vaa * np.nan, vaa * np.nan, vaa * np.nan, vaa * np.nan, \
+    vaa * np.nan, vaa * np.nan, vaa * np.nan, vaa * np.nan, vaa * np.nan, vaa * np.nan, \
+    vaa * np.nan, vaa * np.nan, vaa * np.nan
 
-#%% =========== ozone scattering  ====================================
-BXXX, toa_cor_o3 = sl.molecular_absorption(ozone,tozon, sza, vza,toa)
+alb_sph, rp, refl = toa * np.nan, toa * np.nan, toa * np.nan
+
+# %% =========== ozone scattering  ====================================
+
+BXXX, toa_cor_o3 = sl.ozone_scattering(ozone, tozon, sza, vza, toa)
 
 # Filtering pixels unsuitable for retrieval
-isnow[sza>75] = 100
-isnow[toa_cor_o3[20, :,:] < 0.1] = 102
+isnow[sza > 75] = 100
+isnow[toa_cor_o3[20, :, :] < 0.1] = 102
+
 for i_channel in range(21):
     toa_cor_o3[i_channel, ~np.isnan(isnow)] = np.nan
 
-vaa[ ~np.isnan(isnow)] = np.nan
-saa[ ~np.isnan(isnow)] = np.nan
-sza[ ~np.isnan(isnow)] = np.nan
-vza[ ~np.isnan(isnow)] = np.nan
-height = height.astype(float)
-height[ ~np.isnan(isnow)] = np.nan
+vaa[~np.isnan(isnow)] = np.nan
+saa[~np.isnan(isnow)] = np.nan
+sza[~np.isnan(isnow)] = np.nan
+vza[~np.isnan(isnow)] = np.nan
+
+height[~np.isnan(isnow)] = np.nan
 
 # =========== view geometry and atmosphere propeties  ==============
-raa, cos_sza, cos_vza, ak1, ak2, inv_cos_za, cos_sa = sl.view_geometry(vaa, saa, sza, vza, aot, height)
 
-tau, p, g, gaer,taumol,tauaer = sl.aerosol_properties(aot, height, cos_sa)
-        
+raa, am1, am2, ak1, ak2, amf, co = sl.view_geometry(vaa, saa, sza, vza, aot, height)
+
+tau, p, g, gaer, taumol, tauaer = sl.aerosol_properties(aot, height, co)
+     
 # =========== snow properties  ====================================
+
 D, area, al, r0, bal = sl.snow_properties(toa_cor_o3, ak1, ak2)
 # filtering small D
 D_thresh = 0.1
-isnow[np.logical_and(D<D_thresh, np.isnan(isnow))] = 104
+isnow[D < D_thresh] = 104
 
 for i in range(21):
-    toa_cor_o3[i,D<D_thresh] = np.nan
-area[D<D_thresh] = np.nan
-al[D<D_thresh] = np.nan
-r0[D<D_thresh] = np.nan
-bal[D<D_thresh] = np.nan
-cos_sza[D<D_thresh] = np.nan
-cos_vza[D<D_thresh] = np.nan
-#D[D<D_thresh] = np.nan
+    toa_cor_o3[i, D < D_thresh] = np.nan
+    
+area[D < D_thresh] = np.nan
+al[D < D_thresh] = np.nan
+r0[D < D_thresh] = np.nan
+bal[D < D_thresh] = np.nan
+am1[D < D_thresh] = np.nan
+am2[D < D_thresh] = np.nan
+# D[D<D_thresh] = np.nan
 
 # =========== clean snow  ====================================
+
 # for that we calculate the theoretical reflectance at band 1 of a surface with:
 # r0 = 1, a (albedo) = 1, ak1 = 1, ak2 = 1
 # t1 and t2 are the backscattering fraction
-t1, t2, ratm, r, astra, rms = sl.prepare_coef(tau, g, p, cos_sza, cos_vza, inv_cos_za,gaer,taumol,tauaer)
-rs_1 = sl.alb2rtoa(1, t1[0,:,:], t2[0,:,:], np.ones_like(r0), np.ones_like(ak1), 
-                   np.ones_like(ak2), ratm[0,:,:], r[0,:,:])
+t1, t2, ratm, r, astra, rms = sl.prepare_coef(tau, g, p, am1, am2, amf, gaer, 
+                                              taumol, tauaer)
+rs_1 = sl.alb2rtoa(1, t1[0, :, :], t2[0, :, :], np.ones_like(r0), np.ones_like(ak1), 
+                   np.ones_like(ak2), ratm[0, :, :], r[0, :, :])
 
 # we then compare it to the observed toa[0] value
-ind_clean = toa_cor_o3[0,:,:] >= rs_1
+ind_clean = toa_cor_o3[0, :, :] >= rs_1
 isnow[ind_clean] = 0
       
 # STEP 4a: clean snow retrieval
 # the spherical albedo derivation: alb_sph
-def mult_channel(c,A):
-    tmp = A.T*c
+
+
+def mult_channel(c, A):
+    tmp = A.T * c
     return tmp.T
-alb_sph = np.exp(-np.sqrt(1000.*4.*np.pi* mult_channel(bai/w, np.tile(al,(21,1,1)))))
-alb_sph[alb_sph>0.999]=1
 
-# =========== polluted snow  ====================================
-ind_pol = toa_cor_o3[0,:,:] < rs_1
+
+alb_sph = np.exp(-np.sqrt(1000. * 4. * np.pi
+                          * mult_channel(bai / w, np.tile(al, (21, 1, 1)))))
+alb_sph[alb_sph > 0.999] = 1
+
+# ========== very dirty snow  ====================================
+
+ind_pol = toa_cor_o3[0, :, :] < rs_1
 
 isnow[ind_pol] = 1
 
-#  very dirty snow  
-ind_very_dark = np.logical_and(toa_cor_o3[20,:,:]<0.4, ind_pol)
+ind_very_dark = np.logical_and(toa_cor_o3[20] < 0.4, ind_pol)
 isnow[ind_very_dark] = 6
 
-am11=np.sqrt(1.-cos_sza[ind_very_dark]**2.)
-am12=np.sqrt(1.-cos_vza[ind_very_dark]**2.)
+am11 = np.sqrt(1. - am1[ind_very_dark] ** 2.)
+am12 = np.sqrt(1. - am2[ind_very_dark] ** 2.)
 
-theta = np.arccos(-cos_sza[ind_very_dark] * cos_vza[ind_very_dark] + am11 * am12 * np.cos(raa[ind_very_dark]*3.14159/180.))  *180./np.pi
+tz = np.arccos(-am1[ind_very_dark] * am2[ind_very_dark] + am11 * am12 
+               * np.cos(raa[ind_very_dark] * 3.14159 / 180.)) * 180. / np.pi
              
-pz=11.1*np.exp(-0.087*theta)+1.1*np.exp(-0.014*theta)
+pz = 11.1 * np.exp(-0.087 * tz) + 1.1 * np.exp(-0.014 * tz)
 
-rclean = 1.247 + 1.186 *(cos_sza[ind_very_dark]+cos_vza[ind_very_dark]) +  5.157 * cos_sza[ind_very_dark] * cos_vza[ind_very_dark] + pz
-
-rclean = rclean /4. /(cos_sza[ind_very_dark] + cos_vza[ind_very_dark])
+rclean = 1.247 + 1.186 * (am1[ind_very_dark] + am2[ind_very_dark]) \
+    + 5.157 * am1[ind_very_dark] * am2[ind_very_dark] + pz
 
+rclean = rclean / 4. / (am1[ind_very_dark] + am2[ind_very_dark])
 r0[ind_very_dark] = rclean
-#  end very dirty snow  
 
-ind_pol =  np.logical_or(ind_very_dark, ind_pol)
+# =========== polluted snow  ====================================
+
+ind_pol = np.logical_or(ind_very_dark, ind_pol)
+
 if np.any(ind_pol):
     subs_pol = np.argwhere(ind_pol)
     
     # approximation of the transcendental equation allowing closed-from solution
-    #alb_sph[:,ind_pol] =   (toa_cor_o3[:,ind_pol] - r[:,ind_pol])/(t1[:,ind_pol]*t2[:,ind_pol]*r0[ind_pol] + ratm[:,ind_pol]*(toa_cor_o3[:,ind_pol] - r[:,ind_pol]))
+    # alb_sph[:, ind_pol] = (toa_cor_o3[:, ind_pol] - r[:, ind_pol]) \ 
+    # /(t1[:,ind_pol]*t2[:,ind_pol]*r0[ind_pol] + ratm[:,ind_pol]*(toa_cor_o3[:,ind_pol] - r[:,ind_pol]))
     
     # solving iteratively the transcendental equation
-    alb_sph[:,ind_pol] = 1
+    alb_sph[:, ind_pol] = 1
     
-    def solver_wrapper(toa_cor_o3,tau, t1, t2, r0, ak1, ak2, ratm, r):
+    def solver_wrapper(toa_cor_o3, tau, t1, t2, r0, ak1, ak2, ratm, r):
         def func_solv(albedo):
             return toa_cor_o3 - sl.alb2rtoa(albedo, t1, t2, r0, ak1, ak2, ratm, r)
         # it is assumed that albedo is in the range 0.1-1.0
+        
         return sl.zbrent(func_solv, 0.1, 1, 100, 1.e-6)
     
     solver_wrapper_v = np.vectorize(solver_wrapper)
     # loop over all bands except band 19, 20
     for i_channel in np.append(np.arange(18), [20]):
-        alb_sph[i_channel,ind_pol] = solver_wrapper_v(
-                toa_cor_o3[i_channel,ind_pol],
-                tau[i_channel,ind_pol], 
-                t1[i_channel,ind_pol], 
-                t2[i_channel,ind_pol],
-                r0[ind_pol], ak1[ind_pol], 
-                ak2[ind_pol],
-                ratm[i_channel,ind_pol], 
-                r[i_channel,ind_pol])
         
-        ind_bad = alb_sph[i_channel,:,:]==-999
-        alb_sph[i_channel,ind_bad] = np.nan
-        isnow[ind_bad]= -i_channel
+        alb_sph[i_channel, ind_pol] = solver_wrapper_v(
+            toa_cor_o3[i_channel, ind_pol], tau[i_channel, ind_pol], 
+            t1[i_channel, ind_pol], t2[i_channel, ind_pol],
+            r0[ind_pol], ak1[ind_pol], ak2[ind_pol], ratm[i_channel, ind_pol], 
+            r[i_channel, ind_pol])
+        
+        ind_bad = alb_sph[i_channel, :, :] == -999
+        alb_sph[i_channel, ind_bad] = np.nan
+        isnow[ind_bad] = -i_channel
     
     # INTERNal CHECK FOR CLEAN PIXELS
     # Are reprocessed as clean
-    ind_clear_pol1 = np.logical_and(ind_pol, alb_sph[0,:,:]>0.98)
-    ind_clear_pol2 = np.logical_and(ind_pol, alb_sph[1,:,:]>0.98)
+    ind_clear_pol1 = np.logical_and(ind_pol, alb_sph[0, :, :] > 0.98)
+    ind_clear_pol2 = np.logical_and(ind_pol, alb_sph[1, :, :] > 0.98)
     ind_clear_pol = np.logical_or(ind_clear_pol1, ind_clear_pol2)
-    isnow[ind_clear_pol]= 7
+    isnow[ind_clear_pol] = 7
+    
     for i_channel in range(21):
-        alb_sph[i_channel,ind_clear_pol] = np.exp(-np.sqrt(4.*1000.*al[ind_clear_pol] * np.pi * bai[i_channel] / w[i_channel] )) 
+        alb_sph[i_channel, ind_clear_pol] = np.exp(-np.sqrt(4. * 1000. 
+                                                            * al[ind_clear_pol] 
+                                                            * np.pi * bai[i_channel] 
+                                                            / w[i_channel])) 
         
     # re-defining polluted pixels
-    ind_pol =  np.logical_and(ind_pol, isnow!=7)
+    ind_pol = np.logical_and(ind_pol, isnow != 7)
     
-    #retrieving snow impurities        
+    # retrieving snow impurities        
     ntype, bf, conc = sl.snow_impurities(alb_sph, bal)
 
     # alex   09.06.2019
-    # reprocessing of albedo to remove gaseous absorption using linear polynomial approximation in the range 753-778nm.
-    # Meaning: alb_sph[12],alb_sph[13] and alb_sph[14] are replaced by a linear  interpolation between alb_sph[11] and alb_sph[15]
-    afirn=(alb_sph[15,ind_pol]-alb_sph[11,ind_pol])/(w[15]-w[11])
-    bfirn=alb_sph[15,ind_pol]-afirn*w[15]
-    alb_sph[12,ind_pol] = bfirn + afirn*w[12]
-    alb_sph[13,ind_pol] = bfirn + afirn*w[13]
-    alb_sph[14,ind_pol] = bfirn + afirn*w[14]             
+    # reprocessing of albedo to remove gaseous absorption using linear polynomial 
+    # approximation in the range 753-778nm.
+    # Meaning: alb_sph[12],alb_sph[13] and alb_sph[14] are replaced by a linear 
+    # interpolation between alb_sph[11] and alb_sph[15]
+    afirn = (alb_sph[15, ind_pol] - alb_sph[11, ind_pol]) / (w[15] - w[11])
+    bfirn = alb_sph[15, ind_pol] - afirn * w[15]
+    alb_sph[12, ind_pol] = bfirn + afirn * w[12]
+    alb_sph[13, ind_pol] = bfirn + afirn * w[13]
+    alb_sph[14, ind_pol] = bfirn + afirn * w[14]             
 
     # BAV 09-02-2020: 0.5 to 0.35
-    # pixels that are clean enough in channels 18 19 20 and 21 are not affected by pollution, the analytical equation can then be used
-    ind_ok =  np.logical_and(ind_pol, toa_cor_o3[20,:,:]>0.35)
-    for i_channel in range(17,21):
-        alb_sph[i_channel,ind_ok] = np.exp(-np.sqrt(4.*1000.*al[ind_ok] * np.pi * bai[i_channel] / w[i_channel] ))
+    # pixels that are clean enough in channels 18 19 20 and 21 are not affected 
+    # by pollution, the analytical equation can then be used
+    ind_ok = np.logical_and(ind_pol, toa_cor_o3[20, :, :] > 0.35)
+    
+    for i_channel in range(17, 21):
+        alb_sph[i_channel, ind_ok] = np.exp(-np.sqrt(4. * 1000. * al[ind_ok]
+                                                     * np.pi * bai[i_channel] 
+                                                     / w[i_channel]))
     # Alex, SEPTEMBER 26, 2019
-    # to avoid the influence of gaseous absorption (water vapor) we linearly interpolate in the range 885-1020nm for bare ice cases only (low toa[20])
-    # Meaning: alb_sph[18] and alb_sph[19] are replaced by a linear interpolation between alb_sph[17] and alb_sph[20]
-    delx=w[20]-w[17]
-    bcoef=(alb_sph[20,ind_pol]-alb_sph[17,ind_pol])/delx
-    acoef=alb_sph[20,ind_pol]-bcoef*w[20]
-    alb_sph[18,ind_pol] = acoef + bcoef*w[18]
-    alb_sph[19,ind_pol] = acoef + bcoef*w[19]
+    # to avoid the influence of gaseous absorption (water vapor) we linearly 
+    # interpolate in the range 885-1020nm for bare ice cases only (low toa[20])
+    # Meaning: alb_sph[18] and alb_sph[19] are replaced by a linear interpolation 
+    # between alb_sph[17] and alb_sph[20]
+    delx = w[20] - w[17]
+    bcoef = (alb_sph[20, ind_pol] - alb_sph[17, ind_pol]) / delx
+    acoef = alb_sph[20, ind_pol] - bcoef * w[20]
+    alb_sph[18, ind_pol] = acoef + bcoef * w[18]
+    alb_sph[19, ind_pol] = acoef + bcoef * w[19]
 
 # ========= derivation of plane albedo and reflectance =========== 
-rp = np.power (alb_sph, ak1)
-refl =r0* np.power(alb_sph, (ak1*ak2/r0))
+
+rp = np.power(alb_sph, ak1)
+refl = r0 * np.power(alb_sph, (ak1 * ak2 / r0))
 
 ind_all_clean = np.logical_or(ind_clean, isnow == 7)
 
-## CalCULATION OF BBA of clean snow
+# CalCULATION OF BBA of clean snow
 
 # old method: integrating equation
-#BBA_v = np.vectorize(sl.BBA_calc_clean)
-#p1,p2,s1,s2 = BBA_v(al[ind_all_clean], ak1[ind_all_clean])
-#
-## visible(0.3-0.7micron)
-#rp1[ind_all_clean]=p1/sol1_clean
-#rs1[ind_all_clean]=s1/sol1_clean
-## near-infrared (0.7-2.4micron)
-#rp2[ind_all_clean]=p2/sol2
-#rs2[ind_all_clean]=s2/sol2
-## shortwave(0.3-2.4 micron)
-#rp3[ind_all_clean]=(p1+p2)/sol3_clean
-#rs3[ind_all_clean]=(s1+s2)/sol3_clean
+# BBA_v = np.vectorize(sl.BBA_calc_clean)
+# p1, p2, s1, s2 = BBA_v(al[ind_all_clean], ak1[ind_all_clean])
+
+# visible(0.3-0.7micron)
+# rp1[ind_all_clean] = p1 / sol1_clean
+# rs1[ind_all_clean] = s1 / sol1_clean
+# near-infrared (0.7-2.4micron)
+# rp2[ind_all_clean] = p2 / sol2
+# rs2[ind_all_clean] = s2 / sol2
+# shortwave(0.3-2.4 micron)
+# rp3[ind_all_clean] = (p1 + p2) / sol3_clean
+# rs3[ind_all_clean] = (s1 + s2) / sol3_clean
 
 # approximation
 # planar albedo
-#rp1 and rp2 not derived anymore
-rp3[ind_all_clean]=sl.plane_albedo_sw_approx(D[ind_all_clean],cos_sza[ind_all_clean])
-#     spherical albedo
-#rs1 and rs2 not derived anymore
-rs3[ind_all_clean]= sl.spher_albedo_sw_approx(D[ind_all_clean])
+# rp1 and rp2 not derived anymore
+rp3[ind_all_clean] = sl.plane_albedo_sw_approx(D[ind_all_clean], 
+                                               am1[ind_all_clean])
+# spherical albedo
+# rs1 and rs2 not derived anymore
+rs3[ind_all_clean] = sl.spher_albedo_sw_approx(D[ind_all_clean])
     
 # calculation of the BBA for the polluted snow
 rp1[ind_pol], rp2[ind_pol], rp3[ind_pol] = sl.BBA_calc_pol(
-        rp[:, ind_pol], asol, sol1_pol, sol2, sol3_pol)
+    rp[:, ind_pol], asol, sol1_pol, sol2, sol3_pol)
 rs1[ind_pol], rs2[ind_pol], rs3[ind_pol] = sl.BBA_calc_pol(
-        alb_sph[:, ind_pol], asol, sol1_pol, sol2, sol3_pol)
+    alb_sph[:, ind_pol], asol, sol1_pol, sol2, sol3_pol)
                
-#%% Output
-if os.path.isfile(sys.argv[1]):  
-    
-    print('\nText file output')
-    # data_in = pd.read_csv(sys.argv[1])
-    data_out = data_in
-    data_out['grain_diameter']=D
-    data_out['snow_specific_area']=area
-    data_out['al']=al
-    data_out['r0']=r0
-    data_out['diagnostic_retrieval']=isnow
-    data_out['conc']=conc
-    data_out['albedo_bb_planar_sw']=rp3
-    data_out['albedo_bb_spherical_sw']=rs3
-    
+# %% Output
+
+WriteOutput(BXXX, 'O3_SICE', InputFolder)
+WriteOutput(D, 'grain_diameter', InputFolder)
+WriteOutput(area, 'snow_specific_surface_area', InputFolder)
+WriteOutput(al, 'al', InputFolder)
+WriteOutput(r0, 'r0', InputFolder)
+WriteOutput(isnow, 'diagnostic_retrieval', InputFolder)
+WriteOutput(conc, 'conc', InputFolder)
+WriteOutput(rp3, 'albedo_bb_planar_sw', InputFolder)
+WriteOutput(rs3, 'albedo_bb_spherical_sw', InputFolder)
+
+for i in np.append(np.arange(11), np.arange(15, 21)):
     
-    for i in np.append(np.arange(11), np.arange(15,21)):
     # for i in np.arange(21):
-        data_out['albedo_spectral_spherical_'+str(i+1).zfill(2)]=alb_sph[i,:,:]
-    for i in np.append(np.arange(11), np.arange(15,21)):
-        data_out['rBRR_'+str(i+1).zfill(2)]=rp[i,:,:]
-
-    data_out.to_csv(sys.argv[1][:-4]+'_out.csv')
-# ========= input tif ===============
-elif os.path.isdir(sys.argv[1]):
-    WriteOutput(BXXX,   'O3_SICE',   InputFolder)
-    WriteOutput(D,      'grain_diameter',InputFolder)
-    WriteOutput(area,   'snow_specific_area', InputFolder)
-    WriteOutput(al,   'al',     InputFolder)
-    WriteOutput(r0,   'r0',InputFolder)
-    WriteOutput(isnow,'diagnostic_retrieval',InputFolder)
-    WriteOutput(conc, 'conc',InputFolder)
-    WriteOutput(rp3,  'albedo_bb_planar_sw',InputFolder)
-    WriteOutput(rs3,  'albedo_bb_spherical_sw',InputFolder)
-    
-    for i in np.append(np.arange(11), np.arange(15,21)):
-    # for i in np.arange(21):
-        WriteOutput(alb_sph[i,:,:],    'albedo_spectral_spherical_'+str(i+1).zfill(2), InputFolder)
-        WriteOutput(rp[i,:,:],    'albedo_spectral_planar_'+str(i+1).zfill(2), InputFolder)
-        WriteOutput(refl[i,:,:],   'rBRR_'+str(i+1).zfill(2), InputFolder)
-
-print("End SICE.py %s --- %s CPU seconds ---" % (InputFolder, time.process_time() - start_time))
+    WriteOutput(alb_sph[i, :, :], 'albedo_spectral_spherical_' 
+                + str(i + 1).zfill(2), InputFolder)
+    WriteOutput(rp[i, :, :], 'albedo_spectral_planar_' 
+                + str(i + 1).zfill(2), InputFolder)
+    WriteOutput(refl[i, :, :], 'rBRR_' 
+                + str(i + 1).zfill(2), InputFolder)
+
+print("End SICE.py %s --- %s CPU seconds ---" % 
+      (InputFolder, time.process_time() - start_time))
diff --git a/sice_lib.py b/sice_lib.py
index 83e8408..a8d784c 100644
--- a/sice_lib.py
+++ b/sice_lib.py
@@ -83,7 +83,8 @@
 import numpy as np
 from constants import w, bai, xa, ya, f0, f1, f2, bet, gam, coef1, coef2, coef3, coef4
 
-#%% ================================================
+# %% ================================================
+
 # 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
 
@@ -93,34 +94,42 @@
 # 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
+
+# %% viewing characteristics and aerosol properties
+
 # sza                       solar zenith angle
 # vza                       viewing zenith angle
 # saa                       solar azimuthal angle
@@ -129,195 +138,214 @@ 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
-#%%     
-def aerosol_properties(aot, height, cos_sa):
+    return raa, am1, am2, ak1, ak2, amf, co
+
+# %%     
+
+
+def aerosol_properties(aot, height, co):
     # Atmospheric optical thickness
-    tauaer =aot*(w/0.5)**(-1.3)
-
-    ad =height/7400.
-    ak = height*0+1
-    ak[ad > 1.e-6]=np.exp(-ad[ad > 1.e-6])
-    
-    taumol = np.tile(height*np.nan, (21,1,1))
-    tau = np.tile(height*np.nan, (21,1,1))
-    g = np.tile(height*np.nan, (21,1,1))
-    pa = np.tile(height*np.nan, (21,1,1))
-    p = np.tile(height*np.nan, (21,1,1))
-    g0=0.5263
-    g1=0.4627
-    wave0=0.4685
-    gaer=g0+g1*np.exp(-w/wave0)
-    pr=0.75*(1.+cos_sa**2)
+    tauaer = aot * (w / 0.5) ** (-1.3)
+
+    ad = height / 7400.
+    ak = height * 0 + 1
+    ak[ad > 1.e-6] = np.exp(-ad[ad > 1.e-6])
+    
+    taumol = np.tile(height * np.nan, (21, 1, 1))
+    tau = np.tile(height * np.nan, (21, 1, 1))
+    g = np.tile(height * np.nan, (21, 1, 1))
+    pa = np.tile(height * np.nan, (21, 1, 1))
+    p = np.tile(height * np.nan, (21, 1, 1))
+    g0 = 0.5263
+    g1 = 0.4627
+    wave0 = 0.4685
+    gaer = g0 + g1 * np.exp(-w / wave0)
+    pr = 0.75 * (1. + co ** 2)
     
     for i in range(21):
-        taumol[i,:,:] = ak*0.00877/w[i]**(4.05)
-        tau[i,:,:] = tauaer[i] + taumol[i,:,:]
+        
+        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,:,:]*pr + tauaer[i]*pa[i,:,:])/tau[i,:,:]
+        p[i, :, :] = (taumol[i, :, :] * pr + tauaer[i] * pa[i, :, :]) / tau[i, :, :]
     
-    return tau, p, g, gaer,taumol,tauaer
+    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)                       
     
-    akap2=2.25e-6                    
-    alpha2=4.*np.pi*akap2/1.020                        
+    akap2 = 2.25e-6                    
+    alpha2 = 4. * np.pi * akap2 / 1.020                        
     eps = 1.549559365010611
     
     # 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) * np.log(rr2/r0)/alpha2/(ak1*ak2/r0)**2
-    al = bal/1000.
+    bal = np.log(rr2 / r0) * np.log(rr2 / r0) / alpha2 / (ak1 * ak2 / r0) ** 2
+    al = bal / 1000.
     
     # effective grain size(mm):diameter
-    D=al/16.36              
+    D = al / 16.36              
     # 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):
-    astra=tau*np.nan
-    rms=tau*np.nan
-    t1=tau*np.nan
-    t2=tau*np.nan
+def prepare_coef(tau, g, p, am1, am2, amf, gaer, taumol, tauaer):
+    astra = tau * np.nan
+    rms = tau * np.nan
+    t1 = tau * np.nan
+    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):
-        astra[i,:,:]=(1.-np.exp(-tau[i,:,:]*inv_cos_za))/(cos_sza+cos_vza)/4.
-        rms[i,:,:] = 1.- b1[i,:,:]*b2[i,:,:]/oskar[i,:,:]  \
-        + (3.*(1.+g[i,:,:])*cos_sza*cos_vza - 2.*(cos_sza+cos_vza))*astra[i,:,:]
-        #backscattering fraction
-        # t1[i,:,:] = np.exp(-(1.-g[i,:,:])*tau[i,:,:]/cos_sza/2.)
-        # t2[i,:,:] = np.exp(-(1.-g[i,:,:])*tau[i,:,:]/cos_vza/2.)
-        t1[i,:,:]=np.exp(-(1.-g[i,:,:])*tau[i,:,:]/cos_sza/2./sssss[i,:,:])
-        t2[i,:,:]=np.exp(-(1.-g[i,:,:])*tau[i,:,:]/cos_vza/2./sssss[i,:,:])
+        
+        astra[i, :, :] = (1. - np.exp(-tau[i, :, :] * amf)) / (am1 + am2) / 4.
+        rms[i, :, :] = 1. - b1[i, :, :] * b2[i, :, :] / oskar[i, :, :] \
+            + (3. * (1. + g[i, :, :]) * am1 * am2 - 2. * (am1 + am2)) * astra[i, :, :]
+        # backscattering fraction
+        # t1[i, :, :] = np.exp(-(1. - g[i, :, :]) * tau[i, :, :] / am1 / 2.)
+        # t2[i, :, :] = np.exp(-(1. - g[i, :, :]) * tau[i, :, :] / am2 / 2.)
+        t1[i, :, :] = np.exp(-(1. - g[i, :, :]) * tau[i, :, :] / am1 / 2.
+                             / sssss[i, :, :])
+        t2[i, :, :] = np.exp(-(1. - g[i, :, :]) * tau[i, :, :] / am2 / 2.
+                             / sssss[i, :, :])
       
-    rss = p*astra
+    rss = p * astra
     r = rss + rms
     
     # SALBED
-#    ratm = salbed(tau, g)
-    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_cst =     a_s[0]*g**0  + a_s[1]*g**1 + a_s[2]*g**2 + a_s[3]*g**3
-    b_cst =     bs[0]*g**0   + bs[1]*g**1 + bs[2]*g**2 + bs[3]*g**3
-    c_cst =     cs[0]*g**0   + cs[1]*g**1 + cs[2]*g**2 + cs[3]*g**3
-    al_cst=     als[0]*g**0  + als[1]*g**1 + als[2]*g**2 + als[3]*g**3
-    bet_cst=    bets[0]*g**0 + bets[1]*g**1 + bets[2]*g**2 + bets[3]*g**3        
+    # ratm = salbed(tau, g)
+    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_cst = a_s[0] * g ** 0 + a_s[1] * g ** 1 + a_s[2] * g ** 2 + a_s[3] * g ** 3
+    b_cst = bs[0] * g ** 0 + bs[1] * g ** 1 + bs[2] * g ** 2 + bs[3] * g ** 3
+    c_cst = cs[0] * g ** 0 + cs[1] * g ** 1 + cs[2] * g ** 2 + cs[3] * g ** 3
+    al_cst = als[0] * g ** 0 + als[1] * g ** 1 + als[2] * g ** 2 + als[3] * g ** 3
+    bet_cst = bets[0] * g ** 0 + bets[1] * g ** 1 + bets[2] * g ** 2 + bets[3] * g ** 3        
    
-    ratm = tau*(a_cst*np.exp(-tau/al_cst)+b_cst*np.exp(-tau/bet_cst)+c_cst)
+    ratm = tau * (a_cst * np.exp(-tau / al_cst) + b_cst * np.exp(-tau / bet_cst) 
+                  + c_cst)
+    
     return t1, t2, ratm, r, astra, rms
-#%% snow_imputirities
+
+# %% 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
 
-#%% ===========================================================================
+# %% ===========================================================================
+
+
 def salbed(tau, g):
     # WARNING: NOT USED ANYMORE
     # SPHERICAL ALBEDO OF TERRESTRIAL ATMOSPHERE:      
@@ -327,29 +355,31 @@ 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):
     # Equation solver using Brent's method
     # 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 
@@ -360,10 +390,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):
@@ -375,7 +405,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)
@@ -387,13 +417,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
  
@@ -415,7 +445,9 @@ def zbrent(f, x0, x1, max_iter=100, tolerance=1e-6):
  
     return x1
 
-#%% =====================================================================     
+# %% =====================================================================     
+
+
 def funp(x, al, sph_calc, ak1):
     #     Spectral planar albedo
     # Inputs:
@@ -433,104 +465,125 @@ 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):
-    anka= 0.7389  -0.1783*cos_sza    +0.0484*cos_sza**2.
-    banka=0.0853  +0.0414*cos_sza    -0.0127*cos_sza**2.
-    canka=0.1384  +0.0762*cos_sza    -0.0268*cos_sza**2.
-    diam1=187.89  -69.2636*cos_sza     +40.4821*cos_sza**2.
-    diam2=2687.25 -405.09*cos_sza   +94.5*cos_sza**2.
-    return anka+banka*np.exp(-1000*D/diam1)+canka*np.exp(-1000*D/diam2)
 
 def spher_albedo_sw_approx(D):
-    anka= 0.6420
-    banka=0.1044
-    canka=0.1773
-    diam1=158.62
-    diam2=2448.18
-    return anka+banka*np.exp(-1000*D/diam1)+canka*np.exp(-1000*D/diam2)
-
-#%%   CalCULATION OF BBA for clean pixels
+    anka = 0.6420
+    banka = 0.1044
+    canka = 0.1773
+    diam1 = 158.62
+    diam2 = 2448.18
+    return anka + banka * np.exp(-1000 * D / diam1) + canka \
+        * np.exp(-1000 * D / diam2)
+
+# %%   CalCULATION OF BBA for clean pixels
+
+
 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):
+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)<eps*abs(os)):
+                
+                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) < 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
+# %% Calculation f BBA for polluted snow
+
 
 def BBA_calc_pol(alb, asol, sol1_pol, sol2, sol3_pol):
     # polluted snow
@@ -543,73 +596,76 @@ 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
-    rati=r7/r8
-    alasta = (alam8-alam7)/np.log(rati)
-    an=1./alasta
-    p   = r7 * np.exp(alam7/alasta)
-    
-    aj31=(1./an)*(np.exp(-an*z2)-np.exp(-an*z1))
-    aj32=(1./(bet+an))*(np.exp(-(bet+an)*z2)-np.exp(-(an+bet)*z1))
-    aj33=(1./(gam+an))*(np.exp(-(gam+an)*z2)-np.exp(-(an+gam)*z1))
-    aj3=(-f0*aj31-f1*aj32-f2*aj33)*p
-    
-    BBA_vis = aj1/sol1_pol
-    BBA_nir = (aj2+aj3)/sol2 #here segment 2.1 and 2.2 are summed
-    BBA_sw   = (aj1+aj2+aj3)/sol3_pol 
-
-    return BBA_vis,BBA_nir, BBA_sw
-
-#%% ==========================================================================
-def quad_func(x0,x1,x2,y0,y1,y2):
+    z1 = 0.865
+    z2 = 2.4
+    rati = r7 / r8
+    alasta = (alam8 - alam7) / np.log(rati)
+    an = 1. / alasta
+    p = r7 * np.exp(alam7 / alasta)
+    
+    aj31 = (1. / an) * (np.exp(-an * z2) - np.exp(-an * z1))
+    aj32 = (1. / (bet + an)) * (np.exp(-(bet + an) * z2) - np.exp(-(an + bet) * z1))
+    aj33 = (1. / (gam + an)) * (np.exp(-(gam + an) * z2) - np.exp(-(an + gam) * z1))
+    aj3 = (-f0 * aj31 - f1 * aj32 - f2 * aj33) * p
+    
+    BBA_vis = aj1 / sol1_pol
+    BBA_nir = (aj2 + aj3) / sol2  # here segment 2.1 and 2.2 are summed
+    BBA_sw = (aj1 + aj2 + aj3) / sol3_pol 
+
+    return BBA_vis, BBA_nir, BBA_sw
+
+# %% ==========================================================================
+
+
+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