subroutine getmtmatrix(nopr,iout,&
      wmic,wavenumber,mre,mim,&
      XI,XJ,&
      SIGMA2,&
      VI1,VI2,VI3,&
      VR1,VR2,VR3,&
      m11,m12,m21,m22,&
      m11div,m12div,m21div,m22div,&
      theta86,r62mtm11,r62mtm12,ratiomtm21m11,&
      dmodmt,factor62,factor62mult,shadow87)

! ********
! Purpose of getmtmatrix.f90 is to calculate
! the reflection matrix from scratch from Appendix A of
! Mishchenko and Travis, "Satellite retrieval
! of aerosol properties over the ocean using polarization as
! well as intensity sunlight"
! JGR, vol 102, No. D14, pages 16989-17013, July 27, 1997

! ********
      integer nopr,iout
      DOUBLE PRECISION XI,XJ
      DOUBLE PRECISION SIGMA2,SIGMA
      DOUBLE PRECISION VI1, VI2, VI3, VR1, VR2, VR3

      double precision wmic,wavenumber
      double precision pi,convr

! see equations 1-3)
      double precision n0vecx,n0vecy,n0vecz
      double precision nvecx,nvecy,nvecz
      double precision n0nvecx,n0nvecy,n0nvecz
      double precision kn0nvecx,kn0nvecy,kn0nvecz

      double precision phivecx,phivecy,phivecz
      double precision thetavecx,thetavecy,thetavecz

      double precision phi0vecx,phi0vecy,phi0vecz
      double precision theta0vecx,theta0vecy,theta0vecz

      double precision vec1,dmodmt,factor62mult

! See equation 63
      double precision kdz4,kd4
      double precision nn0x,nn0y,nn0z
      double precision nn04
      double precision term1,term2,factor62
      double precision r62mtm11,r62mtm12,ratiomtm21m11

      double precision zx,zy,zz
      double precision vnorm

      double precision znx,zny,znz
      double precision pnx,pny,pnz

      double precision zn0x,zn0y,zn0z
      double precision pn0x,pn0y,pn0z

! see equations 86 and 87
      double precision theta86,shadow87

! Used in equations 84-86
      double precision mre,mim,mreal2
      double precision sum1,csval,csval2
      double precision a1val,a2val,a3val
      double precision a1,a2,a3,a4,a51,a52
      double precision b1,b2
      double precision c1,c2
      double precision rperp,rpar

! See equations 80-83 for dot products
      double precision dotphi0n,dotphin0,dottheta0n,dotthetan0
! Equations 80-83
      double precision fthetatheta,fthetaphi,fphitheta,fphiphi
      double precision fthetatheta2,fthetaphi2,fphitheta2,fphiphi2
