module polarization
  interface stokes
    module procedure stokes_ce,stokes_ell
  end interface

contains
  subroutine stokes_ce(ext,eyt,qq,uu,vv)
    use const_and_precisions, only : wp_,two
    implicit none
! arguments
    complex(wp_), intent(in) :: ext,eyt
    real(wp_), intent(out) :: qq,uu,vv

    qq = abs(ext)**2 - abs(eyt)**2
    uu = two* dble(ext*dconjg(eyt))
    vv = two*dimag(ext*dconjg(eyt))
  end subroutine stokes_ce


  subroutine stokes_ell(chi,psi,qq,uu,vv)
    use const_and_precisions, only : wp_,two
    implicit none
! arguments
    real(wp_), intent(in) :: chi,psi
    real(wp_), intent(out) :: qq,uu,vv

    qq=cos(two*chi)*cos(two*psi)
    uu=cos(two*chi)*sin(two*psi)
    vv=sin(two*chi)
  end subroutine stokes_ell


  subroutine polellipse(qq,uu,vv,psi,chi)
    use const_and_precisions, only : wp_,half
    implicit none
! arguments
    real(wp_), intent(in) :: qq,uu,vv
    real(wp_), intent(out) :: psi,chi
!    real(wp_) :: ll,aa,bb,ell

!    ll  = sqrt(qq**2 + uu**2)
!    aa  = sqrt(half*(1 + ll))
!    bb  = sqrt(half*(1 - ll))
!    ell = bb/aa
    psi = half*atan2(uu,qq)
    chi = half*asin(vv)
  end subroutine polellipse

  subroutine pol_limit(anv,bv,bres,sox,ext,eyt) !,gam)
    use const_and_precisions, only : wp_,ui=>im,pi,zero,one
    implicit none
! arguments
    real(wp_), dimension(3), intent(in) :: anv,bv
    real(wp_), intent(in) :: bres,sox
    complex(wp_), intent(out) :: ext,eyt
!    real(wp_), optional, intent(out) :: gam
! local variables
    real(wp_), dimension(3) :: bnv
    real(wp_) :: anx,any,anz,an2,an,anpl2,anpl,anpr,anxy, &
       btot,yg,den,dnl,del0,ff,ff2,sngam,csgam
!
    btot = sqrt(bv(1)**2+bv(2)**2+bv(3)**2)
    bnv  = bv/btot
    yg   = btot/bres

    anx  = anv(1)
    any  = anv(2)
    anz  = anv(3)
    an2  = anx**2 + any**2 + anz**2
    an   = sqrt(an2)
    anxy = sqrt(anx**2 + any**2)

    anpl = (anv(1)*bnv(1) + anv(2)*bnv(2) + anv(3)*bnv(3))
    anpl2= anpl**2
    anpr = sqrt(an2 - anpl2)

    dnl  = one - anpl2
    del0 = sqrt(dnl**2 + 4.0_wp_*anpl2/yg**2)

    sngam =  (anz*anpl   - an2*bnv(3))/(an*anxy*anpr)
    csgam = -(any*bnv(1) - anx*bnv(2))/   (anxy*anpr)

    ff  = 0.5_wp_*yg*(dnl - sox*del0)
    ff2 = ff**2
    den = ff2 + anpl2
    if (den>zero) then
      ext  =  (ff*csgam - ui*anpl*sngam)/sqrt(den)
      eyt  = (-ff*sngam - ui*anpl*csgam)/sqrt(den)
    else  ! only for XM (sox=+1) when N//=0
      ext  = -ui*sngam
      eyt  = -ui*csgam
    end if

!    gam = atan2(sngam,csgam)/degree
  end subroutine pol_limit

  subroutine polarcold(anpl,anpr,xg,yg,sox,exf,eyif,ezf,elf,etf)
    use const_and_precisions, only : wp_,zero,one
    implicit none
! arguments
    real(wp_), intent(in) :: anpl,anpr,xg,yg,sox
    real(wp_), intent(out) :: exf,eyif,ezf,elf,etf
! local variables
    real(wp_) :: anpl2,anpr2,an2,yg2,dy2,aa,e3,qq,p

    if(xg <= zero) then
      exf  = zero
      if(sox < zero) then
        ezf  = one
        eyif = zero
      else
        ezf  = zero
        eyif = one
      end if
      elf = zero
      etf = one
    else
      anpl2 = anpl**2
      anpr2 = anpr**2
      an2   = anpl2 + anpr2

      yg2=yg**2
      aa=1.0_wp_-xg-yg2

      dy2  = one - yg2
      qq   = xg*yg/(an2*dy2 - aa)

      if (anpl == zero) then
        if(sox < zero) then
          exf  = zero
          eyif = zero
          ezf  = one
        else
          qq   = -aa/(xg*yg)
          exf  =  one/sqrt(one + qq**2)
          eyif =  qq*exf
          ezf  =  zero
        end if
      else
        e3   = one - xg
        p    = (anpr2 - e3)/(anpl*anpr) ! undef for anpr==0
        exf  =  p*ezf
        eyif = qq*exf
        ezf  = one/sqrt(one + p**2*(one + qq**2))
      end if

      elf = (anpl*ezf + anpr*exf)/sqrt(an2)
      etf = sqrt(one - elf**2)
    end if
  end subroutine polarcold

end module polarization