PROGRAM APOL
  IMPLICIT NONE
  INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(14) 
  REAL(KIND=DP), DIMENSION(80,80) :: LL,MM,NN,II
  REAL(KIND=DP), DIMENSION(80) :: TAU,LAM,BNU,OML,S,P
  REAL(KIND=DP), DIMENSION(160,160) :: GG 
  REAL(KIND=DP), DIMENSION(160) :: RHS 
  REAL(KIND=DP), DIMENSION(159,80) :: EE  
  REAL(KIND=DP), DIMENSION(159) :: MUA,SI,PI,Q,INT,POL 
  INTEGER :: J,K
  READ(1,*) LL    ! read in matrix representation of lambda-transform
  READ(2,*) MM    ! read in M-transform matrix
  READ(3,*) NN    ! read in N-transform matrix
  READ(4,*) TAU   ! array of 80 optical depths
  DO J=1,80
  READ(5,*) BNU(J), LAM(J)          ! lambda = kappa/(kappa+sigma)
  END DO 
  II= 0.0                           ! make 80x80 identity matrix
  DO J=1,80
  II(J,J)=1.0
  END DO 
  OML= 1.0-LAM  ! (1-lambda)
  DO J=1,80
    DO K=1,80
    LL(J,K)= OML(J)*LL(J,K)
    MM(J,K)= OML(J)*MM(J,K)
    NN(J,K)= OML(J)*NN(J,K)
    END DO
  END DO
  ! construct GG, the left side of equations (36) [see Harrington 1970]
  GG(1:80,1:80)= II - LL  
  GG(1:80,81:160)= (-1.0/3.0)*MM
  GG(81:160,1:80)= (-3.0/8.0)*MM
  GG(81:160,81:160)= II - (3.0/8.0)*NN
  RHS=0.0
  RHS(1:80)=LAM(1:80)*BNU(1:80)  ! right-hand side of linear system
  CALL MATINV(GG,RHS,160)        ! solve system of linear equations
  S=RHS(1:80)                    ! RHS is now the solution vector
  P=RHS(81:160)                  ! 1st half is s(tau), 2nd is p(tau)
  DO J=1,80
  WRITE(6,*) TAU(J),BNU(J),LAM(J),S(J),P(J)
  END DO
  write(6,*) '-----------------------------------------------------------------'
  READ(20,*) MUA     ! read in the 159 mu=cos(theta) values
  READ(21,*) EE      ! read in the exp-transform matrix (e_ij of equation 35)
  SI= MATMUL(EE,S)
  PI= MATMUL(EE,P)
  Q=PI*(1.0-MUA*MUA)              ! Q(mu) (negative if plane is perp. to radius)
  INT=SI+PI*((1.0/3.0)-MUA*MUA)   ! I(mu)
  POL=Q/INT                       ! fractional polarization 
  DO J=1,159
  WRITE(6,*) J,MUA(J),INT(J),Q(J),POL(J)
  END DO
  END PROGRAM

 SUBROUTINE MATINV(A,B,N)
 INTEGER, PARAMETER :: DP = SELECTED_REAL_KIND(14)
 REAL(KIND=DP), DIMENSION (160,160) :: A
 REAL(KIND=DP), DIMENSION (160,2) :: INDEX 
 REAL(KIND=DP), DIMENSION (160) :: B,IPIVOT,PIVOT
 REAL(KIND=DP) :: SWAP
 INTEGER :: N,J,K,L
!
!     INITIALIZATION
!
      DETERM=1.0
      DO 20 J=1,N
   20 IPIVOT(J)=0
      DO 550 I=1,N
!
!     SEARCH FOR PIVOT ELEMENT
!
      AMAX=0.0
      DO 105 J=1,N
      IF (IPIVOT(J)-1) 60, 105, 60
   60 DO 100 K=1,N
      IF (IPIVOT(K)-1) 80, 100, 740
   80 IF (ABS (AMAX)-ABS (A(J,K))) 85, 100, 100
   85 IROW=J
      ICOLUM=K
      AMAX=A(J,K)
  100 CONTINUE
  105 CONTINUE
      IPIVOT(ICOLUM)=IPIVOT(ICOLUM)+1
!
!     INTERCHANGE ROWS TO PUT PIVOT ELEMENT ON DIAGONAL
!
      IF (IROW-ICOLUM) 140, 260, 140
  140 DETERM=-DETERM
      DO 200 L=1,N
      SWAP=A(IROW,L)
      A(IROW,L)=A(ICOLUM,L)
  200 A(ICOLUM,L)=SWAP
      SWAP=B(IROW)
      B(IROW)=B(ICOLUM)
      B(ICOLUM)=SWAP
  260 INDEX(I,1)=IROW
      INDEX(I,2)=ICOLUM
      PIVOT(I)=A(ICOLUM,ICOLUM)
      DETERM=DETERM*PIVOT(I)
!
!     DIVIDE PIVOT ROW BY PIVOT ELEMENT
!
      A(ICOLUM,ICOLUM)=1.0
      PINV=1.0/PIVOT(I)
      DO 350 L=1,N
  350 A(ICOLUM,L)=A(ICOLUM,L)*PINV
      B(ICOLUM)=B(ICOLUM)*PINV
!
!     REDUCE NON-PIVOT ROWS
!
      DO 550 L1=1,N
      IF(L1-ICOLUM) 400, 550, 400
  400 T=A(L1,ICOLUM)
      A(L1,ICOLUM)=0.0
      DO 450 L=1,N
  450 A(L1,L)=A(L1,L)-A(ICOLUM,L)*T
      B(L1)=B(L1)-B(ICOLUM)*T
  550 CONTINUE
!
!     INTERCHANGE COLUMNS
!
      DO 710 I=1,N
      L=N+1-I
      IF (INDEX(L,1)-INDEX(L,2)) 630, 710, 630
  630 JROW=INDEX(L,1)
      JCOLUM=INDEX(L,2)
      DO 705 K=1,N
      SWAP=A(K,JROW)
      A(K,JROW)=A(K,JCOLUM)
      A(K,JCOLUM)=SWAP
  705 CONTINUE
  710 CONTINUE
  740 RETURN
      END