;;
;; 横田さんHapkeプログラムをIDLに移植した版
;;
Forward_function func_RADF1984
Forward_function func_HenyeyGreen_1param
Forward_function func_psi
Forward_function func_phase
Forward_function func_E1
Forward_function func_E2
Forward_function func_Cai
Forward_function func_Eta
function hapke_models
;;
;; Wrap関数とする ?
;;
end
;;
;;
;;
function func_RADF1984, i_deg, e_deg, phase_deg, w, g, theta_deg, B0, h
;; 初期化 ;;
ans = i_deg * 0d
ans_pos = where( (i_deg lt 90d and e_deg lt 90d) and (i_deg ge -180 and e_deg ge -180), t_count)
if t_count eq 0 then begin
return, ans ;; 0
endif
i_rad = i_deg[ans_pos] * !dpi/180.
e_rad = e_deg[ans_pos] * !dpi/180.
phase_rad = phase_deg[ans_pos] * !dpi/180.
theta_rad = theta_deg * !dpi/180.
psi = func_psi(i_rad, e_rad, phase_rad)
f = exp(-2d * tan(psi/2d))
Eta_i = func_Eta(theta_rad, i_rad)
Eta_e = func_Eta(theta_rad, e_rad)
Cai = func_Cai(theta_rad)
E1_i = func_E1(theta_rad, i_rad)
E1_e = func_E1(theta_rad, e_rad)
E2_i = func_E2(theta_rad, i_rad)
E2_e = func_E2(theta_rad, e_rad)
;; 初期化 ;;
S = i_rad * 0
Mu_eff = i_rad * 0
Mu0_eff = i_rad * 0
;;;;;;;;;;;
;; Case 1;;
;;;;;;;;;;;
i_rad_pos = where(i_rad le e_rad, i_count)
if i_count gt 0 then begin
tmp1 = cos(psi[i_rad_pos]) * E2_e[i_rad_pos] + (sin(psi[i_rad_pos]/2d))^2. * E2_i[i_rad_pos]
tmp2 = 2d - E1_e[i_rad_pos] - (psi[i_rad_pos]/!dpi) * E1_i[i_rad_pos]
Mu0_eff[i_rad_pos] = Cai * (cos(i_rad[i_rad_pos]) + sin(i_rad[i_rad_pos])*tan(theta_rad) * tmp1/tmp2)
tmp1 = E2_e[i_rad_pos] - (sin(psi[i_rad_pos]/2d))^2. * E2_i[i_rad_pos]
tmp2 = 2d - E1_e[i_rad_pos] - (psi[i_rad_pos]/!dpi) * E1_i[i_rad_pos]
Mu_eff[i_rad_pos] = Cai * (cos(e_rad[i_rad_pos]) + sin(e_rad[i_rad_pos])*tan(theta_rad) * tmp1/tmp2)
S[i_rad_pos] = (Mu_eff[i_rad_pos]/Eta_e[i_rad_pos]) * (cos(i_rad[i_rad_pos])/Eta_i[i_rad_pos]) $
* Cai/(1d - f[i_rad_pos] + f[i_rad_pos] * Cai * (cos(i_rad[i_rad_pos]) / Eta_i[i_rad_pos] ))
endif
;;;;;;;;;;;
;; Case 2;;
;;;;;;;;;;;
e_rad_pos = where(i_rad gt e_rad, e_count)
if e_count gt 0 then begin
tmp1 = E2_i[e_rad_pos] - (sin(psi[e_rad_pos]/2d))^2. * E2_e[e_rad_pos]
tmp2 = 2d - E1_i[e_rad_pos] - psi[e_rad_pos]/!dpi * E1_e[e_rad_pos]
Mu0_eff[e_rad_pos] = Cai * (cos(i_rad[e_rad_pos]) + sin(i_rad[e_rad_pos])*tan(theta_rad) * tmp1/tmp2)
tmp1 = cos(psi[e_rad_pos])*E2_i[e_rad_pos] + (sin(psi[e_rad_pos]/2d))^2. * E2_e[e_rad_pos]
tmp2 = 2d - E1_i[e_rad_pos] - psi[e_rad_pos]/!dpi * E1_e[e_rad_pos]
Mu_eff[e_rad_pos] = Cai * (cos(e_rad[e_rad_pos]) + sin(e_rad[e_rad_pos])*tan(theta_rad) * tmp1/tmp2)
S[e_rad_pos] = (Mu_eff[e_rad_pos]/Eta_e[e_rad_pos]) * (cos(i_rad[e_rad_pos])/Eta_i[e_rad_pos]) $
* Cai/(1d - f[e_rad_pos] + f[e_rad_pos] * Cai * (cos(e_rad[e_rad_pos]) / Eta_e[e_rad_pos]))
endif
;; Wrap up ;;
BSOE = B0/ (1d + tan(phase_rad / 2d) / h)
p = func_HenyeyGreen_1param(phase_rad, g)
H1 = (1d + 2d*Mu0_eff) / (1d + 2d * Mu0_eff * sqrt(1d - w) )
H2 = (1d + 2d*Mu_eff) / (1d + 2d * Mu_eff * sqrt(1d - w) )
ans[ans_pos] = w/4d * Mu0_eff/( Mu0_eff + Mu_eff) * ((1d + BSOE)*p + H1*H2 - 1d) * S
return, ans
end
;;
;;
;;
function func_HenyeyGreen_1param, phase_rad, g
tmp1 = 1d - g*g
tmp2 = 1d + 2d*g*cos(phase_rad) + g*g
ans = tmp1/(tmp2*sqrt(tmp2))
return, ans
end
;;
;;
;;
function func_psi, i_rad, e_rad, phase_rad
tmp1 = cos(phase_rad) - cos(e_rad) * cos(i_rad)
tmp2 = sin(e_rad) * sin(i_rad)
tmp2_pos = where(tmp2 ne 0)
tmp3 = tmp1[tmp2_pos]/tmp2[tmp2_pos]
tmp3 = tmp3 > (-1d) < 1d
;; 初期化 ;;
ans = i_rad * 0d
ans[tmp2_pos] = acos(tmp3)
return, ans
end
;;
;;
;;
function func_phase, i_rad, e_rad, psi
tmp1 = cos(e_rad) * cos(i_rad) + sin(e_rad) * sin(i_rad) * cos(psi)
tmp1 = tmp1 > (-1d) < 1d
ans = acos(tmp1)
return, ans
end
;;
;;
;;
function func_E1, theta_rad, y
;; 初期化 ;;
;ans = dblarr(n_elements(theta_rad))
ans = y * 0d
ans_pos = where(y gt 0)
tmp1 = -2d/!dpi / tan(theta_rad) / tan(y[ans_pos])
ans[ans_pos] = exp(tmp1)
return, ans
end
;;
;;
;;
function func_E2, theta_rad, y
;; 初期化 ;;
ans = y * 0d
ans_pos = where(y gt 0)
tmp1 = -1d/!dpi / (tan(theta_rad) / tan(y[ans_pos]))^2.
ans[ans_pos] = exp(tmp1)
return, ans
end
;;
;;
;;
function func_Cai, theta_rad
tmp1 = 1d + !dpi * tan(theta_rad) * tan(theta_rad)
ans = 1d/sqrt(tmp1)
return, ans
end
;;
;;
;;
function func_Eta, theta_rad, y
tmp1 = cos(y) + sin(y) * tan(theta_rad) * func_E2(theta_rad, y) / (2d - func_E1(theta_rad, y))
ans = func_Cai(theta_rad) * tmp1
return, ans
end