subroutine collett(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)

! hmmm, m12 is of opposite sign as m12 of hafermanmatris.f90
! and otarmatrix.f90

! ******************************
! See Collett "Mueller-Stokes Matrix Formulation of
! Fresnel's Equations
! American Journal of Physics, vol 39, pg 517, 1971

! The full 4x4 Mueller-Stokes matrix (MSmatrix) (equation(7)) is 
!    m11  m12  0   0 
!    m12  m11  0   0
!    0    0    m33 0
!    0    0    0   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 Mueller-Stokes reflection
!              matrix

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

! ******************************
! 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

    thetanew=thetasun-betadeg
     s1=convr*thetanew
    snthetanew=sin(s1)
    csthetanew=cos(s1)

! obtain refreacted angle from Snell's law
    snr=(1.0/nreal)*snthetanew
    rval=asin(snr)
    rangle=rval/convr

    alphaminus=thetanew-rangle
    alphaplus=thetanew+rangle

    csplus=cos(convr*alphaplus)
    csminus=cos(convr*alphaminus)

    csplus2=csplus*csplus
    csminus2=csminus*csminus

    a1=tan(convr*alphaminus)/sin(convr*alphaplus)
    factor=0.5*(a1*a1)

    m11=factor*(csminus2+csplus2)
    m12=factor*(csminus2-csplus2)

    m33=-2.0*factor*csplus*csminus
    m44=m33

! ratio of m12 to m11
    ratiom12m11=m12/m11

! ratio of m33 to m11
    ratiom33m11=m33/m11

! ******************************
   if (nopr .eq. 1) then 
    write(iout,100) thetasun,betadeg,nreal,thetanew,rangle
    100 format(/,2x,'collett: thetasun,betadeg,nreal,thetanew,rangle',/,&
    2x,5(1x,f10.4))
    write(iout,105) alphaplus,alphaminus
    105 format(2x,'collett: alphaplus,alphaminus',/,&
    2x,2(1x,f10.4))
    write(iout,110) csplus,csminus
    110 format(2x,'collett: csplus,csminus',/,&
    2x,2(1x,f10.4))
    write(iout,112) factor
    112 format(2x,'collett: factor',/,&
    2x,1(1x,f10.4))
    write(iout,115) m11,m12,m33,m44
    115 format(2x,'collett: m11,m12,m33,m44',/,&
    2x,4(1x,f10.4))
    write(iout,120) ratiom12m11,ratiom33m11
    120 format(2x,'collett: ratiom12m11,ratiom33m11',/,&
    2x,2(1x,f10.4))
   endif

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