C @(#)necripp1.for	19.1 (ESO-DMD) 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 ECHRP1
C
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.IDENTIFICATION
C
C  program  ECHRIPP1              version 1.0       870820
C
C  Author :   J.D.Ponz            ESO - Garching
C
C.MODIFICATIONS
C
C  version 1.1   P. Ballester   901009   Option to freeze K and Alpha
C
C.KEYWORDS
C
C  ECHELLE, CASPEC, BLAZE FUNCTION
C
C.PURPOSE
C
C  COMPUTE THE ECHELLE CONSTANTS ALPHA AND K FOR EACH ECHELLE ORDER
C  (only if the flag is not set to FREEZE)
C
C.ALGORITHM
C
C  a new one...
C
C.INPUT/OUTPUT
C
C  RIPPLE/ECHELLE INPUT OUTPUT  K0,A0
C
C.VERSION
C 
C 010713		last modif
C 
C-------------------------------------------------------------
C 
       IMPLICIT  NONE
C 
       INTEGER          NPIXA(2),NPTOT(100),NORDER(100)
       INTEGER          I, ISTAT, KUN, KNUL, INIMA, OUTIMA
       INTEGER          MADRID(1), NAXISA
       INTEGER*8        IPNTRA, IPNTRB
C 
       DOUBLE PRECISION WSTART(100),STEPA(2),STARTA(2)
       DOUBLE PRECISION K0, ALPHA0,LAMBDA1,LAMBDA2
C 
       REAL             CUT(4), PARS(4), XMIN, XMAX
C 
       CHARACTER*80     INFRM, OUTFRM, TABLE, FREEZE
       CHARACTER        IDENTA*72,CUNITA*64
C 
       INCLUDE          'MID_INCLUDE:ST_DEF.INC'
       COMMON/VMR/MADRID
       INCLUDE          'MID_INCLUDE:ST_DAT.INC'
C
C ... INITIALIZE SYSTEM
C
       CALL STSPRO('ECHRP1')
       CALL STKRDC('P1',1,1,80,I,INFRM,KUN,KNUL,ISTAT)
       CALL STKRDC('P2',1,1,80,I,OUTFRM,KUN,KNUL,ISTAT)
       CALL CLNFRA(INFRM,INFRM,0)
       CALL CLNFRA(OUTFRM,OUTFRM,0)
       CALL STKRDR('INPUTR',1,4,I,PARS,KUN,KNUL,ISTAT)
       TABLE = 'riptab'
C       CALL STKRDC('P4',1,1,80,I,TABLE,KUN,KNUL,ISTAT)
       CALL STKRDC ('P5',1,1,80,I,FREEZE,KUN,KNUL,ISTAT)
       CALL STTPUT(' Ripple correction. Method FIT',ISTAT)
       LAMBDA1 = PARS(1)
       LAMBDA2 = PARS(2)
       K0     = PARS(3)
       ALPHA0 = PARS(4)
C
C ... MAP INPUT/OUTPUT FRMS
C
       CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .             2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     .             CUNITA,IPNTRA,INIMA,ISTAT)
       CALL STDRDR(INIMA,'LHCUTS',1,4,I,CUT,KUN,KNUL,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)
       CALL STDRDI(INIMA,'NORDER',1,NPIXA(2),I,NORDER,KUN,KNUL,ISTAT)
       CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .             NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     .             CUNITA,IPNTRB,OUTIMA,ISTAT)
C
C ... RIPPLE FUNCTION
C
       CALL RIPPLE1(NPIXA(1),NPIXA(2),MADRID(IPNTRA),MADRID(IPNTRB),
     .      WSTART,STEPA(1),NPTOT,NORDER,K0,ALPHA0,XMIN,XMAX,TABLE,
     .      LAMBDA1,LAMBDA2,FREEZE)
       XMIN = CUT(3)
       XMIN = CUT(4)
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)
       CALL STSEPI
       END

       SUBROUTINE RIPPLE1(NX,NY,X,Y,WS,WSTEP,NP,NO,K0,A0,
     .                        XMIN,XMAX,TABLE,LAMBDA1,LAMBDA2,FREEZE)
