C @(#)necrebi.for	19.1 (ESO-IPG) 02/25/03 14:20:25
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 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 ECHREB
C
C+
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  10:37 - 11 JAN 1988
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C
C.IDENTIFICATION
C
C  program  ECHREBI1
C
C.MODIFICATIONS
C
C  910128   P.Ballester   Define and set variables
C
C.KEYWORDS
C
C  ECHELLE, CASPEC, REBIN
C
C.PURPOSE
C
C  OBTAIN A 1D FRM CALIBRATED IN WAVELENGTHS FROM A 2D FRM
C  OF TYPE ECHELLE WITH THE DISPERSION COEFFS DEFINED AS DESCRIPTORS
C
C.ALGORITHM
C
C  EACH LINE IN THE INPUT FRM 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] [TABLE]
C  OR
C  REBIN/ECHELLE INPUT OUTPUT REFER_IMAGE [METHOD] [TABLE]
C
C METHOD CAN BE EITHER "NONLINEAR" - NON LINEAR REBINNING USING
C                                     THE COEFFS IN THE TABLE
C                      "LINEAR"    - LINEAR REBINNING
C
C.MODIFICATIONS
C [2.2]  DEFINE THE SAMPLING DOMAIN VIA A REFER IMAGE
C 
C 010201		last modif
C 
C 
C-
      IMPLICIT NONE
C
      INTEGER      MADRID, ID, IR
      INTEGER      I,I1,IMNO
      INTEGER*8    IPNTRA
      INTEGER*8    IPNTRB
      INTEGER      ISTAT,NAXISA,NAXISB
      INTEGER      PIXBN, KUN,KNUL, NNNN
      INTEGER      NPIXA(3),NPIXB(3),IORD(2),IORD2(2),TID
      INTEGER      IORDER(100),NCOE(2),IA1(100),NPTOT(100)
      DOUBLE PRECISION WSTART(100),STEPA(3),STARTA(3)
      DOUBLE PRECISION STEPB(3),STARTB(3)
      DOUBLE PRECISION WSTEP, A1(7,100)
      REAL         CUT(4),XMIN,XMAX
      CHARACTER*60 OUTFRM,INFRM,REFER,TABLE
      CHARACTER    METHOD*1, LINE*80, DSCNAM*8
      CHARACTER    IDENTA*72,CUNITA*64,CUNITB*64
      LOGICAL      IREF
      INCLUDE      'MID_INCLUDE:ST_DEF.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE      'MID_INCLUDE:ST_DAT.INC'
      DATA CUNITB/'Flux            Relative wavel.'/
      DATA NPIXA /3*1/, STARTA /3*0.0/, STEPA /3*1.0/
      DATA NPIXB /3*1/, STARTB /3*0.0/, STEPB /3*1.0/
C
C ... INITIALIZE SYSTEM
C
      CALL STSPRO('ECHREB')
      CALL STKRDC('IN_A',1,1,60,I,INFRM,KUN,KNUL,ISTAT)
      CALL STKRDC('OUT_A',1,1,60,I,OUTFRM,KUN,KNUL,ISTAT)
      CALL CLNFRA(INFRM,INFRM,0)
      CALL CLNFRA(OUTFRM,OUTFRM,0)
      CALL STKRDC('P3',1,1,60,I,REFER,KUN,KNUL,ISTAT)
      CALL STKRDC('P4',1,1,1,I,METHOD,KUN,KNUL,ISTAT)
      I1     = INDEX('.+1234567890',REFER(1:1))
      IF (I1.NE.0) THEN
          CALL GENCNV(REFER,4,1,NNNN,XMIN,WSTEP,ISTAT)
          IF (ISTAT.LT.1) CALL STETER(11,'Wrong sampling step')
