C @(#)necford.for 17.1.1.1 (ESO-DMD) 01/25/02 17:51:29 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 ECHFOR C C+ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 21:22 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.IDENTIFICATION C C program ECHSORD version 1.6 830310 C C.MODIFICATIONS C C 910128 P. Ballester Define and set variables for mmake -i C C.KEYWORDS C C ECHELLE, CASPEC C C.PURPOSE C C SEARCH FOR ORDERS IN AN IMAGE. TYPE ECHELLE - FLAT FIELD C OUTPUT IS A TABLE WITH THE FOLLOWING COLUMNS : C COLUMN NO. LABEL UNIT DESCRIPTION C 1 :ORDER - ORDER NUMBER C 2 :X PIXEL SAMPLE NUMBER C 3 :Y PIXEL 1. MOMENT AS SUM(F*Y)/SUM(F) C 4 :YBKG PIXEL POSITION OF THE BACKGROUND C 5 :BKG DN BACKGROUND LEVEL C C.ALGORITHM C C MAXIMA FOLLOWING ALGORITHM IS USED C C CONSTRAINS : MAXIMUM NUMBER OF ORDERS IN THE IMAGE IS 100 C FOR EACH ORDER 19 POINTS ARE FOUND C C.INPUT/OUTPUT C C COMMAND : C SEARCH/ORDER INPUT W,T,S OUTPUT METH C input keywords are C ECHC(1:8) INPUT - IMAGE NAME C ECHC(21:28) OUTPUT - TABLE NAME C W,T,S - PARAMETERS WHERE : C ECHR(2) T THRESHOLD IN THE DETECTION C C 010226 last modif C C---------------------------------------------------------- IMPLICIT NONE INTEGER MADRID INTEGER ACTVAL,I INTEGER IXSTR INTEGER NACOL,NAXIS,NDIM,NORDER,NP,NPOINT INTEGER NTOT INTEGER STATUS,WINDOW,IACOL,ID INTEGER NPIX(3),ICOL(7),KUN,KNUL,IMNI INTEGER UPPER(100),LOWER(100),ORDPOS(100) INTEGER*8 IADD INTEGER*8 PNTRA C REAL RPAR(4),LEVEL C DOUBLE PRECISION START(3),STEP(3) C CHARACTER*8 FRMIN,TABLE CHARACTER*72 IDENA CHARACTER*80 LINE CHARACTER*64 CUNITA CHARACTER*2 METH CHARACTER*16 LABCOL(5),UNIT(5) CHARACTER*6 FORM(5) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA LABCOL(1)/'ORDER'/,FORM(1)/'I6'/ DATA LABCOL(2)/'X'/,FORM(2)/'F7.1'/ DATA LABCOL(3)/'Y'/,FORM(3)/'F7.1'/ DATA LABCOL(4)/'YBKG'/,FORM(4)/'F7.1'/ DATA LABCOL(5)/'BKG'/,FORM(5)/'E12.3'/ DATA UNIT(1)/'UNITLESS'/ DATA UNIT(2)/'PIXEL'/ DATA UNIT(3)/'PIXEL'/ DATA UNIT(4)/'PIXEL'/ DATA UNIT(5)/'DN'/ DATA NACOL/5/ C NDIM = 3 C C ... initialize system and read parameters C CALL STSPRO('ECHFOR') CALL STKRDC('ECHC',1,1,8,ACTVAL,FRMIN,KUN,KNUL,STATUS) CALL STKRDC('ECHC',1,17,4,ACTVAL,METH,KUN,KNUL,STATUS) TABLE = 'ORDER' CALL STKRDR('ECHR',1,4,ACTVAL,RPAR,KUN,KNUL,STATUS) LEVEL = RPAR(2) IF (METH.EQ.'CO') THEN LABCOL(3) = 'YBKG' LABCOL(4) = 'Y' END IF CALL STTPUT(' ECHELLE DEFINITION',STATUS) CALL STTPUT(' ------------------',STATUS) CALL STTPUT(' INPUT IMAGE : '//FRMIN,STATUS) CALL STTPUT(' OUTPUT TABLE : '//TABLE,STATUS) CALL STTPUT(' PARAMETERS ',STATUS) WRITE (LINE,9010) RPAR(2) CALL STTPUT(LINE,STATUS) 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,IMNI,STATUS) C C ... find order positions in the middle of the image C NORDER = 0 IXSTR = NPIX(1)/2 WINDOW = RPAR(1) NP = NPIX(2) - 2*WINDOW CALL FINDM1(MADRID(PNTRA),NPIX(1),NPIX(2),IXSTR,WINDOW,NP,LEVEL, + NORDER,ORDPOS,UPPER,LOWER) IF (NORDER.EQ.0) THEN CALL STTPUT('Order detection failed',STATUS) CALL STTPUT('Use a lower threshold',STATUS) GO TO 20 END IF WRITE (LINE,9030) NORDER CALL STTPUT(LINE,STATUS) C C ... allocate table C NTOT = (NORDER+1)*NPIX(1)/10 IACOL = 12 CALL TBTINI(TABLE,F_TRANS,F_O_MODE,IACOL,NTOT,ID,STATUS) DO 10 I = 1,NACOL CALL TBCINI(ID,D_R4_FORMAT,1,FORM(I), . UNIT(I),LABCOL(I),ICOL(I), + STATUS) 10 CONTINUE CALL TBCMAP(ID,0,IADD,STATUS) C C ... follow orders C CALL FOLLOW(MADRID(PNTRA),NPIX(1),NPIX(2),NORDER,ORDPOS,UPPER, + LOWER,WINDOW,LEVEL,MADRID(IADD),NTOT,NACOL,NPOINT) C C ... update number of elements C NTOT = NPOINT CALL STKWRI('ECHI',NORDER,4,1,KUN,STATUS) CALL STKWRI('ECHI',NPIX,7,2,KUN,STATUS) CALL TBIPUT(ID,NACOL,NTOT,STATUS) IDENA = 'ORDER POSITION' C CALL SXDPUT(TABLEFULL,'IDENT','C',IDENA,1,72,STATUS) C CALL DSCUPD(TABLEFULL,TABLEFULL,' ',STATUS) CALL TBSINI(ID,STATUS) CALL TBTCLO(ID,STATUS) C C ... free data C 20 CALL STSEPI C 9000 FORMAT (' ORDER WIDTH : ',F6.1) 9010 FORMAT (' THRESHOLD : ',F6.1) C 9020 FORMAT (' SLOPE : ',F6.3) 9030 FORMAT (' NUMBER OF DETECTED ORDERS : ',I4) END SUBROUTINE FOLLOW(X,N1,N2,NORDER,ORDPOS,UPPER,LOWER,WINDOW,THRES, + TAB,NTOT,NACOL,NP) C C DEFINE ORDER POSTIONS BY FOLLOWING STRUCTURES ABOVE A GIVEN THRESHOLD C IMPLICIT NONE C INPUT ARGUMENTS C C X(N1, N2) REAL*4 INPUT ARRAY C N1, N2 INTG*4 DIMENSIONS OF X C NORDER INTG*4 NUMBER OF ORDERS C ORDPOS INTG*4 Y POSITION OF THE ORDERS IN PIXELS C LOWER,UPPER INTG*4 LOWER AND UPPER ORDER POSITIONS C WINDOW INTG*4 NUMBER OF PIXELS IN THE EXCLUDING WINDOW C THRES REAL THRESHOLD C NTOT INTG*4 ALLOCATED SPACE C C OUTPUT ARGUMENTS C C TAB REAL ARRAY WITH THE LOCATED POSITIONS IN FORMAT C (NTOT,0:NACOL) 1 - ORDER NUMBER C 2 - SAMPLE C 3 - MOMENT ORDER 1 C 4 - position of back. C 5 - BACKGROUND C NP INTG NUMBER OF LOCATED POSITIONS C C------------------------------------------------------------- C C IMPLICIT NONE INTEGER N1,N2,NORDER,NTOT,NP,NACOL INTEGER NSAMP,ISTAT,IORD,IRL,IL REAL THRES REAL X(N1,N2) INTEGER ORDPOS(NORDER),LOWER(NORDER),UPPER(NORDER) INTEGER ICPOS(20),ILPOS(20),IUPOS(20),WINDOW INTEGER TINULL DOUBLE PRECISION TDNULL REAL TAB(NTOT,0:NACOL),TRNULL CHARACTER*80 LINE LOGICAL NEXT INTEGER NF,IHW,IXSTR,IBACK,ISTEP,IX,IY,I1,IU INTEGER KZ6902,NPTS,NFOUND C CALL TBMNUL(TINULL,TRNULL,TDNULL) NSAMP = 10 IXSTR = N1/2 C C ... DISPLAY RESULTS C CALL STTPUT(' ',ISTAT) WRITE (LINE,9020) CALL STTPUT(LINE,ISTAT) WRITE (LINE,9030) CALL STTPUT(LINE,ISTAT) C C ... ITERATION ON ORDERS C NP = 0 DO 40 IORD = 1,NORDER C C ... INCLUDE START + END POINTS C NF = 0 IF (IORD.NE.NORDER) IHW = (ORDPOS(IORD+1)-ORDPOS(IORD))/2 NP = NP + 1 TAB(NP,1) = IORD TAB(NP,2) = 1. TAB(NP,3) = TRNULL TAB(NP,4) = TRNULL TAB(NP,5) = TRNULL NP = NP + 1 TAB(NP,1) = IORD TAB(NP,2) = N1 TAB(NP,3) = TRNULL TAB(NP,4) = TRNULL TAB(NP,5) = TRNULL NP = NP + 1 TAB(NP,1) = IORD TAB(NP,2) = IXSTR TAB(NP,3) = ORDPOS(IORD) IBACK = ORDPOS(IORD) + IHW IF (IBACK.LT.N2-WINDOW) THEN TAB(NP,4) = IBACK TAB(NP,5) = X(IXSTR,IBACK) ELSE TAB(NP,4) = TRNULL TAB(NP,5) = TRNULL END IF C C ... LEFT/RIGHT PART OF THE IMAGE C DO 30 IRL = 1,2 IF (IRL.EQ.1) THEN ISTEP = 1 ELSE ISTEP = -1 END IF IX = IXSTR IL = LOWER(IORD) IU = UPPER(IORD) NEXT = .TRUE. 99 KZ6902 = 1 IF ( .NOT. (NEXT)) GO TO 20 IX = IX + ISTEP IY = MAX(IL-3,WINDOW) I1 = MIN(N2-WINDOW,IU+3) NPTS = I1 - IY + 1 CALL FINDM1(X,N1,N2,IX,IY,NPTS,THRES,NFOUND,ICPOS, + IUPOS,ILPOS) C C ... CHECK IF ONE ORDER FOUND C IF (NFOUND.EQ.1) THEN C C ... STORE INTO TABLE C IF (MOD(IX,NSAMP).EQ.0) THEN NF = NF + 1 NP = NP + 1 TAB(NP,1) = IORD TAB(NP,2) = IX TAB(NP,3) = ICPOS(1) IBACK = ICPOS(1) + IHW IF (IBACK.LT.N2-WINDOW) THEN TAB(NP,4) = IBACK TAB(NP,5) = X(IXSTR,IBACK) ELSE TAB(NP,4) = TRNULL TAB(NP,5) = TRNULL END IF END IF IL = ILPOS(1) IU = IUPOS(1) NEXT = IL .GT. WINDOW .AND. IU .LT. + (N2-WINDOW) .AND. IX .GT. WINDOW .AND. + IX .LT. (N1-WINDOW) ELSE NEXT = .FALSE. END IF GOTO 99 20 CONTINUE 30 CONTINUE WRITE (LINE,9010) IORD,IXSTR,ORDPOS(IORD),IHW,NF CALL STTPUT(LINE,ISTAT) 40 CONTINUE WRITE (LINE,9040) CALL STTPUT(LINE,ISTAT) RETURN C 9000 FORMAT (I4) 9010 FORMAT (1X,I7,2X,I8,2X,I8,2X,I10,2X,I8) 9020 FORMAT (' SEQ.NO. X CENTER Y CENTER INTERORDER TEMPLA') 9030 FORMAT (' ------- -------- -------- ---------- --------') 9040 FORMAT (' -------------------------------------------------') END SUBROUTINE FINDM1(X,N1,N2,IX,IY,NP,THRES,NORDER,ORDPOS,UPPER, + LOWER) C IMPLICIT NONE C C FIND CENTER OF THE ORDERS IN THE MIDDLE OF THE IMAGE C C INPUT ARGUMENTS C X(N1,N2) REAL INPUT ARRAY C N1,N2 INTG ARRAY DIMENSION C IX, IY INTG STARTING PIXEL POSITIONS C NP INTG NUMBER OF PIXELS TO SEARCH C THRES REAL THRESHOLD C C OUTPUT ARGUMENTS C NORDER INTG NUMBER OD ORDERS C ORDPOS INTG POSITION OF THE ORDER (UPPER-LOWER)/2 C UPPER INTG UPPER ORDER LIMIT C LOWER INTG LOWER ORDER LIMIT C INTEGER N1,N2,IX,IY,NP,NORDER,ISTART,NO,I,NFIRST,NPIX INTEGER ORDPOS(1),UPPER(1),LOWER(1),WIDTH,WIDTH1,STATUS C REAL THRES REAL X(N1,N2) C LOGICAL FIRST C C C ISTART = IX WIDTH = 0 C C ... FIND UPPER LEVEL C FIRST = .TRUE. NO = 0 DO 10, I = IY,IY + NP - 1 IF (X(ISTART,I).GT.THRES) THEN C C ... ABOVE THE THRESHOLD C IF (FIRST) THEN C C ... FIRST VALUE ABOVE THE THRESHOLD C FIRST = .FALSE. NFIRST = I NPIX = 1 END IF ELSE C C ... BELOW THE THRESHOLD C IF ( .NOT. FIRST) THEN C C ... FIRST VALUE BELOW THE THRESHOLD C FIRST = .TRUE. NO = NO + 1 UPPER(NO) = I - 1 LOWER(NO) = NFIRST ORDPOS(NO) = LOWER(NO) + (UPPER(NO)-LOWER(NO))/2 IF (WIDTH.EQ.0) THEN WIDTH = UPPER(NO) - LOWER(NO) ELSE WIDTH1 = UPPER(NO) - LOWER(NO) IF (ABS(WIDTH1-WIDTH).GT.0.1*WIDTH) + CALL STTPUT( + 'Warning: Order width changes',STATUS) WIDTH = WIDTH1 END IF END IF END IF 10 CONTINUE NORDER = NO RETURN END