      SUBROUTINE POL(RN1,RN2,IA,IAP,CINDEX,THICKNESS,N)
C     THIS COMPUTES RAY POLARIZATION INTENSITIES AND PHASES
C
      REAL*8 RN1,RN2,IA,IAP
C
      INTEGER N,NLAY
C
      COMPLEX*16 ROOTONE,CPII,CTWO,INDEXI,INDEXIP
C
      COMPLEX*16 CINDEX(N),CTEMP,COSTHETAL
C
      REAL*8 THICKNESS(N),THICKNESSL
C
      REAL*8 CTEMPR,CTEMPI
C
C
C       COSPOL IS AN IMPORTANT NUMBER. IT IS THE COSINE OF THE ANGLE
C       BETWEEN THE INCIDENT PLANE DIRECTION AND THE Y=PLANE OF THE RAY
C       FOR ON-AXIS IN THE Y-Z PLANE OF A CENTERED SYSTEM, THE ANGLE IS 0
C       FOR ON-AXIS RAYS IN THE X-Z PLANE OF A CENTERED SYSTEM, THE ANGLE IS
C       PII/2
C
C     NOTE THAT IT IS ASSUMED THAT THE IMAGINARY PART OF THE INCIDENT MEDIUM'S
C     REFRACTIVE INDEX IS ALWAYS ZERO.
C     NOTE THAT IT IS ASSUMED THAT THE IMAGINARY PART OF THE SUBSTRATE MEDIUM'S
C     REFRACTIVE INDEX IS ALWAYS ZERO.
C
C       COMPLEX REPRESENTATION OF i, THE ROOT OF -1
      ROOTONE=DCMPLX(0.0D0,1.0D0)
C       COMPLEX REPRESENTATION OF PI
      CPII=DCMPLX(3.14159265358979323846D0,0.0D0)
C       COMPLEX REPRESENTATION OF 2.0
      CTWO=DCMPLX(2.0D0,0.0D0)

C     RN1 IS THE REAL INDEX OF THE INCIDENT MEDIUM
C     RN2 IS THE REAL INDEX OF THE SUBSTRATE
C     IA IS THE ANGLE OF INCIDENCE IN RADIANS
C     IAP IS THE ANGLE OF REFLECTION,REFRACTION OR DIFFRACTION IN RADIANS
C
C     SET COMPLEX INDECIES AND ANGLES
C       INDEX OF THE INCIDENT MEDIUM
                INDEXI=DCMPLX(RN1,0.0D0)
C       INDEX OF THE SUBSTRATE
                INDEXIP=DCMPLX(RN2,0.0D0)
C       ANGLE OF INCIDENCE
                THETAI=DCMPLX(IA,0.0D0)
C       ANGLE OF REFRACTION,REFLECTION OR DIFFRACTION
                THETAIP=DCMPLX(IAP,0.0D0)
C

C     CINDEX(I) IS THE COMPLEX INDEX OF COATING LAYER I
C     THICKNESS(I) IS THE THICKNESS  OF COATING LAYER I IN MICRONS
C
C     A MULTILAYER COATING, ASSUMING NO-COHERENCE
C     IF RN1 AND RN2 HAVE OPPOSITE SIGNS THEN WE HAVE A REFLECTION
C     SO WE USE THE R VALUE
C     IF RN1 AND RN2 HAVE THE SAME SIGNS THEN WE HAVE A REFRACTION
C     SO WE USE THE T VALUE
                        NLAY=N
C
                       DO L=1,NLAY
C     CONSTRUCT THE 2X2 L-LAYER MATRIX
C     INDEXL IS THE COMPLEX INDEX AT THE TRACED WAVELENGTH FOR THE L-TH LAYER
      INDEXL=CINDEX(L)
C       T1 IS THE REAL PART OF THE LAYER INDEX
      T1=DBLE(INDEXL)
C       T2 IS THE IMAGINARY PART OF THE LAYER INDEX
      T2=DIMAG(INDEXL)
C       T3 IS THE LAYER THICKNESS
      T3=THICKNESS(L)
C       T4 IS THE REAL PART OF THE INCIDENT MEDIA INDEX
      T4=DBLE(INDEXI)
C       T5 IS THE REAL PART OF THE SINE OF THE INCIDENT ANGLE
      T5=DSIN(IA)
C
      CTEMPR=(T1**2)-(T2**2)-((T4**2)*(DSIN(IA)**2))
      CTEMPI=2.0D0*T1*T2
      CTEMP=DCMPLX(CTEMPR,CTEMPI)
      CTEMP=CDSQRT(CTEMP)
      COSTHETAL=CTEMP/INDEXL
C     THICKNESS OF L IN MICRONS IS:
      THICKNESSL=DCMPLX(T3)
