C @(#)necexor.for	19.1 (ESO-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 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.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  21:13 - 3 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C.IDENTIFICATION
C  program  ECHEXOR  version 1.2 840114
C
C.KEYWORDS
C
C  ECHELLE, CASPEC, DATA EXTRACTION
C
C.PURPOSE
C
C  EXTRACT ECHELLE ORDERS IN CASPEC IMAGES. ORDER POSITION ARE DEFINED
C  BY COEFFICIENTS OF A REGRESSION. SLIT WIDTH, ANGLE AND OFFSET
C  ARE DEFINED BY PARAMETERS
C            WIDTH,ANGLE,OFFSET
C
C.ALGORITHM
C
C   ORDERS ARE EXTRACTED BY PASSING A NUMERICAL SLIT OF GIVEN WIDTH 
C  AND ANGLE IN POSITIONS DEFINED BY THE REGRESSION COEFFS. THE SLIT
C  WILL SAMPLE AT INCREMENTS OF ONE STEP IN X.
C
C  THERE ARE THREE METHODS :
C   LINEAR   - PIXEL VALUES ARE ADDED, LINEAR INTERPOLATION
C   AVERAGE  - PIXEL VALUES ARE AVERAGED, LINEAR INTERPOLATION
C   WEIGHTED - PIXEL VALUES ARE WEIGHTD PROPORTIONALLY TO THE
C              PROFI PERPENDICULAR TO THE DISPERSION DIRECTION
C.INPUT/OUTPUT
C
C   EXTRACT/ORDER output = input width,angle,offset 
C				 table coeffs [ord1,ord2]
C  the connection file contains
C  IN_A   input FRM
C  IN_B   input table with coeffs
C  OUT_A  output FRM
C  INPUTR  width,angle
C  INPUTC  coeffs
C  INPUTI  ord1,ord2
C  ACTION               T - table option
C                       B - batch order by order
C 
C 
C 010123		last modif
C 
C------------------------------------------------------------------
C
      PROGRAM ECHEXO
C 
      IMPLICIT NONE
C 
      INTEGER NDIM
      PARAMETER (NDIM=2)
      INTEGER MADRID, NAC, NAR, IMNO
      INTEGER NWDTH
      INTEGER I,IAV,INDX,IORD,IORD1
      INTEGER IORD2,ISTAT,NA,NAXISA,NAXISB
      INTEGER NCOL,NPAR,NROW,NSORT
      INTEGER*8 PNTRA
      INTEGER*8 PNTRB
      INTEGER*8 PNTRC
      INTEGER   STATUS,TID,INO,TMPID
      INTEGER NPIXA(3),NPIXB(3),IPAR(10),KUN,KNUL
C
      DOUBLE PRECISION STEPA(3),STEPB(3), STARTA(3),STARTB(3)
      REAL  RPAR(10)
      REAL  CUTS(6)
      REAL XMIN,XMAX,XMIN1,XMAX1
      REAL PARAMS(12)
C
      DOUBLE PRECISION DPAR(100)
C
      CHARACTER*72   IDENA,IDENB
      CHARACTER*64   CUNITA,CUNITB,DESNAM
      CHARACTER*1    METH
      CHARACTER*60   FRMIN,FRMOUT,INA1,INA2
      CHARACTER*60   TABLE,NAME,TYPE
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      DATA CUNITB/
     +     'FLUX            PIXEL                           ORDER'/
C
C  initialize system
C
      CALL STSPRO('ECHEXO')
      CALL STKRDC('IN_A',1,1,60,I,INA1,KUN,KNUL,STATUS)
      CALL STKRDC('IN_B',1,1,60,I,TABLE,KUN,KNUL,STATUS)
      CALL STKRDC('IN_C',1,1,1,I,METH,KUN,KNUL,STATUS)
      CALL STKRDC('OUT_A',1,1,60,I,INA2,KUN,KNUL,STATUS)
      CALL STKRDC('INPUTC',1,1,60,I,NAME,KUN,KNUL,STATUS)
      CALL STKRDR('INPUTR',1,3,I,RPAR,KUN,KNUL,STATUS)
      CALL STKRDI('INPUTI',1,2,I,IPAR,KUN,KNUL,STATUS)
      IORD1  = IPAR(1)
      IORD2  = IPAR(2)
      FRMIN = INA1
      FRMOUT = INA2
      IF (METH(1:1) .EQ. 'l') METH = 'L'
      IF (METH(1:1) .EQ. 'a') METH = 'A'
      IF (METH(1:1) .EQ. 'm') METH = 'M'
C
C  input parameters in table
C
      CALL TBTOPN(TABLE,F_U_MODE,TID,STATUS)
      CALL TBIGET(TID,NCOL,NROW,NSORT,NAC,NAR,STATUS)
C
C  map input FRM
C
      CALL STIGET(FRMIN,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .            NDIM,NAXISA,NPIXA,STARTA,STEPA,IDENA,
     +            CUNITA,PNTRA,IMNO,STATUS)
      PARAMS(1) = RPAR(1)/ABS(STEPA(2))
      PARAMS(2) = RPAR(2)
      PARAMS(3) = RPAR(3)/ABS(STEPA(2))

      IF (METH(1:1) .EQ. 'M') THEN
         NWDTH = PARAMS(1) + 2
         CALL STFCRE('column',D_R4_FORMAT,F_X_MODE,1,NWDTH,TMPID,STATUS)
         CALL STFMAP(TMPID,F_X_MODE,1,NWDTH,IAV,PNTRC,STATUS)
      ENDIF
C
C  read coeffs
C
      NAME   = NAME(1:59)//' '
      INDX   = INDEX(NAME,' ') - 1
      DESNAM = NAME(1:INDX)//'C'
      CALL STDRDC(TID,DESNAM,1,17,4,IAV,TYPE,KUN,KNUL,ISTAT)
      TYPE   = 'MULT'
      IF (TYPE.EQ.'MULT') THEN
         DESNAM = NAME(1:INDX)//'I'
         CALL STDRDI(TID,DESNAM,4,4,IAV,IPAR,KUN,KNUL,ISTAT)
         NA     = (IPAR(3)+1)* (IPAR(4)+1)
         DESNAM = NAME(1:INDX)//'R'
         CALL STDRDR(TID,DESNAM,1,4,IAV,PARAMS(4),KUN,KNUL,ISTAT)
         DESNAM = NAME(1:INDX)//'D'
         CALL STDRDD(TID,DESNAM,1,NA,IAV,DPAR,KUN,KNUL,ISTAT)
         CALL TBTCLO(TID,STATUS)
C
C ... map output FRM
C
         IF (IORD1.EQ.IORD2 .AND. IORD1.EQ.0) THEN
            IORD1  = PARAMS(6)
            IORD2  = PARAMS(7)
         END IF
         STARTB(1) = STARTA(1)
         STARTB(2) = IORD1
         STEPB(1)  = STEPA(1)
         STEPB(2)  = 1.
         IDENB     = IDENA
         NAXISB    = 2
         NPIXB(1)  = NPIXA(1)
         NPIXB(2)  = IORD2 - IORD1 + 1
         CALL STIPUT(FRMOUT,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     +        NAXISB,NPIXB,STARTB,
     +        STEPB,IDENB,CUNITB,PNTRB,INO,STATUS) 
C     
C  extract order
C
         NPAR      = 10
         PARAMS(8) = IPAR(3)
         PARAMS(9) = IPAR(4)
         XMIN      = 1.E30
         XMAX      = -XMIN
         DO 10 I = IORD1,IORD2
            
            PARAMS(10) = I
            IORD       = I - IORD1 + 1
            
            IF (METH(1:1) .EQ. 'M') THEN
               CALL EXTR2M(MADRID(PNTRA),NPIXA(1),NPIXA(2),
     +              MADRID(PNTRB),NPIXB(1),NPIXB(2),MADRID(PNTRC),NWDTH,
     +              IORD,NPAR,PARAMS,NA,DPAR,XMIN1,XMAX1,METH,
     +              STARTA,STEPA)
            ELSE
               CALL EXTR1M(MADRID(PNTRA),NPIXA(1),NPIXA(2),
     +              MADRID(PNTRB),NPIXB(1),NPIXB(2),
     +              IORD,NPAR,PARAMS,NA,DPAR,XMIN1,XMAX1,METH,
     +              STARTA,STEPA)
            ENDIF

            XMIN   = AMIN1(XMIN,XMIN1)
            XMAX   = AMAX1(XMAX,XMAX1)
 10      CONTINUE
      ELSE
         CALL STTPUT(' Error in dispersion coefficients ')
         GO TO 20
      END IF
C
C  descriptor handling
C
      CUTS(1) = XMIN
      CUTS(2) = XMAX
      CUTS(3) = XMIN
      CUTS(4) = XMAX
      CUTS(5) = 0.
      CUTS(6) = 0.
      CALL DSCUPT(IMNO,INO,' ',STATUS)
      CALL STDWRR(INO,'LHCUTS',CUTS,1,6,KUN,STATUS)
C
C  free data
C
C      CALL STFCLO(INO,STATUS)
 20   CALL STSEPI
      END

      SUBROUTINE EXTR1M(IN,NPIXA1,NPIXA2,OUT,NB1,NB2,INDX,NPAR,PAR,ND,
     +     DPAR,XMIN,XMAX,METH,STARTA,STEPA)
C
C       EXTRACTION ROUTINE - BILINEAR INTERPOLATION
C
C   ARGUMENTS
C IN     REAL*4  INPUT IMAGE
C NPIXA1 INTG*4  NO. OF SAMPLES
C NPIXA2 INTG*4  NO. OF LINES
C OUT    REAL*4  EXTRACTED ORDERS
C NB1    INTG*4  NO. OF SAMPLES
C NB2    INTG*4  NO. OF ORDERS
C INDX   INTG*4  SEQUENTIAL NO. OF THE ORDER
C NPAR   INTG*4  NO. OF EXTRACTION PARAMETERS
C PAR    REAL*4  EXTRACTION PARAMETERS
C ND     INTG*4  NO. OF COEFFS.
C DPAR   REAL*8  COEFFS. TO DEFINE ORDER POSITION
C XMIN   REAL*4  MINIMUM VALUE
C XMAX   REAL*4  MAXIMUMVALUE
C
      IMPLICIT NONE
C 
      INTEGER  NB1, NB2, NPAR, NPIXA1, NPIXA2, ND
      INTEGER  LOSLIT, HISLIT
      INTEGER  NWW, IWIDTH, IYLIM, NX, N1, N2, K, L, IX, INDX
      INTEGER  IX1, IP, L1, K1, IW
C 
      REAL OUT(NB1,NB2),PAR(NPAR)
      REAL IN(NPIXA1,NPIXA2)
      REAL XMIN,XMAX
C 
      DOUBLE PRECISION DPAR(ND),STARTA(2),STEPA(2)
      DOUBLE PRECISION XVAL,YVAL,AX
      DOUBLE PRECISION WIDTH,SKEW,HWIDTH,OFFSET,DC
      DOUBLE PRECISION YV,WW(50)
      DOUBLE PRECISION SKWDX,SKWDY,SKWOF,SUM
      DOUBLE PRECISION DEGREE,DD,DSIN,DCOS
      DOUBLE PRECISION YLOW, YHIGH, WSUM
C 
      LOGICAL IF1D
C 
      CHARACTER*1      METH
C 
      DOUBLE PRECISION IX_R8R8, XI_R8R8
C
C     DEGREE = RAD / DEGREE
C
      DATA DEGREE/0.017453292519943295769237D+0/
C

      NWW    = 50
C
C ... define extraction constants
C
      WIDTH  = PAR(1)
      IWIDTH = WIDTH
      HWIDTH = 0.5*WIDTH
      SKEW   = PAR(2)
      SKWDX  = 0.0
      SKWDY  = 1.0
      SKWOF  = 0.0
      XMIN   = 0.0
      XMAX   = XMIN
C      OFFSET = PAR(3) - HWIDTH + 0.5
      OFFSET = PAR(3)
      IF (SKEW.NE.0.0) THEN       
         DD      = SKEW*DEGREE
         SKWDX  = -DSIN(DD)     !   DX PER PIXEL
         SKWDY  = DCOS(DD)      !   DY PER PIXEL
         SKWOF  = -SKWDX*HWIDTH !   2D INTERPOLATION
         IF1D   = .FALSE.
      END IF
      AX     = 1.D0 + SKWOF
      IYLIM  = NPIXA2 - IWIDTH
      NX     = SKWOF
      N1     = NX
      N2     = NPIXA1 - NX
C
C ... get ordinate values
C
      K      = PAR(8)
      L      = PAR(9)
      YV     = PAR(10)
C
C ... loop over the samples
C
      DO 10 IX = 1,N1
         OUT(IX,INDX) = 0.
 10   CONTINUE

      DO 100 IX = N2 + 1,NPIXA1
         OUT(IX,INDX) = 0.
 100  CONTINUE

      DO 50 IX1 = N1 + 1,N2

         IP     = 0
         DC     = 1.D0
         YVAL   = 0.D0
         XVAL   = IX_R8R8(AX,STARTA(1),STEPA(1))

C Computes position of the center of the slit from the coefficients of
C the bivariate polynomial stored in table order.

         DO 30 L1 = 0,L
            IP     = IP + 1
            WW(IP) = DC
            YVAL   = YVAL + WW(IP)*DPAR(IP)
            DO 20 K1 = 1,K
               IP     = IP + 1
               WW(IP) = WW(IP-1)*XVAL
               YVAL   = YVAL + WW(IP)*DPAR(IP)
 20         CONTINUE
            DC     = DC*YV
 30      CONTINUE
C
C Corrects the central position for the offset
C OFFSET does not take into account HWIDTH 
C
         YVAL   = XI_R8R8(YVAL,STARTA(2),STEPA(2)) + OFFSET
C     
         YLOW   = YVAL - HWIDTH
         YHIGH  = YVAL + HWIDTH
         LOSLIT = INT(YLOW)
         HISLIT = INT(YHIGH)

         IF (LOSLIT.LT.1)      LOSLIT = 1
         IF (HISLIT.GT.NPIXA2) HISLIT = NPIXA2

         SUM    = 0.D0
         WSUM   = 0.D0
         IX     = INT(AX)
C
         IF (LOSLIT.EQ.HISLIT) THEN
            SUM = (YHIGH-YLOW)*IN(IX,LOSLIT)
C
         ELSE
C
            DO 40 IW = LOSLIT,HISLIT
               
               IF (IW.EQ.LOSLIT) THEN
                  SUM = SUM + IN(IX,IW)*(LOSLIT+1.-YLOW)
                  
               ELSE IF (IW.EQ.HISLIT) THEN
                  SUM = SUM + IN(IX,IW)*(YHIGH-HISLIT)

               ELSE
                  SUM = SUM + IN(IX,IW)

               ENDIF


 40         CONTINUE
C
         ENDIF
C
         IF (METH(1:1) .EQ. 'L') THEN ! Linear extraction
            OUT(IX1,INDX) = SUM
         ELSE                   ! Average extraction
            OUT(IX1,INDX) = SUM/WIDTH
         ENDIF
         XMIN   = AMIN1(XMIN,OUT(IX1,INDX))
         XMAX   = AMAX1(XMAX,OUT(IX1,INDX))
         AX     = AX + 1.D0
 50   CONTINUE

      RETURN

      END



      SUBROUTINE EXTR2M(IN,NPIXA1,NPIXA2,OUT,NB1,NB2,COL,NWDTH,INDX,NPAR
     $     ,PAR,ND,DPAR,XMIN,XMAX,METH,STARTA,STEPA)
C
C       EXTRACTION ROUTINE - MEDIAN
C
C   ARGUMENTS
C IN     REAL*4  INPUT IMAGE
C NPIXA1 INTG*4  NO. OF SAMPLES
C NPIXA2 INTG*4  NO. OF LINES
C OUT    REAL*4  EXTRACTED ORDERS
C NB1    INTG*4  NO. OF SAMPLES
C NB2    INTG*4  NO. OF ORDERS
C INDX   INTG*4  SEQUENTIAL NO. OF THE ORDER
C NPAR   INTG*4  NO. OF EXTRACTION PARAMETERS
C PAR    REAL*4  EXTRACTION PARAMETERS
C ND     INTG*4  NO. OF COEFFS.
C DPAR   REAL*8  COEFFS. TO DEFINE ORDER POSITION
C XMIN   REAL*4  MINIMUM VALUE
C XMAX   REAL*4  MAXIMUMVALUE
C
      IMPLICIT NONE
C 
      INTEGER  NB1, NB2, NPAR, NPIXA1, NPIXA2, NWDTH, NWDTH2, ND
      INTEGER  NWW, IWIDTH, IYLIM, NX, N1, N2, K, L, IX, INDX
      INTEGER  IX1, IP, L1, K1, IW
      INTEGER  LOSLIT, HISLIT
C 
      REAL OUT(NB1,NB2),PAR(NPAR)
      REAL IN(NPIXA1,NPIXA2),COL(NWDTH)
      REAL XMIN,XMAX,LIMITS(2)
C 
      DOUBLE PRECISION DPAR(ND),STARTA(2),STEPA(2)
      DOUBLE PRECISION XVAL,YVAL,AX
      DOUBLE PRECISION WIDTH,SKEW,HWIDTH,OFFSET,DC
      DOUBLE PRECISION YV,WW(50)
      DOUBLE PRECISION SKWDX,SKWDY,SKWOF,SUM
      DOUBLE PRECISION DEGREE,DD,DSIN,DCOS
      DOUBLE PRECISION YLOW, YHIGH, WSUM
      DOUBLE PRECISION IX_R8R8, XI_R8R8
C 
      LOGICAL IF1D
C 
      CHARACTER*1      METH
      CHARACTER*80   OUTPUT
C
C     DEGREE = RAD / DEGREE
C
      DATA DEGREE/0.017453292519943295769237D+0/
      DATA LIMITS/0.0,0.0/
C     

      NWW    = 50
C
C ... define extraction constants
C
      WIDTH  = PAR(1)
      IWIDTH = WIDTH
      HWIDTH = 0.5*WIDTH
      NWDTH2 = HWIDTH
      SKEW   = PAR(2)
      SKWDX  = 0.0
      SKWDY  = 1.0
      SKWOF  = 0.0
      XMIN   = 0.0
      XMAX   = XMIN
C      OFFSET = PAR(3) - HWIDTH + 0.5
      OFFSET = PAR(3)
      IF (SKEW.NE.0.0) THEN       
          DD      = SKEW*DEGREE
          SKWDX  = -DSIN(DD)      !   DX PER PIXEL
          SKWDY  = DCOS(DD)       !   DY PER PIXEL
          SKWOF  = -SKWDX*HWIDTH  !   2D INTERPOLATION
          IF1D   = .FALSE.
      END IF
      AX     = 1.D0 + SKWOF
      IYLIM  = NPIXA2 - IWIDTH
      NX     = SKWOF
      N1     = NX
      N2     = NPIXA1 - NX
C
C ... get ordinate values
C
      K      = PAR(8)
      L      = PAR(9)
      YV     = PAR(10)
C
C ... loop over the samples
C
      DO 10 IX = 1,N1
         OUT(IX,INDX) = 0.
 10   CONTINUE

      DO 100 IX = N2 + 1,NPIXA1
         OUT(IX,INDX) = 0.
 100  CONTINUE

      DO 50 IX1 = N1 + 1,N2

         IP     = 0
         DC     = 1.D0
         YVAL   = 0.D0
         XVAL   = IX_R8R8(AX,STARTA(1),STEPA(1))

C Computes position of the center of the slit from the coefficients of
C the bivariate polynomial stored in table order.

         DO 30 L1 = 0,L
            IP     = IP + 1
            WW(IP) = DC
            YVAL   = YVAL + WW(IP)*DPAR(IP)
            DO 20 K1 = 1,K
               IP     = IP + 1
               WW(IP) = WW(IP-1)*XVAL
               YVAL   = YVAL + WW(IP)*DPAR(IP)
 20         CONTINUE
            DC     = DC*YV
 30      CONTINUE
C
C Corrects the central position for the offset
C OFFSET does not take into account HWIDTH 
C
         YVAL   = XI_R8R8(YVAL,STARTA(2),STEPA(2)) + OFFSET
C     
         YLOW   = YVAL - HWIDTH
         YHIGH  = YVAL + HWIDTH
         LOSLIT = INT(YLOW)
         HISLIT = INT(YHIGH)

         IF (LOSLIT.LT.1)      LOSLIT = 1
         IF (HISLIT.GT.NPIXA2) HISLIT = NPIXA2

         SUM    = 0.D0
         WSUM   = 0.D0
         IX     = INT(AX)
C
         IF (LOSLIT.EQ.HISLIT) THEN
            OUT(IX1,INDX) = IN(IX,LOSLIT)
         ELSE
            K1 = 0
            DO 40 IW = LOSLIT,HISLIT
               K1 = K1 + 1
               IF (IW.EQ.LOSLIT) THEN
                  COL(K1) = IN(IX,IW)!*(LOSLIT+1.-YLOW)
               ELSE IF (IW.EQ.HISLIT) THEN
                  COL(K1) = IN(IX,IW)!*MIN(1.0,YHIGH-HISLIT)
               ELSE
                  COL(K1) = IN(IX,IW)
               ENDIF
 40         CONTINUE            !determine MEDIAN:
            L1 = K1 / 2
            IF (K1.GT.0) THEN
               CALL GET_MEDIAN(COL,NWDTH,K1,L1,OUT(IX1,INDX),IW)
               IF (IW.LT.0) THEN
                  WRITE(OUTPUT,12345) IX1,INDX
                  CALL STTPUT(OUTPUT)
               ENDIF
            ELSE
               OUT(IX1,INDX) = -1.0
            ENDIF
         ENDIF
C
         XMIN   = AMIN1(XMIN,OUT(IX1,INDX))
         XMAX   = AMAX1(XMAX,OUT(IX1,INDX))
         AX     = AX + 1.D0
 50   CONTINUE
C
      RETURN
C 
12345 FORMAT('ERROR at ',I5,',',I5,' occurred')
C
      END

      SUBROUTINE GET_MEDIAN(RA,N,NUSE,LL,VALU,STATUS)

      IMPLICIT NONE
C 
      INTEGER   N,LL,STATUS
      INTEGER   NUSE
C 
      REAL      RA(N),VALU
C 
      VALU = 0.0
      STATUS = 0
C 
      IF (NUSE.LE.3) THEN
         IF (NUSE.LT.1) THEN
            STATUS = -1
            RETURN
         ELSE IF (NUSE.LE.2) THEN
            STATUS = 2
            VALU = RA(1)
c            VALU = (RA(1) + RA(2))/2.0               !for ndim = 1,2, median is 1st elem.
            RETURN
         ELSE 
            STATUS = 3
            VALU = RA(2)               !for ndim = 3, median is 2nd elem.
            RETURN
         ENDIF
      ENDIF
C 
      CALL NSORT(RA,NUSE,LL,VALU)
      RETURN
C 
      END
