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 ECHMRG
C     +
C     .COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C     all rights reserved
C     
C     .VERSION: 1.0  ESO-FORTRAN Conversion, AA  22:26 - 3 DEC 1987
C     
C     .LANGUAGE: F77+ESOext
C     
C     .AUTHOR: D.PONZ
C     
C     
C     .IDENTIFICATION
C     
C     program  ECHMERG 
C     
C     .MODIFICATIONS
C     
C     910128   P. Ballester   Define and set variables for mmake -i
C     981104   O. Stahl       Modified for FEROS
C     
C     .KEYWORDS
C     
C     ECHELLE, CASPEC, ORDER MERGING
C     
C     .PURPOSE
C     
C     OBTAIN A 1D FRM CALIBRATED IN WAVELENGTHS FROM A 2D FRM
C     SAMPLED IN THE SPACE WAVELENGTH-ORDER.
C     
C     .ALGORITHM
C     
C     OVERLAPPING REGIONS OF CONSECUTIVE ORDERS ARE PROCESSED IN TWO WAYS:
C     - CONCATENATION OF CONSECUTIVE ORDERS, USING THE MIDDLE POINT AS
C     CONCATENATING POSITION (METHOD='CONCATENATE', DEFAULT)
C     - AVERAGE OF OVERLAPPING REGION (METHOD='AVERAGE')
C     - NO MERGING. INDIVIDUAL ORDERS ARE WRITTEN IN DIFFERENT FILES
C     (METHOD = 'NOMERGE')
C     - WEIGHTED AVERAGE. WEIGHTS PROPORTIONAL TO SINC**2 (METHOD='SINC')
C     THE METHOD 'NOMERGE' WILL PRODUCE AS MANY FILES AS DEFINED ORDERS
C     WITH NAMES xxxyyyy, xxxx IS THE OUTPUT FILE NAME, yyyy IS THE 
C     ORDER NUMBER
C     
C     .INPUT/OUTPUT
C     
C     MERGE/ECHELLE INPUT OUTPUT  [CONC]
C     MERGE/ECHELLE INPUT OUTPUT  ORD1,ORD2 NOCONCAT
C     MERGE/ECHELLE INPUT OUTPUT  DEL       AVERAGE
C     OR
C     MERGE/ECHELLE INPUT OUTPUT  TABLE SINC
C 
C.VERSION
C 30714		last modif
C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C     
C     
      IMPLICIT NONE
C 
      INTEGER MADRID
      INTEGER I, II1, II2, IORD1
      INTEGER IORD2
      INTEGER*8 IPNTRA, IPNTRB
      INTEGER ISTAT
      INTEGER NAXISA, NAXISB
C 
      DOUBLE PRECISION    DEL, WEND, WINIT, WSTEP
C     
      INTEGER NPIXA(3),NPIXB(3),NPTOT(100),ICOL(3),KUN,KNUL
      INTEGER INIMA, OUTIMA, TID
      DOUBLE PRECISION  WSTART(100),STEPA(3)
      DOUBLE PRECISION  STARTA(3),STEPB(3),STARTB(3)
      REAL CUT(4),XMIN,XMAX,WRANGE(2),VAL(3)
      REAL WAVE(100),ORDER(100),WAVE1(100),WAVE2(100)
      REAL WAVE3(100),WAVE4(100)
      CHARACTER*60 INFRM,OUTFRM
      CHARACTER OUTFIL*60
      CHARACTER METHOD*1,WS*5
      CHARACTER IDENTA*72,IDENTB*72,CUNITA*64,CUNITB*64
      CHARACTER*16 LABEL(3)
      LOGICAL NULL(3)
C     
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C     
C     DATA CUNITB/'WAVELENGTH      FLUX      '/
      DATA CUNITB/'FLUX            WAVELENGTH'/
      DATA LABEL(1)/'WAVE'/,LABEL(2)/'ORDER'/,LABEL(3)/'LNWAVE'/
C     
C     ... INITIALIZE SYSTEM
C     
      CALL STSPRO('ECHMRG')
      INFRM = ' '
      OUTFRM = ' '
      CALL STKRDC('P1',1,1,60,I,INFRM,KUN,KNUL,ISTAT)
      CALL STKRDC('P2',1,1,60,I,OUTFRM,KUN,KNUL,ISTAT)
      CALL CLNFRA(INFRM,INFRM,0)
      CALL CLNFRA(OUTFRM,OUTFRM,0)
      CALL STKRDC('P4',1,1,1,I,METHOD,KUN,KNUL,ISTAT)
      CALL FORUPC(METHOD,METHOD)
      CALL STKRDR('INPUTR',1,2,I,WRANGE,KUN,KNUL,ISTAT)
C
      IF (METHOD.EQ.'A') THEN
         WINIT  = 0.
         WEND   = WINIT
         DEL    = WRANGE(1)
      ELSE
         WINIT  = WRANGE(1)
         WEND   = WRANGE(2)
      END IF
 
C
      IF (METHOD.EQ.'S') THEN
         WINIT  = 0.
         WEND   = WINIT
         DEL    = WRANGE(1)
      END IF
C     
C     ... MAP INPUT FRM
C     
      CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .     3,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA,
     .     IPNTRA,INIMA,ISTAT)
      CALL STDRDD(INIMA,'WSTART',1,NPIXA(2),I,
     .     WSTART,KUN,KNUL,ISTAT)
      CALL STDRDI(INIMA,'NPTOT',1,NPIXA(2),I,NPTOT,
     .     KUN,KNUL,ISTAT)
C     ... CHECK IF NPTOT(NPIXA(2)) IS 0 
      IF (NPTOT(NPIXA(2)).EQ.0) THEN  
         NPIXA(2) = NPIXA(2) - 1
      END IF
      WSTEP  = STEPA(1)
      IF (METHOD.EQ.'N') THEN
C     
C     ... NO MERGE - PRODUCE ONE FILE PER ORDER
C     
         IF (WINIT.EQ.WEND .AND. WINIT.LT.0.5) THEN
            IORD1  = 1
            IORD2  = NPIXA(2)
         ELSE
            IORD1  = MAX(NINT(WINIT),1)
            IORD2  = MIN(NINT(WEND),NPIXA(2))
         END IF
         OUTFRM = OUTFRM(1:59)//' '
         II1    = INDEX(OUTFRM,' ') - 1
         II1    = MIN(II1,20)
         OUTFIL = OUTFRM(1:II1)
         DO 10 I = IORD1,IORD2
            II2    = 10000 + I
            WRITE (WS,9000) II2
            OUTFIL(II1+1:II1+4) = WS(2:5)
            NAXISB    = 1
            NPIXB(1)  = NPTOT(I)
            NPIXB(2)  = 1
            STARTB(1) = WSTART(I)
            STARTB(2) = 1.
            STEPB(1)  = WSTEP
            STEPB(2)  = 0.
            WRITE (IDENTB,9010) I
            IDENTB(11:72) = IDENTA(1:62)
            CALL STIPUT(OUTFIL,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .           NAXISB,NPIXB,STARTB,STEPB,
     .           IDENTB,CUNITB,IPNTRB,OUTIMA,ISTAT)
            CALL COPY(MADRID(IPNTRA),NPIXA(1),NPIXA(2),MADRID(IPNTRB),
     .           NPIXB(1),I,XMIN,XMAX)
            CUT(1) = XMIN
            CUT(2) = XMAX
            CUT(3) = XMIN
            CUT(4) = XMAX
            CALL DSCUPT(INIMA,OUTIMA,' ',ISTAT)
            CALL STDWRR(OUTIMA,'LHCUTS',CUT,1,4,KUN,ISTAT)
C     
            CALL STFCLO(OUTIMA,ISTAT)
            CALL STTPUT('File '//OUTFIL//' created ...',ISTAT)
 10      CONTINUE

      ELSE
C     
C     ... MAP OUTPUT FRAME
C     
         IF (WINIT.EQ.WEND) THEN
            WINIT  = WSTART(1)
            WEND   = WSTART(NPIXA(2)) + (NPTOT(NPIXA(2))-1)*WSTEP
         END IF
         NAXISB    = 1
         NPIXB(1)  = (WEND-WINIT)/WSTEP + 1
         NPIXB(2)  = 1
         STARTB(1) = WINIT
         STARTB(2) = 1.
         STEPB(1)  = WSTEP
         STEPB(2)  = 0.
         CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .        NAXISB,NPIXB,STARTB,STEPB,
     .        IDENTA,CUNITB,IPNTRB,OUTIMA,ISTAT)
C     
C     ... METHODS OF OVERLAPPING
C     
         IF (METHOD.EQ.'C') CALL ECHMR1(MADRID(IPNTRA),
     +        NPIXA(1),NPIXA(2),STARTA,
     +        STEPA,WSTART,NPTOT,MADRID(IPNTRB),
     +        NPIXB,STARTB,XMIN,XMAX)
         IF (METHOD.EQ.'A') CALL ECHMR2(MADRID(IPNTRA),
     +        NPIXA(1),NPIXA(2),STARTA,
     +        STEPA,WSTART,NPTOT,MADRID(IPNTRB),
     +        NPIXB,STARTB,XMIN,XMAX,DEL)


         IF (METHOD.EQ.'S') THEN
C     
C     ... READ RIPPLE TABLE ------------------------------------------
C     
            CALL TBTOPN('BLAZE',F_I_MODE,TID,ISTAT)
            CALL TBLSER(TID,LABEL(1),ICOL(1),ISTAT)
            CALL TBLSER(TID,LABEL(2),ICOL(2),ISTAT)
            CALL TBLSER(TID,LABEL(3),ICOL(3),ISTAT)
C
            DO 20 I = 1,NPIXA(2)-1
               CALL TBRRDR(TID,I,2,ICOL,VAL,NULL,ISTAT)

               IF(WSTART(1).GT.1000) THEN
C     looks like linear wavelength scale
                  WAVE(I)  = VAL(1)
               ELSE
C     looks like logarithmic wavelength scale
                  WAVE(I)  = VAL(3)
               ENDIF

               ORDER(I)   = VAL(2)
 20         CONTINUE
            
            DO 30 I=1,NPIXA(2)
               WAVE1(I) = WSTART(I)+DEL
               WAVE2(I) = WSTART(I)+(NPTOT(I)-1)*STEPA(1) - DEL
 30         CONTINUE

            DO 60 I=1,NPIXA(2)
               WAVE3(I) = WAVE1(I)
               WAVE4(I) = WAVE2(I)
 60         CONTINUE

            DO 40 I = 1,NPIXA(2)-2
               WAVE3(I+1)  = (WAVE(I)+WAVE(I+1))/2 -
     +              (WAVE(I+1)-WAVE(I))/2*1.5
               WAVE4(I+1)  = (WAVE(I)+WAVE(I+1))/2 +
     +              (WAVE(I+1)-WAVE(I))/2*1.5
               
 40         CONTINUE
            
            DO 50 I=1,NPIXA(2)
               WAVE1(I) = MAX(WAVE1(I),WAVE3(I))
               WAVE2(I) = MIN(WAVE2(I),WAVE4(I))
 50         CONTINUE
            

            CALL ECHMR3(MADRID(IPNTRA),
     +           NPIXA(1),NPIXA(2),STARTA,
     +           STEPA,WSTART,NPTOT,MADRID(IPNTRB),
     +           NPIXB,STARTB,XMIN,XMAX,DEL,WAVE1,WAVE2)
            
         END IF

C     ... DONE WITH RIPPLE TABLE --------------------------------------

C     
C     ... WRITE PROCESS DESCRIPTORS
C     
         CUT(1) = XMIN
         CUT(2) = XMAX
         CUT(3) = XMIN
         CUT(4) = XMAX
         CALL DSCUPT(INIMA,OUTIMA,' ',ISTAT)
         CALL STDWRR(OUTIMA,'LHCUTS',CUT,1,4,KUN,ISTAT)

      END IF
      CALL STSEPI
 9000 FORMAT (I5)
 9010 FORMAT ('ORDER:',I3)
      END

      SUBROUTINE ECHMR1(X,NPIXA1,NPIXA2,STARTA,STEPA,WI,NP,Y,NY,YSTR,
     +     XMIN,XMAX)
C     
C     MERGE THE ORDERS, USING SIMPLE CONCATENATION IN THE MIDDLE POINT
      IMPLICIT NONE
C     
      INTEGER  NY,J2,I,J1,JOFF,J,JJ
      REAL     XMIN,XMAX,WEN1,WST2,WSTART
      INTEGER  NPIXA1,NPIXA2
      INTEGER  NP(NPIXA2)
      REAL     X(NPIXA1,NPIXA2)
      DOUBLE PRECISION STARTA(2),STEPA(2),WI(NPIXA2),YSTR
      REAL     Y(NY)
      DOUBLE PRECISION WEND, YSTEP, YEND, WSTR