C 
          WRITE(LINE,9010) WSTEP
          CALL STTPUT(LINE,ISTAT)
          IREF   = .FALSE.
      ELSE
          IREF   = .TRUE.
          CALL CLNFRA(REFER,REFER,0)
      END IF
      IF (METHOD.EQ.'N' .OR. METHOD.EQ.'n') THEN
          CALL STKRDC('IN_B',1,1,60,I,TABLE,KUN,KNUL,ISTAT)
          CALL CLNTAB(TABLE,TABLE,0)
          CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT)
          CALL STDRDI(TID,'ORDER',1,2,I,IORD,KUN,KNUL,ISTAT)
          CALL STDRDI(TID,'COEFS',1,2,I,NCOE,KUN,KNUL,ISTAT)
          CALL STDRDD(TID,'COEFD',1,7*NCOE(2),I,A1,
     .                KUN,KNUL,ISTAT)
          CALL STDRDI(TID,'COEFI',1,NCOE(2),I,IA1,
     .                KUN,KNUL,ISTAT)
          CALL TBTCLO(TID,ISTAT)
      END IF
C
C ... MAP INPUT IMAGE
C
      CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .            3,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA,
     .            IPNTRA,IMNO,ISTAT)

      IF (NAXISA.EQ.1 .OR. NPIXA(2).EQ.1) THEN
          NPIXA(2) = 1
          NAXISA = 1
      END IF

      IF (NPIXA(2).GT.100)  THEN
          CALL STETER(9,
     +   'Module ECHRE31: Too many orders (limit 100). Buffer overflow')
      ENDIF

      IF (STEPA(2).NE.1.)  THEN
          CALL STETER(9,'Module ECHREBI: Wrong step(2) .ne. 1 ')
      ENDIF

      IORD2(1) = STARTA(2)                  ! Starting relative order number
      IORD2(2) = STARTA(2) + NPIXA(2) - 1   ! Final relative order number
C
C ... MAP OUTPUT FRM
C
      IF (METHOD.EQ.'N' .OR. METHOD.EQ.'n') THEN
          NAXISB = 2
          IF (IREF) THEN
              CALL STFOPN(REFER,D_R4_FORMAT,0,F_IMA_TYPE,
     .                    IR,ISTAT)
              CALL STDRDI(IR,'NPIX',1,2,I,NPIXB,KUN,KNUL,ISTAT)
              CALL STDRDD(IR,'START',1,2,I,STARTB,KUN,KNUL,ISTAT)
              CALL STDRDD(IR,'STEP',1,2,I,STEPB,KUN,KNUL,ISTAT)
              WSTEP  = STEPB(1)
              CALL STDRDI(IR,'NORDER',1,NPIXB(2),I,IORDER,
     .                    KUN,KNUL,ISTAT)
              CALL STDRDI(IR,'NPTOT',1,NPIXB(2),I,NPTOT,
     .                    KUN,KNUL,ISTAT)
              CALL STDRDD(IR,'WSTART',1,NPIXB(2),I,WSTART,
     .                    KUN,KNUL,ISTAT)
              DO 10 I = 1,NPIXB(2)
                  WSTART(I) = WSTART(I)/1000.
   10         CONTINUE
          ELSE
              CALL WRANGE(NCOE(2),A1,IA1,WSTEP,IORD,NPIXA(1),WSTART,
     .                    NPIXB(1),NPTOT,IORDER,STARTA,STEPA,IORD2)
              NPIXB(2)  = NPIXA(2)
              STARTB(1) = 0.D0
              STARTB(2) = STARTA(2)
              STEPB(1)  = WSTEP
              STEPB(2)  = STEPA(2)
          END IF

          CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .                NAXISB,NPIXB,STARTB,STEPB,
     .                IDENTA,CUNITB,IPNTRB,ID,ISTAT)
C
C ... NON LINEAR REBINNING
C
          CALL ECHRE3(MADRID(IPNTRA),NPIXA(1),NPIXA(2),STARTA,STEPA,
     .                NCOE(2),A1,IA1,IORD,MADRID(IPNTRB),NPIXB(1),
     .                NPIXB(2),WSTART,WSTEP,IORD2)
