# Demo entry 6786027

Judy Yao

Submitted by anonymous on Mar 22, 2019 at 16:08
Language: Python 3. Code size: 3.0 kB.

```# This code is for RS HW5.
# Programmer: Zhu Yao
# Date: March 20, 2019

import numpy as np
import matplotlib.pyplot as plt

# ----------------------- Function Definition----------------------------------

# define the initial solar radiance at different wavelengths

def SR(lam):
a = 1.2e14/pow(lam, 5)
b = np.exp(2.5/lam)-1
sr = a/b
return(sr)

# --------------Solar Radiance reaching at the TOA(scale height)---------------
Ltoa1 = SR(0.45)
Ltoa2 = SR(0.55)
Ltoa3 = SR(0.65)

# --------------Solar Radiance after distance D -------------------------------
LD1 = np.zeros(shape=(90, 1000))
LD2 = np.zeros(shape=(90, 1000))
LD3 = np.zeros(shape=(90, 1000))

# parameters
D = np.zeros(shape=(90, 1000))  # distance before the fourth term
S = np.zeros(shape=(90, 1000))  # distance after the fourth term
cossun = np.cos(51.53*np.pi/180)
beta1 = 1e-6/pow(0.45, 4)
beta2 = 1e-6/pow(0.55, 4)
beta3 = 1e-6/pow(0.65, 4)

for i in range(90):
step = (8553/np.cos(i*np.pi/180))/1000  # 8553m: RT/g
for j in range(1000):
S[i, j] = j*step
D[i, j] = (8553 - S[i, j]*np.cos(i*np.pi/180))/cossun
LD1[i, j] = Ltoa1*np.exp(-beta1*D[i, j])
LD2[i, j] = Ltoa2*np.exp(-beta2*D[i, j])
LD3[i, j] = Ltoa3*np.exp(-beta3*D[i, j])

# --------------initial Radiance after scattered -------------------------------
LSbeg1 = np.zeros(shape=(90, 1000))
LSbeg2 = np.zeros(shape=(90, 1000))
LSbeg3 = np.zeros(shape=(90, 1000))
costheta = np.zeros(90)

for i in range(90):
costheta[i] = 0.34*np.sin(i*np.pi/180)+0.62*np.cos(i*np.pi/180)
step = (8553/np.cos(i*np.pi/180))/1000
for j in range(1000):
cons1 = (1./(4*np.pi))*beta1*step*6.1e-5
cons2 = (1./(4*np.pi))*beta2*step*6.1e-5
cons3 = (1./(4*np.pi))*beta3*step*6.1e-5
LSbeg1[i, j] = cons1*LD1[i, j]*(0.75*(1+np.square(costheta[i])))
LSbeg2[i, j] = cons2*LD2[i, j]*(0.75*(1+np.square(costheta[i])))
LSbeg3[i, j] = cons3*LD3[i, j]*(0.75*(1+np.square(costheta[i])))

# --------------Final Radiance after scattered over S---------------------------
LSend1 = np.zeros(shape=(90, 1000))
LSend2 = np.zeros(shape=(90, 1000))
LSend3 = np.zeros(shape=(90, 1000))

for i in range(90):
step = (8553/np.cos(i*np.pi/180))/1000  # 8553m: RT/g
for j in range(1000):
S[i, j] = j*step
LSend1[i, j] = LSbeg1[i, j]*np.exp(-beta1*S[i, j])
LSend2[i, j] = LSbeg2[i, j]*np.exp(-beta2*S[i, j])
LSend3[i, j] = LSbeg3[i, j]*np.exp(-beta3*S[i, j])
print(LSend1[30, :])
LS1 = np.zeros(shape=(90, 1))
LS2 = np.zeros(shape=(90, 1))
LS3 = np.zeros(shape=(90, 1))

LS1 = np.sum(LSend1[:, :], axis=1)
LS2 = np.sum(LSend2[:, :], axis=1)
LS3 = np.sum(LSend3[:, :], axis=1)

fig1 = plt.figure()
x = np.arange(90)
plt.xlim(0, 90)
plt.plot(x, LS1[:]*1e-6, label='lambda=0.45')
plt.plot(x, LS2[:]*1e-6, label='lambda=0.55')
plt.plot(x, LS3[:]*1e-6, label='lambda=0.65')
plt.xlabel('viewer''s angle')