C     
      XMIN   = 0.
      XMAX   = 0.
      DO 10 I = 1,NY
         Y(I)   = 0.
 10   CONTINUE
      YSTEP  = STEPA(1)
      YEND   = YSTR + (NY-1)*YSTEP
C     
C     ... ITERATION ON ORDERS
C     
      WEND   = 0.
      DO 30 I = 1,NPIXA2
         WSTR   = DMAX1(WI(I),WEND+YSTEP)
         IF (I.EQ.NPIXA2) THEN
            WEND   = WI(I) + (NP(I)-1)*YSTEP
         ELSE
            WEN1   = WI(I) + (NP(I)-1)*YSTEP
            WST2   = WI(I+1)
            IF (WST2.LT.WEN1) THEN
               WEND   = 0.5* (WEN1+WST2)
            ELSE
               WEND   = WEN1
            END IF
         END IF

         IF (WSTR.GE.YEND) RETURN
         IF (WEND.GT.YSTR) THEN
C     
C     ... ITERATION ON WAVELENGTHS
C     
            WSTART = DMAX1(YSTR,WSTR)
            WEND   = DMIN1(WEND,YEND)
            J1     = NINT((WSTART-WI(I))/YSTEP) + 1
            J2     = NINT((WEND-WI(I))/YSTEP) + 1
            JOFF   = NINT((WI(I)-YSTR)/YSTEP)
            DO 20 J = J1,J2
               JJ     = J + JOFF
               IF (JJ.GT.0) THEN
                  Y(JJ)  = X(J,I)
                  XMAX   = MAX(XMAX,Y(JJ))
                  XMIN   = MIN(XMIN,Y(JJ))
               END IF

 20         CONTINUE
         END IF

 30   CONTINUE
      RETURN

      END

      SUBROUTINE ECHMR2(X,NPIXA1,NPIXA2,STARTA,STEPA,WI,NP,Y,NY,YSTR,
     +     XMIN,XMAX,DEL)
C     
C     MERGE THE ORDERS, AVERAGE ON THE OVERLAPPING REGION
C     
      IMPLICIT NONE
C     
      INTEGER NPIXA1, NPIXA2, NY
      INTEGER NP(NPIXA2)
      REAL X(NPIXA1,NPIXA2)
      DOUBLE PRECISION STARTA(2),STEPA(2),YSTR,DEL,WI(NPIXA2)
      REAL Y(NY), XMIN, XMAX
C     
      INTEGER I, IORD1, IORD2, IPIX, IPIX1, IPIX2
      DOUBLE PRECISION YSTEP, W0, W1, WL, P1, P2
C     
      XMIN   = 0.
      XMAX   = 0.
      DO 10 I = 1,NY
         Y(I)   = 0.
 10   CONTINUE
      YSTEP  = STEPA(1)
      IORD1  = 1
      IORD2  = 2
      W0     = WI(IORD2) + DEL
      W1     = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL
C     
C     ... ITERATION ON ORDERS
C     
      DO 20 IPIX = 1,NY
         WL     = YSTR + (IPIX-1)*YSTEP
         IF (WL.LT.W0) THEN
            IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
            Y(IPIX) = X(IPIX1,IORD1)
            XMAX   = MAX(Y(IPIX),XMAX)
            XMIN   = MIN(Y(IPIX),XMIN)
         ELSE IF (WL.LT.W1) THEN
            IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
            IPIX2  = NINT((WL-WI(IORD2))/YSTEP) + 1
            P2     = (WL-W0)/ (W1-W0)
            P1     = 1.D0 - P2
            IF (X(IPIX1,IORD1).LE.0.0) THEN
               P2     = 1.D0
               P1     = 0.D0
            END IF
            IF (X(IPIX2,IORD2).LE.0.0) THEN
               P2     = 0.D0
               P1     = 1.D0
            END IF
            Y(IPIX) = X(IPIX1,IORD1)*P1 + X(IPIX2,IORD2)*P2
            XMAX   = MAX(Y(IPIX),XMAX)
            XMIN   = MIN(Y(IPIX),XMIN)
         ELSE
            IORD1  = IORD1 + 1
            IF (IORD1.GT.NPIXA2) RETURN
            IORD2  = IORD2 + 1
            IF (IORD2.GT.NPIXA2) THEN
               W0     = 1.E20
            ELSE
               W0     = WI(IORD2) + DEL
            END IF
            W1     = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL
            IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
            Y(IPIX) = X(IPIX1,IORD1)
            XMAX   = MAX(Y(IPIX),XMAX)
            XMIN   = MIN(Y(IPIX),XMIN)
         END IF
 20   CONTINUE
      RETURN
      END

      SUBROUTINE COPY(A,NA1,NA2,B,NB,I,XMIN,XMAX)
