C @(#)necconv.for	19.1 (ES0-DMD) 02/25/03 14:20:22
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  21:03 - 3 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C.IDENTIFICATION
C
C  program ECHCONV2.FOR
C
C.PURPOSE
C
C  Execute the command
C  CONVERT/ECHELLE input output sampling_domain function parms option
C
C  input            - input image with lines corresponding to 
C			echelle orders
C  output           - output image with lines corresponding to 
C			echelle orders
C  sampling_domain  - start,step or filename defining start,step
C  function         - rebinning function
C  params           - define the transformation coeffs (def 0.,1.)
C  option           - rebinning method as PIX (def)
C                                         LIN
C                                         SPG
C.ALGORITHM
C
C  Linear rebinning
C
C.KEYWORDS
C
C  echelle, flux calibration, rebin
C
C.INPUT/OUTPUT
C
C  P1 - P5    contain input parameters
C
C.MODIFS
C 901024, M Peron , remove list_directed internal I/O
C 901217, P.Ballester, Skip bad orders of extracted standard star spectrum
C-----------------------------------------------------------
C
C
      PROGRAM ECHCV2
      IMPLICIT NONE
C
      INTEGER     NFPAR
      PARAMETER   (NFPAR=12)
      INTEGER     MADRID
      INTEGER     NPTOT1(100),NPIX1(2),NPIXR(2),IDUM(2)
      INTEGER     NPTOT2(100),NPIX2(2),NORDER(100)
      INTEGER*8   PNTR1, PNTR2
      INTEGER     IMNO, NA, KUN, KNUL, ISTAT
      INTEGER     INOP, I, J, NAXIS, NFPACT, NN, NDIM 
      INTEGER     II, NPMAX, STATUS,NFUNC, IMNO1, IMR
      DOUBLE PRECISION STEP1(2),START1(2),WSTR1(100)
      DOUBLE PRECISION STEP2(2),START2(2),WSTR2(100)
      REAL        CUTS(4),RMIN,RMAX,RDUM(2)
      DOUBLE PRECISION RPAR(2)
      DOUBLE PRECISION DPARM(NFPAR),FPARM(NFPAR)
      CHARACTER*8  IMAGE1,IMAGE2,METHOD,OPTION
      CHARACTER*64 IMAGER
      CHARACTER*72 IDENT
      CHARACTER*80 CUNIT
      CHARACTER*3  FUN(9)
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      COMMON/VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      DATA FUN/'LIN','POL','INV','EXP','DEX','LOG','DLG','IPO','U01'/
C
C ... get into MIDAS
C
      CALL STSPRO('ECHCV2')
C
C ... get input parameters + default
C
      CALL STKRDC('P1',1,1,8,NA,IMAGE1,KUN,KNUL,ISTAT)
      CALL STKRDC('P2',1,1,8,NA,IMAGE2,KUN,KNUL,ISTAT)
      CALL STKRDC('P3',1,1,64,NA,IMAGER,KUN,KNUL,ISTAT)
      CALL STKRDC('P4',1,1,8,NA,METHOD,KUN,KNUL,ISTAT)
      CALL STKRDC('P6',1,1,8,NA,OPTION,KUN,KNUL,ISTAT)
      CALL STKRDD('INPUTD',1,NFPAR,NA,FPARM,KUN,KNUL,ISTAT) 
      INOP   = 4
      IF (OPTION(1:3).EQ.'PIX') INOP   = 1
      IF (OPTION(1:3).EQ.'LIN') INOP   = 2
      IF (OPTION(1:3).EQ.'SPG') INOP   = 3
C
C ... translate FUNCTION into function number
C
      NFUNC  = 0
      DO 10 I = 1,9
          IF (METHOD(1:3).EQ.FUN(I)) NFUNC  = I
   10 CONTINUE
      IF (NFUNC.EQ.0) THEN
          CALL STTPUT(' Specified function non-existent...',ISTAT)
          GO TO 50
      END IF
C
C ... search for number of active parameters and DOUBLE
C
      DO 20 J = NFPAR,1,-1
          NFPACT = J
          IF (FPARM(J).NE.0.0D0) GO TO 30
   20 CONTINUE
   30 DO 40 J = 1,NFPACT
          DPARM(J) = FPARM(J)
   40 CONTINUE
C
C
C ... read input image
C
      CALL STIGET(IMAGE1,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .     2,NAXIS,NPIX1,START1,STEP1,IDENT,CUNIT,
     .     PNTR1,IMNO1,ISTAT)
      IF (NAXIS.GT.1 .AND. NPIX1(2).GT.1) THEN
CCCCCCCCCCCCC          CALL SXDFND(IMAGE1,'WSTART',TYPE,NO,NB,ISTAT)
C          IF (TYPE(1:1).EQ.' ') THEN
C              CALL STTPUT(' Wrong input image ',ISTAT)
C              GO TO 50
C          END IF
          CALL STDRDD(IMNO1,'WSTART',1,NPIX1(2),NN,WSTR1,
     .                KUN,KNUL,ISTAT)
          CALL STDRDI(IMNO1,'NPTOT',1,NPIX1(2),NN,NPTOT1,
     .                KUN,KNUL,ISTAT)
          NDIM   = 2
      ELSE
          NPIX1(2) = 1
          NDIM   = 1
          NPTOT1(1) = NPIX1(1)
          WSTR1(1) = START1(1)
      END IF
      START2(1) = START1(1)
      START2(2) = START1(2)
      STEP2(2) = STEP1(2)
C
C ... find start, step, npix for each order in the output image
C
      II     = INDEX('-+.0123456789',IMAGER(1:1))
      IF (II.EQ.0) THEN
          CALL STFOPN(IMAGER,D_R4_FORMAT,0,F_IMA_TYPE,IMR,ISTAT)
          CALL STDRDI(IMR,'NPIX',1,2,NN,NPIXR,KUN,KNUL,ISTAT)
          IF (NPIXR(2).NE.NPIX1(2)) THEN
              CALL STTPUT(' Error in reference image',ISTAT)
              GO TO 50
          END IF
          CALL STDRDD(IMR,'STEP',1,1,NN,STEP2,KUN,KNUL,ISTAT)
          CALL STDRDD(IMR,'WSTART',1,NPIX1(2),NN,WSTR2,
     .                KUN,KNUL,ISTAT)
          CALL STDRDI(IMR,'NPTOT',1,NPIX1(2),NN,NPTOT2,
     .                KUN,KNUL,ISTAT)
          CALL STDRDI(IMR,'NORDER',1,NPIX1(2),NN,NORDER,
     .                KUN,KNUL,ISTAT)
          NPMAX  = NPIXR(1)
      ELSE
          CALL GENCNV(IMAGER,4,2,IDUM,RDUM,RPAR,ISTAT)
          CALL OUTIMA(NPIX1(2),STEP1,WSTR1,NPTOT1,RPAR(1),RPAR(2),
     +            START2,STEP2,WSTR2,NPTOT2,NPMAX)
      END IF
C
C ... map output image
C
      NPIX2(1) = NPMAX
      NPIX2(2) = NPIX1(2)
      CALL STIPUT(IMAGE2,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,NDIM,
     .            NPIX2,START2,STEP2,IDENT,CUNIT,
     +            PNTR2,IMNO,STATUS)
C
C ... rebin image
C
      CALL APPREB(MADRID(PNTR1),NPIX1(1),NPIX1(2),STEP1,WSTR1,NPTOT1,
     +            MADRID(PNTR2),NPIX2(1),NPIX2(2),STEP2,WSTR2,NPTOT2,
     +            NFUNC,NFPAR,NFPACT,DPARM,INOP,RMIN,RMAX)
C
C ... write non standard descriptors
C
CCCCCCCCCCCCCCCCC      CALL SXDCOP(IMAGE1,IMAGE2,3,' ',ISTAT)
      CUTS(1) = RMIN
      CUTS(2) = RMAX
      CUTS(3) = RMIN
      CUTS(4) = RMAX
      CALL STDWRR(IMNO,'LHCUTS',CUTS,1,4,KUN,ISTAT)
      IF (NPIX1(2).GT.1) THEN
          CALL STDWRD(IMNO,'WSTART',WSTR2,1,NPIX2(2),KUN,ISTAT)
          CALL STDWRI(IMNO,'NPTOT',NPTOT2,1,NPIX2(2),KUN,ISTAT)
          CALL STDWRI(IMNO,'NORDER',NORDER,1,NPIX2(2),KUN,ISTAT)
      END IF
C
C ... end
C
   50 IF (ISTAT.NE.0) THEN
          CALL STTPUT(' Error in ECHCONV',STATUS)
      END IF
      CALL STSEPI
      END

      SUBROUTINE OUTIMA(NORDER,STEP1,WSTR1,NPTOT1,STARTR,STEPR,
     +                  START2,STEP2,WSTR2,NPTOT2,NPMAX)
C
C FIND SAMPLING DOMAIN IN THE OUTPUT IMAGE BASED ON THE
C REFERENCE START AND STEP VALUES
C
C INPUT ARGUMENTS
C NORDER  NO OF ORDERS
C STEP1   STEP IN INPUT IMAGE
C WSTR1   START IN INPUT IMAGE
C NPTOT1  NO OF PIXELS/ORDER
C STARTR  START USED AS REFERENCE
C STEPR   STEP USED AS REFERENCE
C OUTPUT ARGUMENTS
C START2  START IN OUTPUT IMAGE IF NORDER = 1
C STEP2   STEP IN OUTPUT IMAGE
C WSTR2   START IN OUTPUT IMAGE IF NORDER > 1
C NPTOT2  NO OF PIXELS/ORDER IF NORDER > 1
C NPMAX   MAX(NPTOT)
C
      IMPLICIT NONE
      DOUBLE PRECISION WSTR1(1),WSTR2(1)
      DOUBLE PRECISION STARTR, STEPR, W, STEP1, STEP2, START2
      INTEGER NPTOT1(1),NPTOT2(1)
      INTEGER I, NPMAX, NORDER, NPOS
C
C ... iterate for each order
C
      NPMAX  = 0
      DO 10 I = 1,NORDER
C
C ...     start, end for each order
C
          NPOS      = NINT((WSTR1(I)-STARTR)/STEPR)
          WSTR2(I)  = STARTR + (NPOS+1)*STEPR
          NPOS      = INT((WSTR1(I)+ (NPTOT1(I)-1)*STEP1-STARTR)/STEPR)
          W         = STARTR + (NPOS-1)*STEPR
          NPTOT2(I) = (W-WSTR2(I))/STEPR + 1
          NPMAX     = MAX(NPMAX,NPTOT2(I))
C   TYPE *,WSTR2(I),NPTOT2(I)
   10 CONTINUE
      STEP2  = STEPR
      IF (NORDER.EQ.1) START2 = WSTR2(1)
      RETURN
      END

      SUBROUTINE APPREB(X,NPIX11,NPIX12,STEP1,WSTR1,NPTOT1,
     +                  Y,NPIX21,NPIX22,STEP2,WSTR2,NPTOT2,
     +                  NFUNC,NFPAR,NFPACT,DPARM,INOP,RMIN,RMAX)
C
C REBIN EACH ORDER
C
      IMPLICIT NONE
      INTEGER NPIX11,NPIX12,NPIX21,NPIX22,NPTOT1(1),NPTOT2(1)
      REAL    X(NPIX11,NPIX12)
      REAL    Y(NPIX21,NPIX22)
      DOUBLE PRECISION WSTR1(1),WSTR2(1), STEP1, STEP2, WX
      REAL    RMIN,RMAX,RMIN1,RMAX1
      DOUBLE PRECISION DPARM(1), X0
      INTEGER   MADRID, NI, NI1, IS
      INTEGER*8 IP1, IP2, IP3
      INTEGER   NO, NO1
      INTEGER*8 JP1, JP2
      INTEGER   I, N1, N2, NINT, NP1, NFUNC,NFPAR
      INTEGER   NFPACT, INOP, J1, J2, J
C
      COMMON/VMR/MADRID(1)
C
C ALLOCATE WORK SPACE
C
      NI     = NPIX11
      NI1    = 8*NI
      CALL TDMGET(NI1,IP1,IS)
      CALL TDMGET(NI1,IP2,IS)
      CALL TDMGET(NI1,IP3,IS)
      NO     = NPIX21
      NO1    = 8*NO
      CALL TDMGET(NO1,JP1,IS)
      CALL TDMGET(NO1,JP2,IS)
      RMIN   = 0.
      RMAX   = 0.
      DO 20 I = 1,NPIX12
C
C   INPUT ARRAYS
C
          X0     = WSTR1(I)
          WX     = STEP1
          N1     = NPTOT1(I)
          CALL IMVAL3(N1,X0,WX,X(1,I),MADRID(IP1),MADRID(IP2),
     +                 MADRID(IP3),NP1)
C
C   OUTPUT ARRAYS
C
          N2     = NPTOT2(I)
          X0     = WSTR2(I)
          WX     = STEP2
          CALL IMVAL2(N2,X0,WX,MADRID(JP1),MADRID(JP2))
C
C         CALL INTERPOLATING (AND INTEGRATING) ROUTINES:
C
C  !!!! Watch out for VECI,VECD arrays in INTERP, 
C	dimension => NINT !!!!!!
C
          NINT   = 8
          CALL REBMET(NP1,MADRID(IP1),MADRID(IP3),MADRID(IP2),N2,
     +                 MADRID(JP1),MADRID(JP2),NFUNC,NFPAR,NFPACT,DPARM,
     +                 INOP,NINT,Y(1,I),RMIN1,RMAX1)
          RMIN   = AMIN1(RMIN1,RMIN)
          RMAX   = AMAX1(RMAX1,RMAX)
          J1     = N2 + 1
          J2     = NPIX21
          DO 10 J = J1,J2
              Y(J,I) = 0.
   10     CONTINUE
   20 CONTINUE
C
C  FREE MEMORY
C
      CALL TDMFRE(NI1,IP1,IS)
      CALL TDMFRE(NI1,IP2,IS)
      CALL TDMFRE(NI1,IP3,IS)
      CALL TDMFRE(NO1,JP1,IS)
      CALL TDMFRE(NO1,JP2,IS)
      RETURN
      END

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

      SUBROUTINE IMVAL3(NP,XS,WS,YY,X,W,Y,NP1)
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C  Subroutine IMVAL3
C
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
      DOUBLE PRECISION X(NP),W(NP),Y(NP),WSS,XSS
      DOUBLE PRECISION WS, XS
      REAL YY(NP)
      INTEGER NP1, I
C
      NP1    = 0
      WSS    = WS
      XSS    = XS
      DO 10 I = 1,NP
           IF (YY(I).GE.0.0) THEN         ! Modif. 901217 P.B (old .GT.)
              NP1    = NP1 + 1
              X(NP1) = XSS + WSS* (I-1)
              W(NP1) = WSS
              IF (YY(I).GT.1.E30) THEN
                 Y(NP1) = 0.D0
              ELSE 
                 Y(NP1) = YY(I)
              ENDIF
           END IF
   10 CONTINUE
      RETURN
      END
