subroutine hafermanmatrix(nopr,iout,&
   thetasun,thetaview,phisun,phiview,nreal,nim,betadeg,&
   m11,m12,m33,m44,ratiom12m11,ratiom33m11)

! Caution: seems your implementation is good only for the 
!          scattering plane (with sun vector and observation vector in same plane)

! ******************************
! See Haferman "A multi-dimensional discrete ordinates method for
! polarized radiative transfer PartI: Validation for randomly
! oriented axisymmetric particles"
! JQSRT, vol 58, No 1, p379-398, 1997

! The full 4x4 bidirectional reflection matrix is equation 17
!    m11  m12  0    0 
!    m12  m11  0    0
!    0    0    m33  -m43
!    0    0    m43  m33

! For input Stokes vector V0=(s0,s1,s2,s3)T 
!  T, do transpose to get column vector
! You get output vector V1=(soR,s1R,s2R,s3R)T
! With    V1 = MSmatrix V0

! ******************************
! Input
! thetasun     solar zenith angle (degrees)
! thetaview    sensor zenith angle (degrees)
! phisun       solar azimuth angle (degrees)
! phiview      ssensor azimuth angle (degrees)
! betadeg      surface tilt angle (degrees)
! nreal,nim    real and imaginary indices

! Output
! m11,m12      the 1,1 and 1,2 matrix elements of the reflection matrix
! m33,m44      the 3,3 and 4,4 matrix elements of the reflection matrix

! ******************************
    real :: nreal,nim
    real :: nreal2
    real :: m11,m12,m33,m44,m43
    real :: ratiom12m11

! ******************************
! Helpful for conversions
    radius=6370.0
    pi=3.14159265
    pi=atan(1.)*4.
    twopi=2.0*pi
    circum=2.00*pi*radius
! radians per degree
    convr=pi/180.0
! km per degree (110.0)
    conv=circum/360.0

! ******************************
! Asuume you can use previous info to calculate the
! reflection matrix elements

! Note that you use betadeg (tilt angle) to offset thetasun
    thetanew=thetasun-betadeg

     s1=convr*thetanew
    csbeta=cos(s1)
    cs2=csbeta*csbeta
    snbeta=sin(s1)
    sn2=sbeta*sbeta

! square of real index
! Original paper uses square of complex index
    nreal2=nreal*nreal

! eauation (18)
    a1=nreal2*csbeta
    a2=sqrt(nreal2+cs2-1.00)
    rv=(a1-a2)/(a1+a2)
    rv2=rv*rv

! eauation (19)
    b1=csbeta
    rh=(b1-a2)/(b1+a2)
    rh2=rh*rh

    m11=0.5*(rv2+rh2)
    m12=0.5*(rv2-rh2)

    m33=rv*rh
    m44=m33

! m43 is the imaginary part of rv rh*  with * being complex conjugate
! Since you just use real part here, m43=0.0
    m43=0.0

   ratiom12m11=m12/m11
   ratiom33m11=m33/m11

! ******************************
   if (nopr .eq. 1) then 
    write(iout,100) thetasun,betadeg,nreal,thetanew
    100 format(/,2x,'hafermanmatrix: thetasun,betadeg,nreal,thetanew',/,&
    2x,5(1x,f12.5))
    write(iout,102) csbeta,snbeta
    102 format(2x,'hafermanmatrix: csbeta,snbeta',/,&
    2x,2(1x,f12.5))
    write(iout,104) rv,rh
    104 format(2x,'hafermanmatrix: rv,rh',/,&
    2x,2(1x,f12.5))
    write(iout,106) rv2,rh2
    106 format(2x,'hafermanmatrix: rv2,rh2',/,&
    2x,2(1x,f12.5))
    write(iout,115) m11,m12,m33,m44
    115 format(2x,'hafermanmatrix: m11,m12,m33,m44',/,&
    2x,4(1x,f12.5))
    write(iout,120) ratiom12m11,ratiom33m11
    120 format(2x,'hafermanmatrix: ratiom12m11,ratiom33m11',/,&
    2x,2(1x,f12.5))
   endif

! ******************************
  return
  end