C
C THIS PROGRAM COMPUTES CHANDRASEKHARS' H FUNCTION IN FORTRAN
C DEDICATION: TO MRS. FRANCES H. BREEN.
C
      REAL*8 XMU(100),WEIGHT(100),XL,XU,H(100),Y(100)
      REAL*8 OMEGA,U, HALF, RHS,SUM,SUM1,SUM2, SUM3, PI
      INTEGER *2 IU
      INTEGER*4 ITERATION
      PARAMETER (NTHETA=8, XL=0.0D0, XU=1.0D0, HALF=0.5D0)
      PARAMETER( OMEGA= 0.999999, ITERATION= 10000)
      PARAMETER( IU=6 )
      DATA PI/3.141592654D0/

C
C THIS IS THE PSI FUNCTION FOR ISOTROPIC SCATTERING (RT P106 # 5)
C
      PSI(U)=HALF*OMEGA
C
      CALL GAUSST(NTHETA,XL,XU,XMU,WEIGHT)
      DO 100 I=1,NTHETA
100   H(I)=1.0D0
C
      DO 200 N=1,ITERATION
      
      DO 300 I=1,NTHETA
              SUM=0.0D0
            DO 400 J=1,NTHETA
              SUM=SUM+(WEIGHT(J)*H(J)*PSI(XMU(J)))/(XMU(I)+XMU(J))
400         CONTINUE
      H(I)=1.0D0+ XMU(I)*H(I)*SUM
300   CONTINUE
C
C CHECK RESULTS AFTER EVERY 200 ITERATIONS
C FOR CONVERGENCE (MANY WAYS OF IMPROVING THIS )
C
      IF(MOD(N,200).NE.0)GOTO 200
      WRITE(IU,*)'HFUNCTIONS FOR ITERATION=',N
      DO 500 L=1,NTHETA
500   WRITE(IU,*)L,XMU(L),H(L)
200   CONTINUE
C
C CALCULATE MOMENTS OF H FUNCTION FOR THE GIVEN PSI
C USED TO CHECK ACCURACY OF H FUNCTIONS
C
       SUM1=0.0D0
       SUM2=0.0D0
       SUM3=0.0D0
       DO 600 J=1,NTHETA
          SUM1=SUM1+WEIGHT(J)*H(J)*PSI(XMU(J))
          SUM2=SUM2+WEIGHT(J)*PSI(XMU(J))
          SUM3=SUM3+WEIGHT(J)*H(J)
600    CONTINUE

C
C CHECK THE CORRECTNESS OF H FUNCTIONS
C
       RHS=1.0-DSQRT(1.0-2.*SUM2)
       WRITE(IU,*)' LHS OF EQ11 P106 RT= ', SUM1
       WRITE(IU,*)' RHS OF EQ11 P106 RT= ', RHS
       WRITE(IU,*)' MEAN VALUE OF H FUNCTION=',SUM3
C 
C Compute Final radiances
      WRITE(6,*)'     MU','            MUO', '              RAD'
      DO I=1, NTHETA
      DO J=1, NTHETA
      RAD=  .25*OMEGA*XMU(J)*H(I)*H(J)/(XMU(I)+XMU(J))
      WRITE(6,*)XMU(I),XMU(J),RAD
      ENDDO
      ENDDO
      STOP
      END
      
      SUBROUTINE GAUSST(N,XL,XU,SS,TT)                                  
      IMPLICIT REAL*8(A-H,O-Z)                                           
      DIMENSION Z(200),PA(200),W(200),R(200),SS(1),TT(1)               
      TOL = 1.0D-16                                                     
      PI = 3.141592653589793D+00                                        
      AA = 2.0D+00/PI**2                                                
      AB = -62.0D+00/(3.0D+00*PI**4)                                    
      AC = 15116.0D+00/(15.0D+00*PI**6)                                
      AD = -12554474.D+00/(105.0D+00*PI**8)                             
      PA(1) = 1.D0                                                      
      EN = N                                                           
      NP1 = N+1                                                         
      U = 1.0D+00-(2.0D+00/PI)**2                                       
      D = 1.0D+00/DSQRT((EN+0.5D+00)**2+U/4.0D+00)                      
      DO 100  I = 1,N                                                   00000150
      SM = I                                                            00000160
      AZ = 4.0D+00*SM-1.0D+00                                           00000170
      AE = AA/AZ                                                        00000180
      AF = AB/AZ**3                                                     00000190
      AG = AC/AZ**5                                                     00000200
      AH = AD/AZ**7                                                     00000210
  100 Z(I) = 0.25D+00*PI*(AZ+AE+AF+AG+AH)                               00000220
      DO 200  K = 1,N                                                   00000230
      X = DCOS(Z(K)*D)                                                  00000240
    1 PA(2) = X                                                         00000250
      DO 201  NN = 3,NP1                                                00000260
      ENN = NN-1                                                        00000270
  201 PA(NN) = ((2.0D+00*ENN-1.0D+00)*X*PA(NN-1)-(ENN-1.0D+00)*PA(NN-2))00000280
     +         /ENN                                                     00000290
      PNP = EN*(PA(N)-X*PA(NP1))/(1.0D+00-X*X)                          00000300
      XI = X-PA(NP1)/PNP                                                00000310
      XD = DABS(XI-X)                                                   00000320
      XDD = XD-TOL                                                      00000330
      IF (XDD) 3,3,2                                                    00000340
    2 X = XI                                                            00000350
      GO TO 1                                                           00000360
    3 R(K) = X                                                          00000370
  200 W(K) = 2.0D+00*(1.0D+00-X*X)/(EN*PA(N))**2                        00000380
      AP = (XU-XL)/2.D0                                                 00000390
      BP = (XU+XL)/2.D0                                                 00000400
      DO 300  I = 1,N                                                   00000410
      M=N-I+1                                                           00000420
      SS(M)=BP+AP*R(I)                                                  00000430
300   TT(M)=AP*W(I)                                                     00000440
      RETURN                                                            00000450
      END                                                               00000460