C
C       CORRECT FOR BLAZE EFFECT USING THE OVERLAPPING REGION
C
C       INPUT ARGUMENTS
C       NX       INTG       MAXIMUM NUMBER OF SAMPLES PER ORDER
C       NY       INTG       NUMBER OF ORDERS
C       X       REAL       X(NX,NY) INPUT IMAGE
C       WS       REAL       STARTING WAVELENGTH PER EACH ORDER
C       WSTEP       REAL       WAVELENGTH STEP
C       NP       INTG       NUMBER OF USEFUL SAMPLES PER ORDER
C       NO       INTG       SPECTRAL ORDER NUMBER
C       K0       REAL       INITIAL GUESS FOR THE CONSTANT K
C       A0       REAL       INITIAL GUESS FOR ALPHA
C       FREEZE   CHARACTER  Flag (NULL/FREEZE) to freeze values of K and Alpha
C
C       OUTPUT ARGUMENTS
C       Y       REAL       Y(NX,NY) CORRECTED IMAGE
C       XMIN       REAL       MINIMUN AND
C       XMAX       REAL       MAXIMUM VALUES
C
       IMPLICIT NONE
C 
       INTEGER           NX, NY
       INTEGER           NP(NY),NO(NY),IERR(2),ICOL(7), TID
       INTEGER           ISTAT, ISEQ, I, IORDER, I1, I2, ILR
C 
       REAL              X(NX,NY), Y(NX,NY), XMIN, XMAX
       REAL              VAL(7)
C 
       DOUBLE PRECISION  WS(NY), WSTEP, K0, A0,LAMBDA1,LAMBDA2
       DOUBLE PRECISION  K(2), KM, AM, A(2)
       DOUBLE PRECISION  DVAL1, DVAL2
C 
       LOGICAL           NULL(2)
C 
       CHARACTER*16      LABEL(7),UNIT(7),FORM
       CHARACTER*80      LINE,FREEZE
       CHARACTER*(*)     TABLE
C 
       INCLUDE          'MID_INCLUDE:ST_DEF.INC'
       INCLUDE          'MID_INCLUDE:ST_DAT.INC'
       DATA              FORM/'G12.4'/
       DATA              LABEL(1)/'ORDER'/,UNIT(1)/'UNITLESS'/
       DATA              LABEL(2)/'K'/,UNIT(2)/'UNITLESS'/
       DATA              LABEL(3)/'ALPHA'/,UNIT(3)/'UNITLESS'/
       DATA              LABEL(4)/'KFIT'/,UNIT(4)/'UNITLESS'/
       DATA              LABEL(5)/'AFIT'/,UNIT(5)/'UNITLESS'/
       DATA              LABEL(6)/'YMAX'/,UNIT(6)/'FLUX'/
       DATA              LABEL(7)/'STATUS'/,UNIT(7)/'UNITLESS'/
       DATA              TID /-1/
C
       ISEQ = 0
       XMIN = 1.E30
       XMAX = -XMIN
       IF (K0.LE.0.0) THEN                     
C
C ... CORRECT RIPPLE IF K0 IS NEGATIVE
C
         CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT)
         CALL TBLSER(TID,LABEL(2),ICOL(1),ISTAT)
         CALL TBLSER(TID,LABEL(3),ICOL(2),ISTAT)
         DO 10 IORDER=1, NY
           CALL TBRRDR(TID,IORDER,2,ICOL,VAL,NULL,ISTAT)
           DVAL1 = VAL(1)
           DVAL2 = VAL(2)
           CALL RIPCOR(WS(IORDER),WSTEP,NP(IORDER),X(1,IORDER),
     .                 Y(1,IORDER),NO(IORDER),DVAL1,DVAL2,NX)
10       CONTINUE
       ELSE
C
C ... FIND BLAZE PARAMETERS
C
         CALL TBTINI(TABLE,F_TRANS,F_O_MODE,8,NY,TID,ISTAT)
         DO 20 I = 1, 6
           CALL TBCINI(TID,D_R4_FORMAT,1,FORM,UNIT(I),LABEL(I),
     .                 ICOL(I),ISTAT)
20       CONTINUE
C
C ... ITERATION ON ALL CONSECUTIVE PAIRS OF ORDERS
C
         CALL STTPUT(' ORDER NUMER  FITTED K  FITTED ALPHA',ISTAT)
         CALL STTPUT(' ----------- ---------- ------------',ISTAT)
         DO 40 IORDER = 2, NY-1

         IF (FREEZE(1:1).EQ.'F'.OR.FREEZE(1:1).EQ.'f') THEN

           KM = K0
           AM = A0

         ELSE

           I1 = IORDER-2
