📄 gomsphere_f.htm
字号:
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<!-- saved from url=(0049)http://www.ees.nmt.edu/zhou/GOMsphere/GOMsphere.f -->
<HTML><HEAD>
<META http-equiv=Content-Type content="text/html; charset=gb2312">
<META content="MSHTML 6.00.2800.1106" name=GENERATOR></HEAD>
<BODY><PRE> subroutine GOMsphere(lamda,radius,theta,
& Qabs,QNscat,QFscat,gN,gF,PN,PF)
c Purpose:
c To calculate absorption coefficienct, scattering efficiency, asymmetry factor and phase
c function for both near and far fields.
c Input variables:
c lamda: wavelength in meter.
c radius: radius of the spherical scatters ( snow grains)in meters.
c theta: scattering angle.
c 1. If only need to calculate Qabs, just specify lamda and radius as inputs.
c 2. If only need to calculate Qscat and g, just specify lamda and radius as inputs.
c 3. theta must be specified only for phase (PN, PF) calculation. For other calculation,theta
c can be given any value, which do not affect the results.
c Output variables:
c Qabs: absorption efficiency. It is the same for both near and far fields.
c QNscat: scattering efficiency for near field scattering.
c QFscat: scattering efficiency for far field scattering.
c gN: asymmetry factor for near-filed scattering.
c gF: asymmetry factor for far-field scattering.
c PN: phase function for near-field scattering.
c PF: phase function for far-field scattering.
c Local variable:
c n: number of roots for Gauss-Legendre n-point quadratures.
c NOTE:
c This program is distributed publically along a paper published in Applied Optics and a report documented
c the details of the algorithm. Citation is suggested as:
c Zhou, X., S. Li, and K. Stamnes, A new geometrical optics code for computing optical properties of large
c dielectric spheres, Applied Optics, 42 (21), 4295-4306, 2003.
c
c Any errors or comments please direct them to Dr. Xiaobing Zhou at xzhou@nmt.edu
c or xzhou_2000_2001@yahoo.com. Thanks. May, 2003.
c
real*8 lamda, radius
real*8 Qabs,QNscat,QFscat,gN,gF,theta,PN,PF
integer*4 n
n=20
c To calcualte the absorption efficiency Qabs.
call abs_efficiency(lamda,radius,n,Qabs)
c To calcualte the scattering efficiency and asymmetry factor for both near and far fields
call QscatAndG(lamda,radius,n,QNscat,QFscat,gN,gF)
c To calculate the phase function for at a scattering angle.
call phasefunction(lamda,radius,theta,QNscat,QFscat,
& PN,PF)
return
end
c **************************************************************
c *********To calculate the absorption efficiency Qabs**********
c **************************************************************
subroutine abs_efficiency(lamda,radius,n,Qabs)
c PURPOSE:
c To calculate absorption efficiency, given wavelength, radius of sphere
c REFERENCE: Davis P J, P Rainowitz, "Methods of numerical integration",Second edition,
c Orlando, FL: Academic Press,1984. p.481-483.
c INPUT:
c lamda: wavelength
c radius: radius of sphere
c n: number of roots for Gaussian quadratures.
c OUTPUT:
c Qabs: absorption efficiency
c LOCAL:
c temp: temperature, used to obtain the complex refractive index; for VNIR,
c m is independent on temp, but for microwave, it is.
c m: refractive index
c beta: absorption coefficient
c eps: tolerance for truncation. = truncation error.
c NN: number of interfaces of reflection for ray inside a sphere.
real*8 lamda,radius,temp,Qabs
integer*4 n,i,NN
real*8 eps
complex*16 m
real*8 re,im,pi,xi
real*8 x(n),w(n)
integer*4 kind
real*8 endpts(2), bb(n)
real*8 a, b, G, absFF
integer*4 mm,k
real*8 xm,xk,dx
eps=1.0d-15
pi = dacos(-1.0d0)
a=0.0d0
b=1.0d0
xm=0.5d0*(b+a)
xk=0.5d0*(b-a)
mm=(n+1)/2
k=2*mm-n+1
Qabs=0.0d0
call GaussLeg(n,x,w)
c To make in_angle lie in 0 to pi/2, x(i) must be selected as positive.
do i=mm+1,n
dx=xk*x(i)
Qabs=Qabs+w(i)*(absFF(xm+dx,lamda,radius)
& +absFF(xm-dx,lamda,radius))
end do
Qabs=Qabs*xk
return
end
real*8 function absFF(x,lamda,radius)
c PURPOSE:
c To calculate the integrand function f(x) in the expression of absorption efficiency
c
c REFERENCE:
c Zhou, X., S. Li, and K. Stamnes, A new geometrical optics code for computing optical properties of large
c dielectric spheres, Applied Optics, 42 (21), 4295-4306, 2003.
c INPUT VARIABLES:
c x: varible, cosine of the incident angle
c lamda: wavelength in meter
c radius: radius of scatter in meter
c OUTPUT VARIABLES:
c absFF: value of the integrand function for absorption efficiency calculation. eq.(21b)
c LOCAL VARIABLES:
c j: loop varible
c beta: absorption coefficient
c xi: path length between two reflection events inside the sphere
c temp: snow temperature in K, only as an input to refice,no effect in VNIR
c m: complex refractive index
c in_angle:incident angle in radian
c TT: transmission for natural light
c RR: reflectance for natural light
c Trp: transmission for parallel-polarized light
c Trv: transmission for vertical-polarized light
c Rep: reflection for parallel-polarized light
c Rev: reflection for vertical-polarized light
c tp: complex, amplitude coefficient of transmission for parallel-polarization
c tv: complex, amplitude coefficient of transmission for vertical-polarization
c rp: complex, amplitude coefficient of reflection for parallel-polarization
c rv: complex, amplitude coefficient of reflection for vertical-polarization
c t: transmission coefficient for natural light t=sqrt(0.5(tp*2+tv**2))
c r: reflection coefficient for natural light r=sqrt(0.5(rp*2+rv**2))
c u: mcos(theta)=u+iv
c v:
real*8 x,lamda,radius
complex*16 m
real*8 beta,xi,temp,in_angle,TT, RR,Trp,Trv
complex*16 Rep,Rev,tp,tv,rp,rv,t,r,u,v
temp=270.0d0
c To obtain refractive index m
call refice( lamda, temp, m )
c To obtain absorption coefficient beta
call abs_coef(lamda, m, beta)
c To obtain a path length xi
in_angle=dacos(x)
call path_length(radius,in_angle,m,xi)
call TransReflec(in_angle,m,TT,RR,Trp,Trv,
& Rep,Rev,tp,tv,rp,rv,t,r,u,v)
absFF=2.0*TT/(1.-RR*dexp(-beta*xi))*(1.0d0-dexp(-beta*xi))*x
return
end
c ************************************************************************************
c *********To calculate the scattering efficiency Qscat and asymmetry factor**********
c ************************************************************************************
subroutine QscatAndG(lamda,radius,n,QNscat,QFscat,gN,gF)
c PURPOSE:
c To calculate scattering efficiency and asymmetry factor for both near filed (QNscat)
c and far filed (QFscat), given wavelength, radius of sphere, number of roots for
c Gaussian quadratures. Eq.(26c)+eq.(41).
c REFERENCE:
c Davis P J, P Rainowitz, "Methods of numerical integration",Second edition,
c Orlando, FL: Academic Press,1984. p.481-483.
c INPUT:
c lamda: wavelength
c radius: radius of sphere
c n: number of roots for Gaussian quadratures.
c OUTPUT:
c QFscat: scattering efficiency for far field
c QNscat: scattering efficiency for near field
c gN: asymmetry factor for near field
c gF: asymmetry factor for far field
c LOCAL:
c temp: temperature, used to obtain the complex refractive index; for VNIR,
c m is independent on temp, but for microwave, it is.
c m: refractive index
c beta: absorption coefficient
c eps: tolerance for truncation. = truncation error. 10^-10
c NN: number of interfaces of reflection for ray inside a sphere.
c theta: scattering angle
c in_angle: real, incident angle
c Trans: transmission
c Reflec: reflectance
c Trp: transmission for parallel-polarized light
c Trv: transmission for vertical-polarized light
c Rep: reflection for parallel-polarized light
c Rev: reflection for vertical-polarized light
c tp: real, amplitude coefficient of transmission for parallel-polarization
c tv: real, amplitude coefficient of transmission for vertical-polarization
c rp: real, amplitude coefficient of reflection for parallel-polarization
c rv: real, amplitude coefficient of reflection for vertical-polarization
c t: transmission coefficient for natural light t=sqrt(0.5(tp*2+tv**2))
c r: reflection coefficient for natural light r=sqrt(0.5(rp*2+rv**2))
c i: loop variable
c j: jth interface
c gNQscat: g*Qscat for near field
c gFQscat: g*Qsact for far field
c scatF1,scatF2: value of integrand function in Qscat expression
c asymNF1,asymNF2: value of integrand function in g*Qscat for near field (see Eq.(48b))
c asymFF1,asymFF2: value of integrand function in g*Qscat for far field (see Eq.(48b))
c u: mcos(theta)=u+iv
c v:
real*8 lamda,radius
real*8 QFscat,QNscat,gNQscat,gFQscat,gN,gF
real*8 Fscat,Nscat
real*8 scatfN1,scatfN2,scatfFB1,scatfFB2
real*8 asymNF1,asymNF2,asymFFB1,asymFFB2
integer*4 n
real*8 temp, theta, in_angle,Trans, Reflec,Trp,Trv
complex*16 Rep,Rev,tp,tv,rp,rv,t,r,thetat
real*8 beta,eps,re,im,pi,xi,u,v
complex*16 m
real*8 x(n),w(n),scatF,xm,xk,dx, xmp,xkp,dxp
integer*4 mm,k,i,NN
real*8 a,b,bp,xacc,x1,x2,epsilon
parameter (maxdiv=10)
real*8 D
integer*4 nbb, N2
real*8 xb1(maxdiv),xb2(maxdiv), rts(maxdiv)
nsegments=maxdiv
nb=maxdiv
pi = dacos(-1.0d0)
a=0.0d0
b=pi
bp=pi/2.0
xacc = 1.0d-15
x1=0.0d0
x2=1.0d0
epsilon = 1.0e-8
temp=270.0d0
eps=1.0d-15
c To obtain refractive index m: lamda(wavelength) and temp(temperature) are inputs.
call refice( lamda, temp, m )
c
c To obtain the abscissas x(n) and weights w(n) of the Gauss-Legendre n-point quadrature
c formula: n (number of roots for Gaussian quadratures) is input.
call GaussLeg(n,x,w)
c
xm=0.5d0*(b+a)
xk=0.5d0*(b-a)
xmp=0.5d0*(bp+a)
xkp=0.5d0*(bp-a)
mm=(n+1)/2
k=2*mm-n+1
QNscat=0.0d0
QFscat=0.0d0
Nscat=0.0d0
Fscat=0.0d0
c Product of g and Q for both near and far fields
gNQscat=0.0d0
gFQscat=0.0d0
c Calculating scattering efficiency:
call scat_efficiency(lamda,radius,m,QNscat,QFscat)
c write(*,*)'QNscat=',QNscat,'QFscat=',QFscat
c To make theta lie in 0 to pi, x(i) must be selected as positive.
do i=mm+1,n
dx=xk*x(i)
dxp=xkp*x(i)
c Calculating scattering asymmetry factor g:
call QtimesG(xm+dx,m,lamda,radius,
& scatfN1,scatfFB1,asymNF1,asymFFB1)
call QtimesG(xm-dx,m,lamda,radius,
& scatfN2,scatfFB2,asymNF2,asymFFB2)
Nscat=Nscat+w(i)*(scatfN1 +scatfN2)
gNQscat=gNQscat+w(i)*(asymNF1 +asymNF2)
call QtimesG(xmp+dxp,m,lamda,radius,
& scatfN1,scatfFB1,asymNF1,asymFFB1)
call QtimesG(xmp-dxp,m,lamda,radius,
& scatfN2,scatfFB2,asymNF2,asymFFB2)
c write(*,*)'scatfFB1=',scatfFB1
c write(*,*)'scatfFB2=',scatfFB2
Fscat=Fscat+w(i)*(scatfFB1 +scatfFB2)
gFQscat=gFQscat+w(i)*(asymFFB1 +asymFFB2)
end do
Nscat=Nscat*xk
gNQscat=gNQscat*xk
gN=gNQscat/Nscat
Fscat=Fscat*xkp
gFQscat=gFQscat*xkp
Fscat=Fscat+Nscat
gFQscat=gFQscat+gNQscat
gF=gFQscat/Fscat
return
end
subroutine scat_efficiency(lamda,radius,m,QNscat,QFscat)
c PURPOSE:
c To calculate scattering efficiency for both near filed (QNscat) and far filed
c (QFscat), given wavelength, radius of sphere, number of roots for Gaussian
c quadratures.
c REFERENCE:
c Davis P J, P Rainowitz, "Methods of numerical integration",Second edition,
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -