C @(#)rebirbr.for	19.1 (ES0-DMD) 02/25/03 14:25:23
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 @(#)rebirbr.for	19.1  (ESO)  02/25/03  14:25:23 
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++
C .COPYRIGHT   (C) 1993 European Southern Observatory    
C .IDENT       .prg                                      
C .AUTHORS     Pascal Ballester (ESO/Garching)           
C              Cristian Levin   (ESO/La Silla)           
C .KEYWORDS    Spectroscopy, Long-Slit                   
C .PURPOSE                                               
C .VERSION     1.0  Package Creation  17-MAR-1993        
C -------------------------------------------------------

C @(#)rebirbr.for	1.1 (ESO-IPG) 2/17/93 16:46:20
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  8:58 - 9 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.IDENTIFICATION
C
C  SUBROUTINES TDRBII MODIFIED TO REBIN A 2D FRAME ROW BY ROW
C                                          (M.Pierre Nov. 1989)
C.KEYWORDS
C
C  IMAGES, NONLINEAR REBIN, INTERPOLATION, INTEGRATION
C
C.PURPOSE
C
C  NONLINEAR REBIN IMAGE TO IMAGE IN 2 DIMENSIONS
C
C.ALGORITHM
C
C  Rebin input data YI sampled at XI with binwidth WI to
C  output data YO sampled at XO with binwidth WO.
C
C  Restriction:  Input and output independent variables are
C                assumed to be monotonic increasing or decreasing
C
C  Several FUNCTal relations between input and output independent
C  variables are foreseen, additional ones may be added in dummy call
C  U01 as follows:
C
C SUBROUTINE REB_U01 (DV,P,NP,NPC)
C   REAL*8 DV,P(NP)
C   DV = ...........
C   RETURN
C END
C
C  where DV is input-output value, P(NP) the parameter array and
C  NPC the nuber of active parameters in p(np).
C
C  According to the definition of IMAGES, pixels have constant
C  bin width (STEP) [WI;WO], are not overlapping and have filling
C  factors of exactly one. World coordinates [XI;XO] and flux values
C  [YI;YO] are pertinent
C  to the center of a pixel. In flux conservation mode (the default)
C  the flux is interpreted as the area under a smooth curve between
C  the boundaries of the pixel (i.e. a histogram representation of
C  a smooth curve).
C
C  Detailed flux conservation obtained by REAL*8 spline interpolation
C  and  integration of input data over the projected
C  bin WI for fractions of input pixels and summation of flux in
C  entire input pixels covered within the valid domain of input data
C  i.e. from x(1)-0.5*step to x(n)+0.5*step. Several options available
C  to handle undefined input at the extremes and for NULLVALUE assigned
C  parts of input array.
C
C   Optimized for long arrays by search for valid range, optimum
C  direction or DO LOOPS and updating of loop boundaries - therefore
C  overdoing it for small arrays.
C
C.INPUT
C
C  Keys:  IN_A/C/1/8         input data array
C         IN_B/C/1/8         reference frame or dim of output
C         CFUNC/C/1/8        FUNCTal relation between indep.
C         INPUTD/D/1/12      parameters or FUNCT
C         COPT/C/1/8         options
C
C.OUTPUT
C
C  Key:   OUT_A/C/1/8        output data array
C
C.MODIFS
C M Peron sep 90 : remove stupid limit to the size of the output frame!!!!!
C--------------------------------------------------------------------------
C
       PROGRAMM REBIRBR
      IMPLICIT NONE
      INTEGER   NFPAR
      PARAMETER (NFPAR = 12)
C
      INTEGER    MADRID
      INTEGER    I,IACT
      INTEGER*8  IP1,IP2,IP3,IPNTR
      INTEGER    IS,ISTAT
      INTEGER*8  JP1,JP2,JP3,JPNTR
      INTEGER    NAXISI,NAXISO,KUN,KNUL
      INTEGER    NFPACT,NFUNC
      INTEGER    NI,NINT,NN,NO,NBY,NBO,NCOL,NTOT
      INTEGER NPIXI(3),NPIXO(3),INOP,IMNO,IMNI,IREF
       INTEGER JJ,TID,NROW,NSC,NAC,NAR,ICOL(12),STATUS
      REAL    RMIN, RMAX, CUTS(4)
      DOUBLE PRECISION STEPI(3),STARTI(3)
        REAL X0O,WXO,X0I,WXI
      DOUBLE PRECISION DPARM(NFPAR)
      DOUBLE PRECISION DSTART(3), DSTEP(3),YPOS
      CHARACTER*3 FUN(9)
      CHARACTER*60 INPIMA,OUTIMA,REFIMA,TBCOEFF
      CHARACTER*80 FUNCT,OPTION
      CHARACTER*72 IDENTO,IDENTI
      CHARACTER*80 CUNITO,CUNITI
       LOGICAL NULL(12)
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:TABLED.INC'
C
      DATA FUN/'LIN','POL','INV','EXP','DEX','LOG','DLG','IPO','U01'/

C
C..... GET INTO MIDAS
       CALL STSPRO('REBLL')
C ... get input parameters
C
      CALL STKRDC('OUT_A',1,1,60,IACT,OUTIMA,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_A', 1,1,60,IACT,INPIMA,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_B', 1,1,60,IACT,REFIMA,KUN,KNUL,ISTAT)
      CALL STKRDC('CFUNC',1,1,8,IACT,FUNCT,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_C', 1,1,60,IACT,TBCOEFF,KUN,KNUL,ISTAT)
      CALL STKRDC('COPT',1,1,8,IACT,OPTION,KUN,KNUL,ISTAT)

C
C ... det. interpolation/integration options
C  !   Not implemented
      INOP   = 4
C      IF (OPTION(1:3).EQ.'PIX') INOP   = 1
C      IF (OPTION(1:3).EQ.'LIN') INOP   = 2
C      IF (OPTION(1:3).EQ.'SPG') INOP   = 3
C .......................................... the actual values are:
      CALL FORUPC(OPTION,OPTION)
      IF (OPTION(1:3).EQ.'PIX') INOP   = 2
      IF (OPTION(1:3).EQ.'LIN') INOP   = 3
      IF (OPTION(1:3).EQ.'SPG') INOP   = 1
C .......................................... Bitte, bitte
C
C ... translate FUNCT into FUNCT number
C
      CALL FORUPC(FUNCT,FUNCT)
      NFUNC  = 0
      DO 10 I = 1,9
          IF (FUNCT(1:3).EQ.FUN(I)) NFUNC  = I
   10 CONTINUE
      IF (NFUNC.EQ.0) THEN
          CALL STTPUT(' Specified FUNCT non-existent...',ISTAT)
          GO TO 60
      END IF
C
C..... init the input table of dispersion coefficients
        CALL TBTOPN(TBCOEFF,F_U_MODE,TID,ISTAT)
       CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
       NFPACT = NCOL-1

       DO 40 I = 1,NFPACT
       ICOL(I) = I+1
40       CONTINUE       

C
C ... get parameters form reference image
C
      CALL STFOPN(REFIMA,D_R4_FORMAT,0,F_IMA_TYPE,IREF,ISTAT)
      CALL STDRDI(IREF,'NAXIS',1,1,NN,NAXISO,KUN,KNUL,ISTAT)
      CALL STDRDI(IREF,'NPIX',1,NAXISO,NN,NPIXO,KUN,KNUL,ISTAT)
      CALL STDRDD(IREF,'START',1,NAXISO,NN,DSTART,KUN,KNUL,ISTAT)
      CALL STDRDD(IREF,'STEP',1,NAXISO,NN,DSTEP,KUN,KNUL,ISTAT)
      CALL STDRDC(IREF,'IDENT',1,1,72,NN,IDENTO,KUN,KNUL,ISTAT)
      CALL STDRDC(IREF,'CUNIT',1,1,80,NN,CUNITO,KUN,KNUL,ISTAT)
      CALL STFCLO(IREF,ISTAT)
       NTOT = NPIXO(1)*NPIXO(2)
C
C ... get input image
C
      CALL STIGET(INPIMA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .2,NAXISI,NPIXI,STARTI,STEPI,IDENTI,CUNITI,IPNTR,IMNI,ISTAT)
 
C..... check dimensions
       IF(NPIXI(2).NE.NROW) THEN
       CALL STTPUT('Input table and image have not 
     +       the same number of rows',ISTAT)
       GO TO 50
       ELSE 
       IF(NPIXI(2).NE.NPIXO(2)) THEN
       CALL STTPUT('Input and output images have not 
     +  the same number of rows',ISTAT)
       GO TO 50
             ENDIF
       ENDIF
       CALL TBRRDD(TID,1,1,1,YPOS,NULL,STATUS)
       IF (ABS(YPOS-DSTART(2)).GT.0.001) THEN
        CALL STTPUT(' WARNING:',ISTAT)
       CALL STTPUT('Y_position recorded in the 1st line
     +        of input table is not the same',ISTAT)
        CALL STTPUT('as Y_start of the  output image',ISTAT)
       ENDIF
C ... map output image
C
      CALL STIPUT(OUTIMA,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .NAXISO,NPIXO,DSTART,DSTEP,IDENTO,CUNITO,JPNTR,IMNO,ISTAT)
C

C ... fill the arrays XI,WI and XO,WO
C
      NI     = NPIXI(1)
      NBY    = 8*NI
      CALL TDMGET(NBY,IP1,IS)
      CALL TDMGET(NBY,IP2,IS)
      CALL TDMGET(NBY,IP3,IS)
      X0O     = STARTI(1)
      WXO     = STEPI(1)
C
      NO     = NPIXO(1)
      NBO    = 8*NO
      CALL TDMGET(NBO,JP1,IS)
      CALL TDMGET(NBO,JP2,IS)
      CALL TDMGET(NBO,JP3,IS)
      X0I     = DSTART(1)
      WXI     = DSTEP(1)


C..... DO LOOP ON THE IMAGE ROWS

       DO 100       JJ = 1,NROW

      CALL IMVAL3(NI,X0O,WXO,MADRID(IPNTR),MADRID(IP1),MADRID(IP2),
     +             MADRID(IP3),JJ,NTOT)
      CALL IMVAL2(NO,X0I,WXI,MADRID(JP1),MADRID(JP2))
C
C...... Read dispersion coefficients for row JJ
       CALL TBRRDD(TID,JJ,NFPACT,ICOL,DPARM,NULL,STATUS)
       WRITE(6,*) 'ROW',JJ
C       TYPE *,'COEFF',DPARM(1),DPARM(2),DPARM(3),DPARM(4),DPARM(5)

C  ...  Now, at last, we call the interpolating (and integrating) routines:
C
C  !!!! Watch out for VECI,VECD arrays in INTERP, dimension => NINT !!!!!!
C
      NINT   = 8
      CALL REBMET(NI,MADRID(IP1),MADRID(IP3),MADRID(IP2),NO,
     +             MADRID(JP1),MADRID(JP2),NFUNC,NFPAR,NFPACT,DPARM,
     +             INOP,NINT,MADRID(JP3),RMIN,RMAX)
C
       CALL IMVAL5(NO,JJ,MADRID(JP3),MADRID(JPNTR),NTOT)


C
100       CONTINUE

      CALL TDMFRE(NBY,IP1,IS)
      CALL TDMFRE(NBY,IP2,IS)
      CALL TDMFRE(NBY,IP3,IS)
      CALL TDMFRE(NBO,JP1,IS)
      CALL TDMFRE(NBO,JP2,IS)

C ... write cuts
C
      CUTS(1) = RMIN
      CUTS(2) = RMAX
      CUTS(3) = RMIN
      CUTS(4) = RMAX
      CALL STDWRR(IMNO,'LHCUTS',CUTS,1,4,KUN,ISTAT)
C
C ... free memory
C
 
C ... end
C
      CALL STFCLO(IMNO,ISTAT)
   50 CALL STFCLO(IMNI,ISTAT)
       CALL TBTCLO(TID,ISTAT)
   60 CALL STSEPI
      END

C**************************************************************
      SUBROUTINE IMVAL2(NP,XS,WS,X,W)
C
C  Fill arrays X,W for image world coordinates and pixel size data
C
      IMPLICIT NONE
      INTEGER  I, NP
      REAL     XS, WS
      DOUBLE PRECISION X(NP),W(NP),WSS,XSS
C
      WSS    = WS
      XSS    = XS
      DO 10 I = 1,NP
          X(I)   = XSS + WSS* (I-1)
          W(I)   = WSS
   10 CONTINUE
      RETURN
      END

C------------------------------------------------------------------
      SUBROUTINE IMVAL3(NP,XS,WS,YY,X,W,Y,NL,NTOT)
C
C  Fill arrays X,W for image world coordinates and pixel size data
C  Convert YY(real*4) into Y(real*8)
C
      IMPLICIT NONE
      INTEGER  NP, I,NTOT
       INTEGER NL,IPT
      DOUBLE PRECISION X(NP),W(NP),Y(NP),WSS,XSS
      REAL     XS, WS
      REAL YY(NTOT)
C
      WSS    = WS
      XSS    = XS
      DO 10 I = 1,NP
       IPT = (NL-1)*NP+I
          X(I)   = XSS + WSS* (I-1)
          W(I)   = WSS
          Y(I)   = YY(IPT)
   10 CONTINUE
      RETURN
      END

C----------------------------------------------------------
       SUBROUTINE IMVAL5(N,NL,X,Y,NTOT)
C       FILLS THE BUFFER
       INTEGER I,N,NL,NTOT,NPT

       REAL*4 Y(NTOT),X(N)
       DO 10 I = 1,N
       NPT = (NL-1)*N+I
       Y(NPT) = X(I)
10        CONTINUE
       
       RETURN
       END

 
