C @(#)necmerge.for	17.1.1.1 (ESO-DMD) 01/25/02 17:51:30
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
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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  991005   S.Wolf         ECHMR2O added: does an optimal merging.
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  - WEIGHTED AVERAGE OF OVERLAPPING REGION (METHOD='OPTIMAL')
C    spmerged = (w1*sp1 + w2*sp2)/(w1+w2)
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  DELTA     AVERAGE
C  MERGE/ECHELLE INPUT OUTPUT  DELTA     OPTIMAL WEIGHTS
C  or
C  MERGE/ECHELLE INPUT OUTPUT  TABLE SINC
C
C 
C 010702		last modif
C 
C------------------------------------------------------------------------
C 
      PROGRAM ECHMRG
C 
      IMPLICIT  NONE
C
      INTEGER   MADRID
      INTEGER   I, II1, II2, IORD1
      INTEGER   IORD2
      INTEGER   ISTAT
      INTEGER   NAXISA, NAXISB, NAXISW
      INTEGER NPIXA(3),NPIXB(3),NPIXW(3),NPTOT(100),ICOL(3),KUN,KNUL
      INTEGER INIMA, OUTIMA, VARIMA, WGTIMA, TID
C 
      INTEGER*8 IPNTRA, IPNTRB, IPNTRV, IPNTRW              !pointers ...
C   
      DOUBLE PRECISION  DEL, WEND, WINIT, WSTEP
      DOUBLE PRECISION  WSTART(100)
      DOUBLE PRECISION  STARTA(3),STEPA(3)
      DOUBLE PRECISION  STARTB(3),STEPB(3)
      DOUBLE PRECISION  STARTW(3),STEPW(3)
C
      REAL CUT(4),XMIN,XMAX,WRANGE(2),VAL(3),K(100),A(100)
      REAL RO(100)
C
      CHARACTER*60 INFRM, OUTFRM, VARFRM, WGTFRM
      CHARACTER OUTFIL*12,TABLE*12
      CHARACTER METHOD*1,WS*5
      CHARACTER IDENTA*72,IDENTB*72,CUNITA*64,CUNITB*64
      CHARACTER IDENTW*72,CUNITW*64
      CHARACTER*16 LABEL(3)
C
      LOGICAL NULL(3)
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      DATA CUNITB/'FLUX            WAVELENGTH'/
      DATA LABEL(1)/'ORDER'/,LABEL(2)/'KFIT'/,LABEL(3)/'AFIT'/
C
C ... INITIALIZE SYSTEM
C
      CALL STSPRO('ECHMRG')
      INFRM = ' '
      OUTFRM = ' '
      VARFRM = ' '
      WGTFRM = ' '
      CALL STKRDC('P1',1,1,60,I,INFRM,KUN,KNUL,ISTAT)
      CALL STKRDC('P2',1,1,60,I,OUTFRM,KUN,KNUL,ISTAT)
      CALL STKRDC('P6',1,1,60,I,VARFRM,KUN,KNUL,ISTAT)
      CALL STKRDC('P5',1,1,60,I,WGTFRM,KUN,KNUL,ISTAT)
      CALL CLNFRA(INFRM,INFRM,0)
      CALL CLNFRA(OUTFRM,OUTFRM,0)
      CALL CLNFRA(VARFRM,VARFRM,0)
      CALL CLNFRA(WGTFRM,WGTFRM,0)
      CALL STKRDC('P4',1,1,1,I,METHOD,KUN,KNUL,ISTAT)
      CALL FORUPC(METHOD,METHOD)
      IF (METHOD.EQ.'S') 
     .    CALL STKRDC('P3',1,1,8,I,TABLE,KUN,KNUL,ISTAT)
      CALL STKRDR('INPUTR',1,2,I,WRANGE,KUN,KNUL,ISTAT)
      IF (METHOD.EQ.'A' .OR. METHOD.EQ.'O') THEN
          WINIT  = 0.
          WEND   = WINIT
          DEL    = WRANGE(1)
      ELSE
          WINIT  = WRANGE(1)
          WEND   = WRANGE(2)
      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)