C
C ... WRITE PROCESS DESCRIPTORS
C
          NNNN = NPIXB(1)*NPIXB(2)
          CALL MNMX(MADRID(IPNTRB),NNNN,XMIN,XMAX)
          CUT(1) = XMIN
          CUT(2) = XMAX
          CUT(3) = XMIN
          CUT(4) = XMAX
          DO 20 I = 1,NPIXB(2)
              WSTART(I) = 1000.*WSTART(I)
   20     CONTINUE
          CALL STDCOP(IMNO,ID,3,DSCNAM,ISTAT)
          CALL STDWRD(ID,'WSTART',WSTART,1,NPIXB(2),KUN,ISTAT)
          CALL STDWRI(ID,'NPTOT',NPTOT,1,NPIXB(2),KUN,ISTAT)
          CALL STDWRI(ID,'NORDER',IORDER,1,NPIXB(2),KUN,ISTAT)
          CALL STDWRR(ID,'LHCUTS',CUT,1,4,KUN,ISTAT)
      ELSE
C ... LINEAR REBINNING
          IF (NPIXA(2).GT.1) THEN
CCCCCCCCCCCCC              CALL STDFND(IMNO,'WSTART',I1,I2,I3,ISTAT)
C              IF (TYPE(1:1).EQ.' ') THEN
C                  CALL STTPUT(' Wrong input image ... ',ISTAT)
C                  GO TO 40
CCCCCCCCCCCCC              END IF
              CALL STDRDD(IMNO,'WSTART',1,NPIXA(2),I,WSTART,
     .                    KUN,KNUL,ISTAT)
              CALL STDRDI(IMNO,'NPTOT',1,NPIXA(2),I,NPTOT,
     .                    KUN,KNUL,ISTAT)
          END IF
          CALL STDRDR(IMNO,'LHCUTS',1,4,I,CUT,KUN,KNUL,ISTAT)
          PIXBN     = STEPA(1)/WSTEP
          NAXISB    = NAXISA
          NPIXB(1)  = NPIXA(1)*PIXBN
          NPIXB(2)  = NPIXA(2)
          STARTB(1) = STARTA(1)
          STARTB(2) = STARTA(2)
          STEPB(1)  = WSTEP
          STEPB(2)  = STEPA(2)
          IF (NPIXB(2).GT.1) THEN
              DO 30 I = 1,NPIXB(2)
                  NPTOT(I) = NPTOT(I)*PIXBN
   30         CONTINUE
          ELSE
              WSTART(1) = STARTB(1)
              NPTOT(1) = NPIXB(1)
          END IF
          CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .                NAXISB,NPIXB,STARTB,STEPB,IDENTA,
     .                CUNITB,IPNTRB,ID,ISTAT)
          CALL ECHR01(MADRID(IPNTRA),NPIXA(1),NPIXA(2),STARTA,STEPA,
     .                NPTOT,WSTART,
     .                MADRID(IPNTRB),NPIXB(1),NPIXB(2),WSTEP)

          CALL STDCOP(IMNO,ID,3,DSCNAM,ISTAT)
          IF (NPIXB(2).GT.1) THEN
            CALL STDWRD(ID,'WSTART',WSTART,1,NPIXB(2),KUN,ISTAT)
            CALL STDWRI(ID,'NPTOT',NPTOT,1,NPIXB(2),KUN,ISTAT)
          END IF
          CALL STDWRR(ID,'LHCUTS',CUT,1,4,KUN,ISTAT)
      END IF
      CALL STSEPI
C10000 FORMAT(G<I1>.0)
 9000 FORMAT (G6.0)
 9010 FORMAT(' Sampling step = ',F10.4)
      END

      SUBROUTINE ECHRE3(X,NPIXA1,NPIXA2,STARTA,STEPA,NO,A,IA,IORD,
     +                  Y,NPIXB1,NPIXB2,YSTR,YSTP,IORD2)
