!XXXXXXXX0XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

! BEGIN polarlib.f

!   ----------------------------------------------------------------
!   |                                                              |
!   |                      polarlib                                |
!   |                                                              |
!   | Bibliotheque de calcul de la polarisation d'une onde plane.  |
!   | Extrait des bibliotheques developpees pour le satellite GEOS |
!   | au CRPE en 1977-1980. Conversion succinte des .f en .f90     |
!   |                                                              |
!   |           P. Robert, CNRS/CETP, Novembre 2001                |
!   |                                                              |
!   ----------------------------------------------------------------
!
!XXXXXXXX0XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      subroutine cpolar(sx,sy,sz,ptot,a,b,rkx,rky,rkz,ax,ay,az,phase)

!     ---------------------------------------------------------------+--
!     calcul de l'ellipse de polarisation d'une onde plane, a partir
!     de ses coefficients de Fourrier sur 3 axes orthogonaux.
!     input : complex sx,sy,sz
!     output: a,b        = axes de l'ellipse
!             rkx,rky,rkz= normale au plan de l'ellipse
!             ax ,ay ,az = direction du grand axe
!             phase      = phase absolue du champ le long de l'ellipse
!
!     P. Robert, CETP, Novembre 2001
!     ---------------------------------------------------------------+--

      complex sx,sy,sz

      real aii(3),fii(3),GXI(3),GZI(3)


! *** calcul des amplitudes et phase sur chaque composante

      call CALAIFI(sx,sy,sz,aii,fii)

! *** puissance totale

      ptot= 0.5*(aii(1)**2 +aii(2)**2 + aii(3)**2)

! *** calcul des parametres de l'ellipse

      call ELLIPSE(aii,fii,phase,a,b,EXC,GXI,GZI)

      rkx=GZI(1)
      rky=GZI(2)
      rkz=GZI(3)

      ax=GXI(1)
      ay=GXI(2)
      az=GXI(3)

      return
      end

!XXXXXXXX0XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      SUBROUTINE CALAIFI(SX,SY,SZ,AII,FII)

      DIMENSION AII(3),FII(3)
      COMPLEX SX,SY,SZ

!     ******************************************************************
!     AUTEUR   : P. ROBERT CRPE-CNET, 1980  (revu Novembre 2001)
!     CATEGORIE: DEPOUILLEMENT SPECIFIQUE UBF/GEOS
!     OBJET    : CALCUL DES COEFFICIENTS AI ET FII DU PROGRAMME POLAR
!     ******************************************************************


      AII(1)= 2.*CABS(SX)
      AII(2)= 2.*CABS(SY)
      AII(3)= 2.*CABS(SZ)

      RR= REAL(SX)
      RI=AIMAG(SX)
      FII(1)=0.
      IF(ABS(RR).GT.1.E-30.OR.ABS(RI).GT.1.E-30) FII(1)=ATAN2(RI,RR)

      RR= REAL(SY)
      RI=AIMAG(SY)
      FII(2)=0.
      IF(ABS(RR).GT.1.E-30.OR.ABS(RI).GT.1.E-30) FII(2)=ATAN2(RI,RR)

      RR= REAL(SZ)
      RI=AIMAG(SZ)
      FII(3)=0.
      IF(ABS(RR).GT.1.E-30.OR.ABS(RI).GT.1.E-30) FII(3)=ATAN2(RI,RR)

      RETURN
      END

!XXXXXXXX0XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

      SUBROUTINE ELLIPSE(AII,FII,PHA,AA,BB,EXC,GXI,GZI)

      DIMENSION AII(3),FII(3),GXI(3),GZI(3)
      double precision SN, CN, dblePHA, dbleAA, dbleBB

!     ******************************************************************
!     AUTEUR   : P. ROBERT CRPE-CNET, 1980
!     CATEGORIE: TRAITEMENT DU SIGNAL
!     OBJET    : CALCUL D UNE ELLIPSE A PARTIR DE SES 3 PROJECTIONS
!     ******************************************************************

!      CALCUL DE L"ELLIPSE DE POLARISATION A PARTIR DE SES 3 PROJECTIONS
!      DONNEES PAR LE MODULE AII ET LA PHASE FII DE CHAQUE COMPOSANTE


! *   CALCUL DE LA PHASE ABSOLUE

      SN=0.d0
      CN=0.d0
      dblePHA=0.d0

      DO 76 J=1,3
         SN=SN+(dble(AII(J))**2)*SIN(2.d0*dble(FII(J)))
         CN=CN+(dble(AII(J))**2)*COS(2.d0*dble(FII(J)))
   76 CONTINUE
      IF(DABS(SN).GT.1.D-30.AND.DABS(CN).GT.1.D-30) &
           dblePHA=DATAN2(SN,CN)/2.D0

! *   CALCUL DU MODULE DU GRAND AXE, DU PETIT AXE ET EXENTRICITE

      dbleAA=1.D-30
      dbleBB=1.D-30

      DO 78 J=1,3
         dbleAA=dbleAA+(dble(AII(J))*DCOS(dblePHA-dble(FII(J))))**2
         dbleBB=dbleBB+(dble(AII(J))*DSIN(dblePHA-dble(FII(J))))**2
   78 CONTINUE

      dbleAA=DSQRT(dbleAA)
      dbleBB=DSQRT(dbleBB)
      EXC=sngl(dbleBB/dbleAA)

! *   CALCUL DE LA DIRECTION DES AXES DE L"ELLIPSE

      DO 82 J=1,3
         GXI(J)=sngl(dble(AII(J))*DCOS(dblePHA-dble(FII(J)))/dbleAA)
   82 CONTINUE
!
      GZI(1)=sngl(dble(AII(2))*dble(AII(3))* &
             DSIN(dble(FII(2))-dble(FII(3)))/(dbleAA*dbleBB))
      GZI(2)=sngl(dble(AII(3))*dble(AII(1))* &
             DSIN(dble(FII(3))-dble(FII(1)))/(dbleAA*dbleBB))
      GZI(3)=sngl(dble(AII(1))*dble(AII(2))* &
             DSIN(dble(FII(1))-dble(FII(2)))/(dbleAA*dbleBB))
     
! *   CALCUL DU GRAND ET PETIT AXE, ET DE LA PHASE ABSOLUE

      PHA=sngl(dblePHA)
      AA=sngl(dbleAA)
      BB=sngl(dbleBB)

      RETURN
      END

!XXXXXXXX0XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX0XX

! END polarlib.f90