C
C ... MAP WEIGHT FRM
C
      IF (METHOD.EQ.'O') THEN
         CALL STIGET(WGTFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .        3,NAXISW,NPIXW,STARTW,STEPW,IDENTW,CUNITW,
     .        IPNTRW,WGTIMA,ISTAT)
         IF (NAXISW.NE.NAXISA.OR.
     .        STARTW(1).NE.STARTA(1).OR.STARTW(2).NE.STARTA(2).OR.
     .        STEPW(1).NE.STEPA(1).OR.STEPW(1).NE.STEPA(1)) THEN
            GOTO 9998
         END IF
      END IF
      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,4)
          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)
              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.'O') THEN
             CALL STIPUT(VARFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .            NAXISB,NPIXB,STARTB,STEPB,
     .            IDENTA,CUNITB,IPNTRV,VARIMA,ISTAT)
             CALL ECHMR2O(MADRID(IPNTRA),MADRID(IPNTRW),
     +            NPIXA(1),NPIXA(2),STARTA,
     +            STEPA,WSTART,NPTOT,MADRID(IPNTRB),MADRID(IPNTRV),
     +            NPIXB,STARTB,XMIN,XMAX,DEL)
             CALL DSCUPT(INIMA,VARIMA,' ',ISTAT)
             
          ENDIF
          IF (METHOD.EQ.'S') THEN
C
C ... READ RIPPLE TABLE
C
              CALL TBTOPN(TABLE,F_IO_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)
              DO 20 I = 1,NPIXA(2)
                  CALL TBRRDR(TID,I,3,ICOL,VAL,NULL,ISTAT)
                  RO(I)  = VAL(1)
                  K(I)   = VAL(2)
                  A(I)   = VAL(3)
   20         CONTINUE
              CALL ECHMR3(MADRID(IPNTRA),
     +                      NPIXA(1),NPIXA(2),STARTA,STEPA,WSTART,
     +                      NPTOT,MADRID(IPNTRB),NPIXB,STARTB,XMIN,XMAX,
     +                      RO,K,A)
          END IF
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
      GOTO 9999
 9998 CALL STTPUT
     +     ('Error: Bad weight map! '/
     +     /'(START, STEP and NAXIS same as in INPUT frame !?)',ISTAT)
 9999 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
      INTEGER  NPIXA1,NPIXA2
      INTEGER  NP(NPIXA2)
C 
      DOUBLE PRECISION STARTA(2),STEPA(2),WI(NPIXA2),YSTR
C 
      REAL     X(NPIXA1,NPIXA2)
      REAL     XMIN,XMAX,RVAL,WEN1,WST2,WSTART
      REAL     Y(NY)
C 
      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
                      RVAL = X(J,I)
                      IF (RVAL .GT. XMAX) XMAX = RVAL
                      IF (RVAL .LT. XMIN) XMIN = RVAL
                      Y(JJ) = RVAL
                  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)
      INTEGER IORD1, IORD2, IPIX, IPIX1, IPIX2
C 
      REAL X(NPIXA1,NPIXA2)
      REAL Y(NY),XMIN,XMAX,RVAL
C 
      DOUBLE PRECISION STARTA(2),STEPA(2),YSTR,DEL,WI(NPIXA2)
      DOUBLE PRECISION YSTEP, W0, W1, WL, P1, P2
