C @(#)neczebord.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_ord.for C C.KEYWORDS C C echelle, order searching, maxima following C C.PURPOSE C C Search for echelle(ZEBRA) orders on any echelle image. C Output is a table with the following columns : C 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 Maxima following algorithm. C C.INPUT/OUTPUT C C COMMAND : C SEARCH/ORDER INPUT W OUTPUT C input keywords are C ECHC(1:8) INPUT - IMAGE NAME C ECHC(21:28) OUTPUT - TABLE NAME C W - PARAMETER WHERE : C ECHR(1) average widht of orders C C OUTPUT PARAMETRES: C C----------------------------------------------------------------- PROGRAM ZEBORD IMPLICIT NONE INTEGER MAXORD PARAMETER (MAXORD = 100) INTEGER MADRID INTEGER ACTVAL,I,ISTAT,IXSTR,COLUMN(3) INTEGER NACOL,NAXIS,NDIM,NORDER C INTEGER NPOINT INTEGER IVEC1, IVEC2, IVEC3, IVEC4, ITBL, NTBL INTEGER NTOT,PNTRA,STATUS,WINDOW,IACOL INTEGER NPIX(3),ICOL(7), KUN, KNUL, IMNO, IJ INTEGER UPPER(MAXORD),LOWER(MAXORD) REAL FORD(MAXORD), RPAR(2), CORD(MAXORD) REAL INCLIN INTEGER ISTEP, IS, ISTART, ISTOP, IDIR, TID INTEGER UP0(MAXORD),LOW0(MAXORD) REAL CORD0(MAXORD) DOUBLE PRECISION START(3),STEP(3) DOUBLE PRECISION ST CHARACTER*80 CA CHARACTER*60 TABLE CHARACTER*72 IDENA CHARACTER*80 LINE CHARACTER*64 CUNITA CHARACTER*16 LABCOL(5),UNIT(5) CHARACTER*6 FORM(5) CHARACTER*8 FRMIN CHARACTER*1 PLOT LOGICAL DBGPLT C WIDTH - HALF 0F WIDTH FOR AVERAGE INTEGER WIDTH, NTRACE, NTR, KX, KY 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/ DATA PLOT/'N'/ CA = 'COUNTS' ST = 1.D0 NDIM = 3 C C ... initialize system and read parameters C CALL STSPRO('ZEBORD') CALL STKRDC('IN_B',1,1,60,ACTVAL,TABLE,KUN,KNUL,STATUS) CALL STKRDC('IN_A',1,1,60,ACTVAL,FRMIN,KUN,KNUL,STATUS) CALL STKRDC('INPUTC',1,1,1,ACTVAL,PLOT,KUN,KNUL,STATUS) CALL STKRDI('INPUTI',1,3,ACTVAL,COLUMN,KUN,KNUL,STATUS) CALL STKRDR('INPUTR',1,2,ACTVAL,RPAR,KUN,KNUL,STATUS) IXSTR = COLUMN(1) WIDTH = COLUMN(2) NTRACE = COLUMN(3) INCLIN = RPAR(1) WINDOW = RPAR(2) + 0.5 DBGPLT = (PLOT.EQ.'Y') .OR. (PLOT.EQ.'y') C CALL STTPUT(' Echelle(ZEBRA) definition',ISTAT) CALL STTPUT(' ------------------',ISTAT) CALL STTPUT(' InputImage : '//FRMIN,ISTAT) CALL STTPUT(' OutputTable : '//TABLE,ISTAT) CALL STTPUT(' Parameters ',ISTAT) WRITE (LINE,9040) IXSTR CALL STTPUT(LINE ,ISTAT) WRITE (LINE,9050) WIDTH CALL STTPUT(LINE ,ISTAT) WRITE (LINE,9020) WINDOW CALL STTPUT(LINE,ISTAT) WRITE (LINE,9030) INCLIN CALL STTPUT(LINE,ISTAT) WRITE (LINE,9010) NTRACE CALL STTPUT(LINE,ISTAT) CALL STTPUT(' PLOT FLAG : '//PLOT,ISTAT) 9010 FORMAT ('Follow_Steps_Number: ',I6) 9020 FORMAT ('Order_Width: ',I6) 9030 FORMAT ('Order_Inclination: ',F8.4) 9040 FORMAT ('Center_Column: ',I6) 9050 FORMAT ('Averagin_Width: ',I6) 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('ZEB_ORD.WRK',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE,1, . NPIX(2)*4 + MAXORD*NTRACE*3, . ST,ST,CA,CA,IVEC1,IJ,STATUS) IVEC2 = IVEC1 + NPIX(2) IVEC3 = IVEC2 + NPIX(2) IVEC4 = IVEC3 + NPIX(2) ITBL = IVEC4 + NPIX(2) C C ... extract and plot "averaging" column CALL EXTRAV(MADRID(PNTRA),NPIX(1),NPIX(2),IXSTR,WIDTH,INCLIN, . MADRID(IVEC1)) IF ( DBGPLT ) THEN CALL PLTARR(MADRID(IVEC1),NPIX(2),100,NPIX(2)-150,FRMIN,0) ENDIF C C ... find order positions in the middle of the image by BGD-thresholding C CALL ORDTHR(MADRID(IVEC1),MADRID(IVEC2),MADRID(IVEC3), . MADRID(IVEC4),NPIX(2),WINDOW,UPPER,LOWER,FORD, . CORD,NORDER,FRMIN,DBGPLT) C CD... simulate order miss ( for ORDCHK() debuging ) CD NORDER = NORDER - 1 CD DO I = 2, NORDER CD UPPER(I) = UPPER(I+1) CD LOWER(I) = LOWER(I+1) CD CORD(I) = CORD(I+1) CD ENDDO C C ... try to find missing orders CALL ORDCHK(NORDER,UPPER,LOWER,CORD,MADRID(ITBL),MADRID(IVEC3) . ,NPIX(2),FRMIN,DBGPLT) C C ... allocate table C NTOT = (NORDER+1)*NTRACE IACOL = 12 CALL TBTINI(TABLE,F_TRANS,F_O_MODE,IACOL,NTOT,TID,STATUS) DO 10 I = 1,NACOL CALL TBCINI(TID,D_R4_FORMAT,1,FORM(I), . UNIT(I),LABCOL(I),ICOL(I),STATUS) 10 CONTINUE C C ... now we need to follow orders ISTEP = NPIX(1)/NTRACE NTBL = 1 NTR = 0 KY = MAX(0, MIN(3, (NORDER+3)/5)) DO IDIR = -1,1,2 DO I = 1,NORDER UP0(I) = UPPER(I) LOW0(I) = LOWER(I) CORD0(I) = CORD(I) ENDDO IF (IDIR.GT.0) THEN CALL REVTBL(MADRID(ITBL),NTBL) ISTART = ((IXSTR/ISTEP) + 1) * ISTEP + 1 ISTOP = NPIX(1) - ISTEP ELSE ISTART = (IXSTR/ISTEP) * ISTEP + 1 ISTOP = ISTEP + 1 ENDIF IS = ISTART - IXSTR DO I = ISTART, ISTOP, ISTEP*IDIR WRITE(*,*) '=========================================' WRITE(*,*) 'AvrColunm=',I NTR = NTR + 1 KX = MAX(0, MIN(2, NTBL/8/(KY+1)-1 )) C C23456789012345678901234567890123456789012345678901234567890123456789012345 C CALL ORDFIT(MADRID(ITBL),NTBL,KX,KY,I,UP0,LOW0,CORD0,NORDER . ,INCLIN,IS,MADRID(IVEC3),NPIX(1),NPIX(2),FRMIN,.FALSE.) IS = ISTEP*IDIR CALL EXTRAV(MADRID(PNTRA),NPIX(1),NPIX(2),I,ISTEP,INCLIN, . MADRID(IVEC1)) CD READ(*,*) III IF ( DBGPLT ) THEN CALL PLTARR(MADRID(IVEC1),NPIX(2),100,NPIX(2)-150,FRMIN,0) ENDIF CALL ORDCOR(MADRID(IVEC1),MADRID(IVEC2),MADRID(IVEC3), . MADRID(IVEC4),NPIX(2),WINDOW,UP0,LOW0,FORD, . CORD0,NORDER,FRMIN,DBGPLT) CALL ORDTBL(I,FORD,CORD0,NORDER,MADRID(ITBL),NTBL) ENDDO ENDDO CALL WRTTBL(TID,NORDER,MADRID(ITBL),NTBL,NPIX(1)) C C ... update number of elements C NTOT = NPOINT CALL STKWRI('ECHORD',NORDER,1,1,KUN,KNUL,STATUS) IDENA = 'ORDER POSITION' C CALL SXDPUT(TABLEFULL,'IDENT','C',IDENA,1,72,STATUS) C CALL DSCUPD(TABLEFULL,TABLEFULL,' ',STATUS) CD CALL TBCSRT(TID,1,1,1,STATUS) CALL TBSINI(TID,STATUS) CALL TBTCLO(TID,STATUS) C C ... that's all ... C ... free data 20 CALL STSEPI END C***************************************************************** SUBROUTINE ORDTHR(VEC1,VEC2,VEC3,IVEC4,N2,WINDOW, . UP,LOW,FORD,CORD,NORD,FRMIN,DBGPLT) C C Search for orders in the image trace (BGD-threshholding) C INTEGER MAXORD REAL SIGMA3 PARAMETER (MAXORD = 100) INTEGER N2 REAL VEC1(N2) ! input trace REAL VEC2(N2),VEC3(N2) ! working space INTEGER IVEC4(N2) ! working space INTEGER WINDOW ! order width (IN/OUT!) INTEGER UP(MAXORD) ! orders beg. (OUT) INTEGER LOW(MAXORD) ! orders end (OUT) REAL FORD(MAXORD) ! avr.flux (OUT) REAL CORD(MAXORD) ! ord. center (OUT) INTEGER NORD ! # of orders (OUT) CHARACTER*(*) FRMIN ! label LOGICAL DBGPLT ! debug plotting flag C INTEGER I, WND, ISTAT, IB, IE, K, L C REAL AM, XC, HWHI REAL S, S0, S1 REAL NOISE0, NOISE1, NOISE2 REAL MAXR, MAXL, MINC, Q CHARACTER*80 LINE C CD... culculate "autocorrelation" CD CALL CONVL(VEC1, 1, N2, 1, 1, VEC2, WINDOW*5) CD K = N2/2-WINDOW*5 CD L = WINDOW*10 CD CALL PLTARR(VEC2, N2, K, K+L, FRMIN, 1) C CD... determine order width CD WINDOW = LWORD(VEC2, N2, N2/2, AM, XC, HWHI) CD WRITE (LINE,9000) WINDOW*2 CD CALL STTPUT(LINE,ISTAT) C C ... estimate noise level NOISE0 = SIGMA3(VEC1, N2) C C ... extract "background" (variable threshold level) WND = WINDOW * 0.75 + 0.5 I = 0 CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE0, 4, I) CD CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1) NOISE1 = SIGMA3(VEC2, N2) CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE1, 4, I) CD CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1) NOISE2 = SIGMA3(VEC2, N2) CD CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE2, 4, I) IF ( DBGPLT ) THEN CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1) CD READ(*,*) III END IF WRITE (LINE,9010) NOISE0, NOISE1, NOISE2 CALL STTPUT(LINE,ISTAT) C C ... use threshold level to form orders mask IVEC4(1) = 0 IVEC4(N2) = 0 DO I = 2, N2-1 VEC2(I) = VEC1(I)-VEC2(I) IF ( VEC1(I) .GT. VEC3(I)+NOISE1) THEN IVEC4(I) = 1 ELSE IVEC4(I) = 0 ENDIF VEC3(I) = IVEC4(I) ENDDO IF ( DBGPLT ) THEN CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1) CD READ(*,*) III ENDIF C C ... try to find orders WRITE (LINE, *) '_No._Begin__End__AvrFlux_Center_______' CALL STTPUT(LINE,ISTAT) NORD = 0 DO I = 1, N2-1 VEC3(I) = 0. IF (IVEC4(I).EQ.0 .AND. IVEC4(I+1).NE.0) THEN IB = I+1 S = 0. S0 = 0. S1 = 0. ELSEIF (IVEC4(I).NE.0 .AND. IVEC4(I+1).NE.0) THEN S = S + VEC2(I) S0 = S0 + VEC1(I) S1 = S1 + VEC1(I)*I ELSEIF (IVEC4(I).NE.0 .AND. IVEC4(I+1).EQ.0) THEN IE = I IF (IE.GT.IB) THEN S1 = S1/S0 S0 = S0/(IE-IB+1) ENDIF IF (IE-IB+1 .GE. WINDOW/2) THEN Q = (IE-IB)/3.0 ! | Q | | Q | K = IB+Q+0.5 ! | | | | L = IE-Q+0.5 ! ..IB.....K.....L.....IE.. CALL MINMAX(VEC2(IB), K-IB, Q, MAXL) CALL MINMAX(VEC2(K), L-K+1, MINC, Q) CALL MINMAX(VEC2(L+1),IE-L, Q, MAXR) IF ((IE-IB+1.0 .GE. WINDOW*1.4) .AND. . (MAXR.GT.MINC .AND. MAXL.GT.MINC)) THEN IVEC4((IE+IB)/2) = 0 I = IB-2 CD WRITE (*,*) IB,K,L,IE CD WRITE (*,*) MAXL,MINC,MAXR WRITE (LINE,9040) IB, IE, S0, S1 CALL STTPUT(LINE,ISTAT) C C23456789012345678901234567890123456789012345678901234567890123456789012345 C WRITE (LINE, *) . '-- Warning! Problems with order detection! -- ' CALL STTPUT(LINE,ISTAT) ELSEIF (S .GT. WINDOW*NOISE1*2.) THEN IF (NORD.LT.MAXORD) THEN NORD = NORD + 1 ELSE GO TO 100 ENDIF UP(NORD) = IB LOW(NORD) = IE FORD(NORD) = S0 CORD(NORD) = S1 WRITE (LINE,9020) NORD, IB, IE, S0, S1 CALL STTPUT(LINE,ISTAT) VEC3( INT(S1+0.5) ) = FORD(NORD) ELSE GO TO 100 ENDIF ELSE 100 CONTINUE WRITE (*,9030) IB, IE, S0, S1 ENDIF ENDIF ENDDO WRITE(*,*) 'NumOrders=',NORD CD READ(*,*) I IF ( DBGPLT ) THEN CALL PLTARR(VEC1, N2, 100, N2-150, FRMIN, 0) CD READ(*,*) III CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1) CD READ(*,*) III ENDIF C RETURN 9000 FORMAT ('OrderWidth: ',I6) 9010 FORMAT ('NoiseLevels: ', 3F8.3) 9020 FORMAT (I4,2I6,2F8.1) 9030 FORMAT (' ---',2I6,2F8.1) 9040 FORMAT (' ***',2I6,2F8.1) END C***************************************************************** SUBROUTINE ORDCHK(NORD,UP,LOW,CORD,TBL,VEC,N2,FRMIN,DBGPLT) C C Check orders positions (to find missing orders) C INTEGER MAXORD REAL EVAPOL PARAMETER (MAXORD = 100) INTEGER NORD ! # of orders (IN/OUT) INTEGER UP(MAXORD) ! orders beg. (IN/OUT) INTEGER LOW(MAXORD) ! orders end (IN/OUT) REAL CORD(MAXORD) ! ord. center (IN/OUT) REAL TBL(3,MAXORD) ! wrk.space INTEGER N2 REAL VEC(N2) ! wrk.space CHARACTER*(*) FRMIN ! label LOGICAL DBGPLT ! debug plotting flag C INTEGER I, ISTAT, J, K, L, LM, NEW REAL COEF(30), ERR, ERR0, ERR1, ERR2, CL, C1, C2 CHARACTER*80 LINE LOGICAL OK C C ... perform polynomial fitting of order positions OK = .FALSE. 10 DO I = 1,NORD TBL(1,I) = CORD(I) TBL(2,I) = I TBL(3,I) = 0. ENDDO C K = MAX(0, MIN(3, (NORDER+3)/5)) IF (NORD .LE. 4) THEN RETURN ELSE IF (NORD .LE. 6) THEN K = 1 ELSE IF (NORD .LE. 11) THEN K = 2 ELSE K = 3 ENDIF CALL POLFIT(TBL, NORD, K, 0, COEF, ERR) ERR0 = SQRT(ERR**2/NORD) C IF ( OK ) GO TO 100 C C ... for debuging plot: CD WRITE (*, 9000) K, 0 CD WRITE (*, 9010) (COEF(I),I=1,K+1) CD WRITE (*, 9020) ERR0 CD READ (*,*) I IF ( DBGPLT ) THEN J = N2/NORD DO I = 1,N2 VEC(I) = CORD( MIN(NORD, (I - 1)/J + 1) ) ENDDO CALL PLTARR(VEC, N2, 1, N2, FRMIN, 0) DO I = 1,N2 VEC(I) = EVAPOL(COEF, K, 0, FLOAT(I-1)/J+0.5, 0.) ENDDO CALL PLTARR(VEC, N2, 1, N2, FRMIN, 1) ENDIF C C ... try to insert new orders between every pair of old orders NEW = 0 DO L = 2,NORD C ... at first - one order DO I = 1,NORD+1 IF (I .LT. L) THEN TBL(1,I) = CORD(I) ELSE IF (I .GT. L) THEN TBL(1,I) = CORD(I-1) END IF TBL(2,I) = I TBL(3,I) = 0. ENDDO TBL(1,L) = (CORD(L) + CORD(L-1)) / 2.0 CALL POLFIT(TBL, NORD+1, K, 0, COEF, ERR) ERR1 = SQRT(ERR**2/(NORD+1)) CL = EVAPOL(COEF, K, 0, FLOAT(L), 0.) C ... then - two (if posible) IF (NORD .GE. 7) THEN DO I = 1,NORD+2 IF (I .LT. L) THEN TBL(1,I) = CORD(I) ELSE IF (I .GT. L+1) THEN TBL(1,I) = CORD(I-2) END IF TBL(2,I) = I TBL(3,I) = 0. ENDDO TBL(1,L) = CORD(L-1) + (CORD(L) - CORD(L-1)) * 0.333 TBL(1,L+1) = CORD(L-1) + (CORD(L) - CORD(L-1)) * 0.667 CALL POLFIT(TBL, NORD+2, K, 0, COEF, ERR) ERR2 = SQRT(ERR**2/(NORD+2)) ELSE ERR2 = ERR1*2. END IF CD WRITE (*, 9030) L-1, L, ERR1, ERR2 C IF (ERR1.LT.ERR0 .OR. ERR2.LT.ERR0) THEN IF (ERR1 .LE. ERR2) THEN ERR0 = ERR1 NEW = 1 C1 = CL ELSE ERR0 = ERR2 NEW = 2 C1 = EVAPOL(COEF, K, 0, FLOAT(L), 0.) C2 = EVAPOL(COEF, K, 0, FLOAT(L+1), 0.) END IF LM = L END IF ENDDO C C ... if there are new orders, fix them IF (NEW .GT. 0) THEN WRITE (LINE, 9040) NEW, LM-1, LM CALL STTPUT(LINE,ISTAT) DO I = NORD+NEW, LM+NEW, -1 CORD(I) = CORD(I-NEW) UP(I) = UP(I-NEW) LOW(I) = LOW(I-NEW) ENDDO CORD(LM) = C1 L = (LOW(LM-1) - UP(LM-1) + LOW(LM) - UP(LM))/2 + 1 LOW(LM) = MIN(INT(CORD(LM) + L/2.0 + 0.5), UP(LM)-1) UP(LM) = MAX(INT(CORD(LM) - L/2.0 + 0.5), LOW(LM-1)+1) IF (NEW .GT. 1) THEN CORD(LM+1) = C2 UP(LM+1) = MAX(INT(CORD(LM+1) - L/2.0 + 0.5), LOW(LM)+1) LOW(LM+1) = MIN(INT(CORD(LM+1) + L/2.0 + 0.5), UP(LM+2)-1) END IF NORD = NORD+NEW ELSE OK = .TRUE. END IF GO TO 10 C 100 DO I = 1, NORD CORD(I) = EVAPOL(COEF, K, 0, FLOAT(I), 0.) ENDDO C RETURN CD9000 FORMAT ('Degrees: ',2I2) CD9010 FORMAT ('Coef-s: ',3F8.3) CD9020 FORMAT ('Error: ',F7.3) CD9030 FORMAT ('Ord:',I2,'-',I2,' Err: ',2F6.3) 9040 FORMAT ('-- Attention! ',I1,'-ord.miss.in:',I2,'-',I2) END C***************************************************************** SUBROUTINE ORDCOR(VEC1,VEC2,VEC3,VEC4,N2,WINDOW, . UP,LOW,FORD,CORD,NORD,FRMIN,DBGPLT) C C Correct orders positions (by correlation with model contour) C INTEGER MXBEAM, MAXORD REAL SIGMA3 PARAMETER (MAXORD = 100) PARAMETER (MXBEAM = 100) INTEGER N2 REAL VEC1(N2) ! input trace REAL VEC2(N2),VEC3(N2),VEC4(N2) ! working space INTEGER WINDOW ! order width (IN) INTEGER UP(MAXORD) ! orders beg. (IN/OUT) INTEGER LOW(MAXORD) ! orders end (IN/OUT) REAL FORD(MAXORD) ! avr.flux (OUT) REAL CORD(MAXORD) ! ord. center (IN/OUT) INTEGER NORD ! # of orders (IN) CHARACTER*(*) FRMIN ! label LOGICAL DBGPLT ! debug plotting flag C INTEGER I, WND, J, NBEAM, K, KMAX, IMAX, L INTEGER IB, IE, PIX(2) REAL S, SS, SX, SB, SMAX, CUTS(2) REAL NOISE0, NOISE1, NOISE2 REAL BEAM(MXBEAM) C DATA BEAM /MXBEAM*0./ C C ... estimate noise level NOISE0 = SIGMA3(VEC1, N2) C C ... extract and subtract "background" WND = WINDOW * 0.75 + 0.5 I = 0 CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE0, 4, I) NOISE1 = SIGMA3(VEC2, N2) CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE1, 4, I) NOISE2 = SIGMA3(VEC2, N2) C CALL BGRND(VEC1, N2, VEC2, VEC3, WND, NOISE2, 4, I) IF ( DBGPLT ) THEN CALL PLTARR(VEC3, N2, 100, N2-150, FRMIN, 1) ENDIF CD WRITE (LINE,9010) NOISE0, NOISE1, NOISE2 CD CALL STTPUT(LINE,ISTAT) DO I = 1, N2 VEC2(I) = VEC1(I)-VEC3(I) ENDDO C C ... construct sample order profile (we use gaussian) NBEAM = MIN(WINDOW,MXBEAM) BEAM(1) = 1.0 SB = 1.0 DO I = 2, NBEAM BEAM(I) = EXP(-FLOAT(I-1)**2 / (WINDOW/2.0)**2) SB = SB + 2*BEAM(I)**2 ENDDO CD WRITE(*,9030) (BEAM(I),I=1,NBEAM) C C ... for debuging plot: IF ( DBGPLT ) THEN DO I = 1,N2 VEC4(I) = 0. ENDDO END IF C C ... for every known order... DO I = 1, NORD C ... search for max.correlation with model order within UP and LOW SMAX = -1. KMAX = CORD(I) C ... try to modify UP and LOW edges if need IF (I .GT. 1) THEN IB = LOW(I-1) + 1 ELSE IB = 2 ENDIF K = UP(I) DO J = UP(I), IB, -1 IF (VEC2(J-1) .GE. VEC2(J)) GO TO 100 K = J-1 ENDDO 100 IB = K IF (I .LT. NORD) THEN IE = UP(I+1) - 1 ELSE IE = N2-1 ENDIF K = LOW(I) DO J = LOW(I), IE IF (VEC2(J+1) .GE. VEC2(J)) GO TO 200 K = J+1 ENDDO 200 IE = K SS = 0. SX = 0. DO K = IB, IE C ... calc.next val. of conv-n with BEAM S = VEC2(K)*BEAM(1) DO J = 2, NBEAM L = K - (J-1) IF (L .GE. IB) THEN S = S + VEC2(L)*BEAM(J) ELSE S = S + VEC2(IB)*BEAM(IB-L+1)*BEAM(J) ENDIF L = K + (J-1) IF (L .LE. IE) THEN S = S + VEC2(L)*BEAM(J) ELSE S = S + VEC2(IE)*BEAM(L-IE+1)*BEAM(J) ENDIF ENDDO S = S/SB SS = SS + S SX = SX + K*S IF (S .GT. SMAX) THEN SMAX = S KMAX = K ENDIF IF ( DBGPLT ) VEC4(K) = S ENDDO C ... check is it significant or not? CALL MNMX(VEC2(IB), IE-IB+1, CUTS, PIX) IMAX = PIX(2)-1 + IB CD WRITE (*,*) IB,IE,IMAX,KMAX,SMAX IF ( (SMAX .GT. 5.*NOISE1) .AND. . (IMAX .GT. IB) .AND. (IMAX .LT. IE) .AND. . (KMAX .GT. IB) .AND. (KMAX .LT. IE) ) THEN CORD(I) = SX / SS FORD(I) = SMAX ELSE FORD(I) = -1. ENDIF ENDDO IF ( DBGPLT ) THEN CALL PLTARR(VEC4, N2, 100, N2-150, FRMIN, 1) ENDIF C C ... for debuging plot IF ( DBGPLT ) THEN DO I = 1,N2 VEC4(I) = 0. ENDDO DO I = 1, NORD IF (FORD(I) .GT. 0.) THEN VEC4( INT(CORD(I)+0.5) ) = FORD(I) ENDIF ENDDO CALL PLTARR(VEC4, N2, 100, N2-150, FRMIN, 1) ENDIF C C ... as the last step correct edges location CD WRITE (LINE, *) '_No._Begin__End__AvrFlux_Center_______' CD CALL STTPUT(LINE,ISTAT) DO I = 1, NORD IF (I .GT. 1) THEN J = (CORD(I-1) + CORD(I))/2.0 + 1.5 ELSE J = 2 ENDIF IF (I .LT. NORD) THEN K = (CORD(I+1) + CORD(I))/2.0 + 0.5 ELSE K = N2-1 ENDIF L = LOW(I) - UP(I) + 1 UP(I) = MAX(INT(CORD(I)+0.5) - L/2, J) LOW(I) = MIN(INT(CORD(I)+0.5) + (L-1)/2, K) CD WRITE (LINE,9020) I, UP(I), LOW(I), FORD(I), CORD(I) CD CALL STTPUT(LINE,ISTAT) ENDDO C RETURN CD9010 FORMAT ('NoiseLevels: ', 3F8.3) CD9020 FORMAT (I4,2I6,2F8.1) CD9030 FORMAT (3F8.3) END C***************************************************************** SUBROUTINE ORDFIT(TBL,NTBL,KX,KY,XPOS,UP,LOW,CORD,NORD,INCL . ,STEP,VEC,N1,N2,FRMIN,DBGPLT) C C 2D-poly.fit of current table C eval. edges¢ers for new X-position C INTEGER MAXORD REAL EVAPOL PARAMETER (MAXORD = 100) INTEGER NTBL ! # of rows in table (IN) REAL TBL(3,NTBL) ! "in core" table (IN) INTEGER KX, KY ! polinomial degrees (IN) INTEGER XPOS ! X-posit. (IN) INTEGER UP(MAXORD) ! orders beg. (IN/OUT) INTEGER LOW(MAXORD) ! orders end (IN/OUT) REAL CORD(MAXORD) ! ord.cent. (IN/OUT) INTEGER NORD ! # of orders (IN) REAL INCL ! inclination (IN/OUT) INTEGER STEP ! step of XPOS(IN) INTEGER N1, N2 REAL VEC(N1) ! wrk.space CHARACTER*(*) FRMIN ! label LOGICAL DBGPLT ! debug plotting flag C REAL COEF(30), ERR, SI INTEGER OVR, I, J C IF (KX.GT.0 .AND. KY.GT.0) THEN CALL POLFIT(TBL, NTBL, KX, KY, COEF, ERR) ERR = SQRT(ERR**2/NTBL) WRITE(*,*) 'TblLength:',NTBL WRITE(*,*) 'FitDegrees:',KX,KY WRITE(*,*) 'FitError:',ERR IF ( DBGPLT ) THEN OVR = 0 VEC(1) = N2-150. VEC(N1) = N2-150. VEC(2) = 100. VEC(N1-1) = 100. ENDIF SI = 0. DO I = 1,NORD CORD(I) = EVAPOL(COEF, KX, KY, FLOAT(XPOS), FLOAT(I)) UP(I) = EVAPOL(COEF, KX, KY, FLOAT(XPOS), I-0.4) + 0.5 LOW(I) = EVAPOL(COEF, KX, KY, FLOAT(XPOS), I+0.4) + 0.5 SI = SI + (EVAPOL(COEF, KX, KY, XPOS+1.0, FLOAT(I)) - . EVAPOL(COEF, KX, KY, XPOS-1.0, FLOAT(I))) / 2.0 IF ( DBGPLT ) THEN DO J = 2, N1-2 VEC(J) = EVAPOL(COEF, KX, KY, FLOAT(J), FLOAT(I)) ENDDO CALL PLTARR(VEC, N1, 1, N1, FRMIN, OVR) OVR = 1 END IF ENDDO INCL = SI / NORD WRITE (*,*) 'NewInclination:',INCL ELSE WRITE (*,*) 'OldInclination:',INCL SI = INCL * STEP DO I = 1,NORD CORD(I) = CORD(I) + SI UP(I) = UP(I) + SI + 0.5 LOW(I) = LOW(I) + SI + 0.5 ENDDO END IF C RETURN END C***************************************************************** SUBROUTINE ORDTBL(XPOS,FORD,CORD,NORD,TBL,NTBL) C C Put next portion of rows to TBL-array C INTEGER MAXORD PARAMETER (MAXORD = 100) INTEGER XPOS ! X-posit.(:X) (IN) REAL FORD(MAXORD) ! avr.flux (IN) REAL CORD(MAXORD) ! ord.cent.(:Y)(IN) INTEGER NORD ! # of orders (IN) REAL TBL(3,1) ! "in core" table (IN/OUT) INTEGER NTBL ! # of rows in table (OUT) C INTEGER NC ! curr.row number INTEGER I C DATA NC/1/ C DO I = 1, NORD IF (FORD(I) .GT. 0.) THEN TBL(1,NC) = CORD(I) TBL(2,NC) = XPOS TBL(3,NC) = I NTBL = NC NC = NC + 1 ENDIF ENDDO C RETURN END C***************************************************************** SUBROUTINE REVTBL(TBL,NTBL) C C Reverse TBL-array C INTEGER NTBL ! # of rows in TBL (IN) REAL TBL(3,NTBL) ! "in core" table (IN/OUT) C REAL RVAL(3) INTEGER I, J, K C IF (NTBL .GT. 1) THEN DO I = 1, NTBL/2 K = NTBL-I+1 DO J = 1, 3 RVAL(J) = TBL(J,I) ENDDO DO J = 1, 3 TBL(J,I) = TBL(J,K) ENDDO DO J = 1, 3 TBL(J,K) = RVAL(J) ENDDO ENDDO END IF C RETURN END C***************************************************************** SUBROUTINE WRTTBL(TID,NORD,TBL,NTBL,N1) C C Form the output table sorted by :ORDER C INTEGER TID ! table id (IN) INTEGER NORD ! # of orders (IN) INTEGER NTBL ! # of rows in table (IN) REAL TBL(3,NTBL) ! "in core" table (IN) INTEGER N1 ! image X-size (IN) C INTEGER NC ! curr.row number REAL RVAL(5) INTEGER I, K, ICOL(5), ISTAT C DATA ICOL/1,2,3,4,5/ DATA NC/1/ C DO I = 1, 5 RVAL(I) = 0. ENDDO C DO I = 1, NORD RVAL(1) = I RVAL(2) = 1.0 CALL TBRWRR(TID,NC,2,ICOL,RVAL,ISTAT) NC = NC + 1 DO K = 1, NTBL IF ( INT(TBL(3,K)+0.1) .EQ. I) THEN RVAL(1) = TBL(3,K) RVAL(2) = TBL(2,K) RVAL(3) = TBL(1,K) CALL TBRWRR(TID,NC,5,ICOL,RVAL,ISTAT) NC = NC + 1 ENDIF ENDDO RVAL(1) = I RVAL(2) = N1 CALL TBRWRR(TID,NC,2,ICOL,RVAL,ISTAT) NC = NC + 1 ENDDO C RETURN END C***************************************************************** SUBROUTINE EXTRAV(X,N1,N2,I,WIDTH,INCL,OUT) C C Extract and average 2*WIDHT+1 vertical columns C from image array X(N1,N2) centered at I-column C ( with order inclination account ) C INTEGER N1 INTEGER N2 REAL X(N1, N2) ! input image INTEGER I ! position of center column INTEGER WIDTH ! +- from I-column REAL INCL ! orders inclination coefficient REAL OUT(N2) ! 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 CC***************************************************************** C SUBROUTINE CONVL(X,N1,N2,I1,I2,OUT,NSTEP) CC CC Perform "convolution" of two vertical columns: I1, I2 CC of image array X(N1,N2) CC C INTEGER N1 C INTEGER N2 C REAL X(N1, N2) C INTEGER I1, I2 C REAL OUT(N2) C INTEGER NSTEP CC C INTEGER J, K, N, N22 C REAL S, S1 CC C N22 = N2/2 C N = MIN0(NSTEP,N22-1) C S1 = 0. C DO 10 K = 1, N2 C OUT(K) = 0. C S1 = S1 + X(I1,K) C10 CONTINUE C DO 30 J = -N, N C S = 0. C DO 20 K = 1, N2 C S = S + X(I1,K)*X(I2,MOD(N2+K+J-1,N2)+1) C20 CONTINUE C OUT(N22+J) = S/S1 C30 CONTINUE C RETURN C END C CC***************************************************************** C INTEGER FUNCTION LWORD(INP,NP,CENTER,AM,XC,HWHI) CC CC Determine new order width CC C INTEGER NP ! IN arr.length C REAL INP(NP) ! IN input array C INTEGER CENTER ! IN center pixel of contour C REAL AM, XC, HWHI ! OUT output params CC C INTEGER I, J, K, WND, IC C REAL A, B, D, A11, A12, A22, B1, B2 C REAL X, Y, Y1, Y2 CC CC ... at first - rough estimation of HWHI C K = NP C DO I = CENTER+1,NP-1 C IF (INP(I) .LE. INP(CENTER)/2.) THEN C K = I C GO TO 10 C ENDIF C IF ((INP(I-1).GE.INP(I)) .AND. (INP(I+1).GE.INP(I))) THEN C K = I-2 C GO TO 10 C ENDIF C ENDDO CC C10 WND = K - CENTER C LWORD = WND CC CC ... and then try to determine parameters with more precision C AM = INP(CENTER) C IC = CENTER CC ... make sure that we stay at maximum C DO I = CENTER-WND+1, CENTER+WND C IF (INP(I) .GT. AM) THEN C AM = INP(I) C IC = I C ENDIF C ENDDO CC ... redefine window of interest C Y1 = AM C Y = (INP(IC+1) + INP(IC-1))/2. C DO I = 1, WND C Y2 = (INP(IC+I+1) + INP(IC-I-1))/2. C IF (Y1-Y .GE. Y-Y2) THEN C WND = I+1 C GO TO 20 C ENDIF C Y1 = Y C Y = Y2 C ENDDO CC ... perform very simple fit C20 A11 = 0. C A12 = 0. C A22 = WND*2. C B1 = 0. C B2 = 0. C DO I = IC-WND+1, IC+WND C AM = MAX(AM, INP(I)) C X = I - 0.5 C Y = INP(I) - INP(I-1) C A11 = A11 + X**2 C A12 = A12 + X C B1 = B1 + X*Y C B2 = B2 + Y C ENDDO C D = (A11*A22 - A12*A12) C A = (B1*A22 - B2*A12)/D C B = (B2*A11 - B1*A12)/D C XC = -B/A C HWHI = SQRT( MAX(0., -AM/A)) C WRITE (*,9000) AM, XC, HWHI*2. CC C RETURN C 9000 FORMAT ('Am=',F7.1,' Xc=',F6.1,' FWHI=',F6.3) C END C C***************************************************************** REAL FUNCTION SIGMA3(INP,NP) C C Rough noise estimation by a "3-point dispersion" algorithm C INTEGER NP REAL INP(NP) C INTEGER K REAL S C S = 0. DO 10 K = 2, NP-1 S = S + (INP(K) - (INP(K-1)+INP(K+1))/2.)**2 10 CONTINUE SIGMA3 = SQRT(S/(NP-3)) RETURN END C***************************************************************** SUBROUTINE BGRND(INP,NP,WRK,OUT,PORT,NOISE,NIT,CIT) C C Background extraction C INTEGER NP ! vectors lehgth REAL INP(NP) ! input vector REAL WRK(NP) ! working array REAL OUT(NP) ! output vector INTEGER PORT ! smoothing window = PORT*2+1 REAL NOISE ! "cutting" level INTEGER NIT ! # of steps to pass INTEGER CIT ! current step (IN/OUT!!!) C INTEGER IT, I, II, IO, N REAL S C C ... copy input values to the working area IF (CIT .EQ. 0) THEN DO 5 I = 1,NP 5 WRK(I) =INP(I) ENDIF C C ... execute NIT steps of iteration: "smoothing" + "cutting" DO 60 IT = 1,NIT CIT = CIT + 1 C C ... first - simlple smooth with a "sliding-window" S = 0. DO 10 I = 1,PORT 10 S = S + WRK(I) ! fill smoothing window N = PORT DO 40 I = 1,NP II = I + PORT IO = I - PORT-1 IF (II .LE. NP) THEN S = S + WRK(II) ! add input element to window N = N + 1 ENDIF IF (IO .GE. 1 ) THEN S = S - WRK(IO) ! del. output elem. from window N = N - 1 ENDIF OUT(I) = S / N 40 CONTINUE C C ... and then supperss strong details with noise level account DO 50 I = 1,NP IF (INP(I) .GT. OUT(I)+NOISE) THEN WRK(I) = OUT(I) ELSE WRK(I) = INP(I) ENDIF 50 CONTINUE C 60 CONTINUE C RETURN END C***************************************************************** SUBROUTINE PLTARR(OUT,N2,NB,NE,NAME,OVR) C C Plot elements 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 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 C ********************************************************************* C INSERTED POLYFIT.FOR C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 15:55 - 19 NOV 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: J.D.PONZ C C.PURPOSE C Evaluation of 2D-polynomial C C.KEYWORDS C REGRESSION C C.ALGORITHM C C C----------------------------------------------------------- C REAL FUNCTION EVAPOL(PAR, K, L, VAL1, VAL2) C INTEGER K, L ! polynomial degrees REAL PAR(30) ! arr. of (K+1)*(L+1) coef-s REAL VAL1, VAL2 ! arguments C INTEGER NA, K1, L1, IP DOUBLE PRECISION WORK(30),RESULT,DC C IP = 0 DC = 1.D0 RESULT = 0.D0 NA = (K+1)*(L+1) IF ( NA.GE.1 .AND. NA.LE.30) THEN DO L1 = 0,L IP = IP + 1 WORK(IP) = DC RESULT = RESULT + WORK(IP)*PAR(IP) IF (K .GT. 0) THEN DO K1 = 1,K IP = IP + 1 WORK(IP) = WORK(IP-1)*VAL1 RESULT = RESULT + WORK(IP)*PAR(IP) END DO END IF DC = DC*VAL2 ENDDO END IF C EVAPOL = RESULT C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 16:02 - 19 NOV 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: J.D.PONZ C C.PURPOSE C C 2D-polynomial fitting C C.KEYWORDS C C POLYNOMIAL FIT C C.ALGORITHM C C fitting algorithm using Householder transformation C C----------------------------------------------------------- C SUBROUTINE POLFIT(VAL,NROW,IDEG1,IDEG2,COEF,ERROR) IMPLICIT NONE C INTEGER NROW ! # data elements: REAL VAL(3,NROW) ! { Y, X1, X2 } INTEGER IDEG1,IDEG2 ! polynomial degrees REAL COEF(30) ! polyn-l coefficients (OUT) REAL ERROR ! fitting error (OUT) C INTEGER IPAR(13) INTEGER I, ISTAT, KUN, KMIN INTEGER NVAR, K, L, N, LL, N1 INTEGER IROW, K1 CHARACTER*20 CPAR REAL RPAR(13) C DOUBLE PRECISION XMIN,XMAX,YMIN,YMAX,V1,V2,V3 DOUBLE PRECISION G(30,30),X(30) C ERROR = -1. DO 15 I = 1, 30 X(I) = 0.D0 15 CONTINUE K = IDEG1 L = IDEG2 IF (K.EQ.0 .AND. L.EQ.0) THEN GO TO 70 ELSE N = (K+1) * (L+1) END IF IF (N.GT.29) THEN GO TO 70 END IF IF (K.EQ.0 .OR. L.EQ.0) THEN NVAR = 2 ELSE NVAR = 3 END IF C C C ... SOLVE LEAST SQUARE PROBLEM C XMIN = 1.D30 YMIN = 1.D30 XMAX = -XMIN YMAX = -YMIN LL = 0 N1 = 0 DO 50 IROW = 1,NROW V1 = VAL(1,IROW) V2 = VAL(2,IROW) V3 = VAL(3,IROW) XMIN = MIN(XMIN,V2) XMAX = MAX(XMAX,V2) YMIN = MIN(YMIN,V3) YMAX = MAX(YMAX,V3) CALL TDSET2(LL+1,V2,V3,V1,K,L,G,X,N,30) IF (LL.NE.0) THEN KMIN = MIN(LL,N+1) DO 40 K1 = 1,KMIN CALL TDHHTR(K1,LL+1,G,X,N,30) 40 CONTINUE END IF LL = MIN(LL+1,N+1) 50 CONTINUE CALL TDSOLV(G,X,N,30) ERROR = ABS(G(N+1,N+1)) C C ... COPY RESULT OF REGRESSION C CPAR(9:16) = 'ORDER ' IPAR(1) = NROW IPAR(2) = NVAR - 1 IPAR(3) = 3 ! :Y CPAR(17:20) = 'MULT' IPAR(4) = 2 ! :X IPAR(5) = 1 ! :ORDER IPAR(6) = IDEG1 IPAR(7) = IDEG2 RPAR(1) = XMIN RPAR(2) = XMAX RPAR(3) = YMIN RPAR(4) = YMAX RPAR(5) = SQRT((ERROR*ERROR)/IPAR(1)) CALL STKWRC('OUTPUTC',1,CPAR,1,20,KUN,ISTAT) CALL STKWRI('OUTPUTI',IPAR,1,13,KUN,ISTAT) CALL STKWRR('OUTPUTR',RPAR,1,5,KUN,ISTAT) CALL STKWRD('OUTPUTD',X,1,N,KUN,ISTAT) 70 DO I = 1, N COEF(I) = X(I) END DO C C ... END C RETURN END C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 18:24 - 11 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: J.D.PONZ C C.KEYWORDS TABLE, APPLICATIONS C.ENVIRONMENT MIDAS C.PURPOSE C LINEAR REGRESSION C.ALGORITHM C C HOUSEHOLDER TRANSFORMATION C REF.: CH.L.LAWSON AND R.J.HANSON C SOLVING LEAST SQUARES PROBLEMS C PRENTICE HALL INC. 1974 C C C------------------------------------------------------------------ C SUBROUTINE TDSET2(LL,X1,X2,Y,K,L,G,X,N,M) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.PURPOSE C C PREPARE DATA IN WORK SPACE FOR 1 OR 2D POLYNOMIAL C C------------------------------------------------------------------ IMPLICIT NONE INTEGER LL DOUBLE PRECISION X1 DOUBLE PRECISION X2 DOUBLE PRECISION Y INTEGER K INTEGER L INTEGER N INTEGER M DOUBLE PRECISION G(M,M) DOUBLE PRECISION X(M) C INTEGER IP,L1,K1 DOUBLE PRECISION DC IP = 0 DC = 1.D0 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 TDHHTR(K,L1,G,X,N,M) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.PURPOSE C C Apply a Householder transformation to the two rows K and L1 of G. C C------------------------------------------------------------------ IMPLICIT NONE INTEGER K INTEGER L1 ! IN : RANK OF THE MATRIX INTEGER N INTEGER M DOUBLE PRECISION G(M,M) DOUBLE PRECISION X(M) C INTEGER KH,J DOUBLE PRECISION S,H,B C C ... CONSTRUCT THE TRANSFORMATION C S = DSQRT(G(K,K)**2+G(L1,K)**2) IF (G(K,K).GE.0.D0) S = -S H = G(K,K) - S G(K,K) = S B = G(K,K)*H C C ... APPLY THE TRANSFORMATION TO THE ROWS K AND L1 OF G. C KH = N + 1 - K IF ((KH.NE.0) .AND. (B.NE.0.D0)) 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 RETURN END SUBROUTINE TDSOLV(G,X,N,M) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.PURPOSE C C C 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------------------------------------------------------------------ IMPLICIT NONE INTEGER N INTEGER M DOUBLE PRECISION G(M,M) DOUBLE PRECISION X(M) C INTEGER I,II,J DOUBLE PRECISION S C 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 RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Author: P. Ballester. Creation 17-AUG-1992 C C Find minimum and maximum in a 1D array of type REAL C SUBROUTINE MINMAX(TAB,NEL,MIN,MAX) INTEGER NEL REAL TAB(NEL),MIN,MAX INTEGER CNT MIN = TAB(1) MAX = TAB(1) IF (NEL .GT. 1) THEN DO 1000 CNT = 2,NEL IF (TAB(CNT) .GT. MAX) MAX = TAB(CNT) IF (TAB(CNT) .LT. MIN) MIN = TAB(CNT) 1000 CONTINUE ENDIF RETURN END