C     
C     COPY THE ORDER NUMBER I
C     
      IMPLICIT NONE
      INTEGER NA1, NA2, NB, I, J
      REAL A(NA1,NA2),B(NB)
      REAL XMIN, XMAX
C     
      XMIN = 0.
      XMAX = 0.
C     
      DO 10 J = 1,NB
         B(J)   = A(J,I)
         XMIN   = MIN(B(J),XMIN)
         XMAX   = MAX(B(J),XMAX)
 10   CONTINUE
      RETURN
      END

      SUBROUTINE ECHMR3(X,NPIXA1,NPIXA2,STARTA,STEPA,WI,NP,Y,NY,YSTR,
     +     XMIN,XMAX,DEL,WAVE1,WAVE2)
C     
C     MERGE THE ORDERS, AVERAGE ON THE OVERLAPPING REGION
C     
      IMPLICIT NONE
C     
      INTEGER NPIXA1, NPIXA2, NY
      INTEGER NP(NPIXA2)
      REAL X(NPIXA1,NPIXA2)
      DOUBLE PRECISION STARTA(2),STEPA(2),YSTR,DEL,WI(NPIXA2)
      REAL Y(NY), XMIN, XMAX, WAVE1(NPIXA2), WAVE2(NPIXA2)
C     
      INTEGER I, IORD1, IORD2, IPIX, IPIX1, IPIX2
      DOUBLE PRECISION YSTEP, W0, W1, WL, P1, P2
C     
      XMIN   = 0.
      XMAX   = 0.
      DO 10 I = 1,NY
         Y(I)   = 0.
 10   CONTINUE
      YSTEP  = STEPA(1)
      IORD1  = 1
      IORD2  = 2

C      W0     = WI(IORD2) + DEL
C      W1     = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL

      W0 = WAVE1(IORD2)
      W1 = WAVE2(IORD1)
C     
C     ... ITERATION ON ORDERS
C     
      DO 20 IPIX = 1,NY
         WL     = YSTR + (IPIX-1)*YSTEP
         IF (WL.LT.W0) THEN
            IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
            Y(IPIX) = X(IPIX1,IORD1)
            XMAX   = MAX(Y(IPIX),XMAX)
            XMIN   = MIN(Y(IPIX),XMIN)
         ELSE IF (WL.LT.W1) THEN
            IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
            IPIX2  = NINT((WL-WI(IORD2))/YSTEP) + 1
            P2     = (WL-W0)/ (W1-W0)
            P1     = 1.D0 - P2
            IF (X(IPIX1,IORD1).LE.0.0) THEN
               P2     = 1.D0
               P1     = 0.D0
            END IF
            IF (X(IPIX2,IORD2).LE.0.0) THEN
               P2     = 0.D0
               P1     = 1.D0
            END IF
            Y(IPIX) = X(IPIX1,IORD1)*P1 + X(IPIX2,IORD2)*P2
            XMAX   = MAX(Y(IPIX),XMAX)
            XMIN   = MIN(Y(IPIX),XMIN)
         ELSE
            IORD1  = IORD1 + 1
            IF (IORD1.GT.NPIXA2) RETURN
            IORD2  = IORD2 + 1
            IF (IORD2.GT.NPIXA2) THEN
               W0     = 1.E20
            ELSE
C               W0     = WI(IORD2) + DEL
               W0     = WAVE1(IORD2)
            END IF
C            W1     = WI(IORD1) + (NP(IORD1)-1)*YSTEP - DEL
            W1     = WAVE2(IORD1)
            IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
            Y(IPIX) = X(IPIX1,IORD1)
            XMAX   = MAX(Y(IPIX),XMAX)
            XMIN   = MIN(Y(IPIX),XMIN)
         END IF
 20   CONTINUE
      RETURN
      END