C
C ... CONSIDER LEFT AND RIGHT OVERLAPPING REGIONS
C
           DO 30 ILR = 0, 1
             I1 = I1 + 1
             I2 = I1 + 1
             CALL FITKA(WS(I1),WSTEP,NP(I1),X(1,I1),NO(I1),     
     .                  WS(I2),NP(I2),X(1,I2),NO(I2),
     .                  K0,A0,K(ILR+1),A(ILR+1),IERR(ILR+1),
     .                  LAMBDA1,LAMBDA2)
30         CONTINUE
C
C ... COMPUTE MEAN VALUE
C
           KM = 0.5*(K(1) + K(2))
           AM = 0.5*(A(1) + A(2))

         ENDIF		! End of test on FREEZE

C           TYPE *,IORDER,KM,AM
C
C ... CORRECT FOR RIPPLE
C
           CALL RIPCOR(WS(IORDER),WSTEP,NP(IORDER),X(1,IORDER),
     .                     Y(1,IORDER),NO(IORDER),KM,AM,NX)
           IF (IORDER.EQ.2) THEN
             CALL RIPCOR(WS(1),WSTEP,NP(1),X(1,1),Y(1,1),
     .                       NO(1),KM,AM,NX)
             VAL(1) = NO(1)
             VAL(2) = KM
             VAL(3) = AM
             VAL(4) = KM
             VAL(5) = AM
             ISEQ = ISEQ + 1
             CALL TBRWRR(TID,ISEQ,6,ICOL,VAL,ISTAT)
           ENDIF
           VAL(1) = NO(IORDER)
           VAL(2) = KM
           VAL(3) = AM
           VAL(4) = KM
           VAL(5) = AM
           ISEQ = ISEQ + 1
           CALL TBRWRR(TID,ISEQ,6,ICOL,VAL,ISTAT)
           IF (IORDER.EQ.NY-1) THEN
             CALL RIPCOR(WS(NY),WSTEP,NP(NY),X(1,NY),Y(1,NY),
     .                   NO(NY),KM,AM,NX)
             VAL(1) = NO(NY)
             VAL(2) = KM
             VAL(3) = AM
             VAL(4) = KM
             VAL(5) = AM
             ISEQ = ISEQ + 1
             CALL TBRWRR(TID,ISEQ,6,ICOL,VAL,ISTAT)
           ENDIF
C
C ... WRITE OUT RESULTS
C
           WRITE(LINE,100) NO(IORDER),KM,AM
           CALL STTPUT(LINE,ISTAT)
40       CONTINUE
         CALL TBSINI(TID,ISTAT)
       ENDIF
       CALL TBTCLO(TID,ISTAT)
       RETURN
100    FORMAT(4X,I8,F12.2,F12.4)
       END

       SUBROUTINE RIPCOR(W1,WS,NP,X,Y,M,K,A,NX)
C
C       CORRECT FOR RIPPLE FUNCTION
C
C       INPUT ARGUMENTS
C       W1      INITIAL WAVELENGTH
C       WS      WAVELENGTH STEP
C       NP      NO. OF POINTS
C       X       INPUT FLUX
C       M       ORDER NUMBER
C       K       K CONSTANT
C       A       ALPHA CONSTANT
C
C       OUTPUT ARGUMENTS
C       Y       CORRECTED FLUX
C
       IMPLICIT NONE
       INTEGER          NP, NX, M
       REAL             X(NP),Y(NX)
       DOUBLE PRECISION W1,WS,K,A
       DOUBLE PRECISION PI,DK,DA,DW,DM,R,DX,PA,DC
       INTEGER          I
       DATA             PI/3.1415926535897932D0/
C
       DK = K
       DA = A
       DM = M
       PA = PI*A
       DC = DM/DK
       DO 10 I = 1, NP
         DW   = W1 + (I-1)*WS
         DX   = (PA*DM*DC)*(DW-1.D0/DC)
         IF (DABS(DX).LT.1.D-10) THEN
           Y(I) = X(I)
         ELSE                  
           R    = DSIN(DX)/DX
           Y(I) = X(I)/(R*R) 
         ENDIF
10     CONTINUE
       DO 20 I = NP+1, NX
         Y(I) = 0.
20     CONTINUE
       RETURN
       END

       SUBROUTINE FITKA(WSTART1,WSTEP,N1,X1,NO1,     
     .                  WSTART2,N2,X2,NO2,K0,A0,K,A,IERR,
     .                  LAMBDA1,LAMBDA2)