! The M matrix (equation 64 -79(
      double precision m11,m12,m21,m22
      double precision m11div,m12div,m21div,m22div

! ********
      pi=atan(1.0d0)*4.0d0
      convr=pi/180.0d0

! ********
! incident n0 vector
    n0vecx=VI1
    n0vecy=VI2
    n0vecz=VI3

! reflection n vector
    nvecx=VR1
    nvecy=VR2
    nvecz=VR3

! vector difference n0 - n
    n0nvecx=n0vecx-nvecx
    n0nvecy=n0vecy-nvecy
    n0nvecz=n0vecz-nvecz

! equation (63) kd vector
    kn0nvecx=wavenumber*n0nvecx
    kn0nvecy=wavenumber*n0nvecy
    kn0nvecz=wavenumber*n0nvecz

! kdz to the 4th power
    kdz4=kn0nvecz**4.0d0

! absolute value of kd to the 4th power
    vnorm=dsqrt( (kn0nvecx*kn0nvecx) + (kn0nvecy*kn0nvecy) +(kn0nvecz*kn0nvecz) )
    kd4=vnorm**4.0d0

! cross product n x n0
    nn0x=(nvecy*n0vecz)-(nvecz*n0vecy)
    nn0y=(nvecz*n0vecx)-(nvecx*n0vecz)
    nn0z=(nvecx*n0vecy)-(nvecy*n0vecx)
    vec1=(nn0x*nn0x) + (nn0y*nn0y) +(nn0z*nn0z)
! This is used to rescale matrix M to values
! which can be compared to breon.f90 matrix elements
    dmodmt=vec1*vec1
    vnorm=dsqrt(vec1)
! length of  n x n0 raised to the 4th power
    nn04=vnorm**4.0d0

! *****************
! see equation (62)
    term1=kd4/(4.0d0*XI*XJ*nn04*kdz4*2.0d0*SIGMA2)

! The exponential term in equation 62
    a1val=( (kn0nvecx*kn0nvecx) + (kn0nvecy*kn0nvecy) )
    a2val=2.0d0*(kn0nvecz*kn0nvecz)*SIGMA2
    a3val=a1val/a2val
    term2=dexp(-a3val)

! This is the scalar prefactor to matrix M in equation 62 of papaer MT
! factor62 multiples the M matrix in equation 62 to get the R ocean 
! reflection matrix
    factor62=term1*term2

! Useful to compare to Breon rspdiv value (should be similar) 
! You calculate m11div=m11/dmodmt so that you can compare
! m11breon and m11divnew
! To retain equation 62 value, need to then multiply factor62 by dmodmt
    factor62mult=factor62*dmodmt

! ********
! z axis unit vector
    zx=0.d0
    zy=0.d0
    zz=1.0d0

! ********
! For reflective ray n

! cross product z x n
    znx=(zy*nvecz)-(zz*nvecy)
    zny=(zz*nvecx)-(zx*nvecz)
    znz=(zx*nvecy)-(zy*nvecx)
    vnorm=dsqrt( (znx*znx) + (zny*zny) + (znz*znz) )
! unit vec phi
    phivecx=znx/vnorm
    phivecy=zny/vnorm
    phivecz=znz/vnorm

! cross product phi x n
    pnx=(phivecy*nvecz)-(phivecz*nvecy)
    pny=(phivecz*nvecx)-(phivecx*nvecz)
    pnz=(phivecx*nvecy)-(phivecy*nvecx)
    vnorm=dsqrt( (pnx*pnx) + (pny*pny) + (pnz*pnz) )
! unit vec phi
    thetavecx=pnx/vnorm
    thetavecy=pny/vnorm
    thetavecz=pnz/vnorm

! ********
! For incident ray n0

! cross product z x n0
    zn0x=(zy*n0vecz)-(zz*n0vecy)
    zn0y=(zz*n0vecx)-(zx*n0vecz)
    zn0z=(zx*n0vecy)-(zy*n0vecx)
    vnorm=dsqrt( (zn0x*zn0x) + (zn0y*zn0y) + (zn0z*zn0z) )
! unit vec phi0
    phi0vecx=zn0x/vnorm
    phi0vecy=zn0y/vnorm
    phi0vecz=zn0z/vnorm

! cross product phi0 x n0
    pn0x=(phi0vecy*n0vecz)-(phi0vecz*n0vecy)
    pn0y=(phi0vecz*n0vecx)-(phi0vecx*n0vecz)
    pn0z=(phi0vecx*n0vecy)-(phi0vecy*n0vecx)
    vnorm=dsqrt( (pn0x*pn0x) + (pn0y*pn0y) + (pn0z*pn0z) )
! unit vec theta0
    theta0vecx=pn0x/vnorm
    theta0vecy=pn0y/vnorm
    theta0vecz=pn0z/vnorm

! ********
   if (nopr .eq. 1) then
      write(iout,fmt=20) VI1,VI2,VI3,VR1,VR2,VR3,&
      n0vecx,n0vecy,n0vecz,nvecx,nvecy,nvecz,&
      n0nvecx,n0nvecy,n0nvecz,&
      kn0nvecx,kn0nvecy,kn0nvecz,wmic,wavenumber,&
      znx,zny,znz,pnx,pny,pnz,&
      phivecx,phivecy,phivecz,thetavecx,thetavecy,thetavecz,&
      zn0x,zn0y,zn0z,pn0x,pn0y,pn0z,&
      phi0vecx,phi0vecy,phi0vecz,theta0vecx,theta0vecy,theta0vecz,&
      a1val,a2val,a3val,term1,term2,factor62,dmodmt,factor62mult
   20 format(/,2x,'getmtmatrix: VI1,VI2,VI3 ',3(1x,f12.5),/, &
      2x,'getmtmatrix: VR1,VR2,VR3 ',3(1x,f12.5),/, &
      2x,'getmtmatrix: n0vecx,n0vecy,n0vecz ',3(1x,f12.5),/, &
      2x,'getmtmatrix: nvecx,nvecy,nvecz ',3(1x,f12.5),/, &
      2x,'getmtmatrix: n0nvecx,n0nvecy,n0nvecz ',3(1x,f12.5),/, &
      2x,'getmtmatrix: kn0nvecx,kn0nvecy,kn0nvecz ',3(1x,f12.5),/, &
      2x,'getmtmatrix: wmic,wavenumber ',2(1x,f12.5),/,&
      2x,'getmtmatrix:   for reflective vector n ',/,&
      2x,'getmtmatrix: znx,zny,znz ',3(1x,f12.5),/, &
      2x,'getmtmatrix: pnx,pny,pnz ',3(1x,f12.5),/, &
      2x,'getmtmatrix: phivecx,phivecy,phivecz ',3(1x,f12.5),/, &
      2x,'getmtmatrix: thetavecx,thetavecy,thetavecz ',3(1x,f12.5),/,&
      2x,'getmtmatrix:   for incident vector n ',/,&
      2x,'getmtmatrix: zn0x,zn0y,zn0z ',3(1x,f12.5),/, &
      2x,'getmtmatrix: pn0x,pn0y,pn0z ',3(1x,f12.5),/, &
      2x,'getmtmatrix: phi0vecx,phi0vecy,phi0vecz ',3(1x,f12.5),/, &
      2x,'getmtmatrix: theta0vecx,theta0vecy,theta0vecz ',3(1x,f12.5),/,&
      2x,'getmtmatrix: a1val,a2val,a3val ',3(1x,f12.5),/,&
      2x,'getmtmatrix: term1,term2,factor62 ',3(1x,f12.5),/,&
      2x,'getmtmatrix: dmodmt,factor62mult ',3(1x,f12.5))
   endif

! ********
! Equation (86)
! use + sign since equation is in terms of (n-no)
      sum1=1.0d0*( (n0vecx*n0nvecx) + (n0vecy*n0nvecy) + (n0vecz*n0nvecz) )
      vnorm=dsqrt( (n0nvecx*n0nvecx) + (n0nvecy*n0nvecy) + (n0nvecz*n0nvecz) )
      csval=sum1/vnorm
      csval2=csval*csval

! The angle theta in degrees from equation 86 of Appendix A of MT
      theta86=(acos(csval))/convr

! Square of the real index
      mreal2=mre*mre

! Equation (84)
      a1=csval
      a2=dsqrt( mreal2 - 1.0d0 + csval2)
      rperp=(a1-a2)/(a1+a2)

! Equation (85)
      b1=mreal2*csval
      b2=a2
      rpar=(b1-b2)/(b1+b2)

! ********
   if (nopr .eq. 1) then
      write(iout,fmt=30) mre,mreal2,csval,csval2,&
      a1,a2,b1,b2,rpar,rperp
   30 format(/,2x,'getmtmatrix: mre,mreal2 ',2(1x,f12.5),/,&
      2x,'getmtmatrix: csval,csval2 ',2(1x,f12.5),/, &
      2x,'getmtmatrix: a1,a2 ',2(1x,f12.5),/, &
      2x,'getmtmatrix: b1,b2 ',2(1x,f12.5),/, &
      2x,'getmtmatrix: rpar,rperp ',2(1x,f12.5))
   endif

! ********
! used in Equations 80-83 
      dotphi0n= (phi0vecx*nvecx) + (phi0vecy*nvecy) + (phi0vecz*nvecz)
      dotphin0= (phivecx*n0vecx) + (phivecy*n0vecy) + (phivecz*n0vecz)

      dottheta0n= (theta0vecx*nvecx) + (theta0vecy*nvecy) + (theta0vecz*nvecz)
      dotthetan0= (thetavecx*n0vecx) + (thetavecy*n0vecy) + (thetavecz*n0vecz)

! ********
! Equation 80 
      fthetatheta= (dotphi0n*dotphin0*rperp) + (dottheta0n*dotthetan0*rpar)

! Equation 81
      fthetaphi= (-dottheta0n*dotphin0*rperp) + (dotphi0n*dotthetan0*rpar)

! Equation 82
      fphitheta= (-dotphi0n*dotthetan0*rperp) + (dottheta0n*dotphin0*rpar)

! Equation 83
      fphiphi= (dottheta0n*dotthetan0*rperp) + (dotphi0n*dotphin0*rpar)

! Takes squares of the terms
      fthetatheta2 = fthetatheta*fthetatheta
      fthetaphi2 = fthetaphi*fthetaphi
      fphitheta2 = fphitheta*fphitheta
      fphiphi2 = fphiphi*fphiphi

! ********
! Equation (64)
      m11=0.5d0*( fthetatheta2 + fthetaphi2 + fphitheta2 + fphiphi2 )

! Equation (65)
      m12=0.5d0*( fthetatheta2 - fthetaphi2 + fphitheta2 - fphiphi2 )

! Equation (66)
      m21=0.5d0*( fthetatheta2 + fthetaphi2 - fphitheta2 - fphiphi2 )

! Equation (67)
      m22=0.5d0*( fthetatheta2 - fthetaphi2 - fphitheta2 + fphiphi2 )

! ********
! ratio of m21 and m11
      ratiom21m11=m21/m11

      ratiomtm21m11=1.0*ratiom21m11

! ********
! For diagnostics
! Can compare m11div with Hafermanmatrix.f90 m11
     m11div=m11/dmodmt
     m12div=m12/dmodmt
     m21div=m21/dmodmt
     m22div=m22/dmodmt

! ********
   if (nopr .eq. 1) then
      write(iout,fmt=40) dotphi0n,dotphin0,dottheta0n,dotthetan0,&
      fthetatheta,fthetaphi,fphitheta,fphiphi,&
      fthetatheta2,fthetaphi2,fphitheta2,fphiphi2,&
      m11,m12,m21,m22,&
      dmodmt,&
      m11div,m12div,m21div,m22div,&
      ratiom21m11,ratiomtm21m11
   40 format(/,2x,'getmtmatrix: dotphi0n,dotphin0,dottheta0n,dotthetan0',/,&
      2x,4(1x,f12.5),/,&
      2x,'getmtmatrix: fthetatheta,fthetaphi,fphitheta,fphiphi',/,&
      2x,4(1x,f12.5),/,&
      2x,'getmtmatrix: fthetatheta2,fthetaphi2,fphitheta2a,fphipha2i',/,&
      2x,4(1x,f12.5),/,&
      2x,'getmtmatrix: m11,m12,m21,m22',/,&
      2x,4(1x,f12.5),/,&
      2x,'getmtmatrix: dmodmt',/,&
      2x,1(1x,f12.5),/,&
      2x,'getmtmatrix: m11div,m12div,m21div,m22div',/,&
      2x,4(1x,f12.5),/,&
      2x,'getmtmatrix: ratiom21m11,ratiomtm21m11',/,&
      2x,1(1x,f12.5),1(1x,f12.5))
   endif

! ********
! Calculate the shadowing term
     SIGMA=dsqrt(SIGMA2)

! Equation 88
! For XI
     b1=1.0d0-(XI*XI)
     a1=XI/dsqrt(2.0d0*b1)
     a2=a1*a1
     a3=(SIGMA/XI)*dsqrt((2.0d0*b1)/pi)
     a4=a2/SIGMA2
     a51=a1/SIGMA
     c1=0.5*( (a3*exp(-a4)) -erfc(a51) )

! Equation 88
! For XI
     b1=1.0d0-(XJ*XJ)
     a1=XJ/dsqrt(2.0d0*b1)
     a2=a1*a1
     a3=(SIGMA/XJ)*dsqrt((2.0d0*b1)/pi)
     a4=a2/SIGMA2
     a52=a1/SIGMA
     c2=0.5*( (a3*exp(-a4)) -erfc(a52) )

! Equation 87
     shadow87=1.0d0/(1.0d0+c1+c2)

! ********
   if (nopr .eq. 1) then
      write(iout,fmt=50) a51,a52,c1,c2,shadow87
   50 format(/,2x,'getmtmatrix: a51,a52,c1,c2,shadow87',/,&
      2x,5(1x,f12.5))
   endif

! ********
! See equation 62
! The equation R11 term
! m11 is the (1,1) term of the matrix M
    r62mtm11=factor62*m11*shadow87

! The equation R12 term
! m11 is the (1,2) term of the matrix M
    r62mtm12=factor62*m21*shadow87

! ********
   if (nopr .eq. 1) then
      write(iout,fmt=60) r62mtm11,r62mtm12
   60 format(/,2x,'getmtmatrix: r62mtm11,r62mtm12 ',/,&
      2x,5(1x,f12.5))
   endif

! ********
      return
      end subroutine getmtmatrix