!+ Program to check the fit to the Planck function.
!
! Method:
!	The integral of the Planck function across a spectral band
!	is evaluated precisely and from the polynomial fit at a 
!	number of temperatures and the results are printed.
!
! 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.
!
!- ---------------------------------------------------------------------
      PROGRAM CHECK_PLANCK
!
!
!
      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/math_cnst.ccf'
      INCLUDE 'incl/error.pcf'
!
!     DECLARTION OF VARIABLES.
!
      LOGICAL
     &     LOCK_CODE
!             LOGICAL TO FORBID INTERACTIVE LOOPING
!
      CHARACTER
     &     FILE_SPECTRAL*80
!             NAME OF SPECTRAL FILE
     &   , FILE_OUT*80
!             NAME OF OUTPUT FILE
      INTEGER
     &     IERR
!             ERROR FLAG
     &   , I_BAND
!             NUMBER OF BAND
     &   , N_TEMP
!             NUMBER OF TEMPERATURES
     &   , N_POINT
!             NUMBER OF POINTS OF INTEGRATION
     &   , I
!             LOOP VARIABLE
     &   , J
!             LOOP VARIABLE
      LOGICAL
     &     L_EXIST
!             EXISTENCE FLAG
      REAL*8
     &     T_LOW
!             LOW TEMPERATURE TO CHECK
     &   , T_HIGH
!             HIGH TEMPERATURE TO CHECK
     &   , T
!             TEMPERATURE
     &   , T_RATIO
!             RATIO OF T TO REFERENCE VALUE
     &   , FITTED_PLANCK
!             FITTED PLANCKIAN FUNCTION
     &   , TRUE_PLANCK
!             TRUE PLANCKIAN FUNCTION
     &   , ERROR
!             ERROR IN FITTED PLANCKIAN
     &   , INTEGRAL_OLD
!             OLD VALUE OF INTEGRAL
     &   , INTEGRAL_NEW
!             NEW VALUE OF INTEGRAL
     &   , U(NPD_THERMAL_ABSCISSA)
!             ABSCISSAE FOR INTEGRATION
     &   , INTEGRAND(NPD_THERMAL_ABSCISSA)
!             VALUES OF INTEGRAND
!
!     DECLARTION OF SPECTRAL VARIABLES.
      INCLUDE 'dec_spec.finc'
!
!     SUBROUTINES CALLED:
      EXTERNAL
     &     IN_SPECTRUM
!     FUNCTIONS CALLED:
      REAL*8
     &     PLANCK
!             PLANCKIAN FUNNCTION
     &   , TRAPEZOID
!             TRAPEZOIDAL INTEGRAL
      EXTERNAL
     &     PLANCK, TRAPEZOID
!
!
      DATA IERR/I_NORMAL/
!
!
!
      WRITE(IU_USER, '(/A)') 'ENTER THE NAME OF THE SPECTRAL FILE.'
100   READ(IU_STDIN, '(A)') FILE_SPECTRAL
      INCLUDE 'read_spec.finc'
      IF (IERR.NE.I_NORMAL) THEN
         IF (IERR.EQ.I_ERR_EXIST) THEN
            WRITE(IU_ERR, '(A)') '*** FILE DOES NOT EXIST.'
            IF (LOCK_CODE(.TRUE.)) THEN
               INCLUDE 'batch_error_main.finc'
            ELSE
               WRITE(IU_USER, '(A)') 'PLEASE RE-ENTER NAME.'
               IERR=I_NORMAL
               GOTO 100
            ENDIF
         ELSE
            STOP
         ENDIF
      ENDIF
!
      WRITE(IU_USER, '(/A)') 'ENTER BAND TO CONSIDER.'
1     READ(IU_STDIN, *) I_BAND
      IF ( (I_BAND.LT.1).OR.(I_BAND.GT.N_BAND) ) THEN
         WRITE(IU_ERR, '(A)') '+++ RESPONSE IS OUT OF RANGE.'
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error_main.finc'
         ELSE
            WRITE(IU_USER, '(A)') 'PLEASE RE-SPECIFY.'
            GOTO 1
         ENDIF
      ENDIF
!
      WRITE(IU_USER, '(/A)') 'ENTER RANGE OF TEMPERATURES.'
      READ(IU_STDIN, *) T_LOW, T_HIGH
      WRITE(IU_USER, '(A)') 'ENTER NUMBER OF TEMPERATURES.'
      READ(IU_STDIN, *) N_TEMP
!
      WRITE(IU_USER, '(/A)') 'ENTER THE NAME OF THE OUTPUT FILE'
2     READ(IU_STDIN, '(A)') FILE_OUT
      INQUIRE(FILE=FILE_OUT, EXIST=L_EXIST)
      IF (L_EXIST) THEN
         WRITE(IU_USER, '(A)') '*** THIS FILE ALREADY EXISTS.'
         IF (LOCK_CODE(.TRUE.)) THEN
            INCLUDE 'batch_error_main.finc'
         ELSE
            WRITE(IU_USER, '(A)') 'PLEASE RE-SPECIFY.'
            GOTO 2
         ENDIF
      ENDIF
      OPEN(UNIT=IU_FILE_OUTPUT, FILE=FILE_OUT, STATUS='NEW'
     &   )
!
      N_TEMP=N_TEMP-1
      DO I=0, N_TEMP
         T=T_LOW+(T_HIGH-T_LOW)*DBLE(I)/DBLE(N_TEMP)
         T_RATIO=T/T_REF_PLANCK
         FITTED_PLANCK=THERMAL_COEFFICIENT(N_DEG_FIT, I_BAND)
         DO J=N_DEG_FIT-1, 0, -1
             FITTED_PLANCK=FITTED_PLANCK*T_RATIO
     &          +THERMAL_COEFFICIENT(J, I_BAND)
         ENDDO
!        EVALUATE THE PLANCKIAN INTEGRAL BY RICHARDSON EXTRAPOLATION.
         U(1)=WAVE_LENGTH_SHORT(I_BAND)
         U(3)=WAVE_LENGTH_LONG(I_BAND)
         U(2)=0.5*(U(1)+U(3))
         INTEGRAL_OLD=0.0
         N_POINT=3
300      DO J=1, N_POINT
            INTEGRAND(J)=PLANCK(T, U(J))
         ENDDO
         INTEGRAL_NEW=TRAPEZOID(N_POINT, U, INTEGRAND)
         IF (ABS(INTEGRAL_NEW-INTEGRAL_OLD).GT.ABS(INTEGRAL_OLD)
     &      *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 PLANCKIAN FUNCTION. '
               WRITE(IU_ERR, '(/A)')
     &            'INCREASE NPD_THERMAL_ABSCISSA AND RE-COMPILE.'
               IERR=I_ERR_FATAL
               STOP
            ENDIF
            DO J=N_POINT, 2, -1
               U(2*J-1)=U(J)
               U(2*J-2)=0.5*(U(J)+U(J-1))
            ENDDO
            N_POINT=2*N_POINT-1
            GOTO 300
         ENDIF
!
!        USE RICHARDSON EXTRAPOLATION TO REFINE 
!        THE VALUS OF THE INTEGRAL.
         TRUE_PLANCK=PI*(4.*INTEGRAL_NEW-INTEGRAL_OLD)/3.
         ERROR=FITTED_PLANCK-TRUE_PLANCK
!
         WRITE(IU_FILE_OUTPUT, '(4(3X, 1PE12.5))')
     &      T, TRUE_PLANCK, FITTED_PLANCK, ERROR
!
      ENDDO
!
      CLOSE(IU_FILE_OUTPUT)
!
!
      STOP
      END