C
C       COMPUTE CONTANT K AND ALPHA
C
C       INPUT ARGUMENTS
C       WSTEP    WAVELENGTH STEP
C       WSTART1  STARTING WAVELENGTH ORDER M+1
C       N1       NUMBER OF POINTS IN ORDER M+1
C       X1       FLUX IN ORDER M+1
C       NO1      ORDER M+1
C       WSTART2  STARTING WAVELENGTH ORDER M
C       N2       NUMBER OF POINTS IN ORDER M
C       X2       FLUX IN ORDER M
C       NO2      ORDER M
C       K0       INITIAL K VALUE
C       A0       ALPHA
C
C       OUTPUT ARGUMENTS
C       K       FINAL K VALUE
C       A       FINAL A VALUE
C       IERR    CONVERGENCE FLAG 0 - GOOD VALUE
C                               OTHERWISE IT DOES NOT CONVERGE
C
C
       IMPLICIT NONE
       INTEGER           N1, N2
       REAL              X1(N1), X2(N2)
       DOUBLE PRECISION  WSTART1, WSTART2, WSTEP
       DOUBLE PRECISION  DW1, DW2, W2, K0, K, A0, A,LAMBDA1,LAMBDA2
       DOUBLE PRECISION  ETA, XTOL, STEPMX, XPAR(2), FSUMSQ
       DOUBLE PRECISION  FVEC(300), FJAC(300,2),S(2), V(2,2), W(1500)
       DOUBLE PRECISION  XX1(300),XX2(300),XSTART,XSTEP
       INTEGER           ORD1, ORD2, IW(1), ISTAT, NPIX, I1, I2
       INTEGER           NO1, NO2, J, J1, IPRINT, NPAR, MAXCAL
       INTEGER           IERR, NF, NITER, IFAIL, LV, LW, LJ, LIW
C       INTEGER           I
C       CHARACTER*80      LINE
C
       DOUBLE PRECISION  DSQRT, X02AAF
       EXTERNAL          LSQFUN, LSQMON
       COMMON/E04PAR/XSTART,XSTEP,ORD1,ORD2,XX1,XX2
C
       DW1  = (WSTART1 + (N1-1)*WSTEP) - WSTART2
       IF (DW1.LE.0.D0) THEN
         CALL STTPUT('Warning: There is no order overlap',ISTAT)
         RETURN
       ENDIF
C       DW2  = DW1*0.8D0
C       W2   = DW1*0.1D0
       DW2 = LAMBDA2
       W2 = LAMBDA1
       IF ((DW2+W2).GT.DW1) THEN 
          CALL STTPUT('Warning: Wrong wavelengths !',ISTAT)
C          TYPE *,DW1
       ENDIF
C       TYPE *,' *** ORDERS :',NO1,NO2,' ***',W2,DW2
       NPIX = DW2/WSTEP
       NPIX = MIN(NPIX,300)
       I1   = (WSTART2+W2-WSTART1)/WSTEP 
       I2   = W2/WSTEP 
C
C ... PREPARE INPUT
C
       ORD1 = NO1
       ORD2 = NO2
       XSTEP  = WSTEP
       XSTART = WSTART2 + W2 
       DO 10 J = 1, NPIX
         J1 = J-1
         XX1(J) = X1(I1 + J1)
         XX2(J) = X2(I2 + J1)
10     CONTINUE
C100    CONTINUE
C
C ... START MINIMIZATION 
C
       IPRINT = 1
       NPAR   = 2
       MAXCAL = 50*NPAR
       ETA    = 0.9
       XTOL   = 10.*DSQRT(X02AAF(XTOL))
       STEPMX = 10.
       XPAR(1)= K0
       XPAR(2)= A0
       LV     = 2
       LJ     = 300
       LIW    = 1
       LW     = 1500
       IFAIL  = 1
       CALL E04GDF(NPIX,NPAR,LSQFUN,LSQMON,IPRINT,MAXCAL,ETA,XTOL,
     .      STEPMX,XPAR,FSUMSQ,FVEC,FJAC,LJ,S,V,LV,NITER,NF,
     .      IW,LIW,W,LW,IFAIL)
C
C ... TYPE RESULTS
C
C       TYPE *,'IFAIL,SUMSQ',IFAIL,FSUMSQ
C       TYPE *,'K,A',XPAR(1),XPAR(2)
       K = XPAR(1)
       A = XPAR(2)
       RETURN
       END

       SUBROUTINE LSQFUN(IFLAG,M,N,XC,FVECC,FJACC,LJC,IW,LIW,W,LW)
