📄 gomspheremain_f.htm
字号:
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
<!-- saved from url=(0053)http://www.ees.nmt.edu/zhou/GOMsphere/GOMsphereMAIN.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> program testGOM
c This main program is an example to show how to call GOMsphere subroutine to calculate the
c spectral variation of single scattering albedo and asymmetry factors for both near- and far-
c field.
c Input variables:
c lamda: wavelength in meter. (For this program, lamda=0.3-4.0 micrometer)
c radius: radius of the spherical scatters ( snow grains)in meters.
c theta: scattering angle.
c 1. If only need to calculate absorption efficiency Qabs, just specify lamda and radius as inputs.
c 2. If only need to calculate the scattering efficiency Qscat and asymmetry factor g, just specify
c lamda and radius as inputs.
c 3. theta must be specified only for phase function (PN, PF) calculations. For other calculations,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 a0: single scattering albedo=scattering efficiency/(absorption+scattering) efficiency
c a0N: single scattering albedo for near-field
c a0F: single scattering albedo for far-field
c Local variable:
c n: number of roots for Gauss-Legendre n-point quadratures.
c pi: =acos(-1.0)=3.1415926
c
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.
real*8 lamda, radius
integer*4 n
real*8 Qabs,QNscat,QFscat,gN,gF,theta,PN,PF,a0
pi=dacos(-1.0d0)
open (10,file='QFNscat-size-0.5.dat', status='unknown')
open (11,file='QFNscat-size-2.0.dat', status='unknown')
open (12,file='spectral-SSalbedo&g1.0cm.dat', status='unknown')
write (10,'(a4,1x,a4,1x,a6,1x,a6,1x,a6,1x,a5)')
& 'size','Qabs','QNscat','QFscat','g4near','g4far'
write (11,'(a4,1x,a4,1x,a6,1x,a6,1x,a6,1x,a5)')
& 'size','Qabs','QNscat','QFscat','g4near','g4far'
write (12,'(a10,1x,3(a14,1x),1x,a14)')' lamda ',
& 'SS-albedo-near','SS-albedo-far', 'ASYM-fact-near',
& 'ASYM-fact-far'
c Specify a scattering angle just as an input, any value of theta will be OK of we only
c concern the quantities except the phase functions.
theta=60./180.*pi
c Task1 to calculate spectral single scattering albedo a0 and asymmetry factor g for
c lamda from 0.3-4.0um and snow grains with radius a=1cm=1.0d-2.
write (*,*) 'Working on Task 1...'
lamda = 0.3d-6
radius=1.0d-2
do while (lamda .lt.4.1d-6)
call GOMsphere(lamda,radius,theta,
& Qabs,QNscat,QFscat,gN,gF,PN,PF)
a0N=QNscat/(Qabs+QNscat)
a0F=QFscat/(Qabs+QFscat)
write (12,'(d10.4,1x,3(d14.5,1x),d14.5,)')
& lamda*1.d6,a0N,a0F,gN,gF
write (*,*) lamda*1.d6,a0N,a0F,gN,gF
lamda = lamda + 0.01d-6
end do
write (*,*) 'End of Task 1.'
c Task2 to calculate far-field scattering efficiency QFscat versus size parameter x for lamda=0.5micrometer
write (*,*) 'Working on Task 2...'
lamda=0.5d-6
do kk =0,10000,10
radius=lamda*float(kk)/2.0d0/pi
call GOMsphere(lamda,radius,theta,
& Qabs,QNscat,QFscat,gN,gF,PN,PF)
write (10,'(I10,3(1x,e10.3),2(1x,f10.5))')kk,Qabs,QNscat,
& QFscat,gN,gF
end do
write (*,*) 'End of Task 2.'
c Task3 to calculate far-field scattering efficiency QFscat versus size parameter x for lamda=2.0micrometer
write (*,*) 'Working on Task 3...'
lamda=2.0d-6
do kk =0,10000,10
radius=lamda*float(kk)/2.0d0/pi
call GOMsphere(lamda,radius,theta,
& Qabs,QNscat,QFscat,gN,gF,PN,PF)
write (11,'(I10,3(1x,e10.3),2(1x,f10.5))')kk,Qabs,QNscat,
& QFscat,gN,gF
end do
write (*,*) 'End of Task 3.'
stop
end
subroutine GOMsphere(lamda,radius,theta,
& Qabs,QNscat,QFscat,gN,gF,PN,PF)
c Purpose:
c To calculate absorption efficiency, 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 to be distributed publically along a paper to be submitted to Applied Optics.
c Any errors or comments please direct them to Xiaobing Zhou at xzhou@gi.alaska.edu or xzhou@nmt.edu
c or xzhou_2000_2001@yahoo.com. Thanks. March, 2002; May 2003.
real*8 lamda, radius
real*8 Qabs,QNscat,QFscat,gN,gF,theta,PN,PF
integer*4 n
n=20
c To calculate 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
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. Eqs.(43)-(44) and (47-49).
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.(44b))
c asymFF1,asymFF2: value of integrand function in g*Qscat for far field (see Eq.(44b))
c u: mcos(theta)=u+iv
c v:
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -