C===========================================================================
C Copyright (C) 1995,2004 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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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: (c) 1987,2004 European Southern Observatory,
C                                         all rights reserved
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/80         name of input image
C NFTAB/C/1/80        name of table with central wavelengths, bin 
C                       widths, will also receive fluxes derived
C
C.VERSION
C 040312	last modif
C 
C --------------------------------------------------------------------
C
      IMPLICIT NONE
C 
      INTEGER NPIX,ISTAT,IAV,ICOL(3),IRCOL(2)
      INTEGER NCOL,NROW,LL,I,NSC1,NAXIS,KUN,KNUL
      INTEGER IMNO,TID,NAC,NAR,DTYPE
      INTEGER MADRID
C 
      INTEGER*8 PNTR
C 
      CHARACTER*60 NFIMA,NFTAB
      CHARACTER    FORM*6,IDENT*72,CUNIT*72
C 
      REAL    DATA(2),DAT2H,BINFLX
C 
      DOUBLE PRECISION DEND, DSTEP, DSTART, LAMLO, LAMHI
C 
      LOGICAL CURRY,NULL(3)
C
      INCLUDE  'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE  'MID_INCLUDE:TABLED.INC'
      DATA  NFIMA /' '/, NFTAB /' '/
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)
      CALL STIGET(NFIMA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,1,
     +            NAXIS,NPIX,DSTART,DSTEP,IDENT,CUNIT,PNTR,IMNO,ISTAT)
      IF (NAXIS.NE.1) 
     +   CALL STETER(22,'input spectrum must be 1-dim frame')
C 
C  so, this command only works for linear wc systems!
      DEND = DSTART + (NPIX-1)*DSTEP
C
C ... get name of table, open it + get info
C
      CALL STKRDC('NFTABLE',1,1,60,IAV,NFTAB,KUN,KNUL,ISTAT)
      CALL TBTOPN(NFTAB,F_IO_MODE,TID,ISTAT)
      CALL TBIGET(TID,NCOL,NROW,NSC1,NAC,NAR,ISTAT)
C
C ... find column 'X_AXIS' (wavelength)
C
      CALL TBLSER(TID,'X_AXIS',ICOL(1),ISTAT)
      IF (ICOL(1).EQ.-1) THEN
         CALL STETER(10,'column not found')
      END IF
      CALL TBFGET(TID,ICOL(1),FORM,LL,DTYPE,ISTAT)
      IF (DTYPE.EQ.D_C_FORMAT) THEN
         CALL STETER(11,' wrong column format',ISTAT)
      ENDIF
C
C ... find column 'BIN_WIDTH'
C
      CALL TBLSER(TID,'BIN_WIDTH',ICOL(2),ISTAT)
      IF (ICOL(2).EQ.-1) THEN
         CALL STETER(12,'column not found')
      ENDIF
      CALL TBFGET(TID,ICOL(2),FORM,LL,DTYPE,ISTAT)
      IF (DTYPE.EQ.D_C_FORMAT) THEN
         CALL STETER(13,' wrong column format',ISTAT)
      ENDIF
C
C ... find column 'Y_AXIS' (flux)
C
      CALL TBLSER(TID,'Y_AXIS',ICOL(3),ISTAT)
      IF (ICOL(3).EQ.-1) THEN
         CALL STETER(14,'column not found')
      ENDIF
      CALL TBFGET(TID,ICOL(3),FORM,LL,DTYPE,ISTAT)
      IF (DTYPE.EQ.D_C_FORMAT) THEN
         CALL STETER(15,' wrong column format',ISTAT)
      ENDIF
      IRCOL(1) = ICOL(1)
      IRCOL(2) = ICOL(2)
C
C ... read table + calculate (integrated) flux
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        
            DAT2H = DATA(2)/2.
            LAMLO = DATA(1) - DAT2H       !lower limit in lambda of bin
            LAMHI = DATA(1) + DAT2H       !upper limit in lambda of bin
            IF ((LAMLO.GE.DSTART) .AND. (LAMHI.LE.DEND)) THEN
               CALL DOINTG(MADRID(PNTR),DSTART,DSTEP,NPIX,
     +                     LAMLO,LAMHI,BINFLX)
               CALL TBEWRR(TID,I,ICOL(3),BINFLX,ISTAT)
            ENDIF
         ENDIF
10    CONTINUE
C
C the end
C
      CALL TBTCLO(TID,ISTAT)
      CALL STSEPI
      STOP 
      END

      SUBROUTINE DOINTG(FLUX,DSTART,DSTEP,NPIX,LAMLO,LAMHI,BINFLX)
C
      IMPLICIT NONE
C
      INTEGER NPIX,IPIXHI,IPIXLO,I
C
      REAL FLUX(1),PIXLO,PIXHI,FPIXHI,FPIXLO,BINFLX
C
      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