C     COS(PHIL)*THE COMPLEX INDEX OF THE LAYER IS JUST:
C     U SUB L PARALLEL IS:
      ULPAR=INDEXL/COSTHETAL
C     U SUB L PERPENDICULAR IS:
      ULPER=INDEXL*COSTHETAL
C     DELTAL IS:
      DELTAL=(CTWO*CPII/CACTIVE_WAVE)*(INDEXL*THICKNESSL*COSTHETAL)
C     THE LAYER MATRIX FOR PARALLEL POLARIZATION IS:
      MPARL(1,1)=CDCOS(DELTAL)
      MPARL(1,2)=-(ROOTONE/ULPAR)*CDSIN(DELTAL)
      MPARL(2,1)=-ROOTONE*ULPAR*CDSIN(DELTAL)
      MPARL(2,2)=MPARL(1,1)
C     THE LAYER MATRIX FOR PERPENDICULAR POLARIZATION IS:
      MPERL(1,1)=CDCOS(DELTAL)
      MPERL(1,2)=-(ROOTONE/ULPER)*CDSIN(DELTAL)
      MPERL(2,1)=-ROOTONE*ULPER*CDSIN(DELTAL)
      MPERL(2,2)=MPERL(1,1)
C     MULTIPLY THE CURRENT L-LAYER M-MATRIX INTO THE FINAL M-MATRIX
                       IF(L.EQ.1) THEN
        MPAR(1,1)=MPARL(1,1)
        MPAR(1,2)=MPARL(1,2)
        MPAR(2,1)=MPARL(2,1)
        MPAR(2,2)=MPARL(2,2)
        MPER(1,1)=MPERL(1,1)
        MPER(1,2)=MPERL(1,2)
        MPER(2,1)=MPERL(2,1)
        MPER(2,2)=MPERL(2,2)
                       ELSE
C     L GREATER THAN 1
        MPAR=MATMUL(MPARL,MPAR)
        MPER=MATMUL(MPERL,MPER)
                       END IF
                       END DO
C
C       NOW COMPUTE THE AMPLITUDE TRANSMITTANCE AND REFLECTANCE
C       COEFFICIENTS SMALLTPAR, SMALLTPER, SMALLRPAR AND SMALLRPER
        SMALLRPAR=((MPAR(1,1)*INDEXI)
     1  +(MPAR(1,2)*INDEXI*INDEXIP)
     2  -(MPAR(2,1))
     3  -(MPAR(2,2)*INDEXIP))/
     4   ((MPAR(1,1)*INDEXI)
     5  +(MPAR(1,2)*INDEXI*INDEXIP)
     6  +(MPAR(2,1))
     7  +(MPAR(2,2)*INDEXIP))
        SMALLRPER=((MPER(1,1)*INDEXI)
     1  +(MPER(1,2)*INDEXI*INDEXIP)
     2  -(MPER(2,1))
     3  -(MPER(2,2)*INDEXIP))/
     4   ((MPER(1,1)*INDEXI)
     5  +(MPER(1,2)*INDEXI*INDEXIP)
     6  +(MPER(2,1))
     7  +(MPER(2,2)*INDEXIP))
        SMALLTPAR=(CTWO*INDEXI)/
     1   ((MPAR(1,1)*INDEXI)
     2  +(MPAR(1,2)*INDEXI*INDEXIP)
     3  +(MPAR(2,1))
     4  +(MPAR(2,2)*INDEXIP))
        SMALLTPER=(CTWO*INDEXI)/
     1   ((MPER(1,1)*INDEXI)
     2  +(MPER(1,2)*INDEXI*INDEXIP)
     3  +(MPER(2,1))
     4  +(MPER(2,2)*INDEXIP))
C
C       THE INTENSITY TRANSMISSION AND REFLECTION COEFFICIENTS ARE:
        BIGTPAR=CDABS(INDEXIP/INDEXI)*((CDABS(SMALLTPAR))**2)
        BIGTPER=CDABS(INDEXIP/INDEXI)*((CDABS(SMALLTPER))**2)
        BIGRPAR=CDABS(SMALLRPAR)**2
        BIGRPER=CDABS(SMALLRPER)**2
        PHASETPAR=DATAN2(DBLE(SMALLTPAR),DIMAG(SMALLTPAR))
        PHASETPER=DATAN2(DBLE(SMALLTPER),DIMAG(SMALLTPER))
        PHASERPAR=DATAN2(DBLE(SMALLRPAR),DIMAG(SMALLRPAR))
        PHASERPER=DATAN2(DBLE(SMALLRPER),DIMAG(SMALLRPER))
