📄 gomsphere_f.htm
字号:
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)
c METHOD : Tabular interpolation, assuming
c (1) real index is linear in log(wavelength)
c and linear in temperature
c (2) log(imag. index) is linear in log(wavelength)
c and linear in temperature
c AUTHORS OF subroutine refice( WAVMET, TEMP, ref_ice ):
c Stephen Warren, Univ. of Washington (1983)
c (sgw@cloudy.atmos.washington.edu)
c Bo-Cai Gao, JCESS, Univ. of Maryland (1995)
c (gao@imagecube.gsfc.nasa.gov)
c Warren Wiscombe, NASA Goddard (1995)
c (wiscombe@climate.gsfc.nasa.gov)
c MODIFICATIONS IN 1995 :
c Gao, Warren, and (to a small extent) Wiscombe modified the
c original Warren refice program from 1984 to change values of
c imaginary refractive index in the 0.161-0.410 and 1.445-2.50
c micron regions. The values in 0.161-0.410 were incorrect and
c the values in 1.445-2.50 were among the most uncertain in 1984.
c New measurements have made it possible to improve both regions.
c No changes were made to real refractive indices (by re-doing a
c Kramers-Kronig analysis), because the values of imaginary
c index MIM involved are so small (below 0.001) that the
c resulting changes in real index MRE would be in the third
c decimal place at most. (MIM has negligible effect on MRE
c when MIM << MRE.)
c The 0.161-0.410 micron region was changed using data provided
c by Warren, which correct his misinterpretation of Minton's
c measurements for 0.181-0.185 micron, and incorporate new
c measurements of Perovich and Govoni (1991) for 0.250-0.400
c micron. Warren (1984) correctly represented UV measurements
c of Seki et al. and visible measurements of Grenfell/Perovich,
c but he plotted Minton's measurements a factor of 2.3 too low
c because he misinterpreted base-10 as base-e. (The UV
c measurements of Dressler/Schnepp and Shibaguchi et al are also
c probably expressed as absorption coefficients on base-10;
c therefore those values also were probably plotted a factor of
c 2.3 too low in Warren's (1984) Figure 2.)
c The details of how the present imaginary index data for
c 0.161-0.410 micron is obtained are as follows. Point A in
c Warren's Figure 2 at 161 nm is joined with a straight line to
c Minton's corrected point B at 181 nm. Minton's reported
c values for 181-185 nm have been smoothed within his stated
c uncertainty. Now a smooth curve is drawn to join Minton at
c 185 nm to Perovich/Govoni (PG) at 250 nm. PG's values from
c their Table 1 show some unrealistic wiggles that are smaller
c than their error bars, so a smooth curve was fitted through
c them and values were taken from the smoothed curve at 10-nm
c intervals. PG ends at 400 nm, where Grenfell/Perovich (GP)
c starts. At 400 nm we take imaginary index=2.82E-9, the
c average of PG (2.93E-9) and GP (2.71E-9).
c The Warren (1984) values of imaginary index in the 1.445-2.50
c micron region were replaced by those of Kou et al.(1994). In
c order to remove the resulting discontinuities near 1.445 and
c 2.5 micron, the Warren values at 1.43 and 1.44 micron were
c changed to 0.9E-04 and 1.3E-04 respectively, and his values at
c 2.52, 2.55, and 2.565 micron were changed to 8.255E-04,
c 8.578E-04, and 8.739E-04, respectively. The latter change
c eliminated a small local maximum at 2.5 micron which was not
c realistic and has never been seen in spectra of snow bracketing
c that wavelength.
c REFERENCES :
c Warren, S., 1984: Optical Constants of Ice from the Ultraviolet
c to the Microwave, Appl. Opt. 23, 1206-1225
c Kou, L., D. Labrie, and P. Chylek, 1994: Refractive indices
c of water and ice in the 0.65- to 2.5-micron spectral range,
c Appl. Opt. 32, 3531-3540
c Perovich, D., and J. Govoni, 1991: Absorption Coefficients
c of Ice from 250 to 400 nm, Geophys. Res. Lett. 18, 1233-1235
c ======================================================================
IMPLICIT NONE
c .. Parameters ..
INTEGER NWL, NWLT
PARAMETER ( NWL = 574, NWLT = 62)
c ..
c .. Scalar Arguments ..
REAL*8 TEMP, WAVLEN,WAVMET
complex*16 ref_ice
c ..
c .. Local Scalars ..
CHARACTER MESSAG*40
LOGICAL PASS1
INTEGER I, L
REAL*8 FRAC, MIM, MRE, YHI, YLO
c ..
c .. Intrinsic Functions ..
INTRINSIC LOG, CMPLX, EXP
c ..
c ** Refractive index table
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -