C @(#)integbin.for	19.1 (ES0-DMD) 02/25/03 14:29:41
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
      PROGRAM INGBIN
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  17:13 - 3 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C.IDENTIFICATION
C  program    INTEGBIN              version 1.0          851028
C
C.Keywords
C  flux integration over fixed bins in 1-D images
C
C.PURPOSE
C   to be used for normalization of spectra or construction of 
C   response curves
C
C.ALGORITHM
C  get central positions of bins and their widths from table (in keyword
C  NFTAB), integrate flux in image (keyword NFIMA) over bins (where
C  necessary do a linear interpolation), write results to table
C
C.INPUT/OUTPUT
C the following keywords are used:
C NFIMA/C/1/60         name of input image
C NFTAB/C/1/60        name of table with central wavelengths, bin 
C                       widths, will also receive fluxes derived
C
C. modif MPeron 130291 add TBTCLO
C --------------------------------------------------------------------
C
      IMPLICIT NONE
      INTEGER MADRID
      CHARACTER*60 NFIMA,NFTAB
      CHARACTER    FORM*6,IDENT*72,CUNIT*72
      INTEGER   NPIX,ISTAT,IAV,ERROR(2),ICOL(3),IRCOL(2)
      INTEGER   NCOL,NROW,LL,I,NSC1
      INTEGER*8 PNTR
      INTEGER   NAXIS,KUN,KNUL
      INTEGER   IMNO,TID,NAC,NAR,DTYPE
      REAL      DATA(2),BINFLX
      DOUBLE PRECISION DEND, DSTEP, DSTART, LAMLO, LAMHI
      LOGICAL CURRY,NULL(3)
C
      INCLUDE  'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE  'MID_INCLUDE:TABLED.INC'
C
C ... set up MIDAS environment
C
      CALL STSPRO('INTEGBIN')
C
C ... get name of input frame, map it
C
      CALL STKRDC('NFIMAGE',1,1,60,IAV,NFIMA,KUN,KNUL,ISTAT)
      IF (ISTAT.NE.0) GO TO 20
      CALL STIGET(NFIMA,10,0,1,1,NAXIS,NPIX,DSTART,DSTEP,
     +IDENT,CUNIT,PNTR,IMNO,ISTAT)
      IF (ISTAT.NE.0) GO TO 20
      IF (NAXIS.NE.1) THEN
          CALL STTPUT('ERROR: input frame must be one-dimensional',
     +                ISTAT)
          ERROR(1) = 1
          GO TO 20
      END IF
      DEND    = DSTART + (NPIX-1)*DSTEP
C
C ... get name of table
C
      CALL STKRDC('NFTABLE',1,1,60,IAV,NFTAB,KUN,KNUL,ISTAT)
      IF (ISTAT.NE.0) GO TO 20
C
C ... intialize table
C
      CALL TBTOPN(NFTAB,2,TID,ISTAT)
      IF (ISTAT.NE.0) GO TO 20
      CALL TBIGET(TID,NCOL,NROW,NSC1,NAC,NAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 20
C
C ... find column 'X_AXIS' (wavelength)
C
      CALL TBLSER(TID,'X_AXIS',ICOL(1),ISTAT)
      IF (ICOL(1).EQ.-1) THEN
          ISTAT  = 10
          CALL STTPUT('ERROR: column not found',ISTAT)
          GO TO 20
      END IF
      CALL TBFGET(TID,ICOL(1),FORM,LL,DTYPE,ISTAT)
      IF (DTYPE.EQ.D_C_FORMAT) THEN
          ISTAT  = 11
          CALL STTPUT('ERROR: wrong column format',ISTAT)
          GO TO 20
      END IF
C
C ... find column 'BIN_WIDTH'
C
      CALL TBLSER(TID,'BIN_WIDTH',ICOL(2),ISTAT)
      IF (ICOL(2).EQ.-1) THEN
          ISTAT  = 10
          CALL STTPUT('ERROR: column not found',ISTAT)
          GO TO 20
      END IF
      CALL TBFGET(TID,ICOL(2),FORM,LL,DTYPE,ISTAT)
      IF (DTYPE.EQ.D_C_FORMAT) THEN
          ISTAT  = 11
          CALL STTPUT('ERROR: wrong column format',ISTAT)
          GO TO 20
      END IF
C
C ... find column 'Y_AXIS' (flux)
C
      CALL TBLSER(TID,'Y_AXIS',ICOL(3),ISTAT)
      IF (ICOL(3).EQ.-1) THEN
          ISTAT  = 10
          CALL STTPUT('ERROR: column not found',ISTAT)
          GO TO 20
      END IF
      CALL TBFGET(TID,ICOL(3),FORM,LL,DTYPE,ISTAT)
      IF (DTYPE.EQ.D_C_FORMAT) THEN
          ISTAT  = 11
          CALL STTPUT('ERROR: wrong column format',ISTAT)
          GO TO 20
      END IF
      IRCOL(1) = ICOL(1)
      IRCOL(2) = ICOL(2)
C
C ... read table
C
      DO 10 I = 1,NROW
          CALL TBRRDR(TID,I,2,IRCOL,DATA,NULL,ISTAT)
          CURRY = ( .NOT. NULL(1)) .AND. ( .NOT. NULL(2))
          IF (CURRY) THEN  !   lower limit in lambda of bin
              LAMLO  = DATA(1) - DATA(2)/2.  !   upper limit in lambda of bin
              LAMHI  = DATA(1) + DATA(2)/2.
              IF ((LAMLO.GE.DSTART) .AND. (LAMHI.LE.DEND)) THEN
C
C ... let subroutine DOINTG do the job
C
                CALL DOINTG(MADRID(PNTR),DSTART,DSTEP,NPIX,LAMLO,LAMHI,
     +                         BINFLX)
                CALL TBEWRR(TID,I,ICOL(3),BINFLX,ISTAT)
              END IF
          END IF
   10 CONTINUE
C
C the end
C
      CALL TBTCLO(TID,ISTAT)
      CALL STSEPI
      STOP 
C
C ... direct trouble makers here:
C
   20 ERROR(1) = ISTAT
      CALL STKWRI('PROGSTAT',ERROR,1,2,KUN,ISTAT)
      CALL STSEPI
      STOP 
C
      END

      SUBROUTINE DOINTG(FLUX,DSTART,DSTEP,NPIX,LAMLO,LAMHI,BINFLX)
C
C
      IMPLICIT NONE
      INTEGER NPIX,IPIXHI,IPIXLO,I
      REAL FLUX(1),PIXLO,PIXHI,FPIXHI,FPIXLO,BINFLX
      DOUBLE PRECISION DSTART, DSTEP, LAMLO, LAMHI
C
C ... all operations below assume that the flux is evenly distributed
C ... ver the pixel and that its coordinates refer to the pixel's 
C ... center !
C  !   first pixel hit (only partly)
      PIXLO  = 0.5 + (LAMLO-DSTART)/DSTEP  !   last pixel hit (also only partly)
      PIXHI  = 0.5 + (LAMHI-DSTART)/DSTEP  !   nmubers of first ...
      IPIXLO = INT(PIXLO)  !   ... and last pixel
      IPIXHI = INT(PIXHI)  !   fraction of first and ...
      FPIXLO = IPIXLO + 1 - PIXLO  !   ... last pixel to be considered
      FPIXHI = PIXHI - IPIXHI
C
C ... do a linear interpolation across the first and last pixel,
C ... add up the pixels in between, return mean flux per pixel
C
      BINFLX = FPIXLO*FLUX(IPIXLO) + FPIXHI*FLUX(IPIXHI)
      DO 10 I = IPIXLO + 1,IPIXHI - 1
          BINFLX = BINFLX + FLUX(I)
   10 CONTINUE
      BINFLX = BINFLX/ (IPIXHI-IPIXLO-1+FPIXLO+FPIXHI)
C
      RETURN
      END
