      PROGRAM RNGL
      IMPLICIT COMPLEX (C)
      REAL AL(2),OL,PHI,PS,DE,DS,AN,PI,PSI,DLT,PHIR,S,SI2,ANI,EP,
     +F0,F1,G0,G1,G2,FK,AL2,AL4,AK2,AK4,BK4,DK4,T1,T2,RS,RS2,RP
      DIMENSION CZ(2)
      DATA OL,PHI,PS,DE,DS,AN,EP/632.8,70.,39.44,78.24,.001,1.5,100./
      DATA POL,BIR,ROT,ANA,REF/90.,90.,0.,-45.,0./
      DATA CN2/(3.866,-0.028)/
      DATA H0P,H2P,H2D,H4P,H4D/2.819,2.266,92.776,0.417,45.034/
      PI=3.141592654/180.
C FORTRAN 77 
C Calcul des harmoniques d'un ellipsometre a birefringent tournant :
C Psi et Delta pour une couche mince transparente sur substrat
C Indice de substrat a partir de Psi et Delta
C Epaisseur d'une couche mince a indice fixe sur substrat
C Epaisseur et indice d'une couche mince sur substrat
C Coefficients de reflexion Rp et Rs
C Harmoniques H0, H2 et H4 pour une couche mince donnee sur substrat
C Psi et Delta en fonction des harmoniques H2 et H4
C auteur : aime.vareille@wanadoo.fr 
C logiciel libre licence GPL (cf. http://www.aful.org)
C                                                    1:18:07 21/4/1988
      K=1
1     WRITE(6,200)K
200   FORMAT(1H0,' 1: Psi et Delta.'/' 2: Indice d un substrat.'/
     +' 3: Epaisseur a indice fixe.'/' 4: Epaisseur et indice.'/
     +' 5: Rp, Rs.'/' 6: H0, H2, H4.'/' 7: Psi, Delta(H0,H2,H4).'/
     +1X,I2,' ? ')
      READ(5,*)K
        IF(K.EQ.1)THEN
        WRITE(6,100)OL
100   FORMAT(1H+,' Longueur d onde : ',F8.2,' nm ? ')
        READ(5,*)OL
        WRITE(6,101)PHI
101   FORMAT(1H+,' Angle d incidence : ',F7.3,' degres ? ')
        READ (5,*)PHI
        WRITE(6,102)CN2
102   FORMAT(1H+,' Indice du substrat : ',2(F9.5,1X),' ? ')
        READ(5,*)CN2
        WRITE(6,107)AN
107   FORMAT(1H+,' Indice de la couche mince : ',F9.5,' ? ')
        READ(5,*)AN
        WRITE(6,201)EP
201   FORMAT(1H+,' Epaisseur de la couche mince : ',F9.2,' nm ? ')
        READ(5,*)EP
        CALL RNG1(OL,PHI,CN2,K,PS,DE,AN,EP,PSI,DLT,P)
        PS=PSI
        DE=DLT 
        WRITE(6,108)OL,AN,EP,P,PS,DE
108   FORMAT(1H+,' ',F8.1,' nm, Indice : ',F9.5,' Epaisseur : ',F8.2,
     +' nm'/
     +' Periode : ',F8.2,' nm'/
     +' Psi, Delta : ',F9.3,1X,F9.3,' degres'/)
        ENDIF
        IF(K.EQ.2)THEN
        WRITE(6,100)OL
C100   FORMAT(1H+,' Longueur d onde : ',F8.2,' nm ? ')
        READ(5,*)OL
        WRITE(6,101)PHI
C101   FORMAT(1H+,' Angle d incidence : ',F7.3,' degres ? ')
        READ (5,*)PHI
        WRITE(6,103)PS
103   FORMAT(1H+,' Psi : ',F8.3,' degres ? ')
        READ(5,*)PS
        WRITE(6,104)DE
104   FORMAT(1H+,' Delta : ',F8.3,' degres ? ')
        READ(5,*)DE
        PHIR=PHI*PI
        SI2=SIN(PHIR)**2
        CRO=TAN(PS*PI)*EXP(CMPLX(0.,PI*DE))
        CN2=CSQRT(TAN(PHIR)**2*(1.-4.*SI2*CRO/((1.+CRO)**2)))
        WRITE(6,202)CN2
202   FORMAT(1H+,' Indice du substrat : ',2(F9.5,1X),' . '/)
        ENDIF
        IF(K.EQ.3)THEN
        WRITE(6,100)OL
C100   FORMAT(1H+,' Longueur d onde : ',F8.2,' nm ? ')
        READ(5,*)OL
        WRITE(6,101)PHI
C101   FORMAT(1H+,' Angle d incidence : ',F7.3,' degres ? ')
        READ (5,*)PHI
        WRITE(6,102)CN2
C102   FORMAT(1H+,' Indice du substrat : ',2(F9.5,1X),' ? ')
        READ(5,*)CN2
        WRITE(6,107)AN
C107   FORMAT(1H+,' Indice de la couche mince : ',F9.5,' ? ')
        READ(5,*)AN
        WRITE(6,103)PS
C103   FORMAT(1H+,' Psi : ',F8.3,' degres ? ')
        READ(5,*)PS
        WRITE(6,104)DE
C104   FORMAT(1H+,' Delta : ',F8.3,' degres ? ')
        READ(5,*)DE
        CALL RNG1(OL,PHI,CN2,K,PS,DE,AN,EP,PSI,DLT,P)
        WRITE(6,106)OL,AN,EP,P,PS,DE,PSI,DLT
106   FORMAT(1H+,' ',F9.2,' nm, Indice : ',F9.5,' Epaisseur : ',F9.3,
     +' nm'/
     +' Periode : ',F8.2,' nm'/
     +' Psi, Delta avant : ',F9.3,1X,F9.3,' degres'/
     +' Psi, Delta apres : ',E14.7,1X,E14.7,' degres'/)
        ENDIF
        IF(K.EQ.4)THEN
        WRITE(6,100)OL
C100   FORMAT(1H+,' Longueur d onde : ',F9.3,' nm ? ')
        READ(5,*)OL
        WRITE(6,101)PHI
C101   FORMAT(1H+,' Angle d incidence : ',F7.3,' degres ? ')
        READ (5,*)PHI
        WRITE(6,102)CN2
C102   FORMAT(1H+,' Indice du substrat : ',2(F9.5,1X),' ? ')
        READ(5,*)CN2
        WRITE(6,107)AN
C107   FORMAT(1H+,' Indice de la couche mince : ',F9.5,' ? ')
        READ(5,*)AN
        WRITE(6,103)PS
C103   FORMAT(1H+,' Psi : ',F8.3,' degres ? ')
        READ(5,*)PS
        WRITE(6,104)DE
C104   FORMAT(1H+,' Delta : ',F8.3,' degres ? ')
        READ(5,*)DE
        CALL RNG1(OL,PHI,CN2,K,PS,DE,AN,EP,PSI,DLT,P)
        WRITE(6,106)OL,AN,EP,P,PS,DE,PSI,DLT
C106   FORMAT(1H+,' ',F9.2,' nm, Indice : ',F7.3,' Epaisseur : ',F8.2,' nm'/
C     +' Periode : ',F8.2,' nm'/
C     +' Psi, Delta avant : ',F8.2,1X,F8.2,' degres'/
C     +' Psi, Delta apres : ',E14.7,1X,E14.7,' degres')
        ENDIF
        IF(K.EQ.5)THEN
        WRITE(6,100)OL
C100   FORMAT(1H+,' Longueur d onde : ',F8.2,' nm ? ')
        READ(5,*)OL
        WRITE(6,101)PHI
C101   FORMAT(1H+,' Angle d incidence : ',F7.3,' degres ? ')
        READ (5,*)PHI
        WRITE(6,102)CN2
C102   FORMAT(1H+,' Indice du substrat : ',2(F9.5,1X),' ? ')
        READ(5,*)CN2
        WRITE(6,107)AN
C107   FORMAT(1H+,' Indice de la couche mince : ',F9.5,' ? ')
        READ(5,*)AN
        WRITE(6,201)EP
C201   FORMAT(1H+,' Epaisseur de la couche mince : ',F9.2,' nm ? ')
        READ(5,*)EP
      PHIR=PHI*PI
      CS0=COS(PHIR)
      CP0=1/CS0
      SI2=SIN(PHIR)**2
      CS2=CSQRT(CN2*CN2-SI2)
      CP2=CN2*CN2/CS2
      S=SQRT(AN*AN-SI2)
      ANI=0.
      P=OL/(2.*S)
        IF(EP.LT.0.)EP=EP+P
      CAN=CMPLX(AN,ANI)
      CS1=SQRT(CAN*CAN-SI2)
      CP1=CAN*CAN/CS1
      CB=EXP(CMPLX(0.,-4.*3.141592654*EP*REAL(CS1)/OL))
      CR1P=CF(CP1,CP0)
      CR1S=CF(CS0,CS1)
      CR2P=CF(CP2,CP1)
      CR2S=CF(CS1,CS2)
      CH=CG(CR1P,CR2P*CB)
      CI=CG(CR1S,CR2S*CB)
      CRO=CH/CI
      PS=ATAN(ABS(CRO))/PI
      DE=ATAN2(AIMAG(CRO),REAL(CRO))/PI
      RP=ABS(CH)**2
      RS=ABS(CI)**2
        WRITE(6,208)OL,AN,EP,P,PS,DE,CH,RP,CI,RS
208   FORMAT(1H+,' ',F9.2,' nm, Indice : ',F9.5,' Epaisseur : ',F9.3,
     +' nm'/
     +' Periode : ',F9.3,' nm'/
     +' Psi, Delta : ',F9.3,1X,F9.3,' degres'/
     +' rp : ',E14.7,1X,E14.7,'; Rp : ',F8.6/
     +' rs : ',E14.7,1X,E14.7,'; Rs : ',F8.6/)
        ENDIF
        IF(K.EQ.6)THEN
        WRITE(6,100)OL
C100   FORMAT(' Longueur d onde : ',F8.2,' nm ? ')
        READ(5,*)OL
        WRITE(6,101)PHI
C101   FORMAT(' Angle d incidence : ',F7.3,' degres ? ')
        READ (5,*)PHI
        WRITE(6,102)CN2
C102   FORMAT(' Indice du substrat : ',2(F9.5,1X),' ? ')
        READ(5,*)CN2
        WRITE(6,107)AN
C107   FORMAT(' Indice de la couche mince : ',F9.5,' ? ')
        READ(5,*)AN
        WRITE(6,201)EP
C201   FORMAT(' Epaisseur de la couche mince : ',F9.2,' nm ? ')
        READ(5,*)EP
        WRITE(6,400)POL
400   FORMAT(1H+,' Angle du polariseur : ',F8.3,' degres ? ')
        READ (5,*)POL
        WRITE(6,402)BIR
402   FORMAT(1H+,' Birefringence de Q : ',F8.3,' degres ? ')
        READ (5,*)BIR
        WRITE(6,403)ROT
403   FORMAT(1H+,' Pouvoir rotatoire de Q : ',F8.3,' degres ? ')
        READ (5,*)ROT
        WRITE(6,404)REF
404   FORMAT(1H+,' Angle du top zero : ',F8.3,' degres ? ')
        READ (5,*)REF
        WRITE(6,401)ANA
401   FORMAT(1H+,' Angle de l analyseur : ',F8.3,' degres ? ')
        READ (5,*)ANA
      PHIR=PHI*PI
      POLR=POL*PI
      BIRR=BIR*PI
      ROTR=ROT*PI
      REFR=REF*PI
      ANAR=ANA*PI
      CS0=COS(PHIR)
      CP0=1/CS0
      SI2=SIN(PHIR)**2
      CS2=CSQRT(CN2*CN2-SI2)
      CP2=CN2*CN2/CS2
      S=SQRT(AN*AN-SI2)
      ANI=0.
      P=OL/(2.*S)
        IF(EP.LT.0.)EP=EP+P
      CAN=CMPLX(AN,ANI)
      CS1=SQRT(CAN*CAN-SI2)
      CP1=CAN*CAN/CS1
      CB=EXP(CMPLX(0.,-4.*3.141592654*EP*REAL(CS1)/OL))
      CR1P=CF(CP1,CP0)
      CR1S=CF(CS0,CS1)
      CR2P=CF(CP2,CP1)
      CR2S=CF(CS1,CS2)
      CH=CG(CR1P,CR2P*CB)
      CI=CG(CR1S,CR2S*CB)
      CRO=CH/CI
      PS=ATAN(ABS(CRO))/PI
      DE=ATAN2(AIMAG(CRO),REAL(CRO))/PI
      RP=ABS(CH)**2
      RS=ABS(CI)**2
      RS2=RS/2
      T1=TAN(PS*PI)
      T2=T1**2
      H0=RS2*(T2+1+COS(2*ANAR)*(T2-1)+(1+COS(BIRR))/2*
     +(COS(2*(POLR+ROTR))*(T2-1+COS(2*ANAR)*(T2+1))+2*
     +SIN(2*(POLR+ROTR))*SIN(2*ANAR)*T1*COS(DE*PI)))
      H2C=RS*SIN(2*ANAR)*T1*SIN(DE*PI)*SIN(BIRR)*SIN(2*(REFR-ROTR-POLR))
      H2S=RS*SIN(2*ANAR)*T1*SIN(DE*PI)*SIN(BIRR)*COS(2*(REFR-ROTR-POLR))
      H4C=RS2*(1-COS(BIRR))/2*((T2-1+COS(2*ANAR)*(T2+1))*
     +COS(4*REFR-2*(POLR+ROTR))+2*SIN(2*ANAR)*T1*COS(DE*PI)*
     +SIN(4*REFR-2*(POLR+ROTR)))
      H4S=RS2*(1-COS(BIRR))/2*((T2-1+COS(2*ANAR)*(T2+1))*
     +(-SIN(4*REFR-2*(POLR+ROTR)))+2*SIN(2*ANAR)*T1*COS(DE*PI)*
     +COS(4*REFR-2*(POLR+ROTR)))
      HA2=SQRT(H2C**2+H2S**2)
      HP2=ATAN2(H2S,H2C)
      HA4=SQRT(H4C**2+H4S**2)
      HP4=ATAN2(H4S,H4C)
C CALCUL DE H4/H2
      AL2=2*(REFR-POLR-ROTR)
      AL4=2*(2*REFR-POLR-ROTR)
      AK2=RS*SIN(2*ANAR)*SIN(BIRR)
      AK4=RS2*(1+COS(2*ANAR))*(1-COS(BIRR))/2
      BK4=RS*SIN(2*ANAR)*(1-COS(BIRR))/2
      DK4=RS2*(COS(2*ANAR)-1)*(1-COS(BIRR))/2
      TG=SQRT(ABS((DK4+SIN(AL4)*H4S-COS(AL4)*H4C)/AK4))
      PSI=ATAN(TG)/PI
       SDE1=H2S/AK2/TG/COS(AL2)
       SDE2=H2C/AK2/TG/SIN(AL2)
        IF(ABS(SIN(AL2)).LT.0.05)THEN
        SDE=SIGN(SQRT(H2S**2+H2C**2)/AK2/TG,SDE1)
        ELSE
        SDE=SIGN(SQRT(H2S**2+H2C**2)/AK2/TG,SDE2)
        ENDIF
      ODE1=(H4S+(AK4*TG**2+DK4)*SIN(AL4))/BK4/TG/COS(AL4)
      ODE2=(H4C-(AK4*TG**2+DK4)*COS(AL4))/BK4/TG/SIN(AL4)
       IF(ABS(COS(AL4)).LT.0.05)THEN
       ODE=SIGN(SQRT((H4S+(AK4*TG**2+DK4)*SIN(AL4))**2+
     +(H4C-(AK4*TG**2+DK4)*COS(AL4))**2)/BK4/TG,ODE2)
       ELSE
       ODE=SIGN(SQRT((H4S+(AK4*TG**2+DK4)*SIN(AL4))**2+
     +(H4C-(AK4*TG**2+DK4)*COS(AL4))**2)/BK4/TG,ODE1)
       ENDIF
      DLT=ATAN2(SDE,ODE)/PI
      WRITE(6,338)POL,BIR,ROT,REF,ANA
338   FORMAT(1H+,' POL,BIR,ROT,REF,ANA =',5(1X,F8.3),' degres')
      WRITE(6,308)OL,AN,EP,P,CH,RP,CI,RS,H0,H2C,H2S,HA2,HP2/PI,
     +H4C,H4S,HA4,HP4/PI,PS,DE,PSI,DLT
308   FORMAT(1H+,' ',F9.2,' nm, Indice : ',F9.5,' Epaisseur : ',F9.3,
     +' nm'/
     +' Periode : ',F8.2,' nm'/
     +' rp : ',E14.7,1X,E14.7,'; Rp : ',F8.6/
     +' rs : ',E14.7,1X,E14.7,'; Rs : ',F8.6/
     +' H0 : ',F8.6/
     +' H2 C, S, A, P :',3(1X,F8.6),1X,F8.3/
     +' H4 C, S, A, P :',3(1X,F8.6),1X,F8.3/
     +' Psi, Delta avant : ',F9.3,1X,F9.3,' degres'/
     +' Psi, Delta apres : ',E14.7,1X,E14.7,' degres'/)
        ENDIF
        IF(K.EQ.7)THEN
        WRITE(6,100)OL
C100   FORMAT(' Longueur d onde : ',F8.2,' nm ? ')
        READ(5,*)OL
        WRITE(6,101)PHI
C101   FORMAT(' Angle d incidence : ',F7.3,' degres ? ')
        READ (5,*)PHI
        WRITE(6,102)CN2
C102   FORMAT(' Indice du substrat : ',2(F9.5,1X),' ? ')
        READ(5,*)CN2
C        WRITE(6,107)AN
C107   FORMAT(' Indice de la couche mince : ',F9.5,' ? ')
C        READ(5,*)AN
C        WRITE(6,201)EP
C201   FORMAT(' Epaisseur de la couche mince : ',F9.2,' nm ? ')
C        READ(5,*)EP
        WRITE(6,400)POL
C400   FORMAT(' Angle du polariseur : ',F7.3,' degres ? ')
        READ (5,*)POL
        WRITE(6,402)BIR
C402   FORMAT(' Birefringence de Q : ',F7.3,' degres ? ')
        READ (5,*)BIR
        WRITE(6,403)ROT
C403   FORMAT(' Pouvoir rotatoire de Q : ',F7.3,' degres ? ')
        READ (5,*)ROT
        WRITE(6,404)REF
C404   FORMAT(' Angle du top zero : ',F7.3,' degres ? ')
        READ (5,*)REF
        WRITE(6,401)ANA
C401   FORMAT(' Angle de l analyseur : ',F7.3,' degres ? ')
        READ (5,*)ANA
        WRITE(6,500)H0P,H2P,H2D,H4P,H4D
500   FORMAT(1H+,' H0,HA2,HP2,HA4,HP4 : ',F8.6,2(1X,F8.6,1X,F8.3),
     +' degres ? ')
        READ (5,*)H0P,H2P,H2D,H4P,H4D
      H2R=H2D*PI
      H4R=H4D*PI
      PHIR=PHI*PI
      POLR=POL*PI
      BIRR=BIR*PI
      ROTR=ROT*PI
      REFR=REF*PI
      ANAR=ANA*PI
C CALCUL DE PSI ET DELTA A PARTIR DE H2 ET H4
      RS=1.
      RS2=RS/2.
      AL2=2*(REFR-POLR-ROTR)
      AL4=2*(2*REFR-POLR-ROTR)
      H2S=H2P*SIN(H2R)
      H2C=H2P*COS(H2R)
      H4S=H4P*SIN(H4R)
      H4C=H4P*COS(H4R)
      AK2=RS*SIN(2*ANAR)*SIN(BIRR)
      AK4=RS2*(1+COS(2*ANAR))*(1-COS(BIRR))/2
      BK4=RS*SIN(2*ANAR)*(1-COS(BIRR))/2
      DK4=RS2*(COS(2*ANAR)-1)*(1-COS(BIRR))/2
      F1=AK4/(COS(AL4)*H4C-SIN(AL4)*H4S)
      F0=DK4/(COS(AL4)*H4C-SIN(AL4)*H4S)
      G2=((H4C*F1-AK4*COS(AL4))**2+(H4S*F1+AK4*SIN(AL4))**2)/(BK4**2)+
     +(H2C**2+H2S**2)*((F1/AK2)**2)
      G1=2.*(F0*F1*(H2C**2+H2S**2)/(AK2**2)+((H4C*F1-AK4*COS(AL4))*
     +(H4C*F0-DK4*COS(AL4))+(H4S*F1+AK4*SIN(AL4))*
     +(H4S*F0+DK4*SIN(AL4)))/(BK4**2))-1.
      G0=((H4C*F0-DK4*COS(AL4))**2+(H4S*F0+DK4*SIN(AL4))**2)/(BK4**2)+
     +(H2C**2+H2S**2)*((F0/AK2)**2)
      DTA=G1*G1-4.*G0*G2
        IF(DTA.LT.0.)THEN
        WRITE(6,700)DTA
700   FORMAT(1H+,' Discriminant negatif: ',E14.7)
        ELSE
        TG2M=(-G1-SQRT(DTA))/2./G2
        TG2P=(-G1+SQRT(DTA))/2./G2
C        TGM=SQRT(TG2M)
C        TGP=SQRT(TG2P)
        FKM=F1*TG2M+F0
C        PSM=ATAN(TGM)/PI
        FKP=F1*TG2P+F0
C        PSP=ATAN(TGP)/PI
         IF(FKP.GT.0.)THEN
         TG2=TG2P
         FK=FKP
         ELSE
         TG2=TG2M
         FK=FKM
         ENDIF
        TG=SQRT(TG2)
        PS=ATAN(TG)/PI
         IF(ABS(SIN(AL2)).LT.0.05)THEN
         SD1=FK*H2S/AK2/TG/COS(AL2)
         SDE=SIGN(FK*SQRT(H2S**2+H2C**2)/AK2/TG,SD1)
C         SD1M=FKM*H2S/AK2/TGM/COS(AL2)
C         SDEM=SIGN(FKM*SQRT(H2S**2+H2C**2)/AK2/TGM,SD1M)
         ELSE
         SD2=FK*H2C/AK2/TG/SIN(AL2)
         SDE=SIGN(FK*SQRT(H2S**2+H2C**2)/AK2/TG,SD2)
C         SD2M=FKM*H2C/AK2/TGM/SIN(AL2)
C         SDEM=SIGN(FKM*SQRT(H2S**2+H2C**2)/AK2/TGM,SD2M)
         ENDIF
C         IF(ABS(COS(AL4)).LT.0.05)THEN
C         OD2P=(FKP*H4C-(AK4*TG2P+DK4)*COS(AL4))/BK4/TGP/SIN(AL4)
C         ODEP=SIGN(SQRT((FKP*H4S+(AK4*TG2P+DK4)*SIN(AL4))**2+
C     +(FKP*H4C-(AK4*TG2P+DK4)*COS(AL4))**2)/BK4/TGP,OD2P)
C         OD2M=(FKM*H4C-(AK4*TG2M+DK4)*COS(AL4))/BK4/TGM/SIN(AL4)
C         ODEM=SIGN(SQRT((FKM*H4S+(AK4*TG2M+DK4)*SIN(AL4))**2+
C     +(FKM*H4C-(AK4*TG2M+DK4)*COS(AL4))**2)/BK4/TGM,OD2M)
C         ELSE
C         OD1P=(FKP*H4S+(AK4*TG2P+DK4)*SIN(AL4))/BK4/TGP/COS(AL4)
C         ODEP=SIGN(SQRT((FKP*H4S+(AK4*TG2P+DK4)*SIN(AL4))**2+
C     +(FKP*H4C-(AK4*TG2P+DK4)*COS(AL4))**2)/BK4/TGP,OD1P)
C         OD1M=(FKM*H4S+(AK4*TG2M+DK4)*SIN(AL4))/BK4/TGM/COS(AL4)
C         ODEM=SIGN(SQRT((FKM*H4S+(AK4*TG2M+DK4)*SIN(AL4))**2+
C     +(FKM*H4C-(AK4*TG2M+DK4)*COS(AL4))**2)/BK4/TGM,OD1M)
C       ENDIF
        ODE=FK*(H4C*SIN(AL4)+H4S*COS(AL4))/BK4/TG
C        ODEM=FKM*(H4C*SIN(AL4)+H4S*COS(AL4))/BK4/TGM
C         IF(ABS(SDEM*SDEM+ODEM*ODEM-1.).LT.
C     +ABS(SDEP*SDEP+ODEP*ODEP-1.))THEN
C         DE=ATAN2(SDEM,ODEM)/PI
C         PS=PSM
C         ELSE
C         DE=ATAN2(SDEP,ODEP)/PI
C         PS=PSP
C         ENDIF
         DE=ATAN2(SDE,ODE)/PI
        ENDIF
C      WRITE(6,871)ABS(SDEP*SDEP+ODEP*ODEP-1.),
C     +ABS(SDEM*SDEM+ODEM*ODEM-1.),FKP,FKM
C871   FORMAT(' RES+,RES-, FK+,FK- =',2(2X,E14.7,1X,E14.7))
C      WRITE(6,870)PSP,ATAN2(SDEP,ODEP)/PI,PSM,ATAN2(SDEM,ODEM)/PI
C870     FORMAT(' PS+,DE+, PS-,DE- =',2(2X,E14.7,1X,E14.7))
C        WRITE(6,781)FK
C781     FORMAT(' FK = ',E14.7)
C        WRITE(6,782)AK2
C782     FORMAT(' AK2 = ',E14.7)
C        WRITE(6,783)TG
C783     FORMAT(' TG = ',E14.7)
C        WRITE(6,777)SDE1,FK*SQRT(H2S**2+H2C**2)/AK2/TG,SDE
C777     FORMAT(' SDE1,SDEP,SDE = ',3(1X,E14.7))
C        WRITE(6,778)SDE2,FK*SQRT(H2S**2+H2C**2)/AK2/TG,SDE
C778     FORMAT(' SDE2,SDEP,SDE = ',3(1X,E14.7))
C        WRITE(6,780)ODE1,SQRT((FK*H4S+(AK4*TG2+DK4)*SIN(AL4))**2+
C     +(FK*H4C-(AK4*TG2+DK4)*COS(AL4))**2)/BK4/TG,ODE
C780     FORMAT(' ODE1,ODEP,ODE = ',3(1X,E14.7))
CC      SDE=FK*SQRT(H2S**2+H2C**2)/AK2/TG
CC      ODE=SQRT((FK*H4S+(AK4*TG2+DK4)*SIN(AL4))**2+
CC     +(FK*H4C-(AK4*TG2+DK4)*COS(AL4))**2)/BK4/TG
C CALCUL DE L'INDICE ET DE L'EPAISSEUR
C        K1=4
C        AN=1.5
C CALCUL DE L'EPAISSEUR
        K1=3
C        AN=1.5
        CALL RNG1(OL,PHI,CN2,K1,PS,DE,AN,EP,PSI,DLT,P)
C CALCUL DE H0, H2 ET H4 A PARTIR DE PSI ET DELTA
      CS0=COS(PHIR)
      CP0=1/CS0
      SI2=SIN(PHIR)**2
      CS2=CSQRT(CN2*CN2-SI2)
      CP2=CN2*CN2/CS2
      S=SQRT(AN*AN-SI2)
      ANI=0.
      P=OL/(2.*S)
        IF(EP.LT.0.)EP=EP+P
      CAN=CMPLX(AN,ANI)
      CS1=SQRT(CAN*CAN-SI2)
      CP1=CAN*CAN/CS1
      CB=EXP(CMPLX(0.,-4.*3.141592654*EP*REAL(CS1)/OL))
      CR1P=CF(CP1,CP0)
      CR1S=CF(CS0,CS1)
      CR2P=CF(CP2,CP1)
      CR2S=CF(CS1,CS2)
      CH=CG(CR1P,CR2P*CB)
      CI=CG(CR1S,CR2S*CB)
      CRO=CH/CI
      PS1=ATAN(ABS(CRO))/PI
      DE1=ATAN2(AIMAG(CRO),REAL(CRO))/PI
      RP=ABS(CH)**2
      RS=ABS(CI)**2
      RS2=RS/2
      T1=TAN(PS1*PI)
      T2=T1**2
      H0=RS2*(T2+1+COS(2*ANAR)*(T2-1)+(1+COS(BIRR))/2*
     +(COS(2*(POLR+ROTR))*(T2-1+COS(2*ANAR)*(T2+1))+2*
     +SIN(2*(POLR+ROTR))*SIN(2*ANAR)*T1*COS(DE1*PI)))
      H2C=RS*SIN(2*ANAR)*T1*SIN(DE1*PI)*SIN(BIRR)*
     +SIN(2*(REFR-ROTR-POLR))
      H2S=RS*SIN(2*ANAR)*T1*SIN(DE1*PI)*SIN(BIRR)*
     +COS(2*(REFR-ROTR-POLR))
      H4C=RS2*(1-COS(BIRR))/2*((T2-1+COS(2*ANAR)*(T2+1))*
     +COS(4*REFR-2*(POLR+ROTR))+2*SIN(2*ANAR)*T1*COS(DE1*PI)*
     +SIN(4*REFR-2*(POLR+ROTR)))
      H4S=RS2*(1-COS(BIRR))/2*((T2-1+COS(2*ANAR)*(T2+1))*
     +(-SIN(4*REFR-2*(POLR+ROTR)))+2*SIN(2*ANAR)*T1*COS(DE1*PI)*
     +COS(4*REFR-2*(POLR+ROTR)))
      HA2=SQRT(H2C**2+H2S**2)
      HP2=ATAN2(H2S,H2C)
      HA4=SQRT(H4C**2+H4S**2)
      HP4=ATAN2(H4S,H4C)
      WRITE(6,309)H0P,H2P,H2D,H4P,H4D,
     +AN,EP,P,CH,RP,CI,RS,H0,H2C,H2S,HA2,HP2/PI,
     +H4C,H4S,HA4,HP4/PI,PS,DE,PSI,DLT
309   FORMAT(1H+,' H0 : ',F8.6,
     +' H2 A, P :',1X,F8.6,1X,F8.3,
     +' H4 A, P :',1X,F8.6,1X,F8.3/
     +' Indice : ',F9.5,' Epaisseur : ',F9.3,' nm',
     +' Periode : ',F9.3,' nm'/
     +' rp : ',E14.7,1X,E14.7,'; Rp : ',F8.6/
     +' rs : ',E14.7,1X,E14.7,'; Rs : ',F8.6/
     +' H0 : ',F8.6/
     +' H2 C, S, A, P :',3(1X,F8.6),1X,F8.3/
     +' H4 C, S, A, P :',3(1X,F8.6),1X,F8.3/
     +' Psi, Delta avant : ',F9.3,1X,F9.3,' degres'/
     +' Psi, Delta apres : ',E14.7,1X,E14.7,' degres'/)
        ENDIF
      GO TO 1
      END
      SUBROUTINE RNG1(OL,PHI,CN2,K,PS,DE,AN,EP,PSI,DLT,P)
      IMPLICIT COMPLEX (C)
      REAL AL(2),OL,PHI,PS,DE,EP,P,DS,AN,PI,PSI,DLT,PHIR,S,SI2,ANI
      DIMENSION CZ(2)
      DS=0.001
      PI=3.141592654/180.
      PHIR=PHI*PI
      CS0=COS(PHIR)
      CP0=1/CS0
      SI2=SIN(PHIR)**2
      CS2=CSQRT(CN2*CN2-SI2)
      CP2=CN2*CN2/CS2
      S=SQRT(AN*AN-SI2)
        IF(K.EQ.1) GO TO 1
      CRO=TAN(PS*PI)*EXP(CMPLX(0.,PI*DE))
        IF(K.EQ.3) THEN
	CS1=S
        CP1=CS1+SI2/CS1
        CALL CRES(CS0,CS1,CS2,CP0,CP1,CP2,CRO,CZ(1))
        CZ(2)=CZ(1)
        GO TO 7
        ENDIF
C DEBUT NEWTON
	DO 5 I=1,15
	  DO 6 J=1,2
	  CS1=S+DS*(2*J-3)
	  CP1=CS1+SI2/CS1
	  CALL CRES(CS0,CS1,CS2,CP0,CP1,CP2,CRO,CZ(J))
	  AL(J)=ALOG(CABS(CZ(J)))
6	  CONTINUE
	DSS=DS*(AL(1)+AL(2))/(AL(1)-AL(2))
	S=S+DSS
	  IF(ABS(DSS/S).LT..0005)GO TO 7
5	CONTINUE
      WRITE(6,105)S,DSS
105   FORMAT(' Conv. lente ! '/' S = ',F7.3,' DSS = ',F7.4)
7     AN=SQRT(S*S+SI2)
      CZZ=(CZ(1)+CZ(2))/2.
      EP=OL/(4.*3.141592654*S)*ATAN2(-AIMAG(CZZ),REAL(CZZ))
C CALCUL DE PSI DELTA APRES
1     ANI=0.
      P=OL/(2.*S)
        IF(EP.LT.0.)EP=EP+P
      CALL PSD(AN,ANI,EP,OL,CS0,CS2,CP0,CP2,SI2,PSI,DLT)
      RETURN
      END
      FUNCTION CF(CA,CB)
      IMPLICIT COMPLEX (C)
      CF=(CA-CB)/(CA+CB)
      RETURN
      END
      FUNCTION CG(CA,CB)
      IMPLICIT COMPLEX (C)
      CG=(CA+CB)/(1+CA*CB)
      RETURN
      END
      SUBROUTINE PSD(AN,ANI,EP,OL,CS0,CS2,CP0,CP2,SI2,PSI,DLT)
      IMPLICIT COMPLEX (C)
      REAL PI,PSI,DLT,AN,ANI
      PI=3.141592654/180.
      CAN=CMPLX(AN,ANI)
      CS1=SQRT(CAN*CAN-SI2)
      CP1=CAN*CAN/CS1
      CB=EXP(CMPLX(0.,-4.*3.141592654*EP*REAL(CS1)/OL))
      CR1P=CF(CP1,CP0)
      CR1S=CF(CS0,CS1)
      CR2P=CF(CP2,CP1)
      CR2S=CF(CS1,CS2)
      CH=CG(CR1P,CR2P*CB)
      CI=CG(CR1S,CR2S*CB)
      CRO=CH/CI
      PSI=ATAN(ABS(CRO))/PI
      DLT=ATAN2(AIMAG(CRO),REAL(CRO))/PI
      RETURN
      END
      SUBROUTINE CRES(CS0,CS1,CS2,CP0,CP1,CP2,CRO,CZ)
      IMPLICIT COMPLEX (C)
      CR1P=CF(CP1,CP0)
      CR1S=CF(CS0,CS1)
      CR2P=CF(CP2,CP1)
      CR2S=CF(CS1,CS2)
      CA=CR1S*CR2S*CR2P-CRO*CR1P*CR2P*CR2S
      CB=CR2P+CR1P*CR1S*CR2S-CRO*(CR2S+CR1S*CR1P*CR2P)
      CC=CR1P-CRO*CR1S
      CD=SQRT(CB*CB-4.*CA*CC)
      CZ1=(-CB+CD)/(2.*CA)
      CZ2=(-CB-CD)/(2.*CA)
      CZ=CZ1
        IF((ABS(CZ1)-1.).GT.(ABS(CZ2)-1.))CZ=CZ2
      RETURN
      END
C
C
C