C
C REBIN THE ORDERS
C
      IMPLICIT NONE
C
      INTEGER NO, NPIXB1, NPIXB2
      INTEGER NPIXA1,NPIXA2,IORD(2),IA(NO),IORD2(2)
      REAL X(NPIXA1,NPIXA2)
      REAL Y(NPIXB1,NPIXB2)
      DOUBLE PRECISION YSTR(1), YSTP, STARTA(2), STEPA(2)
      DOUBLE PRECISION A(7,NO),B1(3)
C
      INTEGER KMIN, KMAX, I0, I2, I1, K1, II, MIN, MAX, I3
      DOUBLE PRECISION WLSTP, WLST2
C
C ... GET LIMITS FOR THE FIRST ORDER (LOWER WAVELENGTHS)
C
      KMIN   = MIN(IORD(1),IORD(2))
      KMAX   = MAX(IORD(1),IORD(2))
      IF (IORD(1).LT.IORD(2)) THEN
          I0     = NO
          I2     = -1
      ELSE
          I0     = 1  
          I2     = 1
      END IF
      I1     = I0
      WLSTP  = YSTP/1000.
C
C ... ITERATION ON ORDERS
C
      DO 20 K1 = IORD2(2), IORD2(1), -1
          I3 = I1 + IORD2(1) - 1
          IF (IA(I3).LT.0) THEN
              B1(1)  = A(1,I3)
              B1(2)  = -A(2,I3)
              B1(3)  = A(3,I3)
          END IF
          WLST2  = YSTR(I1)
          IF (IA(I3).LT.0) THEN
              CALL LREBIN(X(1,I1),NPIXA1,STARTA,STEPA,B1,3,Y(1,I1),
     +                    NPIXB1,WLST2,WLSTP)
C
C ... SIGN CHANGE
              DO 10 II = 1,NPIXB1
                  Y(II,I1) = -Y(II,I1)
   10         CONTINUE
          ELSE
              CALL LREBI1(X(1,I1),NPIXA1,STARTA,STEPA,A(1,I3),IA(I3),
     +                     Y(1,I1),NPIXB1,WLST2,WLSTP)
          END IF
          I1     = I1 + I2
          IF (I1.GT.NO) RETURN
   20 CONTINUE
      RETURN
      END

      DOUBLE PRECISION FUNCTION POLEV(IA,A,XX)
C
C
C
      IMPLICIT NONE
C
      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
      REAL     X(N), XMIN, XMAX
C
      INTEGER I
C
      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(NCOE,A,IA,WSTEP,IORD,NPIXA,WSTART,NPIXB,NPTOT,
     +                  IORDER,STARTA,STEPA,IORD2)
C
C GET STARTING POSITIONS FOR EACH ORDER AND NO. OF PIXELS PER ORDER
C
C
C ... GET LIMITS FOR THE FIRST ORDER (LOWER WAVELENGTHS)
C
      IMPLICIT NONE
C
      INTEGER NCOE, NPIXA, NPIXB
      DOUBLE PRECISION A(7,NCOE),POLEV,B(3), STARTA(1), STEPA(1)
      DOUBLE PRECISION WSTART(1), WSTEP, WLST1, WLEN1
      INTEGER IA(NCOE),IORD(2),NPTOT(1),IORDER(1), IORD2(2)
C
      INTEGER KMIN, KMAX, I0, I1, I2, I, ISTEP, K, N1, I3
c      INTEGER*8 N1
      DOUBLE PRECISION WS1, WW1, IX_IR8