C
C
      XMIN   = 0.
      XMAX   = 0.
      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
          Y(IPIX) = 0.0
          WL     = YSTR + (IPIX-1)*YSTEP
          IF (WL.LT.W0) THEN
              IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
              RVAL = X(IPIX1,IORD1)
              IF (RVAL .GT. XMAX) XMAX = RVAL
              IF (RVAL .LT. XMIN) XMIN = RVAL
              Y(IPIX) = RVAL
          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
              RVAL = X(IPIX1,IORD1)*P1 + X(IPIX2,IORD2)*P2
              IF (RVAL .GT. XMAX) XMAX = RVAL
              IF (RVAL .LT. XMIN) XMIN = RVAL
              Y(IPIX) = RVAL
          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
              RVAL = X(IPIX1,IORD1)
              IF (RVAL .GT. XMAX) XMAX = RVAL
              IF (RVAL .LT. XMIN) XMIN = RVAL
              Y(IPIX) = RVAL
          END IF
   20 CONTINUE
      RETURN
      END

      SUBROUTINE ECHMR2O(X,WGT,NPIXA1,NPIXA2,STARTA,STEPA,WI,NP,Y,VAR,NY
     $     ,YSTR,XMIN,XMAX,DEL)
C
C MERGE THE ORDERS, AVERAGE ON THE OVERLAPPING REGION USING WEIGHTS
C
      IMPLICIT NONE
C
      INTEGER NPIXA1, NPIXA2, NY
      INTEGER NP(NPIXA2)
      REAL X(NPIXA1,NPIXA2), WGT(NPIXA1,NPIXA2)
      DOUBLE PRECISION STARTA(2),STEPA(2),YSTR,DEL,WI(NPIXA2)
      REAL Y(NY), VAR(NY),XMIN,XMAX,RVAL
C
      INTEGER IORD1, IORD2, IPIX, IPIX1, IPIX2, IBAD
      DOUBLE PRECISION YSTEP, WG1, WG2, W0, W1, WL
      CHARACTER*80 MES
C
      IBAD   = 0
      XMIN   = 0.
      XMAX   = 0.
      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
          Y(IPIX) = 0.0
          VAR(IPIX) = 0.0
          WL     = YSTR + (IPIX-1)*YSTEP
          IF (WL.LT.W0) THEN
              IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
              RVAL = X(IPIX1,IORD1)
              IF (RVAL .GT. XMAX) XMAX = RVAL
              IF (RVAL .LT. XMIN) XMIN = RVAL
              Y(IPIX) = RVAL
              VAR(IPIX) = WGT(IPIX1,IORD1)
              IF (VAR(IPIX).NE.0.0) VAR(IPIX) = 1.0/VAR(IPIX)
          ELSE IF (WL.LT.W1) THEN
              IPIX1  = NINT((WL-WI(IORD1))/YSTEP) + 1
              IPIX2  = NINT((WL-WI(IORD2))/YSTEP) + 1
              WG1 = WGT(IPIX1,IORD1)
              WG2 = WGT(IPIX2,IORD2)

              IF (WG1.LT.1e-10 .AND. WG2.LT.1e-10) THEN
                 Y(IPIX) = 0.0
                 VAR(IPIX) = 0.0
                 IBAD = IBAD + 1
              ELSE
                 VAR(IPIX) = 1.0/(WG1 + WG2)
                 Y(IPIX) = X(IPIX1,IORD1)*WG1 + X(IPIX2,IORD2)*WG2
                 RVAL = Y(IPIX) * VAR(IPIX)
                 IF (RVAL .GT. XMAX) XMAX = RVAL
                 IF (RVAL .LT. XMIN) XMIN = RVAL
                 Y(IPIX) = RVAL
              ENDIF
          ELSE
              IORD1  = IORD1 + 1
              IF (IORD1.GT.NPIXA2) GOTO 99
              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
              RVAL = X(IPIX1,IORD1)
              IF (RVAL .GT. XMAX) XMAX = RVAL
              IF (RVAL .LT. XMIN) XMIN = RVAL
              Y(IPIX) = RVAL
              VAR(IPIX) = WGT(IPIX1,IORD1)
              IF (VAR(IPIX).NE.0.0) VAR(IPIX) = 1.0/VAR(IPIX)
          END IF
 20   CONTINUE
 99   IF (IBAD.GT.0) THEN
         WRITE(MES,*) IBAD,' undefined pixels ... set to 0.0!'
         CALL STTPUT(MES,IBAD)
      END IF
      RETURN
         
      END

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

      SUBROUTINE ECHMR3(X,NPIXA1,NPIXA2,STARTA,STEPA,WI,NP,Y,NY,YSTR,
     +                    XMIN,XMAX,RORDER,K,A)
