C @(#)neczebpar.for 17.1.1.1 (ES0-DMD) 01/25/02 17:51:32 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 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1991 Special Astrophisical Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, 26 Aug 1991 C C.LANGUAGE: F77+ESOext C C.AUTHOR: A.Yu.Kniazev, V.Shergin C C.IDENTIFICATION C ZEB_PAR.FOR C C.KEYWORDS C echelle (ZEBRA), parameters C C.PURPOSE C Determination of the parameters for echelle images (ZEBRA) C in single observation set. C C C.ALGORITHM C C.INPUT/OUTPUT C input keywords are: C ZEBRAC(1:8) INPUT - IMAGE NAME C ZEBRAC(20) INPUT - PLOT KEY C ZEBRAI(1) INPUT - the number of column for begin work C ZEBRAI(2) INPUT - the half of zone for begin work C C Program output parameters in keywords: C ZEBRAR(1) - the mean slope of echelle spectra C ZEBRAR(2) - the mean width of orders C C COMMAND : C DETERMIN/ZEB INPUT C C----------------------------------------------------------------- PROGRAM ZEBPAR IMPLICIT NONE INTEGER ISLOST, LWORD REAL ACCUR PARAMETER (ISLOST = 12) PARAMETER (ACCUR = 0.005) INTEGER MADRID(1) INTEGER ACTVAL,I,K,ISTAT INTEGER NAXIS,NDIM INTEGER IVEC1, IVEC2 INTEGER PNTRA,STATUS,WINDOW INTEGER NPIX(3), KUN, KNUL, IMNO, IJ DOUBLE PRECISION START(3),STEP(3) DOUBLE PRECISION ST CHARACTER*80 CA CHARACTER*72 IDENA CHARACTER*80 LINE CHARACTER*64 CUNITA CHARACTER*8 FRMIN CHARACTER*1 PLOT C REAL INCL, XINCL ! orders inclination coefficient REAL SINCL ! step of change INCL INTEGER ORDWID ! average width of ORDER INTEGER COLUMN(2) ! number of column for begin work INTEGER IDIR, FLAG REAL AM ! amplitude REAL XC ! REAL HWHI ! average width of ORDER (It's main) C REAL CUTS(2) ! minimum, maximum C INTEGER PIX(2) ! index of minimum and maximum pixel REAL TMP(-1:1) ! output value width of order C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C CA = 'COUNTS' ST = 1.D0 NDIM = 3 C C ... initialize system and read parameters C CALL STSPRO('ZEBPAR') CALL STKRDC('IN_A',1,1,60,ACTVAL,FRMIN,KUN,KNUL,STATUS) CALL STKRDI('INPUTI',1,2,ACTVAL,COLUMN,KUN,KNUL,STATUS) CALL STKRDC('INPUTC',1,1,1,ACTVAL,PLOT,KUN,KNUL,STATUS) C CALL STTPUT(' ZEBRA PARAMETERS DETERMINATION',ISTAT) CALL STTPUT(' ------------------------------',ISTAT) CALL STTPUT(' INPUT IMAGE : '//FRMIN,ISTAT) WRITE (LINE,9040) COLUMN(1) CALL STTPUT(LINE ,ISTAT) WRITE (LINE,9050) COLUMN(2) CALL STTPUT(LINE ,ISTAT) CALL STTPUT(' PLOT FLAG : '//PLOT,ISTAT) C C ... read frame C CALL STIGET(FRMIN,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, .NDIM,NAXIS,NPIX,START,STEP,IDENA,CUNITA,PNTRA,IMNO,STATUS) C C ... generate working space C CALL STIPUT('ZEBRA.TMP',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,1, . NPIX(2)*4,ST,ST,CA,CA,IVEC1,IJ,STATUS) IVEC2 = IVEC1 + NPIX(2) C C ... start step-by-step fitting process to find optimal inclination C ... coefficient corresponding to minimum order width WINDOW = NPIX(2)/5 FLAG = 0 IDIR = 0 INCL = 0. XINCL = 0. SINCL = 0.1 C C ... next step (or 3 steps when IDIR = 0) 10 DO I = -1, 1, 1 IF (IDIR .EQ. I .OR. IDIR .EQ. 0) THEN IF (I .NE. 0 .OR. FLAG .EQ. 0) THEN INCL = XINCL + SINCL*I C ... extract and plot "averaging" column C23456789012345678901234567890123456789012345678901234567890123456789012345678 CALL EXTRAV(MADRID(PNTRA),NPIX(1),NPIX(2),COLUMN(1), . COLUMN(2),INCL,MADRID(IVEC1)) IF (PLOT(1:1) .EQ. 'Y' .OR. PLOT(1:1) .EQ. 'y') . CALL PLTARR(MADRID(IVEC1),NPIX(2),100,(NPIX(2)-150),FRMIN,0) C ... calculate "autocorrelation" CALL CONVL(MADRID(IVEC1),1,NPIX(2),1,1,MADRID(IVEC2),WINDOW) IF (PLOT(1:1) .EQ. 'Y' .OR. PLOT(1:1) .EQ. 'y') THEN K = NPIX(2)/2-WINDOW CALL PLTARR(MADRID(IVEC2),NPIX(2),K,K+WINDOW*2,FRMIN,1) ENDIF C ... determine order width ORDWID = LWORD(MADRID(IVEC2),NPIX(2),NPIX(2)/2,AM,XC,HWHI) TMP(I) = HWHI ENDIF ENDIF ENDDO C IF ((TMP(-1) .LT. TMP(0)) .AND. (TMP(-1) .LT. TMP(1))) THEN IDIR = -1 TMP(1) = TMP (0) TMP(0) = TMP (-1) XINCL = INCL - 2*SINCL ELSE IF ((TMP(1) .LT. TMP(0)) .AND. (TMP(1) .LT. TMP(-1))) THEN IDIR = 1 TMP(-1)= TMP (0) TMP(0) = TMP (1) XINCL = INCL ELSE IDIR = 0 FLAG = 1 SINCL = SINCL/2 WRITE (*,*) 'New step =',SINCL ENDIF C C WRITE (LINE,9000) TMP(0)*2 C CALL STTPUT(LINE,ISTAT) WRITE (LINE,9010) XINCL CALL STTPUT(LINE,ISTAT) CALL STTPUT(' ------------------------------',ISTAT) IF (SINCL .GT. ACCUR) GOTO 10 C CALL STKWRR('OUTPUTR',XINCL,1,1,KUN,STATUS) CALL STKWRR('OUTPUTR',TMP(0)*2,2,1,KUN,STATUS) CALL STTPUT('--------------------',ISTAT) WRITE (LINE,9020) TMP(0)*2 CALL STTPUT(LINE,ISTAT) WRITE (LINE,9030) XINCL CALL STTPUT(LINE,ISTAT) C ... that's all ... C ... free data CALL STSEPI 9000 FORMAT ('OrderWidth: ',F8.4) 9010 FORMAT ('Inclination_ of_Order: ',F8.4) 9020 FORMAT ('Min_Order_Width: ',F8.4) 9030 FORMAT ('Optimum_Inclination_of_Order: ',F8.4) 9040 FORMAT ('Center_Column: ',I6) 9050 FORMAT ('Half_of_width: ',I6) END C***************************************************************** SUBROUTINE EXTRAV(X,N1,N2,I,WIDTH,INCL,OUT) C C Extract (with averaging) 2*WIDHT+1 vertical columns C from image array X(N1,N2) centered at I-column C INTEGER N1 ! IN INTEGER N2 ! IN REAL X(N1, N2) ! IN - input image INTEGER I ! IN - position of center column INTEGER WIDTH ! IN - +- from I-column REAL INCL ! IN - orders inclination coefficient REAL OUT(N2) ! OUT- output trace C INTEGER J, J1, K, K1 C DO 10 J = 1, N2 OUT(J) = 0. DO 20 K = I-WIDTH,I+WIDTH J1 = MAX(1, MIN(N2, INT(J + (K-I)*INCL + 0.5))) IF (K .LT. 1) THEN K1 = -(K-1) ELSE IF (K .GT. N1) THEN K1 = N1-(K-N1) ELSE K1 = K ENDIF OUT(J) = OUT(J) + X(K1,J1) 20 CONTINUE OUT(J) = OUT(J)/(2*WIDTH+1) 10 CONTINUE RETURN END C***************************************************************** SUBROUTINE CONVL(X,N1,N2,I1,I2,OUT,NSTEP) C C Perform "convolution" of two vertical columns: I1, I2 C of image array X(N1,N2) C INTEGER N1 ! IN INTEGER N2 ! IN REAL X(N1, N2) ! IN INTEGER I1, I2 ! IN - NUMERICS OF COLUMN REAL OUT(N2) ! OUT INTEGER NSTEP ! IN C INTEGER J, K, N, N22 REAL S, S1 C N22 = N2/2 N = MIN0(NSTEP,N22-1) S1 = 0. DO 10 K = 1, N2 OUT(K) = 0. S1 = S1 + X(I1,K) 10 CONTINUE DO 30 J = -N, N S = 0. DO 20 K = 1, N2 S = S + X(I1,K)*X(I2,MOD(N2+K+J-1,N2)+1) 20 CONTINUE OUT(N22+J) = S/S1 30 CONTINUE RETURN END C***************************************************************** INTEGER FUNCTION LWORD(INP,NP,CENTER,AM,XC,HWHI) C C Determine new order width C INTEGER NP ! IN arr.length REAL INP(NP) ! IN input array INTEGER CENTER ! IN center pixel of contour REAL AM, XC, HWHI ! OUT output params C INTEGER I, K, WND, IC REAL A, B, D, A11, A12, A22, B1, B2 REAL X, Y, Y1, Y2 C C ... at first - rough estimation of HWHI K = NP DO I = CENTER+1,NP-1 IF (INP(I) .LE. INP(CENTER)/2.) THEN K = I GO TO 10 ENDIF IF ((INP(I-1).GE.INP(I)) .AND. (INP(I+1).GE.INP(I))) THEN K = I-2 GO TO 10 ENDIF ENDDO C 10 WND = K - CENTER LWORD = WND C C ... and then try to determine parameters with more precision AM = INP(CENTER) IC = CENTER C ... make sure that we stay at maximum DO I = CENTER-WND+1, CENTER+WND IF (INP(I) .GT. AM) THEN AM = INP(I) IC = I ENDIF ENDDO C ... redefine window of interest Y1 = AM Y = (INP(IC+1) + INP(IC-1))/2. DO I = 1, WND Y2 = (INP(IC+I+1) + INP(IC-I-1))/2. IF (Y1-Y .GE. Y-Y2) THEN WND = I+1 GO TO 20 ENDIF Y1 = Y Y = Y2 ENDDO C ... perform very simple fit 20 A11 = 0. A12 = 0. A22 = WND*2. B1 = 0. B2 = 0. DO I = IC-WND+1, IC+WND AM = MAX(AM, INP(I)) X = I - 0.5 Y = INP(I) - INP(I-1) A11 = A11 + X**2 A12 = A12 + X B1 = B1 + X*Y B2 = B2 + Y ENDDO D = (A11*A22 - A12*A12) A = (B1*A22 - B2*A12)/D B = (B2*A11 - B1*A12)/D XC = -B/A HWHI = SQRT( MAX(0., -AM/A)) WRITE (*,9000) AM, XC, HWHI*2. C RETURN 9000 FORMAT ('Am=',F7.1,' Xc=',F6.1,' FWHI=',F6.3) END C***************************************************************** SUBROUTINE PLTARR(OUT,N2,NB,NE,NAME,OVR) C C Plot elemets from NB to NE of array OUT(N2) with identifier NAME C (for debuging) C INTEGER N2 ! total length of data array REAL OUT(N2) ! array with data to plot INTEGER NB, NE ! edges of plotting part of array INTEGER OVR ! overplot? CHARACTER*(*) NAME ! label C INTEGER ISTAT, IAC, KUN, KNUL INTEGER FMODE, IMAGE(4), NPL REAL XS, SX, YOFF REAL CUTS(2) REAL FRAME(8) REAL MUDAKI(2) CHARACTER LABEL1*80, LABEL2*80, LABEL3*80, LABEL4*80 CHARACTER XFRAME*4, YFRAME*4 C INCLUDE 'MID_INCLUDE:PLTDEC.INC/NOLIST' INCLUDE 'MID_INCLUDE:PLTDAT.INC/NOLIST' C DATA PMODE/-1/ DATA MUDAKI/0.,0./ C CALL STKWRR('PLRSTAT',MUDAKI,19,2,KUN,ISTAT) C CALL STKRDI('PLISTAT',1,1,IAC,FMODE,KUN,KNUL,ISTAT) FMODE = MIN0(FMODE,1) ! we need not too complicated legenda CALL STKRDC('PLCSTAT',1,5,4,IAC,XFRAME,KUN,KNUL,ISTAT) CALL STKRDC('PLCSTAT',1,9,4,IAC,YFRAME,KUN,KNUL,ISTAT) IF (PMODE.LT.0 .OR. OVR.EQ.0) THEN PMODE = 0 ELSE PMODE = 2 ENDIF C IMAGE(1) = MAX( 1, NB) IMAGE(2) = MIN(N2, NE) IMAGE(3) = 1 IMAGE(4) = 1 NPL = MAX( 2, IMAGE(2) - IMAGE(1) + 1 ) XS = 1. SX = 1. C C *** only for plot command IF (PMODE.EQ.0) THEN YOFF = 0.0 C C *** find min and max of the data values C CALL MINMAX(OUT,N2,CUTS(1),CUTS(2)) C C*** get the frame parameters IF (XFRAME(1:4).EQ.'MANU') THEN CALL STKRDR('PLRSTAT',11,4,IAC,FRAME(1),KUN,KNUL,ISTAT) C CALL BOXWTP(FRAME(1),NPL,XS,SX,IMAGE(1)) ELSE XFRAME = 'AUTO' FRAME(1) = FLOAT(IMAGE(1)) FRAME(2) = FLOAT(IMAGE(2)) END IF C IF (YFRAME(1:4).EQ.'MANU') THEN CALL STKRDR('PLRSTAT',15,4,IAC,FRAME(5),KUN,KNUL,ISTAT) ELSE YFRAME = 'AUTO' FRAME(5) = CUTS(1) FRAME(6) = CUTS(2) END IF C CALL GETFRM(FRAME(1),FRAME(2),FRAME(3),FRAME(4),XFRAME) C CALL GETFRM(FRAME(5),FRAME(6),FRAME(7),FRAME(8),YFRAME) C C WRITE(*,*) 'XFRAME=', XFRAME C WRITE(*,*) 'YFRAME=', YFRAME C WRITE(*,*) 'FRAME=', FRAME(1), FRAME(2), FRAME(3), FRAME(4) C WRITE(*,*) 'FRAME=', FRAME(5), FRAME(6), FRAME(7), FRAME(8) C CALL STKWRR('PLRSTAT',FRAME,11,8,KUN,ISTAT) C C ELSE ! overplot mode C YOFF = YOFF + ABS(CUTS(1) - CUTS(2))/10. C CALL STKRDR('INPUTR',1,1,IAC,YOFF,KUN,KNUL,ISTAT) ENDIF C C *** do the plot setup C CALL PLOPN(NAME,PMODE,FMODE) C C *** do the work C CALL PLKEY(OUT(IMAGE(1)),NPL,IMAGE,XS,SX,YOFF) C C *** plot the legenda IF (FMODE.EQ.1 .AND. PMODE.EQ.0) THEN C CALL PLFRAM(FRAME(1),FRAME(2),FRAME(3),FRAME(4), C 2 FRAME(5),FRAME(6),FRAME(7),FRAME(8)) LABEL1 = 'Sequence Number' LABEL2 = 'Value' C IF (FMODE.EQ.1) THEN LABEL3 = 'Array: '//NAME LABEL4 = ' ' C CALL PLLABL(LABEL1,LABEL2,LABEL3,LABEL4) END IF END IF C C C *** good bye and return C CALL PLCLS C RETURN END