C
C       SUBROUTINE TO COMPUTE THE RESIDUALS AND THEIR FIRST DERIVATIVE
C
C
       IMPLICIT NONE
       INTEGER              IFLAG, M, N, LJC, LIW, LW
C
       DOUBLE PRECISION     FJACC(LJC,N), FVECC(M), W(LW), XC(N)       
C
       DOUBLE PRECISION     XX1(300),XX2(300),XSTART,XSTEP
       INTEGER              ORD1, ORD2, IW(LIW), I
C
       DOUBLE PRECISION     PA, PI, WL, PA2, P2
       DOUBLE PRECISION     X10, X1, X13, S1, C1, R1, DKR1, DAR1
       DOUBLE PRECISION     X20, X2, X23, S2, C2, R2, DKR2, DAR2
C
       COMMON/E04PAR/XSTART,XSTEP,ORD1,ORD2,XX1,XX2
       DATA              PI/3.1415926535897932D0/
C
       PA = PI*XC(2)
       PA2 = 2.D0*PA
       P2  = 2.D0*PI
       DO 10 I = 1, M
         WL  = XSTART + (I-1)*XSTEP
         X10 = ORD1 - XC(1)/WL
         X1  = PA*X10
         X13 = X1*X1*X1
         S1  = DSIN(X1)
         C1  = DCOS(X1)
         R1  = S1/X1
         X20 = ORD2 - XC(1)/WL
         X2  = PA*X20
         X23 = X2*X2*X2
         S2  = DSIN(X2)
         C2  = DCOS(X2)
         R2  = S2/X2
         IF (IFLAG.EQ.2) FVECC(I) = R1*R1/XX1(I) - R2*R2/XX2(I)
         DKR1 = (PA2/(WL*X13))*(S1*S1-X1*S1*C1) 
         DKR2 = (PA2/(WL*X23))*(S2*S2-X2*S2*C2)
         DAR1 = (X1*S1*C1-S1*S1)*P2*X10/X13
         DAR2 = (X2*S2*C2-S2*S2)*P2*X20/X23
         FJACC(I,1) = DKR1/XX1(I) - DKR2/XX2(I)
         FJACC(I,2) = DAR1/XX1(I) - DAR2/XX2(I)
10     CONTINUE
       RETURN
       END

       SUBROUTINE LSQMON(M,N,XC,FVECC,FJACC,LJC,S,IGRADE,NITER,NF,
     .                       IW,LIW,W,LW)
C
C       MONITORING SUBROUTINE
C
       IMPLICIT NONE
       INTEGER          M, N, LJC, IGRADE, NITER, NF, LIW, LW
       DOUBLE PRECISION FJACC(LJC,N), FVECC(M), W(LW), XC(N), S(N)
       INTEGER          IW(LIW)
C
C       DOUBLE PRECISION FSUMSQ, GTG, G(50), F01DEF
C       INTEGER          J
C
C       FSUMSQ = F01DEF(FVECC,FVECC, M)
C       CALL LSQGRD(M, N, FVECC, FJACC, LJC, G)
C       GTG    = F01DEF(G,G,N)
C       WRITE(6, 99999) NITER, NF, FSUMSQ, GTG, IGRADE
C       WRITE(6, 99998)
C       DO 20 J = 1,  N
C          WRITE(6, 99997) XC(J), G(J), S(J)
C20     CONTINUE
       RETURN
C99999  FORMAT(///'    ITNS',4X,'F  EVALS',10X, ' SUMSQ',13X,'GTG',
C     . 8X,'GRADE'/1X I4,6X,I5,6X,1PE13.5,6X,1PE9.1,6X,I3)
C99998  FORMAT(/8X,'X',20X,'G',11X,'SINGULAR VALUES')
C99997  FORMAT(1X,1PE13.5,10X,1PE9.1,10X,1PE9.1)
       END

       SUBROUTINE LSQGRD(M, N, FVECC, FJACC, LJC, G)
C
       IMPLICIT NONE
       INTEGER LJC, M, N
       DOUBLE PRECISION FVECC(M), FJACC(LJC, N), G(N)
       DOUBLE PRECISION SUM
       INTEGER          J, I
C
       DO 40 J = 1, N
          SUM = 0.D0
          DO 20 I =1, M
             SUM = SUM + FJACC(I, J)*FVECC(I)
20        CONTINUE
          G(J) = SUM + SUM
40     CONTINUE
C
       RETURN
       END
