📄 gomspheremain_f.htm
字号:
mm=0
scatfN=0.0d0
do while (j.le.Nmax)
j=j+1
c To find nb and rts(incident angles)
call roots(j,m,theta,x1,x2,xacc, nsegments, nb, rts)
if (nb.eq.0) go to 10
do 9 i = 1, nb
thetai=2.0*datan(rts(i))
c subroutine TransReflec(..) is in Trans_Ref.f. To calculate tp,tv,rp,rv,t,r
call TransReflec(thetai,m,Trans,Reflec,Trp,Trv,
& Rep,Rev,tp,tv,rp,rv,t,r,u,v)
c subroutine (..) ScatFieldAmp(...) is in Trans_Ref.f. To calculate epspjp,epspjv
c To obtain path-length xi
call path_length(radius,thetai,m,xi)
c To obtain absorption coefficient beta
call abs_coef(lamda, m, beta)
call ScatFieldAmp(thetai,tp,tv,rp,rv,beta,xi,j,u,v,
& epspjp,epspjv)
c To calculate the geometrical divergence factor Eq. (5.24)
D=dsin(2.*thetai)/(4.0*dsin(theta))
& /dabs(dfloat(j-1)*dcos(thetai)/dsqrt(u*u+v*v)-1.0)
scatfN=scatfN
& +(cdabs(epspjp)**2+cdabs(epspjv)**2)*D
mm=mm+1
9 continue
10 continue
end do
cccc scatfN,scatfFB,asymNF,asymFFB
c See eq.(26b)-(26c)
scatfN=scatfN*dsin(theta)
c refraction part of integrand function of far-filed scattering efficiency
scatfFB=2.0*BessJ1(x*dsin(theta))*BessJ1(x*dsin(theta))
& /dsin(theta)
c integrand function for near-field scattering
asymNF=scatfN*dcos(theta)
c refraction part of integrand function of far-field asymmetry facotr (eq.50)
asymFFB=scatfFB*dcos(theta)
return
end
subroutine ScatFieldAmp(thetai,tp,tv,rp,rv,beta,xi,j,u,v,
& epspjp,epspjv)
c PURPOSE:
c To evaluate amplitude of scattered electric field in term of incident amplitude.
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 thetai: real, the incident angle (in radians)
c m: complex, refractive index
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 beta: absorption coefficient
c xi: path length between two reflection events inside the sphere
c j: integer, jth interface of a ray inteacting with the sphere interface
c OUTPUT VARIABLES:
c epspjp: real, amplitude of scattered parallel-polarized electric field
c epspjv: real, amplitude of scattered vertically-polarized electric field
c LOCAL VARIABLES:
c thetat: real, refraction angle (in radians)
real*8 thetai,beta,xi,u,v
complex*16 tp,tv,rp,rv,thetat
complex*16 epspjp,epspjv
integer*4 j
if (j.eq.1) then
epspjp = rp
epspjv = rv
else
epspjp = complex(u,v)/dcos(thetai)*tp*tp*(-rp)**(j-2)
& *dexp(-beta*xi*(j-1)/2.0d0)
epspjv = complex(u,v)/dcos(thetai)*tv*tv*(-rv)**(j-2)
& *dexp(-beta*xi*(j-1)/2.0d0)
end if
return
end
real*8 function BessJ1(x)
c Purpose:
c Returns the value of the first order Bessel function J1(x) for any real x
c Reference:
c Press W H , S A Teukolsky, W T Vetterling and B P Flannery,"Numerical
c recipes in C: the art of scientific computing", Second edition, Cambridge
c University Press,1992. p.233.
real*8 ax,z
real*8 x,y,ans1,ans2,xx
ax=dabs(x)
if(ax.lt.8.0) then
y=x*x
ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1
& +y*(-2972611.439+y*(15704.48260+y*(-30.16036606))))))
ans2=144725228442.0+y*(2300535178.0+y*(18583304.74
& +y*(99447.43394+y*(376.9991397+y*1.0))))
BessJ1=ans1/ans2
else
z=8.0/ax
y=z*z
xx=ax-2.356194491
ans1=1.0+y*(0.183105d-2+y*(-0.3516396496d-4
& +y*(0.2457520174d-5+y*(-0.240337019d-6))))
ans2=0.04687499995+y*(-0.2002690873d-3
& +y*(0.8449199096d-5+y*(-0.88228987d-6+y*0.105787412d-6)))
BessJ1=dsqrt(0.636619772/ax)*(dcos(xx)*ans1-z*dsin(xx)*ans2)
if(x.lt.0.0d0) BessJ1=-BessJ1
endif
return
end
c *************************************************************************
c *********To calculate the phase function for a scattering angle**********
c *************************************************************************
subroutine phasefunction(lamda,radius,theta,QNscat,QFscat,
& Nphase,Fphase)
c Purpose:
c To obtain phase function for given wavelength and radius of a sphere
c Input:
c lamda: wavelength, in meters
c radius: radius of the scatter, in meters
c theta: scattering angle
c QFscat: scattering efficiency for far field used in scat_efficiency(...)in Q_scat2.f
c QNscat: scattering efficiency for near field used in scat_efficiency(...)in Q_scat2.f
c Output:
c Nphase: phase function for near field scattering
c Fphase: phase function for far field scattering
c Local:
c a: real, lower limit of integration, used in scat_efficiency(...)in Q_scat.f
c b: real, higher limit of integration, used in scat_efficiency(...)in Q_scat.f
c TEMP: = temperature (K) ( for WAVLEN.GT.167 microns only )
c (range: 213.16 to 272.16) used in refice(...) in ref_ice.f
c theta: outgoing (scattering) angle
c xacc: accuracy that is used to find the root.Used in roots(...) in Pphase.f
c ......in subroutine refice( WAVMET, TEMP, ref_ice ) in ref_ice.f
c . epsilon: real, tolerance for cutoff in the summation series.
c . N: integer, maximum number of summation of for absorption efficiency
c . r: real, amplitude coefficient of reflection
c . t: real, amplitude coefficient of transmission for unpolarized wave
c ......in subroutine refice( WAVMET, TEMP, ref_ice ) in ref_ice.f
c
c ......in subroutine roots(j,m,theta,x1,x2,xacc, nsegments, nb, rts)
c . j: jth interface.
c . m: complex refractive index.
c . theta: outgoing angle (scattering angle).
c . x1,x2: limits that all the roots are to be found.
c . xacc: accuracy that is used to find the root.
c . nb: as input, it is the maximum number of roots sought
c . nsegments: number that the inteval [x1,x2] is equally subdivided in order to bracket = maximum
c . possible roots
c . nb: as output, it is the number of different roots that are found.
c . rts: the different roots found that lies within [x1,x2]
c ......in subroutine roots(j,m,theta,x1,x2,xacc, nsegments, nb, rts)
c
c ......in subroutine TransReflec(thetai,m,Trans,Reflec,Trp,Trv,
c . & Rep,Rev,tp,tv,rp,rv,t,r) in Trans_Ref.f
c . thetai: incident angle
c . m: complex refractive index
c . Trans: transmission for natural light
c . Reflec: 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: 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 ......in subroutine TransReflec(thetai,m,Trans,Reflec,Trp,Trv,
c . & Rep,Rev,tp,tv,rp,rv,t,r) in Trans_Ref.f
c
c ......in subroutine ScatFieldAmp(thetai,m,tp,tv,rp,rv,beta,xi,j,
c . & epspjp,epspjv) in Q_scat.f
c . thetai: real, the incident angle (in radians)
c . m: complex, refractive index
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 . beta: absorption coefficient
c . xi: path length between two reflection events inside the sphere
c . j: integer, jth interface of a ray inteacting with the sphere interface
c . epspjp: real, amplitude of scattered parallel-polarized electric field
c . epspjv: real, amplitude of scattered vertically-polarized electric field
c ......in subroutine ScatFieldAmp(thetai,m,tp,tv,rp,rv,beta,xi,j,
c . & epspjp,epspjv) in Q_scat.f
c re: real part of complex refractive index
c D: divergence factor
c x: size parameter. = 2*pi*radius/wavelength
c PQN: P*Q for near field scattering
c PQF: P*Q for far field scattering
real*8 lamda,radius
real*8 theta,Nphase,Fphase,PQN,PQF
real *8 a,b
integer*4 N,j
real*8 QFscat,QNscat,TEMP,xacc,epsilon
real*8 x1,x2
integer*4 nb,nsegments,maxdiv,Nmax
real*8 thetai,Trans,Reflec,Trp,Trv
complex*16 Rep,Rev,tp,tv,rp,rv,t,r,thetat
real*8 beta,xi,pi,x,u,v,BessJ1
complex*16 m, epspjp,epspjv
parameter (maxdiv=10)
c parameter (Nmax=1000)
parameter (Nmax=100)
real*8 re,D
integer*4 nbb,i,mm, N2
real*8 xb1(maxdiv),xb2(maxdiv), rts(maxdiv)
nsegments=maxdiv
nb=maxdiv
xacc = 1.0d-15
x1=0.0d0
x2=1.0d0
pi=dacos(-1.0d0)
epsilon = 1.0e-8
c To obtain complex refractive index (m). subroutine refice( ... ) is in ref_ice.f
TEMP=270.0
call refice( lamda, TEMP, m )
re=dreal(m)
Nphase=0.0d0
Fphase=0.0d0
PQN=0.0d0
PQF=0.0d0
mm=0
j=0
do while (j.le.Nmax)
j=j+1
call roots(j,m,theta,x1,x2,xacc, nsegments, nb, rts)
if (nb.eq.0) go to 10
do 9 i = 1, nb
thetai=2.0*datan(rts(i))
c To obtain absorption coefficient
call abs_coef(lamda, m, beta)
c To obtain path length between two consecutive incident interfaces
call path_length(radius,thetai,m,xi)
c To obtain interface number
call sumNumber3(beta,xi,t,Reflec,Trans,epsilon,N2)
c subroutine TransReflec(..) is in Trans_Ref.f. To calculate tp,tv,rp,rv,t,r
call TransReflec(thetai,m,Trans,Reflec,Trp,Trv,
& Rep,Rev,tp,tv,rp,rv,t,r,u,v)
c subroutine (..) ScatFieldAmp(...) is in Trans_Ref.f. To calculate epspjp,epspjv
call ScatFieldAmp(thetai,tp,tv,rp,rv,beta,xi,j,u,v,
& epspjp,epspjv)
D=dsin(2.*thetai)/(4.0*dsin(theta))
& /dabs(dfloat(j-1)*dcos(thetai)/dsqrt(u*u+v*v)-1.0)
PQN=PQN+2.*(cdabs(epspjp)**2+cdabs(epspjv)**2)*D
mm=mm+1
9 continue
10 continue
end do
x=2.*pi*radius/lamda
if (theta.lt.pi/2.) then
PQF=PQN+4.*BessJ1(x*dsin(theta))*BessJ1(x*dsin(theta))
& /dsin(theta)/dsin(theta)
else
PQF=PQN
endif
Nphase=PQN/QNscat
Fphase=PQF/QFscat
return
end
c ******************************************************************************
c *********To calculate the complex refractive index for a wavelength **********
c ******************************************************************************
subroutine refice( WAVMET, TEMP, ref_ice )
c Calculates complex refractive index of Ice 1H for wavelengths
c between 45 nm and 8.6 m. For wavelengths above 167 microns,
c temperature dependence is included for temperatures between
c 213 and 272K. Mainly intended for applications in Earth ice
c clouds and snow, not other planets or interstellar space;
c the temperature dependence or crystalline form of ice may be
c incorrect for these latter applications.
c I N P U T :
c WAVMET = wavelength (meters)
c WAVLEN = wavelength (microns)
c (range: 0.0443 to 8.600E+06)
c TEMP = temperature (K) ( for WAVLEN.GT.167 only )
c (range: 213.16 to 272.16)
c O U T P U T : ref_ice = complex refractive index
c ( with positive imaginary part )
c (WARNING: input out of range will print a warning message and
c return ref_ice=(0.,0.) in order not to unnecessarily
c halt the calling program; the calling program should
c test the real part of ref_ice to catch these errors)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -