!+ Subroutine to calculate the matrices for a fit to the Plack function.
!
! Method:
!	A polynomial fit in the temperature to the Planck function 
!	integrated across a spectral band is required. This is obtained
!	from a least-squares fit. Minimization of the error in the fit
!	leads to a set of linear equations. This routine determines the
!	matrix and source term in these equations as A and B.
!
! Current Owner of Code: J. M. Edwards
!
! History:
!	Version		Date			Comment
!	1.0		12-04-95		First Version under RCS
!						(J. M. Edwards)
!
! Description of Code:
!   FORTRAN 77  with extensions listed in documentation.
!
!- ---------------------------------------------------------------------
      SUBROUTINE CALC_THERMAL_COEFF(IERR
     &   , N_DEG_FIT
     &   , WAVELENGTH_SHORT, WAVELENGTH_LONG
     &   , THETA_LOW, THETA_HIGH
     &   , ALPHA, Q, A, B)
!
!
!
      IMPLICIT NONE
!
!
!     INCLUDE HEADER FILES.
      INCLUDE 'incl/dimensions_spec.ucf'
      INCLUDE 'incl/dimensions_pp.ucf'
      INCLUDE 'incl/def_std_io.icf'
      INCLUDE 'incl/thermal_pp.tcf'
      INCLUDE 'incl/error.pcf'
!
!     DUMMY ARGUMENTS.
      INTEGER	!, INTENT(OUT)
     &     IERR
!             ERROR FLAG
      INTEGER	!, INTENT(IN)
     &     N_DEG_FIT
!             DEGREE OF FIT TO FUNCTION
      REAL*8	!, INTENT(IN)
     &     WAVELENGTH_SHORT
!             SHORT WAVELENGTH OF BAND
     &   , WAVELENGTH_LONG
!             LONG WAVELENGTH OF BAND
     &   , THETA_LOW
!             LOW SCALED TEMPERATURE
     &   , THETA_HIGH
!             HIGH SCALED TEMPERATURE
     &   , ALPHA
!             NUMERICAL CONSTANT
     &   , Q
!             NUMERICAL CONSTANT
      REAL*8	!, INTENT(OUT)
     &     A(0: NPD_THERMAL_COEFF-1, 0: NPD_THERMAL_COEFF-1)
!             MATRIX OF INNER PRODUCTS
     &   , B(0: NPD_THERMAL_COEFF-1)
!             RHS FOR LEAST SQUARES EQN.
!
!     LOCAL ARGUMENTS.
      INTEGER
     &     N
!             DUMMY FOR STATEMENT FUNCTION
     &   , N_POWER
!             POWER OF TEMPERATURE
     &   , N_POINT
!             NUMBER OF INTEGRATION POINTS
     &   , I
!             LOOP VARIABLE
     &   , J
!             LOOP VARIABLE
      REAL*8
     &     U_1
!             LIMIT OF INTEGRAL
     &   , U_2
!             LIMIT OF INTEGRAL
     &   , U_3
!             LIMIT OF INTEGRAL
     &   , U_4
!             LIMIT OF INTEGRAL
     &   , X
!             DUMMY VARIABLE FOR STATEMENT FUNC.
     &   , INTEGRAL_OLD
!             OLD ESTIMATE OF INTEGRAL
     &   , INTEGRAL_NEW
!             OLD ESTIMATE OF INTEGRAL
     &   , U(NPD_THERMAL_ABSCISSA)
!             ABSCISSAE FOR U-INTEGRATION
     &   , INTEGRAND(NPD_THERMAL_ABSCISSA)
!             INTEGRANDS FOR THIS INTEGRATION
      REAL*8
     &     FNC_11
!             NAME OF STATEMENT FUNCTION
     &   , FNC_12
!             NAME OF STATEMENT FUNCTION
     &   , FNC_13
!             NAME OF STATEMENT FUNCTION
     &   , FNC_21
!             NAME OF STATEMENT FUNCTION
     &   , FNC_22
!             NAME OF STATEMENT FUNCTION
     &   , FNC_23
!             NAME OF STATEMENT FUNCTION
!
!     FUNCTIONS CALLED:
      REAL*8
     &     TRAPEZOID
!             FUNCTION FOR TRAPEZOIDAL INTEGRATION
      EXTERNAL 
     &     TRAPEZOID
!
!
!     STATEMENT FUNCTIONS FOR EVALUATION OF THE PLANCKIAN INTEGRALS
!     IN THE THERMAL SOURCE FUNCTIONS.
      FNC_11(X, N)=Q*X**3*(THETA_HIGH**(N+5)-1.0D+00/(WAVELENGTH_LONG
     &   *ALPHA*X)**(N+5))/(DBLE(N+5)*(EXP(X)-1.0D+00))
      FNC_12(X, N)=Q*X**3*(THETA_HIGH**(N+5)-THETA_LOW**(N+5))
     &   /(DBLE(N+5)*(EXP(X)-1.0D+00))
      FNC_13(X, N)=Q*X**3*(1./(WAVELENGTH_SHORT*ALPHA*X)**(N+5)
     &   -THETA_LOW**(N+5))/(DBLE(N+5)*(EXP(X)-1.0D+00))
      FNC_21(X, N)=Q*X**3*(THETA_HIGH**(N+5)-1.0D+00/(WAVELENGTH_LONG
     &   *ALPHA*X)**(N+5))/(DBLE(N+5)*(EXP(X)-1.0D+00))
      FNC_22(X, N)=Q*X**3*(1.0D+00/(WAVELENGTH_SHORT*ALPHA*X)**(N+5)
     &   -1.0D+00/(WAVELENGTH_LONG*ALPHA*X)**(N+5))
     &   /(DBLE(N+5)*(EXP(X)-1.0D+00))
      FNC_23(X, N)=Q*X**3*(1.0D+00/(WAVELENGTH_SHORT*ALPHA*X)**(N+5)
     &   -THETA_LOW**(N+5))/(DBLE(N+5)*(EXP(X)-1.0D+00))
!
!
!
!     CALCULATE THE TERMS OF THE MATRIX A FIRST OF ALL. THESE TERMS
!     REPRESENT THE INNER PRODUCTS OF THE BASIS FUNCTIONS 
!     OVER TEMPERATURE.
      DO I=0, N_DEG_FIT
         DO J=0, I
            N_POWER=I+J+1
            A(I, J)=(THETA_HIGH**N_POWER-THETA_LOW**N_POWER)
     &          /DBLE(N_POWER)
            A(J, I)=A(I, J)
         ENDDO
      ENDDO
!
!     THE TERMS TO THE RIGHT ARE MORE COMPLICATED, INVOLVING THE 
!     INTEGRAL OF THE PLANCKIAN FUNCTION OVER WAVELENGTH 
!     AND AN INNER PRODUCT OVER TEMPERATURE. THE ORDER OF 
!     INTEGRATION HAS BEEN INVERTED AND
!     THE TEMPERATURE INTEGRATION PERFORMED ANALYTICALLY, LEAVING AN
!     INTEGRATION OVER THE TRANSFORMED VARIABLE
!              U=(H_PLANCK*NU)/(K_BOLTZMANN*T_REF_THERMAL*THETA)
!     TO BE PERFORMED NUMERICALLY. THERE ARE TWO CASES.
!
      IF (WAVELENGTH_LONG*THETA_LOW.GE.WAVELENGTH_SHORT*THETA_HIGH) THEN
         U_1=1.0D+00/(ALPHA*THETA_HIGH*WAVELENGTH_LONG)
         U_2=1.0D+00/(ALPHA*THETA_LOW*WAVELENGTH_LONG)
         U_3=1.0D+00/(ALPHA*THETA_HIGH*WAVELENGTH_SHORT)
         U_4=1.0D+00/(ALPHA*THETA_LOW*WAVELENGTH_SHORT)
!
         DO I=0, N_DEG_FIT
!
            B(I)=0.0D+00
!
!           FIRST INTEGRAL:
            U(1)=U_1
            U(2)=0.5D+00*(U_1+U_2)
            U(3)=U_2
            INTEGRAL_OLD=0.0D+00
            N_POINT=3
1           DO J=1, N_POINT
               INTEGRAND(J)=FNC_11(U(J), I)
            ENDDO
            INTEGRAL_NEW=TRAPEZOID(N_POINT, U, INTEGRAND)
            IF (ABS(INTEGRAL_NEW-INTEGRAL_OLD).GT.ABS(INTEGRAL_NEW)
     &         *TOL_THERMAL) THEN
               INTEGRAL_OLD=INTEGRAL_NEW
               IF (2*N_POINT.GT.NPD_THERMAL_ABSCISSA) THEN
                  WRITE(IU_ERR, '(/A)')
     &               '*** ERROR: TOO MANY POINTS ARE REQUIRED TO '
     &               //'EVALUATE THE SOURCE FUNCTION. '
                  WRITE(IU_ERR, '(A)')
     &               'INCREASE NPD_THERMAL_ABSCISSA AND RE-COMPILE.'
                  IERR=I_ERR_FATAL
                  RETURN
               ENDIF
               DO J=N_POINT, 2, -1
                  U(2*J-1)=U(J)
                  U(2*J-2)=0.5D+00*(U(J)+U(J-1))
               ENDDO
               N_POINT=2*N_POINT-1
               GOTO 1
            ENDIF
!           USE RICHARDSON EXTRAPOLATION TO REMOVE 
!           THE LEADING ERROR TERM:
            INTEGRAL_NEW=(4.0D+00*INTEGRAL_NEW-INTEGRAL_OLD)/3.0D+00
            B(I)=B(I)+INTEGRAL_NEW
!
!           SECOND INTEGRAL:
            U(1)=U_2
            U(2)=0.5D+00*(U_2+U_3)
            U(3)=U_3
            INTEGRAL_OLD=0.0D+00
            N_POINT=3
2           DO J=1, N_POINT
               INTEGRAND(J)=FNC_12(U(J), I)
            ENDDO
            INTEGRAL_NEW=TRAPEZOID(N_POINT, U, INTEGRAND)
            IF (ABS(INTEGRAL_NEW-INTEGRAL_OLD).GT.TOL_THERMAL
     &         *ABS(INTEGRAL_NEW)) THEN
               INTEGRAL_OLD=INTEGRAL_NEW
               IF (2*N_POINT.GT.NPD_THERMAL_ABSCISSA) THEN
                  WRITE(IU_ERR, '(/A)')
     &               '*** ERROR: TOO MANY POINTS ARE REQUIRED TO '
     &               //'EVALUATE THE SOURCE FUNCTION. '
                  WRITE(IU_ERR, '(A)')
     &               'INCREASE NPD_THERMAL_ABSCISSA AND RE-COMPILE.'
                  IERR=I_ERR_FATAL
                  RETURN
               ENDIF
               DO J=N_POINT, 2, -1
                  U(2*J-1)=U(J)
                  U(2*J-2)=0.5D+00*(U(J)+U(J-1))
               ENDDO
               N_POINT=2*N_POINT-1
               GOTO 2
            ENDIF
!           USE RICHARDSON EXTRAPOLATION TO REMOVE 
!           THE LEADING ERROR TERM:
            INTEGRAL_NEW=(4.0D+00*INTEGRAL_NEW-INTEGRAL_OLD)/3.0D+00
            B(I)=B(I)+INTEGRAL_NEW
!
!           THIRD INTEGRAL:
            U(1)=U_3
            U(2)=0.5D+00*(U_3+U_4)
            U(3)=U_4
            INTEGRAL_OLD=0.0D+00
            N_POINT=3
3           DO J=1, N_POINT
               INTEGRAND(J)=FNC_13(U(J), I)
            ENDDO
            INTEGRAL_NEW=TRAPEZOID(N_POINT, U, INTEGRAND)
            IF (ABS(INTEGRAL_NEW-INTEGRAL_OLD).GT.TOL_THERMAL
     &         *ABS(INTEGRAL_NEW)) THEN
               INTEGRAL_OLD=INTEGRAL_NEW
               IF (2*N_POINT.GT.NPD_THERMAL_ABSCISSA) THEN
                  WRITE(IU_ERR, '(/A)')
     &               '** ERROR: TOO MANY POINTS ARE REQUIRED TO '
     &               //'EVALUATE THE SOURCE FUNCTION. '
                  WRITE(IU_ERR, '(A)')
     &               'INCREASE NPD_THERMAL_ABSCISSA AND RE-COMPILE.'
                  IERR=I_ERR_FATAL
                  RETURN
               ENDIF
               DO J=N_POINT, 2, -1
                  U(2*J-1)=U(J)
                  U(2*J-2)=0.5D+00*(U(J)+U(J-1))
               ENDDO
               N_POINT=2*N_POINT-1
               GOTO 3
            ENDIF
!           USE RICHARDSON EXTRAPOLATION TO REMOVE 
!           THE LEADING ERROR TERM:
            INTEGRAL_NEW=(4.0D+00*INTEGRAL_NEW-INTEGRAL_OLD)/3.0D+00
            B(I)=B(I)+INTEGRAL_NEW
!
         ENDDO
!
      ELSE
         U_1=1.0D+00/(ALPHA*THETA_HIGH*WAVELENGTH_LONG)
         U_2=1.0D+00/(ALPHA*THETA_HIGH*WAVELENGTH_SHORT)
         U_3=1.0D+00/(ALPHA*THETA_LOW*WAVELENGTH_LONG)
         U_4=1.0D+00/(ALPHA*THETA_LOW*WAVELENGTH_SHORT)
!
         DO I=0, N_DEG_FIT
!
            B(I)=0.0D+00
!
!           FIRST INTEGRAL:
            U(1)=U_1
            U(2)=0.5D+00*(U_1+U_2)
            U(3)=U_2
            INTEGRAL_OLD=0.0D+00
            N_POINT=3
4           DO J=1, N_POINT
               INTEGRAND(J)=FNC_21(U(J), I)
            ENDDO
            INTEGRAL_NEW=TRAPEZOID(N_POINT, U, INTEGRAND)
            IF (ABS(INTEGRAL_NEW-INTEGRAL_OLD).GT.TOL_THERMAL
     &         *ABS(INTEGRAL_NEW)) THEN
               INTEGRAL_OLD=INTEGRAL_NEW
               IF (2*N_POINT.GT.NPD_THERMAL_ABSCISSA) THEN
                  WRITE(IU_ERR, '(/A)')
     &               '** ERROR: TOO MANY POINTS ARE REQUIRED TO '
     &               //'EVALUATE THE SOURCE FUNCTION. '
                  WRITE(IU_ERR, '(A)')
     &               'INCREASE NPD_THERMAL_ABSCISSA AND RE-COMPILE.'
                  IERR=I_ERR_FATAL
                  RETURN
               ENDIF
               DO J=N_POINT, 2, -1
                  U(2*J-1)=U(J)
                  U(2*J-2)=0.5D+00*(U(J)+U(J-1))
               ENDDO
               N_POINT=2*N_POINT-1
               GOTO 4
            ENDIF
!           USE RICHARDSON EXTRAPOLATION TO REMOVE 
!           THE LEADING ERROR TERM:
            INTEGRAL_NEW=(4.0D+00*INTEGRAL_NEW-INTEGRAL_OLD)/3.0D+00
            B(I)=B(I)+INTEGRAL_NEW
!
!           SECOND INTEGRAL:
            U(1)=U_2
            U(2)=0.5D+00*(U_2+U_3)
            U(3)=U_3
            INTEGRAL_OLD=0.0D+00
            N_POINT=3
5           DO J=1, N_POINT
               INTEGRAND(J)=FNC_22(U(J), I)
            ENDDO
            INTEGRAL_NEW=TRAPEZOID(N_POINT, U, INTEGRAND)
            IF (ABS(INTEGRAL_NEW-INTEGRAL_OLD).GT.TOL_THERMAL
     &         *ABS(INTEGRAL_NEW)) THEN
               INTEGRAL_OLD=INTEGRAL_NEW
               IF (2*N_POINT.GT.NPD_THERMAL_ABSCISSA) THEN
                  WRITE(IU_ERR, '(/A)')
     &               '*** ERROR: TOO MANY POINTS ARE REQUIRED TO '
     &               //'EVALUATE THE SOURCE FUNCTION. '
                  WRITE(IU_ERR, '(A)')
     &               'INCREASE NPD_THERMAL_ABSCISSA AND RE-COMPILE.'
                  IERR=I_ERR_FATAL
                  RETURN
               ENDIF
               DO J=N_POINT, 2, -1
                  U(2*J-1)=U(J)
                  U(2*J-2)=0.5D+00*(U(J)+U(J-1))
               ENDDO
               N_POINT=2*N_POINT-1
               GOTO 5
            ENDIF
!           USE RICHARDSON EXTRAPOLATION TO REMOVE 
!           THE LEADING ERROR TERM:
            INTEGRAL_NEW=(4.0D+00*INTEGRAL_NEW-INTEGRAL_OLD)/3.0D+00
            B(I)=B(I)+INTEGRAL_NEW
!
!           THIRD INTEGRAL:
            U(1)=U_3
            U(2)=0.5D+00*(U_3+U_4)
            U(3)=U_4
            INTEGRAL_OLD=0.0D+00
            N_POINT=3
6           DO J=1, N_POINT
               INTEGRAND(J)=FNC_23(U(J), I)
            ENDDO
            INTEGRAL_NEW=TRAPEZOID(N_POINT, U, INTEGRAND)
            IF (ABS(INTEGRAL_NEW-INTEGRAL_OLD).GT.TOL_THERMAL
     &         *ABS(INTEGRAL_NEW)) THEN
               INTEGRAL_OLD=INTEGRAL_NEW
               IF (2*N_POINT.GT.NPD_THERMAL_ABSCISSA) THEN
                  WRITE(IU_ERR, '(/A)')
     &               '*** ERROR: TOO MANY POINTS ARE REQUIRED TO '
     &               //'EVALUATE THE SOURCE FUNCTION. '
                  WRITE(IU_ERR, '(A)')
     &               'INCREASE NPD_THERMAL_ABSCISSA AND RE-COMPILE.'
                  IERR=I_ERR_FATAL
                  RETURN
               ENDIF
               DO J=N_POINT, 2, -1
                  U(2*J-1)=U(J)
                  U(2*J-2)=0.5D+00*(U(J)+U(J-1))
               ENDDO
               N_POINT=2*N_POINT-1
               GOTO 6
            ENDIF
!           USE RICHARDSON EXTRAPOLATION TO REMOVE 
!           THE LEADING ERROR TERM:
            INTEGRAL_NEW=(4.0D+00*INTEGRAL_NEW-INTEGRAL_OLD)/3.0D+00
            B(I)=B(I)+INTEGRAL_NEW
!
         ENDDO
!
      ENDIF
!
!
      RETURN
      END
