C @(#)necripp1.for 17.1.1.1 (ESO-DMD) 01/25/02 17:51:31 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