#This code simulate TOA reflectance and downwelling irradiance, for typical atmosphere, typical ocean, and geometry, as a function of wavelength, from 350 to 900 nm by step of 5 nm.

import numpy as np
import os
import matplotlib.pyplot as plt

def write_6S_input(filename,solzen,solaz,vzen,vaz,month,day,watvap,o3,aot550,ws,winf,wsup):
    f=open(filename,'w')
    f.write('0 (User defined)\n')
    f.write('%f %f %f %f %d %d\n' % (solzen,solaz,vzen,vaz,month,day))
    f.write('8\n')
    f.write('%f %f\n' % (watvap,o3))
    f.write('2 Maritime Model\n')
    f.write('0\n')
    f.write('%f value\n' % aot550)
    f.write('0 (target level)\n')
    f.write('-1000 (sensor level)\n')
    f.write('0\n')
    f.write('%f %f\n' % (winf/1000,wsup/1000))
    f.write('0 Homogeneous surface\n')
    f.write('1 (dirctional effects)\n')
    f.write('6 Ocean\n')
    f.write('%f 45 35 %f\n' % (ws,chl))
    f.write('-2 No atm. corrections selected\n')
    f.close()

def read_6S_output(filename):
    f=open(filename,'r')
    lines=f.readlines()
    f.close()
    for i in range(len(lines)):
        line=lines[i]
        if "apparent reflectance" in line:
            rho_toa=float(line.split()[3])
        if "app. polarized refl." in line:
            rho_pol_toa=float(line.split()[4])
    return rho_toa,rho_pol_toa


if not os.path.exists('input2'):
    os.system('mkdir input2') # create a folder 'input' to save all the 6S input files

if not os.path.exists('output2'):
    os.system('mkdir output2') # create a folder 'output' to save all the 6S output files

month=1
day=1
sza=30.01 # solar zenith
vza=np.arange(0,76,5) # view zenith
saz=0. # solar azimuth
vaz=np.array([0,180]) # view azimuth
watvap=2.  # water vapor, g/cm2
o3=0.3 # ozone, cm-atm
aot550=0.1 # aerosol optical thickness @ 550 nm
ws=np.array([5,10]) # wind speed, m/s
chl=0.1 # chlorophyll concentration, mg/m3

#wl=np.arange(400,901,5) # from 400 to 900 nm by step of 5 nm
wl=450.

rho_toa=np.zeros((vza.size,vaz.size,ws.size))+np.nan
rho_pol_toa=np.zeros((vza.size,vaz.size,ws.size))+np.nan


for i in range(vza.size):
    for j in range(vaz.size):
        for k in range(ws.size):
            #write_6S_input('input2/sim_wl{}_vza{}_vaz{}_ws{}.txt'.format(wl,vza[i],vaz[j],ws[k]),\
                           #sza,saz,vza[i],vaz[j],month,day,watvap,o3,aot550,ws[k],wl-2.5,wl+2.5)
    
            #os.system('6SV1.1/sixsV1.1 < input2/sim_wl{}_vza{}_vaz{}_ws{}.txt > output2/sim_wl{}_vza{}_vaz{}_ws{}.txt'.format(wl,vza[i],vaz[j],ws[k],wl,vza[i],vaz[j],ws[k]))

            rho_toa[i,j,k],rho_pol_toa[i,j,k]=read_6S_output('output2/sim_wl{}_vza{}_vaz{}_ws{}.txt'.format(wl,vza[i],vaz[j],ws[k]))



fig,ax=plt.subplots(figsize=(8,3))
# Move the left and bottom spines to x = 0 and y = 0, respectively.
ax.spines[["left", "bottom"]].set_position(("data", 0))
# Hide the top and right spines.
ax.spines[["top", "right"]].set_visible(False)
# Draw arrows (as black triangles: ">k"/"^k") at the end of the axes.  In each
# case, one of the coordinates (0) is a data coordinate (i.e., y = 0 or x = 0,
# respectively) and the other one (1) is an axes coordinate (i.e., at the very
# right/top of the axes).  Also, disable clipping (clip_on=False) as the marker
# actually spills out of the Axes.
ax.plot(1, 0, ">k", transform=ax.get_yaxis_transform(), clip_on=False)
ax.plot(0, 1, "^k", transform=ax.get_xaxis_transform(), clip_on=False)
ax.plot(0, 0, "<k", transform=ax.get_yaxis_transform(), clip_on=False)   # Left arrow (new)
ax.plot(vza,rho_toa[:,0,0],'b',label='total reflectance')
ax.plot(vza,rho_toa[:,0,1],'b--')
ax.plot(vza,rho_pol_toa[:,0,0],'r',label='polarized reflectance')
ax.plot(vza,rho_pol_toa[:,0,1],'r--')
ax.plot(-vza,rho_toa[:,1,0],'b')
ax.plot(-vza,rho_toa[:,1,1],'b--')
ax.plot(-vza,rho_pol_toa[:,1,0],'r')
ax.plot(-vza,rho_pol_toa[:,1,1],'r--')
plt.xticks(np.arange(-80,81,20),labels=['80','60','40','20','0','20','40','60','80'])
plt.yticks(np.linspace(0,0.25,6),labels=['','0.05','0.10','0.15','0.20','0.25'])
plt.text(35,0.26,'RAZ=0$^{\circ}$')
plt.text(-55,0.26,'RAZ=180$^{\circ}$')
plt.text(5,0.23,'SZA=30$^{\circ}$')
plt.text(5,0.21,'Chl=0.1 mg/m$^{3}$')
plt.text(5,0.19,'AOT550=0.1, maritime')
plt.text(5,0.17,'UO3 300 Dobson')
plt.text(5,0.15,'UH2O 2 g/cm$^{2}$')
plt.text(20,0.08,'solid line: U = 5 m/s')
plt.text(20,0.06,'dashed line: U = 10 m/s')
plt.title('TOA reflectance')
plt.text(-65,0.22,'450 nm',fontweight='bold')
plt.xlabel('VZA, degree')
plt.ylabel('')
plt.legend(frameon=False)
plt.savefig('6S_presentation.png',dpi=300)
plt.close()