gomspheremain_f.htm

sphere scattering program and statment

      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  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)
          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
c   Up to now, gF is actually the asymmetry factor due only to diffraction       
       gF=gFQscat/Fscat
c   To obtain the asymmetry factor for far-field scattering (Eq.(49))
       gF = ((QNscat * gN) + gF) / (1. + QNscat)
       
       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. eq.(44) & eq.(48)
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     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)
       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.(43)
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                 Eq.(5.47)
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 between 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=g*Qscat for near field.
c      asymFF: value of the integrand function for asymmetry factor calculation.
c              asymFF=g*Qscat 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