C
      KMIN   = MIN(IORD(1),IORD(2))
      KMAX   = MAX(IORD(1),IORD(2))
      IF (IORD(1).LT.IORD(2)) THEN
          I0     = NCOE 
          I2     = -1
      ELSE
          I0     = 1 
          I2     = 1
      END IF
      I1     = I0
      NPIXB  = 0
      WS1    = WSTEP/1000.
C
C ... get order number
C
      ISTEP  = 1
      IF (IORD(1).GT.IORD(2)) ISTEP  = -1 ! Default case in present versions
      IORDER(1) = IORD(1) - IORD2(1) + 1
      DO 10 I = 2, IORD2(2) - IORD2(1) + 1
          IORDER(I) = IORDER(I-1) + ISTEP
   10 CONTINUE
C
C ... iteration on orders to get starting wavelength and no. of pixels
C
      DO 20 K = IORD2(2), IORD2(1), -1
          I3 = I1 + IORD2(1) - 1
          IF (IA(I3).LT.0) THEN
              B(1)   = A(1,I3)
              B(2)   = -A(2,I3)
              B(3)   = A(3,I3)
              CALL ECHORD(B,3,NPIXA,WLST1,WLEN1,STARTA,STEPA)
          ELSE
              WW1    = POLEV(IA(I3),A(1,I3),IX_IR8(1,STARTA,STEPA))
              WLST1  = WW1
              WW1    = POLEV(IA(I3),A(1,I3),IX_IR8(NPIXA,STARTA,STEPA))
              WLEN1  = WW1
          END IF
          WLST1  = WLST1*10.
          WLEN1  = WLEN1*10.
          N1     = (WLEN1-WLST1)/WS1
          IF (N1.GT.65534) THEN !added swolf@eso.org
             N1 = 0             !if WLEN1 and WLST1 are NaN N1 gets wrong!
          ENDIF                 !--> NPIXB gets wrong, too!
          NPTOT(I1) = N1
          NPIXB  = MAX(NPIXB,N1)
          IF (I1.EQ.I0) THEN
              WSTART(I1) = WS1*NINT(WLST1/WS1) !before: WLST1 (swolf@eso.org)
          ELSE
              WSTART(I1) = WSTART(I0) + WS1*NINT((WLST1-WSTART(I0))/WS1)
          END IF
          I1     = I1 + I2
          IF (I1.GT.NCOE) GO TO 30
   20 CONTINUE
   30 NPIXB  = NPIXB + 3
      RETURN
      END

      SUBROUTINE ECHR01(A,NPIXA1,NPIXA2,STARTA,STEPA,NPTOT,WSTART,
     +                  B,NPIXB1,NPIXB2,WSTEP)
C
C REBIN ECHELLE ORDERS
C
      IMPLICIT NONE
C
      INTEGER NPIXA1,NPIXA2,NPIXB1,NPIXB2
      DOUBLE PRECISION  STEPA(2),STARTA(2)
      REAL A(NPIXA1,NPIXA2),B(NPIXB1,NPIXB2)
      DOUBLE PRECISION  WSTART(NPIXA2), WSTEP
      INTEGER NPTOT(NPIXA2)
C
      INTEGER NORD, IORD, I, NPREB
      DOUBLE PRECISION XSTART, XP, ADEX, PIXBN
      REAL   VBIN
C
      PIXBN  = WSTEP/STEPA(1)
      NORD   = NPIXA2
      DO 30 IORD = 1,NORD
          XSTART = WSTART(IORD)
          NPREB  = NPTOT(IORD)
          DO 10 I = 1,NPREB
              XP     = XSTART + (I-1)*WSTEP
              ADEX   = (XP-XSTART)/STEPA(1) + 1.
              B(I,IORD) = VBIN(A(1,IORD),PIXBN,ADEX)
   10     CONTINUE
          DO 20 I = NPREB + 1,NPIXB1
              B(I,IORD) = 0.
   20     CONTINUE
   30 CONTINUE
      RETURN
      END

      REAL FUNCTION VBIN(ARRAY,PIXPBN,ADEX)
C
C FUNCTION TO RETURN A VALUE FOR A REBINNED BIN
C
C INPUT ARGUMENTS
C ARRAY REAL*4 ARRAY WITH REAL ORIGINAL PIXELS
C PIXPBN REAL*4 PIXELS PER BIN, NUMBER OF ORIGINAL
C   PIXELS IN THE NEW BIN, WITH FRACTIONS
C ADEX REAL*4 CENTER INDEX IN THE ARRAY FOR THE NEW BIN,
C   WITH FRACTIONS
C
      IMPLICIT NONE
      INTEGER          IBOTDX,ITOPDX,IENLP,I,ISTLP
      REAL             HBIN,BOTDX,TOPDX,SUM,BOTFRA,TOPFRA
      REAL             ARRAY(1)
      DOUBLE PRECISION PIXPBN, ADEX
C  !   HALF BIN
      HBIN   = PIXPBN*0.5  !   BOTTOM POINTER, FRACTION
      BOTDX  = ADEX - HBIN  !   TOP POINTER, FRACTION
      TOPDX  = ADEX + HBIN  !   BOTTOM POINTER, WITHOUT FRACTION
      IBOTDX = BOTDX  !   TOP POINTER, WITHOUT FRACTION
      ITOPDX = TOPDX
C
      IF (IBOTDX.EQ.ITOPDX) THEN
C
C   BIN SMALLER THAN ORIGINAL PIXELS
C   AND NEW BIN WITHIN AN ORIGINAL PIXEL
C
          SUM    = PIXPBN*ARRAY(IBOTDX+1)
      ELSE
C
C   BIN LARGER THAN ORIGINAL PIXELS, OR
C   SMALL BIN OVER TWO PIXELS
C
          BOTFRA = 1. - (BOTDX-FLOAT(IBOTDX))
          TOPFRA = TOPDX - FLOAT(ITOPDX)  !   START LOOP INDEX
          ISTLP  = IBOTDX + 2
          IENLP  = ITOPDX
          IBOTDX = IBOTDX + 1
          ITOPDX = ITOPDX + 1
C
          SUM    = 0.
          IF (IENLP.GE.ISTLP) THEN
              DO 10 I = ISTLP,IENLP
                  SUM    = SUM + ARRAY(I)
   10         CONTINUE
          END IF  !   LOWER FRACTION
          SUM    = SUM + BOTFRA*ARRAY(IBOTDX)  !   UPPER FRACTION
          SUM    = SUM + TOPFRA*ARRAY(ITOPDX)
      END IF
C
C AREA CORRECTION
C
      VBIN   = SUM/PIXPBN
      RETURN
      END

      SUBROUTINE ECHORD(AC,NC,NPIX,W1,W2,STARTA,STEPA)
C
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  22:34 - 3 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C.KEYWORDS:
C  SPECTROSCOPY, CASPEC, ORDER DEFINITION
C
C.PURPOSE:
C  GET WAVELENGTH RANGE OF A GIVEN ORDER
C
C.ALGORITHM:
C
C.INPUT/OUTPUT
C  INPUT ARGUMENTS
C
C  AC(NC) REAL*8 COEFFICIENTS OF THE TRANSFORMATION
C  NC  REAL*4 NO. OF COEFFICIENTS
C  NPIX  INTG*4 NO OF SAMPLES
C
C  OUTPUT ARGUMENTS
C
C  W1  REAL*4  LOWER LIMIT
C  W2  REAL*4  UPPER LIMIT
C
C-------------------------------------------------------------
C
C
      IMPLICIT NONE
      INTEGER NPIX, NC
      DOUBLE PRECISION AC(NC), STARTA(1), STEPA(1)
      DOUBLE PRECISION W1, W2, IX_IR8
C
      DOUBLE PRECISION A0,A1,A2,F1,FN
C
C
C ... ASSIGN COEFFS
C
C      F1     = 1.D0
C      FN     = NPIX

      F1     = IX_IR8(1,STARTA(1),STEPA(1))
      FN     = IX_IR8(NPIX,STARTA(1),STEPA(1))

      A0     = AC(1)
      A1     = AC(2)
      IF (NC.GE.3) THEN
          A2     = AC(3)
          W1     = (-A1+DSQRT(A1*A1-4.D0*A2* (A0-F1)))/ (2.D0*A2)
          W2     = (-A1+DSQRT(A1*A1-4.D0*A2* (A0-FN)))/ (2.D0*A2)

      ELSE
          W1     = (F1-A0)/A1
          W2     = (FN-A0)/A1
      END IF
      RETURN
      END


      SUBROUTINE LREBI1(X,NX,XSTR,XSTP,AC,NC,Y,NY,YSTR,YSTP)
C
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  0:22 - 4 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C.IDENTIFICATION
C  SUBROUTINE  LREBIN1 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
C
      IMPLICIT NONE
      INTEGER  NX, NY, NC
      REAL X(NX),Y(NY)
      DOUBLE PRECISION YSTR, YSTP, XSTR, XSTP
C
      DOUBLE PRECISION AC(NC),A0,A1,A2,A3,A4,A5,A6,FN,FL,WPOLY
      DOUBLE PRECISION FACTOR,DPOLY,CNORM, IX_R8R8, XN, XNP1
      DOUBLE PRECISION PN,PNP1,TM,TMP1,C,COV,VAL, COVX
      DOUBLE PRECISION B2, B3, B4, B5, B6, DELT, STARTL
      INTEGER          M1, MST, N, K, MPT
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
      CNORM  = 10.D0/YSTP

C
C ... CHECK FOR OUTPUT START
C
      STARTL = YSTR - 0.5*YSTP
      FN     = 0.5D0
      XN     = IX_R8R8(FN,XSTR,XSTP)
      PN     = 10.D0*WPOLY(XN)

      DELT   = PN - STARTL
      IF (DELT.GT.0.) THEN
          M1     = DELT/YSTP + 1.
          TMP1   = PN + YSTP
          TM     = PN
      ELSE
          M1     = 1
          TMP1   = STARTL + YSTP
          TM     = STARTL
      END IF
C
C ... CHECK FOR INPUT START
C
      DO 10 N = 1,NX
          FN     = FN + 1.D0
          XNP1   = IX_R8R8(FN,XSTR,XSTP)
          PNP1   = 10.D0*WPOLY(XNP1)

          IF (PNP1.GT.TM) GO TO 20
          PN     = PNP1
   10 CONTINUE
      N      = NX + 1
      GO TO 30
   20 C      = X(N)
   30 VAL    = 0.

      COVX   = XNP1 - XN
      COV    = (PNP1 - PN)/COVX
C
C ... REBINNING LOOP
C
      MPT    = 1
      MST    = 1
      K      = 0
   35 CONTINUE
          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
          XN     = XNP1

          FN     = FN + 1.D0
          XNP1   = IX_R8R8(FN,XSTR,XSTP)
          PNP1   = 10.D0*WPOLY(XNP1)

          C      = X(N)
          COV    = (PNP1 - PN)/(XNP1-XN)

          GO TO 40
   50     TM     = TMP1
          TMP1   = TMP1 + YSTP
   60     FACTOR = CNORM*DPOLY(IX_R8R8(FN,XSTR,XSTP))
C          FACTOR = FACTOR*10.D0*DABS(DPOLY(IX_R8R8(FN,XSTR,XSTP)))/YSTP
          Y(MPT) = VAL*FACTOR
          MPT    = MPT + 1
          VAL    = 0.
          GOTO 35
   80 CONTINUE
   90 CONTINUE
      DO 100 K = MPT + 1,NY
          Y(K)   = 0.
  100 CONTINUE
      RETURN
      END

      SUBROUTINE LREBIN(X,NX,XSTR,XSTP,AC,NC,Y,NY,YSTR,YSTP)
