C @(#)irsbbfit.for	17.1.1.1 (ES0-DMD) 01/25/02 17:53:05
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.COPYRIGHT   (C) 1992 European Southern Observatory
C.IDENT       .for
C.AUTHOR      E. Oliva,  Firenze-Arcetri
C.KEYWORDS    Spectroscopy, IRSPEC
C
C.PURPOSE     Execute command  ...
C
C.ALGORITHM
C
C
C.INPUT/OUTPUT
C
C
C.VERSION     1.0    Creation     02.09.1992   E. Oliva
C
C-------------------------------------------------------
C
      PROGRAM BBFIT

      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      CHARACTER*60 TABLE
      INTEGER UNITS
      LOGICAL NULL
      DIMENSION WL(10000),FLUX(10000)
      COMMON /VMR/ MADRID(1)

      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'

      DATA MAXDIM/2/

      IRET=1
      CALL STSPRO('BBFIT')
C
C Get units flag
C
      CALL STKRDI('units',1,1,IRET,UNITS,KUNIT,KNUL,ISTAT)
C
C Get name of input table
C
      CALL STKRDC('table',1,1,60,IRET,TABLE,KUNIT,KNUL,ISTAT)
      CALL CLNTAB(TABLE,TABLE,0)
C
C OPEN TABLE FOR READING (IOPEN=0)
C
      IOPEN=0
      CALL TBTOPN(TABLE,IOPEN,IDETAB,ISTAT)
C
C GET NUMBERS CORRESPONDING TO wl AND flux
C
      CALL TBLSER(IDETAB,'wl',NXCOL,ISTAT)
      CALL TBLSER(IDETAB,'flux',NYCOL,ISTAT)
C
C GET INFORMATION UPON THE TABLE SIZE (N.B. NROW = # OF WL points)
C
      CALL TBIGET(IDETAB,NCOL,NROW,NSORT,NCALL,NRALL,ISTAT)
      NWL=NROW
C
C COPY TABLE ELEMENTS INTO VECTORS WL,FLUX
C
      DO IROW=1,NWL
        CALL TBERDR(IDETAB,IROW,NXCOL,RVAL,NULL,ISTAT)
        IF(RVAL.LE.0) 
     +    CALL STETER(1,'One of the wavelengths is <=0')
        WL(IROW)=RVAL
        IF(UNITS.EQ.1) WL(IROW)=WL(IROW)/1.E4
        CALL TBERDR(IDETAB,IROW,NYCOL,RVAL,NULL,ISTAT)
        IF(RVAL.LE.0) 
     +    CALL STETER(1,'One of the fluxes is <=0')
        FLUX(IROW)=RVAL
      ENDDO
      CALL TBTCLO(IDETAB,ISTAT)
C
C CALL ROUTINE WHICH fits the black-body
C
      CALL FIT_BB(WL,FLUX,NWL,TEMP,FACT0)
C
C Write TEMP, FACT0 in OUTPUTR keywords
C
      CALL STKWRR('outputr',TEMP,1,1,KUN,ISTAT)
      CALL STKWRR('outputr',FACT0,2,1,KUN,ISTAT)
C
C RELEASE FILES, UPDATE KEYWORDS AND EXIT
C
      CALL STSEPI
      END
C
C
C
      SUBROUTINE FIT_BB(WL,FLUX,NWL,T,F0)
      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      DIMENSION WL(NWL),FLUX(NWL),F(3)

      T=5000.
      DT=1000.
      FMIN=1E30
      NITER=0
10    CONTINUE
      DO I=1,3
        F(I)=FUNCT(T+(I-1)*DT,F0,WL,FLUX,NWL)
      ENDDO
      IF(F(1).GT.F(2) .AND. F(2).GT.F(3)) THEN
        T=T+2.*DT
        GO TO 10
      ENDIF
      IF(F(1).LT.F(2) .AND. F(2).LT.F(3)) THEN
        T=T-DT
        GO TO 10
      ENDIF
      IF(F(1).EQ.F(2) .AND. F(2).EQ.F(3)) THEN
        T=T+DT
      ELSE
        T=T+2*DT-DT*((F(3)-F(2))/(F(3)-2*F(2)+F(1))+0.50)
      ENDIF
      FNEW=FUNCT(T,F0,WL,FLUX,NWL)
C	WRITE(*,*) 'T,F0,FNEW',T,F0,FNEW
      ERR=ABS(FNEW/FMIN-1)
      IF(ERR.GT.1e-4) THEN
        DT=DT/2.
        DT=MAX(DT,.1)
        FMIN=FNEW
        NITER=NITER+1
        IF(NITER.GT.100) THEN
          CALL STETER(1,'No convergency after 100 iterations')
        ELSE
          GO TO 10
        ENDIF
      ENDIF
      RETURN
      END

      FUNCTION FUNCT(T,F0,WL,FLUX,NWL)
      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
      DIMENSION WL(NWL),FLUX(NWL)

      SUMB2=0.0
      SUMFB=0.0
      DO I=1,NWL
        BBI=1/WL(I)**5/(EXP(14388./WL(I)/T)-1)
        SUMB2=SUMB2+BBI*BBI
        SUMFB=SUMFB+BBI*FLUX(I)
      ENDDO
      F0=SUMFB/SUMB2
      FUNCT=0.0
      DO I=1,NWL
        BBI=1/WL(I)**5/(EXP(14388./WL(I)/T)-1)
        CONTR=F0*BBI-FLUX(I)
        FUNCT=FUNCT+CONTR*CONTR
      ENDDO
      RETURN
      END