C @(#)fitell3.for 17.1.1.1 (ESO-DMD) 01/25/02 17:19:20 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 ELLIPS C***************************************************************************** C* C*... Author : M.Schloetelburg based on formulas of C* R.Bender & C.Moellenhof published in C* Astron. Astrophys. 177,71 (1987) C* C*... Version : 11-December-1990 C* C*... Purpose : fit ellipses to isophotes C* C C 010201 last modif C C****************************************************************************** IMPLICIT NONE C INTEGER NPIX(2),MADRID(1) INTEGER KUN, KNUL, NX, NY INTEGER ISTAT, IAV INTEGER N, NID, NP, NL, NO, NAXIS INTEGER*8 IO,IP REAL CSTEP REAL CENTER(2),CUTS(2),RANGE(2) REAL A(500),B(500),C(500),X(500),Y(500) DOUBLE PRECISION START(2),STEP(2) LOGICAL OPT CHARACTER INIMAG*60,RESULT*60 CHARACTER IDENT*72,CUNIT*64,OPTION*1 CHARACTER*80 TEXT INCLUDE 'MID_INCLUDE:ST_DEF.INC/NOLIST' COMMON /PARA/A,B,C,X,Y COMMON /VMR/MADRID INCLUDE 'MID_INCLUDE:ST_DAT.INC/NOLIST' C 777 FORMAT('*** FATAL: invalid cuts ...') C-- C-- initialize MIDAS environment CALL STSPRO('ELLIPSES') C-- C-- get keys CALL STKRDC('P1',1,1,60,IAV,INIMAG,KUN,KNUL,ISTAT) CALL STKRDC('P2',1,1,60,IAV,RESULT,KUN,KNUL,ISTAT) CALL STKRDR('INPUTR',1,1,IAV,CSTEP,KUN,KNUL,ISTAT) CALL STKRDR('INPUTR',2,2,IAV,CUTS,KUN,KNUL,ISTAT) CALL STKRDR('INPUTR',4,2,IAV,RANGE,KUN,KNUL,ISTAT) CALL STKRDR('INPUTR',6,2,IAV,CENTER,KUN,KNUL,ISTAT) CALL STKRDC('P7',1,1,1,IAV,OPTION,KUN,KNUL,ISTAT) C-- C-- map input frame CALL STIGET(INIMAG,10,0,1,2,NAXIS,NPIX,START,STEP, 1 IDENT,CUNIT,IP,NP,ISTAT) NX = NPIX(1) NY = NPIX(2) C-- C-- get secton to exclude in fit IF (RANGE(1).LT.0..OR.RANGE(2).LT.0.) THEN RANGE(1) = RANGE(1) + 360. RANGE(2) = RANGE(2) + 360. ENDIF C-- C-- get lowest and highest isophote levels IF (CUTS(2).GT.CUTS(1)) GOTO 888 CALL STDRDR(NP,'LHCUTS',1,2,IAV,CUTS,KUN,KNUL,ISTAT) IF (CUTS(2).GT.CUTS(1)) GOTO 888 CALL STDRDR(NP,'LHCUTS',3,2,IAV,CUTS,KUN,KNUL,ISTAT) IF (CUTS(2).GT.CUTS(1)) GOTO 888 WRITE (TEXT,777) CALL STTPUT(TEXT,ISTAT) C-- C-- convert center coordinates from world to pixel frame 888 CENTER(1) = (CENTER(1) - START(1))/STEP(1) + 1. CENTER(2) = (CENTER(2) - START(2))/STEP(2) + 1. C-- C-- create output table 'RESULT'.TBL CALL TBTINI(RESULT,1,F_O_MODE,7,500,NID,ISTAT) CALL TBCINI(NID,10,1,'G15.7','counts','Z',N,ISTAT) CALL TBCINI(NID,10,1,'G15.7','pixels','A',N,ISTAT) CALL TBCINI(NID,10,1,'G15.7','pixels','B',N,ISTAT) CALL TBCINI(NID,10,1,'G15.7','degree','PA',N,ISTAT) CALL TBCINI(NID,10,1,'G15.7','pixels','X',N,ISTAT) CALL TBCINI(NID,10,1,'G15.7','pixels','Y',N,ISTAT) CALL TBCINI(NID,10,1,'G15.7','number','NOISE',N,ISTAT) C-- C-- fit contours and store parameters in 'RESULT'.TBL OPT = .false. IF (OPTION(1:1).EQ.'Y'.OR.OPTION(1:1).EQ.'y') OPT = .true. C-- C-- fit contours and store parameters in 'RESULT'.TBL CALL CFIT(MADRID(IP),NX,NY,CUTS,CENTER,CSTEP, 2 RANGE,NL,NID,OPT) C-- C-- map output frame CALL STIPUT(RESULT,10,1,1,2,NPIX,START,STEP, 1 IDENT,CUNIT,IO,NO,ISTAT) C-- C-- copy descriptors and close input frame CALL STDCOP(NP,NO,1,' ',ISTAT) CALL STFCLO(NP,ISTAT) C-- C-- reconstruct model image from fitted parameters WRITE (TEXT,100) 100 FORMAT ('*** INFO: reconstructing output image now ...') CALL STTPUT(TEXT,ISTAT) CALL FILLBDF(MADRID(IO),NX,NY,NL,CUTS,CSTEP) CALL STFCLO(NO,ISTAT) CALL STSEPI END SUBROUTINE CFIT (M,NPX,NPY,EZ,CENTER,CSTEP,RANGE,NL,NID,OPT) C***************************************************************************** C* C* fit isophotes C* C***************************************************************************** IMPLICIT NONE C INTEGER NPX, NPY REAL M(NPX,NPY) REAL EZ(2) INTEGER HG, VG REAL CENTER(2) REAL CSTEP REAL RANGE(2) INTEGER NL INTEGER NID LOGICAL OPT C INTEGER I, L, ITURN, IMAX, IW, ISTAT REAL HABS(71),VABS(71),Z(71),G(2) REAL A(500),B(500),C(500),X(500),Y(500) REAL WZ, VV, EV REAL AFIT0, BFIT0 REAL AFIT, BFIT, HFIT, VFIT, WFIT REAL HSIN, HCOS, VSIN, VCOS REAL CEP, SEP, EPS REAL T, P, Q, O, U, R REAL H, V, HB, VB REAL BW, BW0, BX REAL USIGN, PI CHARACTER TEXT*80 C COMMON /PARA/A,B,C,X,Y C DATA IW/71/ ! points per isophote DATA G/0.1,0.1/ ! quality limits of fit DATA IMAX/25/ ! maximal # iterations DATA PI/3.1415926535897932/ C 112 FORMAT('*** WARNING: Fit parameters outside quality limits !') 113 FORMAT('*** WARNING: maximum 500 isophotes !') 400 FORMAT(17X,'parameters of each contour') 401 FORMAT(17X,26('-')) 500 FORMAT(1X,'No',1X,'level ',5X,'1.axis',5X,'2.axis', + 6X,'angle',7X,'x-cen',6X,'y-cen',6x,'sigma') 600 FORMAT(1X,I3,3(1X,G10.4),1X,G8.2,3(1X,G10.4)) C EZ(2) = EZ(2) - EZ(1) ! subtract background from maximum WZ = REAL(IW) EPS = 0.0 VV = 0. DO 1 L = 1,IW ! reset array Z(L) = 0. 1 CONTINUE WRITE(TEXT,400) CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,401) CALL STTPUT(TEXT,ISTAT) WRITE(TEXT,500) CALL STTPUT(TEXT,ISTAT) AFIT = 1. ! start with circle BFIT = 1. ! approximately DO 10 L = 1,500 ! isophote loop (not more then 500) ITURN = 0 EV = EZ(2) * 10.**(-CSTEP*REAL(L)) IF (EV.LT.EZ(1)) THEN NL = L - 1 CALL TBTCLO (NID,ISTAT) RETURN ENDIF 100 AFIT0 = AFIT BFIT0 = BFIT HFIT = 0. VFIT = 0. HSIN = 0. HCOS = 0. VSIN = 0. VCOS = 0. DO 20 I = 1,IW ! angle loop T = 360.*REAL(I)/WZ - 0.1 ! polar angle of reference points P = T - EPS Q = EPS C IF (ABS(COSD(P)).GT.0.000001) THEN ! to avoid tan(0.) C Q = EPS + ATAN2D((BFIT0*TAND(P)),AFIT0) C ENDIF IF (ABS(COS(P)/180*PI).GT.0.000001) THEN ! to avoid tan(0.) Q = EPS + ATAN2((BFIT0*TAN(P/180*PI)),AFIT0)*180/PI ENDIF IF (P.LT.-90.) Q = Q - 180. IF (P.GT.+90.) Q = Q + 180. IF (P.GT.270.) Q = Q + 180. IF (Q.LT.0.) Q = Q + 360. IF (Q.GT.360.) Q = Q - 360. O = Q C P = COSD(O) C Q = SIND(O) P = COS(O/180*PI) Q = SIN(O/180*PI) U = AMAX1(Z(I),0.) BW = EZ(2) 200 U = U + 1. BW0 = BW USIGN = 1. IF (O.GT.RANGE(1).AND.O.LT.RANGE(2)) THEN USIGN = -USIGN ENDIF H = CENTER(1) + USIGN*U*P V = CENTER(2) + USIGN*U*Q HG = INT(H) VG = INT(V) IF (HG.GE.NPX .OR. HG.LT.1) THEN USIGN = -USIGN ENDIF H = CENTER(1) + USIGN*U*P V = CENTER(2) + USIGN*U*Q HG = INT(H) VG = INT(V) IF (VG.GE.NPY .OR. VG.LT.1) THEN USIGN = -USIGN ENDIF H = CENTER(1) + USIGN*U*P V = CENTER(2) + USIGN*U*Q HG = INT(H) VG = INT(V) IF (HG.GE.NPX .OR. HG.LT.1) THEN NL = L - 1 CALL TBTCLO(NID,ISTAT) RETURN ENDIF IF (VG.GE.NPY .OR. VG.LT.1) THEN NL = L - 1 CALL TBTCLO(NID,ISTAT) RETURN ENDIF HB = H - REAL(HG) VB = V - REAL(VG) BW = (1.-VB)*((1.-HB)*M(HG,VG) + HB*M(HG+1,VG)) + + VB *((1.-HB)*M(HG,VG+1) + HB*M(HG+1,VG+1)) BW = BW - EZ(1) IF (OPT) THEN USIGN = -USIGN H = CENTER(1) + USIGN*U*P V = CENTER(2) + USIGN*U*Q HG = INT(H) VG = INT(V) IF (HG.GE.NPX .OR. HG.LT.1) THEN USIGN = -USIGN ENDIF H = CENTER(1) + USIGN*U*P V = CENTER(2) + USIGN*U*Q HG = INT(H) VG = INT(V) IF (VG.GE.NPY .OR. VG.LT.1) THEN USIGN = -USIGN ENDIF H = CENTER(1) + USIGN*U*P V = CENTER(2) + USIGN*U*Q HG = INT(H) VG = INT(V) IF (HG.GE.NPX .OR. HG.LT.1) THEN NL = L - 1 CALL TBTCLO(NID,ISTAT) RETURN ENDIF IF (VG.GE.NPY .OR. VG.LT.1) THEN NL = L - 1 CALL TBTCLO(NID,ISTAT) RETURN ENDIF HB = H - REAL(HG) VB = V - REAL(VG) BX = (1.-VB)*((1.-HB)*M(HG,VG) + HB*M(HG+1,VG)) + + VB *((1.-HB)*M(HG,VG+1) + HB*M(HG+1,VG+1)) BX = BX - EZ(1) BW = AMIN1 (BW,BX) ENDIF IF (BW.GT.EV) GO TO 200 R = U - (EV - BW)/(BW0 - BW) HABS(I) = R*P + CENTER(1) VABS(I) = R*Q + CENTER(2) C HSIN = HSIN + HABS(I)*SIND(T) C HCOS = HCOS + HABS(I)*COSD(T) C VSIN = VSIN + VABS(I)*SIND(T) C VCOS = VCOS + VABS(I)*COSD(T) HSIN = HSIN + HABS(I)*SIN(T/180*PI) HCOS = HCOS + HABS(I)*COS(T/180*PI) VSIN = VSIN + VABS(I)*SIN(T/180*PI) VCOS = VCOS + VABS(I)*COS(T/180*PI) HFIT = HFIT + HABS(I)/WZ VFIT = VFIT + VABS(I)/WZ Z(I) = U - 2. - 4.*VV 20 CONTINUE C EPS = 0.5*MOD(ATAN2D((HSIN+VCOS),(HCOS-VSIN))+180.,180.) C CEP = COSD(EPS) C SEP = SIND(EPS) EPS = 0.5*180.0/PI* 2 MOD(ATAN2((HSIN+VCOS),(HCOS-VSIN))+PI,PI) CEP = COS(EPS/180*PI) SEP = SIN(EPS/180*PI) P = HCOS*CEP*CEP + (HSIN+VCOS)*CEP*SEP + VSIN*SEP*SEP Q = HCOS*SEP*SEP - (VCOS+HSIN)*SEP*CEP + VSIN*CEP*CEP P = 2.*P/WZ Q = 2.*Q/WZ AFIT = AMAX1(P,Q) BFIT = AMIN1(P,Q) IF (P.LT.Q) EPS = EPS + 90. P = AFIT / BFIT Q = AFIT0/BFIT0 ITURN = ITURN + 1 IF (ITURN.GT.IMAX) GOTO 111 IF (ABS(Q/P - 1.).GT.G(1)) GO TO 100 P = ABS(HFIT - CENTER(1)) Q = ABS(VFIT - CENTER(2)) CENTER(1) = HFIT CENTER(2) = VFIT IF (P.GT.G(2) .OR. Q.GT.G(2)) GO TO 100 111 IF (ITURN.GT.IMAX) THEN WRITE(TEXT,112) CALL STTPUT(TEXT,ISTAT) ENDIF C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C+ C+ compute standart deviation C+ C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ V = 0. C CEP = COSD(EPS) C SEP = SIND(EPS) CEP = COS(EPS/180*PI) SEP = SIN(EPS/180*PI) DO 30 I = 1,IW H = 360.*REAL(I)/WZ - EPS C P = CEP*COSD(H)*AFIT - SEP*SIND(H)*BFIT C Q = SEP*COSD(H)*AFIT + CEP*SIND(H)*BFIT P = CEP*COS(H/180*PI)*AFIT - SEP*SIN(H/180*PI)*BFIT Q = SEP*COS(H/180*PI)*AFIT + CEP*SIN(H/180*PI)*BFIT R = HABS(I) - HFIT T = VABS(I) - VFIT P = (P - R)*(P - R) Q = (Q - T)*(Q - T) V = V + (P + Q)/(WZ - 1.) 30 CONTINUE VV = SQRT(V/(AFIT*BFIT)) WFIT = EPS - 90. A(L) = AFIT B(L) = BFIT C(L) = WFIT X(L) = HFIT Y(L) = VFIT CALL TBEWRR(NID,L,1,EV,ISTAT) CALL TBEWRR(NID,L,2,AFIT,ISTAT) CALL TBEWRR(NID,L,3,BFIT,ISTAT) CALL TBEWRR(NID,L,4,WFIT,ISTAT) CALL TBEWRR(NID,L,5,HFIT,ISTAT) CALL TBEWRR(NID,L,6,VFIT,ISTAT) CALL TBEWRR(NID,L,7,VV,ISTAT) WRITE(TEXT,600) L,EV,AFIT,BFIT,WFIT,HFIT,VFIT,VV CALL STTPUT(TEXT,ISTAT) 10 CONTINUE C WRITE(TEXT,113) CALL STTPUT(TEXT,ISTAT) NL = L - 1 CALL TBTCLO(NID,ISTAT) RETURN END SUBROUTINE FILLBDF (K,NPX,NPY,NL,EZ,CSTEP) C***************************************************************************** C* C* fill output image C* C***************************************************************************** IMPLICIT NONE INTEGER NPX,NPY REAL K(NPX,NPY) INTEGER NL REAL EZ(2) REAL CSTEP REAL A(500),B(500),C(500),X(500),Y(500) C INTEGER IMIN, IMAX, JMIN, JMAX INTEGER J, I, L REAL U, S, T, V, O, P, Q, D, F, H, W REAL PI C COMMON /PARA/A,B,C,X,Y C DATA PI/3.1415926535897932/ IMIN = NINT(X(1) - A(NL) - 1.) ! there is clearly no IMAX = NINT(X(1) + A(NL) + 1.) ! reconstruction JMIN = NINT(Y(1) - A(NL) - 1.) ! outside of the JMAX = NINT(Y(1) + A(NL) + 1.) ! largest [A(NL)] isophote IMIN = MAX(IMIN,1) IMAX = MIN(IMAX,NPX) JMIN = MAX(JMIN,1) JMAX = MIN(JMAX,NPY) DO J = 1,NPY DO I = 1,NPX K(I,J) = 0. ! initialize output image ENDDO ENDDO L = NL DO J = JMIN,JMAX DO I = IMIN,IMAX U = Y(L) - REAL(J) ! shift to S = X(L) - REAL(I) ! ellipses frame C T = S*COSD(C(L)) + U*SIND(C(L)) ! rotate to C V = U*COSD(C(L)) - S*SIND(C(L)) ! ellipses frame T = S*COS(C(L)/180*PI) + U*SIN(C(L)/180*PI) ! rotate to V = U*COS(C(L)/180*PI) - S*SIN(C(L)/180*PI) ! ellipses frame O = V/A(L) ! scale to P = T/B(L) ! ellipses frame 900 Q = SQRT(O*O + P*P) ! 'radius' of ellipses D = REAL(L) ! find the next L = L + NINT(SIGN(1.,(Q-1.))) ! neigbouring L = MAX(L,1) ! ellipses, if L = MIN(L,NL) ! there is F = REAL(L) ! one. IF (ABS(D-F).GT.0.5) THEN ! switch if R is outside U = Y(L) - REAL(J) ! shift to S = X(L) - REAL(I) ! ellipses frame C T = S*COSD(C(L)) + U*SIND(C(L)) ! rotate to C V = U*COSD(C(L)) - S*SIND(C(L)) ! ellipses frame T = S*COS(C(L)/180*PI) + U*SIN(C(L)/180*PI) ! rotate to V = U*COS(C(L)/180*PI) - S*SIN(C(L)/180*PI) ! ellipses frame O = V/A(L) ! scale to P = T/B(L) ! ellipses frame H = SQRT(O*O + P*P) ! 'radius' of ellipses W = (1.-H)/(Q-H) ! interpolate IF (W.LT.0. .OR. W.GT.1.) GO TO 900 ! not the right intervall F = F + (D-F)*W K(I,J) = EZ(2)*10.**(-CSTEP*F)+EZ(1) ENDIF ENDDO ENDDO RETURN END