C
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  19:08 - 6 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: 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 WAVELENGTH INTO
C  SAMPLE NUMBER IS DEFINED BY
C  S = A0+A1*T+A2*T**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
C
      IMPLICIT NONE
      INTEGER NX, NY, NC, N, M1
      REAL X(NX),Y(NY)
      DOUBLE PRECISION AC(NC),A0,A1,A2,A3,A4,A5,FN
      DOUBLE PRECISION FACTOR,CNORM, START, PN, DELT, STARTL
      DOUBLE PRECISION TM, TMP1
      DOUBLE PRECISION IX_R8R8
      DOUBLE PRECISION XSTR, XSTP, YSTR, YSTP
C
      INTEGER K, MPT, MST
      DOUBLE PRECISION PNP1, C, VAL, COV, DSQRT, DMIN1, DMAX1
      DOUBLE PRECISION XN, XNP1, COVX
C
C WPOLY(FL) = FL*(A1+FL*(A2+FL*(A3+FL*(A4+FL*A5))))
C
C ... ASSIGN COEFFS
C
      A0     = AC(1)
      A1     = -AC(2)
      A2     = AC(3)
      A3     = A1*A1
      A4     = 4.D0*A2
      A5     = 2.D0*A2
      CNORM  = 10.D0/ (A5*YSTP)
C A0 = AC(1)
C A1 = AC(2)
C A2 = AC(3)
C
C ... CHECK FOR OUTPUT START
C
      START  = YSTR - 0.5*YSTP
C STARTL= START - A0
      STARTL = START
      M1     = 1
      FN     = 0.5D0
      XN     = IX_R8R8(FN,XSTR,XSTP)

C PN    = WPOLY(IX_R8R8(FN,XSTR,XSTP))
      PN     = 10.D0* (A1+DSQRT(A3-A4* (A0-XN)))/A5
      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
          XNP1   = IX_R8R8(FN,XSTR,XSTP)

C   PNP1 = WPOLY(IX_R8R8(FN,XSTR,XSTP))
          PNP1   = 10.D0* (A1+DSQRT(A3-A4* (A0-XNP1)))/A5
          IF (PNP1.GT.TM) GO TO 20
          PN     = PNP1
          XN     = XNP1
   10 CONTINUE
      N      = NX + 1
      GO TO 30
   20 C      = X(N)
   30 VAL    = 0.
      COVX   = XNP1 - XN
      COV    = (PNP1 - PN)/COVX
C
C ... REBINNING LOOP
C
      MPT    = 1
      MST    = 1
      K      = 0
C   35 CONTINUE
          IF ( .NOT. (MPT.LE.NY)) GO TO 80
          K      = K + 1
          IF ((K.LT.M1) .OR. (N.GT.NX)) GO TO 60
   40     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
          XN     = XNP1
          FN     = FN + 1.D0
          XNP1   = IX_R8R8(FN,XSTR,XSTP)
C         PNP1 = WPOLY(IX_R8R8(FN,XSTR,XSTP))
          PNP1   = 10.D0* (A1+DSQRT(A3-A4* (A0-XNP1)))/A5

          C      = X(N)

          COVX   = XNP1 - XN
          COV    = (PNP1 - PN)/COVX

          GO TO 40
   50     TM     = TMP1
          TMP1   = TMP1 + YSTP
   60     FACTOR = CNORM*(DSQRT(A3-A4* (A0-XNP1))-
     +             DSQRT(A3-A4* (A0-XNP1-1.D0)))
          Y(MPT) = VAL*FACTOR
          MPT    = MPT + 1
          VAL    = 0.
   80 CONTINUE
      DO 90 K = MPT + 1,NY
          Y(K)   = 0.
   90 CONTINUE
      RETURN
      END

