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')
plt.ylabel('Radiance(W/m^2/sr/micrometer)')
plt.legend()

plt.show()

This snippet took 0.01 seconds to highlight.

Back to the Entry List or Home.

Delete this entry (admin only).