C       NOW POLANG IS ANGLE BETWEEN THE PLANE OF INCIDENCE AND THE
C       PLANE FORMED BY THE "Y" RAY AND THE "Z" RAY VECTORS
C       THE T AND R COMPONENTS IN THE "Y" DIRECTION ARE:
      IF(RN1.GT.0.0D0.AND.RN2.GT.0.0D0.OR.RN1.LT.0.0D0.AND.RN2.LT.
     10.0D0) THEN
C     REFRACTION
      OLDABSCOEF(1:10)=ABSCOEF(1:10)
      ENERGY_FACTOR=(BIGTPAR+BIGTPER)/2.0D0
      ENERGY_FACTOR=ENERGY_FACTOR*
     1DEXP(-DABS(OLDABSCOEF(INT(WA3))*PATHL))
        FACT_PAR=BIGTPAR
        FACT_PER=BIGTPER
        PHASE_PAR=PHASETPAR
        PHASE_PER=PHASETPER
C     CHECK FOR ABSORB
C     SET DEFAULT
      ABSCOEF(1:10)=0.0D0
      READ(UNIT=34,FMT=*,ERR=9994,END=9994) AB
      READ(UNIT=34,FMT=*,ERR=9994,END=9994) ABSCOEF(1)
      READ(UNIT=34,FMT=*,ERR=9994,END=9994) ABSCOEF(2)
      READ(UNIT=34,FMT=*,ERR=9994,END=9994) ABSCOEF(3)
      READ(UNIT=34,FMT=*,ERR=9994,END=9994) ABSCOEF(4)
      READ(UNIT=34,FMT=*,ERR=9994,END=9994) ABSCOEF(5)
      READ(UNIT=34,FMT=*,ERR=9994,END=9994) ABSCOEF(6)
      READ(UNIT=34,FMT=*,ERR=9994,END=9994) ABSCOEF(7)
      READ(UNIT=34,FMT=*,ERR=9994,END=9994) ABSCOEF(8)
      READ(UNIT=34,FMT=*,ERR=9994,END=9994) ABSCOEF(9)
      READ(UNIT=34,FMT=*,ERR=9994,END=9994) ABSCOEF(10)
9994  CALL CLOSE_FILE(34,1)
                       ELSE

C     REFLECTION
      ENERGY_FACTOR=(BIGRPAR+BIGRPER)/2.0D0
      OLDABSCOEF(1:10)=ABSCOEF(1:10)
      ENERGY_FACTOR=ENERGY_FACTOR*
     1DEXP(-DABS(OLDABSCOEF(INT(WA3))*PATHL))
        FACT_PAR=BIGRPAR
        FACT_PER=BIGRPER
        PHASE_PAR=PHASERPAR
        PHASE_PER=PHASERPER
C     CHECK FOR ABSORB
C     SET DEFAULT
      ABSCOEF(1:10)=0.0D0
      READ(UNIT=34,FMT=*,ERR=9995,END=9995) AB
      READ(UNIT=34,FMT=*,ERR=9995,END=9995) ABSCOEF(1)
      READ(UNIT=34,FMT=*,ERR=9995,END=9995) ABSCOEF(2)
      READ(UNIT=34,FMT=*,ERR=9995,END=9995) ABSCOEF(3)
      READ(UNIT=34,FMT=*,ERR=9995,END=9995) ABSCOEF(4)
      READ(UNIT=34,FMT=*,ERR=9995,END=9995) ABSCOEF(5)
      READ(UNIT=34,FMT=*,ERR=9995,END=9995) ABSCOEF(6)
      READ(UNIT=34,FMT=*,ERR=9995,END=9995) ABSCOEF(7)
      READ(UNIT=34,FMT=*,ERR=9995,END=9995) ABSCOEF(8)
      READ(UNIT=34,FMT=*,ERR=9995,END=9995) ABSCOEF(9)
      READ(UNIT=34,FMT=*,ERR=9995,END=9995) ABSCOEF(10)
9995  CALL CLOSE_FILE(34,1)
      IF(TIR) ENERGY_FACTOR=1.0D0
        IF(TIR) THEN
      PHASE_PAR=-PII/2.0D0
      PHASE_PER=-PII/2.0D0
      FACT_PAR=1.0D0
      FACT_PER=1.0D0
      POLANG=0.0D0
                        END IF
                       END IF
      DEALLOCATE (CINDEX,LAYERNAME,THICKNESS,STAT=ALLOERR)
                       POLEXT=.TRUE.
                       RETURN
                       END IF
C**********************************************************************
 999                   ENERGY_FACTOR=1.0D0
      POLANG=0.0D0
                       CALL CLOSE_FILE(34,1)
      DEALLOCATE (CINDEX,LAYERNAME,THICKNESS,STAT=ALLOERR)
C       FILE READ ERROR
C**********************************************************************
                       RETURN
                       END
