C @(#)sperebin.for	19.1 (ES0-DMD) 02/25/03 14:29:43
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 SPEREB
C
C+
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  17:23 - 6 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.IDENTIFICATION
C
C  program  SPEREBIN  version 2.0 840314
C  J.D.Ponz   ESO - Graching
C
C.KEYWORDS
C
C  SPECTRA, REBIN
C
C.PURPOSE
C
C  CONVERT FROM PIXEL INTO WAVELENGTH SAMPLING SPACE
C
C.ALGORITHM
C
C  EACH LINE IN THE INPUT FRAME IS REBINNED IN LINEAR STEP IN WAVELENGTH
C  ACCORDING TO THE DISPERSION COEFFS.
C
C.INPUT/OUTPUT
C
C  REBIN/ECHELLE INPUT OUTPUT STEP [METHOD]
C
C METHOD CAN BE EITHER "NONLINEAR" - NON LINEAR REBINNING USING
C                                     THE COEFFS IN THE TABLE
C-
      IMPLICIT NONE
C
      INTEGER   MADRID,KUN,KNUL
      INTEGER   IA1,IAV,II,II1
      INTEGER*8 IPNTRA,IPNTRB
      INTEGER   ISTAT
      INTEGER   JJ,NAXISA,NAXISB,NPP,IMIN,IMOUT
      INTEGER   NPIXA(2),NPIXB(2),NCOE(2)
      INTEGER   IMF,TID
C
      REAL WSTART,WSTEP
      REAL STEPA(2),STARTA(2),STEPB(2),STARTB(2)
      REAL CUT(4),XMIN,XMAX
C
      DOUBLE PRECISION A1(10),DSTEP(2),DSTART(2)
C
      CHARACTER*80 OUTFR,INFR
      CHARACTER*8  REDSTA
      CHARACTER    HISTOR*80
      CHARACTER    IDENTA*72,CUNITA*64,CUNITB*64
      CHARACTER*60 TABLE
      CHARACTER*80 REFER,TABLEF
C
      INCLUDE    'MID_INCLUDE:ST_DEF.INC/NOLIST'
      COMMON     /VMR/MADRID(1)
      INCLUDE    'MID_INCLUDE:ST_DAT.INC/NOLIST'
C
      DATA CUNITB/'FLUX            WAVELENGTH'/
C
C ... INITIALIZE SYSTEM
C
      CALL STSPRO('SPEREBIN')
C
C ... READ PARAMETERS
C
      CALL STKRDC('P1',1,1,80,IAV,INFR,KUN,KNUL,ISTAT)
      CALL STKRDC('P2',1,1,80,IAV,OUTFR,KUN,KNUL,ISTAT)
      CALL STKRDC('P3',1,1,80,IAV,REFER,KUN,KNUL,ISTAT)
      II     = INDEX('.0123456789',REFER(1:1))
      IF (II.EQ.0) THEN
          CALL STFOPN(REFER,D_R4_FORMAT,0,F_IMA_TYPE,IMF,ISTAT)
          CALL STDRDD(IMF,'START',1,1,IAV,DSTART,KUN,KNUL,ISTAT)
          CALL STDRDD(IMF,'STEP',1,1,IAV,DSTEP,KUN,KNUL,ISTAT)
          CALL STDRDI(IMF,'NPIX',1,1,IAV,NPP,KUn,KNUL,ISTAT)
          WSTART = DSTART(1)
          WSTEP  = DSTEP(1)
      ELSE
          READ (REFER,9000) WSTEP
          WSTART = -1000.
      END IF
C
C     CALL STKRDR('INPUTR',1,1,IAV,WSTEP,KUN,KNUL,STAT)
      CALL STKRDC('MID$CMND',1,1,72,IAV,HISTOR,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_A',1,1,60,IAV,TABLE,KUN,KNUL,ISTAT)
      II1    = INDEX(TABLE,' ') - 1
      TABLEF = TABLE(1:II1)//'.TBL'

      CALL TBTOPN(TABLEF,F_I_MODE,TID,ISTAT)
      CALL STDRDI(TID,'COEFS',1,2,IAV,NCOE,KUN,KNUL,ISTAT)
      CALL STDRDD(TID,'COEFD',1,10,IAV,A1,KUN,KNUL,ISTAT)
      CALL STDRDI(TID,'COEFI',1,1,IAV,IA1,KUN,KNUL,ISTAT)
      DO 10 JJ = IA1 + 1,10
          A1(JJ) = 0.D0
   10 CONTINUE
C
C ... MAP INPUT FRAME
C
      CALL STIGET(INFR,10,0,1,2,NAXISA,NPIXA,DSTART,DSTEP,
     +     IDENTA,CUNITA,IPNTRA,IMIN,ISTAT)
      STARTA(1) = DSTART(1)
      STARTA(2) = DSTART(2)
      STEPA(1) = DSTEP(1)
      STEPA(2) = DSTEP(2)
      IF (NAXISA.EQ.1 .OR. NPIXA(2).EQ.1) THEN
          NPIXA(2) = 1
          NAXISA = 1
      ELSE
          NAXISA = 2
      END IF
C
C ... MAP OUTPUT FRAME
C
      IF (WSTART.LE.0.0) THEN
          CALL WRANGE(A1,IA1,WSTEP,NPIXA(1),STARTB(1),NPIXB(1))
      ELSE
          STARTB(1) = WSTART
          NPIXB(1) = NPP
      END IF
      NAXISB    = NAXISA
      NPIXB(2)  = NPIXA(2)
      STARTB(2) = STARTA(2)
      STEPB(1)  = WSTEP
      STEPB(2)  = STEPA(2)
      REDSTA = 'EXTRACT'
      DSTART(1) = STARTB(1)
      DSTART(2) = STARTB(2)
      DSTEP(1) = STEPB(1)
      DSTEP(2) = STEPB(2)
      CALL STIPUT(OUTFR,10,1,1,NAXISB,NPIXB,DSTART,DSTEP,
     +     IDENTA,CUNITB,IPNTRB,IMOUT,ISTAT)
C
C ... NON LINEAR REBINNIG
C
      CALL REBINB(MADRID(IPNTRA),NAXISA,NPIXA(1),NPIXA(2),
     +     STARTA,STEPA,A1,IA1,MADRID(IPNTRB),NPIXB(1),
     +     NPIXB(2),STARTB,STEPB)
C
C ... WRITE PROCESS DESCRIPTORS
C
      CALL MNMX(MADRID(IPNTRB),NPIXB(1)*NPIXB(2),XMIN,XMAX)
      CUT(1) = XMIN
      CUT(2) = XMAX
      CUT(3) = XMIN
      CUT(4) = XMAX
      CALL STDCOP(IMIN,IMOUT,3,' ',ISTAT)
      CALL STDWRC(IMOUT,'REDSTAT',1,REDSTA,1,8,KUN,ISTAT)
      CALL STDWRC(IMOUT,'HISTORY',1,HISTOR,-1,72,KUN,ISTAT)
      CALL STDWRR(IMOUT,'LHCUTS',CUT,1,4,KUN,ISTAT)
      DSTART(1) = STARTB(1)
      DSTART(2) = STARTB(2)
      CALL STDWRD(IMOUT,'START',DSTART,1,2,KUN,ISTAT)
C      WRITE (*,*) STARTB(1)
      CALL STSEPI
      STOP 
 9000 FORMAT (F8.4)
      END

      DOUBLE PRECISION FUNCTION POLEV(IA,A,XX)
C
C
C
      IMPLICIT NONE
      INTEGER IA
      DOUBLE PRECISION A(IA),XX,YY

      YY     = A(1) + XX* (A(2)+XX* (A(3)+XX* (A(4)+XX* (A(5)+XX* (A(6)+
     +         XX*A(7))))))
      POLEV  = YY
      RETURN
      END

      SUBROUTINE MNMX(X,N,XMIN,XMAX)
C
C FIND MIN MAX VALS
C
      IMPLICIT NONE
      INTEGER  N,I
      REAL     XMIN,XMAX
      REAL X(N)

      XMIN   = X(1)
      XMAX   = XMIN
      DO 10 I = 2,N
          XMIN   = AMIN1(XMIN,X(I))
          XMAX   = AMAX1(XMAX,X(I))
   10 CONTINUE
      RETURN

      END

      SUBROUTINE WRANGE(A,IA,WSTEP,NPIXA,WSTART,NPIXB)
C
C GET STARTING POSITIONS AND NO. OF PIXELS
C
      IMPLICIT NONE
      INTEGER    IA,NPIXA,NPIXB
      REAL       WSTEP
      DOUBLE PRECISION WW,A(1),POLEV
      REAL WSTART
C
      WW     = POLEV(IA,A,1.D0)
      WSTART = WW
      WW     = NPIXA
      WW     = POLEV(IA,A,WW)
      NPIXB  = (WW-WSTART)/WSTEP + 1
      RETURN

      END

      SUBROUTINE LREBI2(X,NX,XSTR,XSTP,AC,NC,Y,NY,YSTR,YSTP)
C
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C  SUBROUTINE  LREBI2 VERSION 1   22 MAR 1983
C  J.D.PONZ
C
C.KEYWORDS:
C  SPECTROSCOPY, CASPEC, REBINNING
C
C.PURPOSE:
C  DATA REBINNING
C  CONVERSION FROM SAMPLES INTO WAVELENGTHS
C
C.ALGORITHM:
C  LINEAR INTERPOLATION. THE TRANSFORMATION SAMPLE NUMBER INTO
C  WAVELENGTH IS DEFINED BY
C  W = A0 + A1*S + A2*S**2 + ...
C
C.INPUT/OUTPUT
C  INPUT ARGUMENTS
C
C  X(NX) REAL*4 INPUT ARRAY
C  NX  INTG*4 NO OF SAMPLES
C  XSTR  REAL*4 STARTING VALUE OF INPUT
C  XSTP  REAL*4 INPUT STEP
C  AC(NC) REAL*8 COEFFICIENTS OF THE TRANSFORMATION
C  NC  REAL*4 NO. OF COEFFICIENTS
C  NY  INTG*4 NO. OF OUTPUT VALUES
C  YSTR  REAL*4 STARTING VALUE OF THE INPUT
C  YSTP  REAL*4 STEP OF THE OUTPUT
C
C  OUTPUT ARGUMENTS
C
C  Y(NY) REAL*4 OUTPUT ARRAY
C
C
C
      IMPLICIT NONE
      INTEGER      NX,NY,NC,M1,N,MPT,MST,K,KZ6904
      REAL         XSTR,XSTP,YSTR,YSTP,DELT
      REAL         B2,B3,B4,B5,B6,STARTL,START
      REAL X(NX),Y(NY)
      DOUBLE PRECISION AC(NC),A0,A1,A2,A3,A4,A5,A6,FN,FL,WPOLY
      DOUBLE PRECISION FACTOR,DPOLY,CNORM
      DOUBLE PRECISION PN,PNP1,TM,TMP1,C,COV,VAL
C
      WPOLY(FL) = A0 + FL* (A1+FL* (A2+FL* (A3+FL* (A4+FL* (A5+FL*
     +            A6)))))
      DPOLY(FL) = A1 + FL* (B2+FL* (B3+FL* (B4+FL* (B5+FL*B6))))
C
C ... ASSIGN COEFFS
C
      A0     = AC(1)
      A1     = AC(2)
      A2     = AC(3)
      A3     = AC(4)
      A4     = AC(5)
      A5     = AC(6)
      A6     = AC(7)
      B2     = 2.D0*A2
      B3     = 3.D0*A3
      B4     = 4.D0*A4
      B5     = 5.D0*A5
      B6     = 6.D0*A6
C CNORM = 10.D0/YSTP
      CNORM  = 1.D0/YSTP
C
C ... CHECK FOR OUTPUT START
C
      START  = YSTR - 0.5*YSTP
      STARTL = START
      M1     = 1
      FN     = 0.5D0
C PN    = 10.D0*WPOLY(FN)
      PN     = WPOLY(FN)
      DELT   = PN - STARTL
      IF (DELT.GT.0.) THEN
          M1     = DELT/YSTP + 1.
          TMP1   = STARTL + M1*YSTP
          TM     = TMP1 - YSTP

      ELSE
          TMP1   = STARTL + YSTP
          TM     = STARTL
      END IF
C
C ... CHECK FOR INPUT START
C
      DO 10 N = 1,NX
          FN     = FN + 1.D0
          PNP1   = WPOLY(FN)
C   PNP1 = 10.D0*WPOLY(FN)
          IF (PNP1.GT.TM) GO TO 20
          PN     = PNP1
   10 CONTINUE
      N      = NX + 1
      GO TO 30

   20 C      = X(N)
   30 VAL    = 0.D0
      COV    = PNP1 - PN
C
C ... REBINNING LOOP
C
C      MMSP   = NSP
      MPT    = 1
      MST    = 1
      K      = 0
99    KZ6904 = 1
          IF ( .NOT. (MPT.LE.NY)) GO TO 80
          K      = K + 1
          IF ((K.LT.M1) .OR. (N.GT.NX)) GO TO 60
C ?
   40     IF (PNP1.EQ.PN) GO TO 90
          VAL    = VAL + C* (DMIN1(TMP1,PNP1)-DMAX1(TM,PN))/COV
          IF (TMP1.LT.PNP1) GO TO 50
          N      = N + 1
          IF (N.GT.NX) GO TO 60
          PN     = PNP1
          FN     = FN + 1.D0
C         PNP1 = 10.D0*WPOLY(FN)
          PNP1   = WPOLY(FN)
          C      = X(N)
          COV    = PNP1 - PN
          GO TO 40

   50     TM     = TMP1
          TMP1   = TMP1 + YSTP
   60     FACTOR = CNORM*DPOLY(FN)
          Y(MPT) = VAL*FACTOR
          MPT    = MPT + 1
          VAL    = 0.
      GOTO 99
   80 CONTINUE
   90 CONTINUE
      DO 100 K = MPT + 1,NY
          Y(K)   = 0.
  100 CONTINUE
      RETURN

      END

      SUBROUTINE REBINB(X,NAXISA,NPIXA1,NPIXA2,STARTA,STEPA,
     +           A1,IA1,Y,NPIXB1,NPIXB2,STARTB,STEPB)
C
C intermediate routine to rebin each line
C
      IMPLICIT NONE
      INTEGER      NAXISA,NPIXA1,NPIXA2,IA1,NPIXB1,NPIXB2,I
      REAL STARTA(2),STEPA(2)
      REAL STARTB(2),STEPB(2)
      DOUBLE PRECISION A1(1)
      REAL X(NPIXA1,NPIXA2),Y(NPIXB1,NPIXB2)

      DO 10 I = 1,NPIXA2
          CALL LREBI2(X(1,I),NPIXA1,STARTA,STEPA,A1,IA1,
     +         Y(1,I),NPIXB1,STARTB,STEPB)
   10 CONTINUE
      RETURN
      END
