gomsphere_f.htm

sphere scattering program and statment

c     Orlando, FL: Academic Press,1984. p.481-483.
c INPUT:
c     lamda: wavelength
c     radius: radius of sphere
c     m: refractive index
c OUTPUT:
c     QFscat: scattering efficiency for far field
c     QNscat: scattering efficiency for near field
c LOCAL:
c     a: real, lower limit of integration
c     b: real, higher limit of integration
c     n: number of roots for Gaussian quadratures.
c     temp: temperature, used to obtain the complex refractive index; for VNIR,
c           m is independent on temp, but for microwave, it is.
c     beta: scatorption coefficient
c     eps: tolerance for truncation. = truncation error. 10^-10
c     NN: number of interfaces of reflection for ray inside a sphere.
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

      real*8 a,b,lamda,radius
      complex*16 m
      real*8 QFscat,QNscat
      integer*4 n
      real*8 pi
      real*8 x(20),w(20),scatFNF,xm,xk,dx
      integer*4 mm,k,i
      n=20
      pi = dacos(-1.0d0)
      a=0.0d0
      b=1.0d0

c  To obtain the abscissas and weights of the Gauss-Legendre n-point quadrature
c     formula.
       call GaussLeg(n,x,w)
c       
       xm=0.5d0*(b+a)
       xk=0.5d0*(b-a)
       mm=(n+1)/2
       k=2*mm-n+1
       QNscat=0.0d0
       QFscat=0.0d0
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)
          
c  Calculating scattering efficiency:
          QNscat=QNscat+w(i)*(scatFNF(xm+dx,lamda,radius)
     &                   +scatFNF(xm-dx,lamda,radius))
       end do
       
       QNscat=QNscat*xk
       QFscat=QNscat+1.0d0
       
       return
       end

       real*8 function scatFNF(x,lamda,radius)
c  PURPOSE:
c      To calculate the integrand function f(x) in the expression of scattering efficiency
C      for near field.
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, also the abscissas of the Gauss-Legendre
c         n-point quadrature.
c      lamda: wavelength of incident wave in unit of meter
c      radius: radius of snow grain as a sphere
c  OUTPUT VARIABLES:
c      scatFNF: value of the integrand function for near filed scattering efficiency 
c               calculation. See eq.(23)
c  LOCAL VARIABLES:
c      NN: the number of interfaces when truncation is carried out with truncation 
c         error = 10^-10
c      m: complex, refractive index
c      beta: absorption coefficients
c      xi: real, path-length bwtween two reflecting interfaces
c      temp: snow temperature in K, only as an input to refice,no effect in VNIR
c      epspjp: real, amplitude of scattered parallel-polarized electric field
c      epspjv: real, amplitude of scattered vertically-polarized electric field
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      j: loop variable, jth interface
c      
      real*8 x,lamda,radius
      integer*4 NN
      complex*16 m
      real*8 temp,beta,thetai,xi,eps
      complex*16 Rep,Rev,tp,tv,rp,rv,t,r
      real*8 Trans, Reflec, Trp,Trv,u,v
      integer*4 j
      scatFNF=0.0d0
      eps=1.0d-15
      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
      thetai = dacos(x)
      call path_length(radius,thetai,m,xi)

       call TransReflec(thetai,m,Trans,Reflec,Trp,Trv,
     &                                     Rep,Rev,tp,tv,rp,rv,t,r,u,v)
c  To obtain truncation number NN
      call sumNumber3(beta,xi,t,Reflec,Trans,eps,NN)

       do j=2,NN
          scatFNF=scatFNF+(Trp*Trp*(cdabs(Rep))**(j-2)
     &                   +Trv*Trv*(cdabs(Rev))**(j-2))
     &                   *dexp(-beta*xi*dfloat(j-1))
       end do
       scatFNF=scatFNF+ 2.*Reflec
       scatFNF = scatFNF * x
       
       return
       end

       subroutine QtimesG(theta,m,lamda,radius,
     &                            scatfN,scatfFB,asymNF,asymFFB)
c  PURPOSE:
c      To calculate the integrand function f(x) in the expression of scattering
c      efficiency (eq.(26b)) and the integrand function f(x)in the expression of 
c      asymmetry factor for near field scattering and far field. 
c      For Far-field, Qsca^far=Qsca^near+1
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      theta: varible, scattering angle in radians
c      m: complex, refractive index
c      beta: absorption coefficients
c      xi: real, path-length bwtween two interfaces
c      NN: the number of interfaces when truncation is carried out with truncation 
c         error = 10^-10
c  OUTPUT VARIABLES:
c      scatfN: value of the integrand function for near-field scattering efficiency calculation
c      asymNF: value of the integrand function for asymmetry factor calculation.
c              asymNF=gQscat for near field.
c      asymFF: value of the integrand function for asymmetry factor calculation.
c              asymFF=gQscat for far field.

c  LOCAL VARIABLES:
c      epspjp: real, amplitude of scattered parallel-polarized electric field
c      epspjv: real, amplitude of scattered vertically-polarized electric field
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      j: loop variable, jth interface
c      thetai: incident angle
c      x: size parameter = 2*pi*radius/wavelength
c      
       real*8 scatfN,scatfFB,asymNF,asymFFB
       real*8 theta,beta,xi,lamda,radius
       complex*16 m,thetat
       integer*4 NN
       integer*4 nb,nsegments
       real*8 thetai,Trans, Reflec,Trp,Trv
       complex*16 Rep,Rev,tp,tv,rp,rv,t,r
       complex*16 epspjp,epspjv
       parameter (maxdiv=10)
       parameter (Nmax=1000)
       real*8 re,D
       integer*4 nbb,i,j,mm, N2
       real*8 xb1(maxdiv),xb2(maxdiv), rts(maxdiv)
       real*8 pi,u,v,a,b,xacc,x1,x2,epsilon,x
       nsegments=maxdiv
       nb=maxdiv

       a=0.0d0
       b=1.0d0
       xacc = 1.0d-15
       x1=0.0d0
       x2=1.0d0
       pi=dacos(-1.0d0)
       epsilon = 1.0e-16
       x=2.0*pi*radius/lamda
       re=dreal(m)
       j=0
       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)
              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 bear-field scattering
         asymNF=scatfN*dcos(theta)
c refraction part of integrand function of far-filed 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