C
C      IMPLICIT NONE
C MERGE THE ORDERS, AVERAGE ON THE OVERLAPPING REGION WITH
C WEIGHTS PROPORTIONAL TO THE SINC**2
C
      IMPLICIT NONE
C 
      INTEGER NY,IORD1,IPIX,NF,IORD,IPIX1,I,IFL,IORD2
      INTEGER NPIXA1,NPIXA2
      INTEGER NP(NPIXA2)
C 
      REAL    XMIN,XMAX,YSTEP,PI,DW,WL
      REAL    X(NPIXA1,NPIXA2)
      REAL    Y(NY),RORDER(NPIXA2)
      REAL    K(NPIXA2),A(NPIXA2)
C 
      DOUBLE PRECISION DK,DA,PA,DM,DC,DX,WEIGHT(3),FLUX(3),SW,SF
      DOUBLE PRECISION YSTR,STARTA(2),STEPA(2),WI(NPIXA2)
C
      XMIN   = 0.
      XMAX   = 0.
      PI     = 0.     !  To be updated (ULTRIX installation)
      DW     = 0.     !  To be updated (ULTRIX installation)
      YSTEP  = STEPA(1)
      IORD1  = 1
      IORD2  = MIN(IORD1+2,NPIXA2)
C
C ... ITERATION ON OUTPUT WAVELENGTHS
C
      DO 30 IPIX = 1,NY
          Y(IPIX) = 0.
          FLUX(1) = 0.D0
          FLUX(2) = 0.D0
          FLUX(3) = 0.D0
          WEIGHT(1) = 0.D0
          WEIGHT(2) = 0.D0
          WEIGHT(3) = 0.D0
          NF     = 0
          WL     = YSTR + (IPIX-1)*YSTEP
C
C ... ITERATION ON ORDERS
C
          DO 10 IORD = IORD1,IORD2
              IPIX1  = (WL-WI(IORD))/YSTEP + 1
              IF (IPIX1.GT.5 .AND. IPIX1.LT. (NP(IORD)-5)) THEN
C
C ... PIXEL IN THE ORDER RANGE
C
                  NF     = NF + 1
                  FLUX(NF) = X(IPIX1,IORD)
                  DK     = K(IORD)
                  DA     = A(IORD)
                  PA     = PI*DA
                  DM     = RORDER(IORD)
                  DC     = DM/DK
                  DX     = (PA*DM*DC)* (DW-1.D0/DC)
                  IF (DABS(DX).LT.1.D-10) THEN
                      WEIGHT(NF) = 1.D0

                  ELSE
                      WEIGHT(NF) = 1.D0/ (DSIN(DX)/DX)**4
                  END IF

              END IF

   10     CONTINUE
C
C ... AVERAGE OVER FLUXES
C
          IF (NF.GE.1) THEN
              SW     = 0.D0
              SF     = 0.D0
              DO 20 IFL = 1,NF
                  SW     = SW + WEIGHT(IFL)
                  SF     = SF + WEIGHT(IFL)*FLUX(IFL)
   20         CONTINUE
              Y(IPIX) = SF/SW
              XMIN   = MIN(XMIN,Y(IPIX))
              XMAX   = MAX(XMAX,Y(IPIX))
          END IF
C
C ... CHECK ON ORDER LIMITS
C
          IPIX1  = (WL-WI(IORD1))/YSTEP + 1
          IF (IPIX1.GE.NP(IORD1)-5) THEN
              IORD1  = IORD1 + 1
              IORD2  = MIN(IORD1+2,NPIXA2)
          END IF

          IF (IORD1.GT.NPIXA2) RETURN
   30 CONTINUE
      DO 40 I = IPIX + 1,NY
          Y(I)   = 0.
   40 CONTINUE
      RETURN

      END