C @(#)necnorm.for 17.1.1.1 (ES0-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 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding 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 ECHNRM C C+ C C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 12:14 - 11 JAN 1988 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.IDENTIFICATION C C program ECHNORM version 1.5 830610 C C.MODIFICATIONS C C 910128 P. Ballester Define and set variables (test mmake -i) C C.KEYWORDS C C ECHELLE, CASPEC, RIPPLE NORMALIZATION C C.PURPOSE C C FIT A POLYNOMIAL TO THE IMAGE ROWS C C.ALGORITHM C C C.INPUT/OUTPUT C C NORMALIZE/ECHELLE INPUT OUTPUT TABLE METHOD PARS C if METHOD='POLY' then PARS=IDEG C- C IMPLICIT NONE INTEGER MADRID, TID INTEGER I,IAV,IDEG INTEGER*8 IPNTRA,IPNTRB INTEGER*8 IS,ISTAT INTEGER J,NAXISA,NCOL,NITER,IMNO,IMNI REAL CONST,ER,RELAX,RIDEN,STEPA1 REAL WSTR1,XMAX1,XMIN1 C INTEGER NPIXA(2),IC(20),KUN,KNUL REAL RPAR(5),VALUE(20),P0(3) REAL XMIN,XMAX,LHCUTS(4) DOUBLE PRECISION STEPA(2),STARTA(2) DOUBLE PRECISION A(10),DSTP(2),DSTA(2) CHARACTER*8 OUTFRM,INFRM,TABLE CHARACTER*16 UNIT(5),LABEL(5),UNIT1(12),LABEL1(12),LABELT,UNITT CHARACTER*6 FORM(5),FORM1(12),FORMT CHARACTER METH*1 CHARACTER IDENTA*72,CUNITA*64 C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA NCOL/5/ DATA LABEL(1)/'ORDER'/,UNIT(1)/'UNITLESS'/,FORM(1)/'I5'/ DATA LABEL(2)/'Y0'/,UNIT(2)/'DN'/,FORM(2)/'E13.3'/ DATA LABEL(3)/'ALFA'/,UNIT(3)/'UNITLESS'/,FORM(3)/'E13.3'/ DATA LABEL(4)/'BETA'/,UNIT(4)/'UNITLESS'/,FORM(4)/'E13.4'/ DATA LABEL(5)/'ERROR'/,UNIT(5)/'UNITLESS'/,FORM(5)/'E13.4'/ DATA LABEL1(1)/'ORDER'/,UNIT1(1)/'UNITLESS'/,FORM1(1)/'I5'/ DATA LABEL1(2)/'DEGREE'/,UNIT1(2)/'UNITLESS'/,FORM1(2)/'I4'/ DATA LABEL1(3)/'A0'/,UNIT1(3)/'UNITLESS'/,FORM1(3)/'E13.3'/ DATA LABEL1(4)/'A1'/,UNIT1(4)/'UNITLESS'/,FORM1(4)/'E13.3'/ DATA LABEL1(5)/'A2'/,UNIT1(5)/'UNITLESS'/,FORM1(5)/'E13.3'/ DATA LABEL1(6)/'A3'/,UNIT1(6)/'UNITLESS'/,FORM1(6)/'E13.3'/ DATA LABEL1(7)/'A4'/,UNIT1(7)/'UNITLESS'/,FORM1(7)/'E13.3'/ DATA LABEL1(8)/'A5'/,UNIT1(8)/'UNITLESS'/,FORM1(8)/'E13.3'/ DATA LABEL1(9)/'A6'/,UNIT1(9)/'UNITLESS'/,FORM1(9)/'E13.3'/ DATA LABEL1(10)/'A7'/,UNIT1(10)/'UNITLESS'/,FORM1(10)/'E13.3'/ DATA LABEL1(11)/'A8'/,UNIT1(11)/'UNITLESS'/,FORM1(11)/'E13.3'/ DATA LABEL1(12)/'A9'/,UNIT1(12)/'UNITLESS'/,FORM1(12)/'E13.3'/ C C ... INITIALIZE SYSTEM C CALL STSPRO('ECHNRM') CALL STKRDC('P1',1,1,8,IAV,INFRM,KUN,KNUL,ISTAT) CALL STKRDC('P2',1,1,8,IAV,OUTFRM,KUN,KNUL,ISTAT) CALL STKRDC('P3',1,1,8,IAV,TABLE,KUN,KNUL,ISTAT) CALL STKRDR('INPUTR',1,3,IAV,RPAR,KUN,KNUL,ISTAT) METH = 'P' IF (METH.EQ.'S') THEN NITER = NINT(RPAR(3)) ELSE IDEG = NINT(RPAR(1)) IDEG = MIN(IDEG,9) NCOL = IDEG + 3 VALUE(2) = IDEG END IF C C ... MAP INPUT FRM C CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, . 2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA, . IPNTRA,IMNI,ISTAT) C C ... MAP OUTPUT FRM C IDENTA = 'FITTED VALUES :'//INFRM DSTA(1) = STARTA(1) DSTA(2) = STARTA(2) DSTP(1) = STEPA(1) DSTP(2) = STEPA(2) CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, . NAXISA,NPIXA,DSTA,DSTP,IDENTA, . CUNITA,IPNTRB,IMNO,ISTAT) C CALL SXIPUT(OUTFRM,'O',NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA, C + IPNTRB,ISTAT) C C ... INITIALIZE TABLE C CALL TBTINI(TABLE,F_TRANS,F_O_MODE,15,NPIXA(2),TID,ISTAT) DO 10 I = 1,NCOL IF (METH.EQ.'S') THEN UNITT = UNIT(I) LABELT = LABEL(I) FORMT = FORM(I) ELSE UNITT = UNIT1(I) LABELT = LABEL1(I) FORMT = FORM1(I) END IF CALL TBCINI(TID,D_R4_FORMAT,1, . FORMT,UNITT,LABELT,IC(I),ISTAT) 10 CONTINUE C C ... ITERATION ON EACH ECHELLE ORDER C IS = 0 ! DUMPING FACTOR RELAX = RPAR(1) ! ERROR ER = RPAR(2) RIDEN = 1.0 STEPA1 = STEPA(1) XMIN = 0. XMAX = XMIN DO 30 I = 1,NPIXA(2) WSTR1 = WSTR1 ! ORDER NUMBER VALUE(1) = STARTA(2) + (I-1)*STEPA(2) IF (METH.EQ.'S') THEN ! INITAL GUESSES P0(2) = 1. C P0(3) = 1./(WSTR(I)+NPIXA(1)*STEPA(1)*0.5) P0(3) = 1./ (WSTR1+NPIXA(1)*STEPA1*0.5) CONST = VALUE(1)*3.1415926535897932/10800.0 CALL FIT1(NPIXA(1),NPIXA(2),MADRID(IPNTRA),MADRID(IPNTRB), + I,CONST,P0,ER,RELAX,NITER,WSTR1,STEPA1, + VALUE(2),VALUE(4)) ELSE CALL FIT2(NPIXA(1),NPIXA(2),MADRID(IPNTRA),MADRID(IPNTRB), + I,STARTA(1),STEPA(1),IDEG,A,XMAX1,XMIN1) DO 20 J = 1,IDEG + 1 VALUE(J+2) = A(J) 20 CONTINUE END IF XMIN = AMIN1(XMIN,XMIN1) XMAX = AMAX1(XMAX,XMAX1) IS = IS + 1 CALL TBRWRR(TID,IS,NCOL,IC,VALUE,ISTAT) 30 CONTINUE LHCUTS(1) = XMIN LHCUTS(2) = XMAX LHCUTS(3) = XMIN LHCUTS(4) = XMAX CALL STDWRR(IMNO,'LHCUTS',LHCUTS,1,4,KUN,ISTAT) CALL TBTCLO(TID,ISTAT) C C ... END C CALL STSEPI END SUBROUTINE FIT1(NX,NY,IN,OUT,I,CONST,P0,ER,RELAX,NITER,X0,DX,P, + ERROR) C C NONLINEAR FIT ROUTINE INTERFACE C IMPLICIT NONE C C INPUT ARGUMENTS C C NX INTG*4 NUMBER OF SAMPLES C NY INTG*4 NUMBER OF LINES (ORDERS) C IN REAL*4 INPUT IMAGE C OUT REAL*4 FITTED VALUES C I INTG*4 NUMBER OF THE LINE C M REAL*4 SPECTRAL ORDER C P0(3) REAL*4 INITIAL GUESSES C ER REAL*4 ERROR BOUND C RELAX REAL*4 DUMPING FACTOR C NITER INTG*4 NUMBER OF ITERATIONS C X0 REAL*4 STARTING WAVELENGTH C DX REAL*4 WAVELENGTH STEP C C OUTPUT PARAMETERS C C P(3) REAL*4 ESTIMATED PARAMETER C ERROR REAL*4 FINAL ERROR C IMPLICIT NONE INTEGER NX,NY,I,NITER,ISTAT REAL CONST,ER,RELAX,X0,DX,ERROR REAL IN(NX,NY),OUT(NX,NY),P0(3),P(3) CHARACTER*80 LINE P0(1) = IN(NX/2,I) WRITE (LINE,9000) I,P0(1),P0(2),P0(3),ER CALL STTPUT(LINE,ISTAT) CALL LSQSIN(NX,X0,DX,IN(1,I),OUT(1,I),RELAX,NITER,ER,CONST,P0,P, + ERROR) RETURN 9000 FORMAT (' Line :',I4,' Init.pars :',3E10.4,' Error :',E10.4) END SUBROUTINE LSQSIN(NP,X0,DX,Y,YFIT,RELAX,NITER,ER,CONST,P0,P, + ERROR) C IMPLICIT NONE C Modified version of the program LSQFIT for the problem : C C find A and B in the expression C Y = (sinc(pi*M*A*(X*B-1)))**2 C INTEGER NP,NITER,ITER,INDEP,KZ6901,IPAR,ISTAT,I,J,K,NF REAL X0,DX,RELAX,ER,CONST,ERROR,X,Y2,W,DET,ERR REAL Y(NP),YFIT(NP),P(3),P0(3),D(3) DOUBLE PRECISION CHI,PP,Q DOUBLE PRECISION AL(3),B(3),A(3,3) CHARACTER*80 LINE C C ... initialize relevant parameters C ITER = 0 C IWEIGH,RELAX ! one independent variable INDEP = 1 ! three parameters IPAR = 3 ! initial guesses P(1) = P0(1) P(2) = P0(2) P(3) = P0(3) NF = NP - 3 C C ... iteration loop C DO 80 KZ6901 = 1,2147483647 IF( .NOT. (ITER.LT.NITER)) GO TO 90 ITER = ITER + 1 WRITE (LINE,9000) ITER,P(1),P(2),P(3) CALL STTPUT(LINE,ISTAT) C C ... initalize normal equation matrix C DO 20 I = 1,IPAR B(I) = 0.0D0 DO 10 J = 1,IPAR A(I,J) = 0.0D0 10 CONTINUE 20 CONTINUE CHI = 0.0D0 DO 50 I = 1,NP C C ... compute function and derivatives C X = X0 + (I-1)*DX CALL FSINC(X,YFIT(I),CONST,P,D) C C ... compute residuals C Y2 = Y(I) - YFIT(I) W = 1.0 PP = DBLE(W) Q = DBLE(Y2) CHI = CHI + Q*Q*PP C C ... enter derivatives into normal equation matrix C DO 40 J = 1,IPAR AL(1) = DBLE(D(J)) DO 30 K = 1,IPAR A(J,K) = A(J,K) + AL(1)*DBLE(D(K))*PP 30 CONTINUE B(J) = B(J) + PP*Q*AL(1) 40 CONTINUE 50 CONTINUE C C ... invert normal equation C WRITE (*,*) A(1,1),A(1,2),A(1,3),B(1) WRITE (*,*) A(2,1),A(2,2),A(2,3),B(2) WRITE (*,*) A(3,1),A(3,2),A(3,3),B(3) C DO K=1,IPAR C KK=IPAR+1-K C PP=A(1,1) C DO I=2,IPAR C Q=A(I,1) C AL(I)=Q/PP C IF (I.LE.KK) AL(I)=-AL(I) C DO J=2,I C IF (I.EQ.2.AND.J.EQ.2)WRITE(*,*) A(2,2) C A(I-1,J-1)=A(I,J)+Q*AL(J) C IF (I.EQ.2.AND.J.EQ.2)WRITE(*,*) K,Q,AL(J),A(1,1) C ENDDO C ENDDO C A(IPAR,IPAR)=1.0D0/PP C DO I=2,IPAR C A(IPAR,I-1)=AL(I) C ENDDO C ENDDO C C ... fill matrix C C DO I=2,IPAR C DO J=1,I C A(J,I)=A(I,J) C ENDDO C ENDDO CALL MATINV(A,3,DET) IF (DET.EQ.0.0) THEN WRITE (*,*) 'ERROR ------------------------------' END IF WRITE (*,*) DET WRITE (*,*) A(1,1),A(1,2),A(1,3) WRITE (*,*) A(2,1),A(2,2),A(2,3) WRITE (*,*) A(3,1),A(3,2),A(3,3) C C ... compute reduced chi**2 and corrections for unknowns C CHI = CHI/DBLE(NF) DO 70 I = 1,IPAR Q = 0.0D0 DO 60 J = 1,IPAR Q = Q + A(I,J)*B(J) 60 CONTINUE C C ... compute errors of the corrections C D(I) = SNGL(DSQRT(DMAX1(A(I,I)*CHI,0.0D0))) P(I) = P(I) + RELAX*SNGL(Q) ! DEBUG <------------- ! WRITE (*,*) D(I),P(I) 70 CONTINUE ERROR = SNGL(DSQRT(DMAX1(CHI,0.0D0))) LINE = ' ' WRITE (LINE,9010) ERROR CALL STTPUT(LINE,ISTAT) IF (ABS(ERROR).LE.ERR) RETURN 80 CONTINUE 90 CONTINUE RETURN 9000 FORMAT (' Iteration :',I4,' Params :',3E13.4) 9010 FORMAT (' Error :',E13.4) END SUBROUTINE FSINC(X,Y1,CONST,P,D) C IMPLICIT NONE C INTERNAL FUNCTION SINC, SINGULARITY AT (X*P(3)->1.) C REAL CONST,X,Y1,A0,A,A2,SINCA,CA,SA,P,D DIMENSION P(3),D(3) A0 = X*P(3) - 1. IF (ABS(A0).LT.0.1E-20) THEN Y1 = 1. D(1) = 1. D(2) = 0. D(3) = 0. ELSE A = CONST*P(2)*A0 A2 = A*A SINCA = SIN(A)/A CA = COS(A) - SA D(1) = SINCA*SINCA Y1 = P(1)*D(1) D(2) = (2.*P(1)*SINCA*CA)/P(2) D(3) = (2.*P(1)*SINCA*CA*X)/A0 END IF RETURN END SUBROUTINE FIT2(NX,NY,IN,OUT,I,X0,DX,IDEG,A,XMIN,XMAX) C IMPLICIT NONE C C LINEAR FIT ROUTINE INTERFACE C C INPUT ARGUMENTS C C NX INTG*4 NUMBER OF SAMPLES C NY INTG*4 NUMBER OF LINES (ORDERS) C IN REAL*4 INPUT IMAGE C OUT REAL*4 FITTED VALUES C I INTG*4 NUMBER OF THE LINE C IDEG INTG*4 POLYNOMIAL DEGREE C C OUTPUT PARAMETERS C C A REAL*8 COEFFS C IMPLICIT NONE INTEGER NX,NY,I,IDEG,NT,J REAL X0,DX,XMIN,XMAX REAL IN(NX,NY),OUT(NX,NY) DOUBLE PRECISION A(1) NT = NX - NX/50 CALL LSQFIT(NT,X0,DX,IN(1,I),OUT(1,I),IDEG,A,XMIN,XMAX) DO 10 J = NT + 1,NX OUT(J,I) = 0. 10 CONTINUE RETURN END SUBROUTINE LSQFIT(NX,X0,DX,X,FIT,K,A,XMIN,XMAX) IMPLICIT NONE C C SOLVE LEAST SQUARE PROBLEM C INTEGER N,NX,K,L,LL,N1,I,KMIN,K1,J REAL X0,DX,XMIN,XMAX,X1,X2,Y,ERROR REAL X(NX),FIT(NX) DOUBLE PRECISION G,XC,A(1) COMMON /LS/G(50,50),XC(50),N C C ITERATION C L = 0 LL = 0 N1 = 0 N = (K+1)* (L+1) DO 20 I = 10,NX - 20 X1 = X0 + (I-1)*DX X2 = X1 Y = X(I) N1 = N1 + 1 CALL SETUP(LL+1,X1,X2,Y,K,L) IF (LL.NE.0) THEN KMIN = MIN(LL,N+1) DO 10 K1 = 1,KMIN CALL HT(K1,LL+1) 10 CONTINUE END IF LL = MIN(LL+1,N+1) 20 CONTINUE CALL SOLVE ERROR = ABS(G(N+1,N+1)) CALL COMP(NX,X0,DX,FIT,K,L,XMIN,XMAX) DO 30 J = 1,K + 1 A(J) = XC(J) 30 CONTINUE RETURN END SUBROUTINE SETUP(LL,X1,X2,Y,K,L) IMPLICIT NONE C C PREPARE DATA IN WORK SPACE C INTEGER N,LL,K,L,IP,L1,K1 REAL X1,X2,Y DOUBLE PRECISION G,X,DC COMMON /LS/G(50,50),X(50),N IP = 0 DC = 1.D0 C DX1 = DBLE(X1)*0.001D0 C DX2 = DBLE(X2)*0.001D0 C DY = DBLE(Y)*0.001D0 DO 20 L1 = 0,L IP = IP + 1 G(LL,IP) = DC DO 10 K1 = 1,K IP = IP + 1 G(LL,IP) = G(LL,IP-1)*X1 10 CONTINUE DC = DC*X2 20 CONTINUE G(LL,N+1) = Y RETURN END SUBROUTINE HT(K,L1) IMPLICIT NONE C FUNCTION : Apply a Householder transformation to the C two rows K and L1 of G. C C INPUT : * matrix G(1:L1,1:n+1) C * numbers L1 and K of the two rows C * n : rank of the matrix C C OUTPUT : modified rows K and L1 of the matrix G. DOUBLE PRECISION G,X,S,H,B COMMON /LS/G(50,50),X(50),N C Construct the transformation INTEGER N,K,L1,KH,J S = DSQRT(G(K,K)**2+G(L1,K)**2) IF (G(K,K).GE.0) THEN S = -S END IF H = G(K,K) - S G(K,K) = S B = G(K,K)*H C Apply the transformation to the rows K and L1 of G. KH = N + 1 - K IF ((KH.NE.0) .AND. (B.NE.0)) THEN DO 10 J = 1,KH S = (G(K,K+J)*H+G(L1,K+J)*G(L1,K))/B G(K,K+J) = G(K,K+J) + S*H G(L1,K+J) = G(L1,K+J) + S*G(L1,K) 10 CONTINUE END IF END SUBROUTINE SOLVE IMPLICIT NONE C FUNCTION : Solve the set of linear equations which are C represented by the upper triangular part C of the matrix G(1:n,1:n+1). C C INPUT : * matrix G(1:n,1:n+1) C * n : rank of the matrix C C OUTPUT : solution vector x(1:n) DOUBLE PRECISION G,X,S COMMON /LS/G(50,50),X(50),N INTEGER N,I,II,J DO 10 I = 1,N X(I) = G(I,N+1) 10 CONTINUE DO 30 II = 1,N I = N - II + 1 S = 0.D0 IF (I.NE.N) THEN DO 20 J = 2,N + 1 - I S = S + G(I,I+J-1)*X(I+J-1) 20 CONTINUE END IF X(I) = (X(I)-S)/G(I,I) 30 CONTINUE END SUBROUTINE COMP(NPIX1,START,STEP,OUT,K,L,XMIN,XMAX) IMPLICIT NONE C INTEGER N,NPIX1,K,L,IX,IP,L1,K1 REAL START,STEP,XMIN,XMAX REAL OUT(NPIX1) DOUBLE PRECISION WORK(50),DC,YVAL,G,A DOUBLE PRECISION X1,X2,XSTART,XSTEP COMMON /LS/G(50,50),A(50),N C C initialize variables XMIN = 0. XMAX = XMIN XSTART = START XSTEP = STEP C C start loop DO 30 IX = 1,NPIX1 X1 = (IX-1)*XSTEP + XSTART IP = 0 DC = 1.D0 YVAL = 0.D0 DO 20 L1 = 0,L IP = IP + 1 WORK(IP) = DC YVAL = YVAL + WORK(IP)*A(IP) DO 10 K1 = 1,K IP = IP + 1 WORK(IP) = WORK(IP-1)*X1 YVAL = YVAL + WORK(IP)*A(IP) 10 CONTINUE DC = DC*X2 20 CONTINUE OUT(IX) = YVAL XMIN = AMIN1(XMIN,OUT(IX)) XMAX = AMAX1(XMAX,OUT(IX)) 30 CONTINUE RETURN END