C @(#)fitell2.for 17.1.1.1 (ES0-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 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) 1987 European Southern Observatory, C all rights reserved C.IDENTIFICATION: ELLIPS.FOR C.PURPOSE: Fit ellipses to 2-dim images,retrieve parameters relevant C to isophotal twisti C.LANGUAGE: F77+ESOext C.AUTHOR: Edwin A. Valentijn C.KEYWORDS: Ellipse fitting,isophotal twisting, display fit C.ALGORITHM: Make polar coord frames, search for points selected on radius and C intensity, least square fit to the polar coord of thee points, C retrieve ellips parameters, display points + fit C.VERSION: 821101 EAV creation C.VERSION: 871002 RHW conversion to standard interfaces C.VERSION: 871015 RHW conversion to standard fortran C --------------------------------------------------------------------- PROGRAM FITELL2 C IMPLICIT NONE C INTEGER MADRID INTEGER I,J,K,NN,IAV,ICAL,IPOL,ISTAT INTEGER MAXLEV,NAXIS,NGAL,NP,NX,NY INTEGER PNTRA,PNTRF,PNTRM,PNTRR INTEGER STAT,XSTEP,YSTEP INTEGER NA,NF,NM,NR INTEGER NPIX(2),N(10) INTEGER KUN,KNUL C DOUBLE PRECISION START(2),STEP(2) REAL IRLIM,XSTRT,YSTRT,CENTER(2),SKY REAL RAD(500),ANG(500),YV(500),X1(500),X2(500),X3(500) REAL CUTS(4),ISOP1(10),ISOP REAL CO1,CO2,CON,CC,B(10),A(10),E(10),EX(10),ANGU(10),PI,EE C CHARACTER IDENT*72,CUNIT*48,TEXT*80,CINPUT*72,INSTRM*72 CHARACTER INFR*60,RADII*8,ANGLI*8,MASK*8,CTYPE C INCLUDE 'MID_INCLUDE:ST_DEF.INC/NOLIST' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC/NOLIST' C DATA PI/3.1415926/ DATA MAXLEV/10/ DATA CTYPE/'R'/ C C *** start code C CALL STSPRO('FITELL2 ') ! set up Midas environment C C ** get keywords C CALL STKRDC('IN_A',1,1,60,IAV,INFR,KUN,KNUL,STAT) ! get frame name CALL STKRDI('INPUTI',1,1,IAV,IPOL,KUN,KNUL,STAT) ! polar frame ? CALL STKRDI('INPUTI',2,1,IAV,ICAL,KUN,KNUL,STAT) ! isophotal window CALL STKRDC('INPUTC',1,1,72,IAV,CINPUT,KUN,KNUL,STAT) ! contours CALL GETINP(CTYPE,CINPUT,MAXLEV,ISOP1,NGAL,STAT) IF (NGAL.GT.1) THEN CALL SORT1(ISOP1,NGAL) ENDIF CALL STKRDR('INPUTR',1,2,IAV,CENTER,KUN,KNUL,STAT) ! center CALL STKRDR('INPUTR',3,1,IAV,IRLIM,KUN,KNUL,STAT) ! radius around center CALL STKRDR('INPUTR',4,1,IAV,SKY,KUN,KNUL,STAT) ! sky background C C *** get the image C CALL STIGET(INFR,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 2,NAXIS,NPIX,START,STEP,IDENT, + CUNIT,PNTRF,NF,STAT) C C initialize coordinates, parameters C NX = NPIX(1) NY = NPIX(2) NP = NX*NY XSTRT = SNGL(START(1)) YSTRT = SNGL(START(2)) XSTEP = SNGL(STEP(1)) YSTEP = SNGL(STEP(2)) C C *** polar coordinates C create two images with polar coord: radii and angles C RADII = 'SCRRADII' ANGLI = 'SCRANGLI' C C *** but enable run of program with preexisting polar frames C IF (IPOL.EQ.1) THEN CALL STIGET(RADII,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 1,NAXIS,NPIX,START,STEP, + INSTRM,CUNIT,PNTRR,NR,STAT) CALL STIGET(ANGLI,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, + 1,NAXIS,NPIX,START,STEP, + INSTRM,CUNIT,PNTRA,NA,STAT) ELSE CALL STIPUT(RADII,10,1,1,NAXIS,NPIX,START,STEP, + 'POLAR-RAD',CUNIT,PNTRR,NR,STAT) CALL STIPUT(ANGLI,10,1,1,NAXIS,NPIX,START,STEP, + 'POL-ANG',CUNIT,PNTRA,NA,STAT) CALL POLAR(NP,NPIX,START,STEP,CENTER, + MADRID(PNTRR),MADRID(PNTRA)) END IF C C create image for retrieved points and ellipses C MASK = 'ELLIPS ' CUTS(1) = 0. CUTS(3) = 0. CUTS(2) = NGAL*10. CUTS(4) = NGAL*10. CALL STIPUT(MASK,10,1,1,NAXIS,NPIX,START,STEP, + 'ELLIPS',CUNIT,PNTRM,NM,STAT) CALL STDWRR(NM,'LHCUTS',CUTS,1,4,KUN,STAT) CALL STDWRC(NM,'INSTRUM',1,INSTRM,1,72,KUN,STAT) C C *** Enter main loop C DO 30 K = 1,NGAL C C *** Search for points wthin radius and isophotal window C For found points retrieve only their polar coord in RAD and ANG C NN = N(K) ISOP = ISOP1(K) CALL FITEL(NP,MADRID(PNTRF),MADRID(PNTRR), + MADRID(PNTRA),RAD,ANG,ISOP, + ICAL,SKY,IRLIM,NN) N(K) = NN ISOP1(K) = ISOP C C *** Check on the actual number of points found C IF (N(K).GT.499) THEN CALL STTPUT('*** WARNING: More then 500 points found '// + ' results unreliable',ISTAT) CALL STTPUT('*** ADVICE: make window smaller',ISTAT) N(K) = 500 END IF C C *** Continue the proces when no or few points were found C IF (N(K).LT.6) THEN A(K) = 0. B(K) = 0. E(K) = 0. EX(K) = 0. ANGU(K) = 0. GO TO 20 END IF C C *** set up for least square fit C DO 10 I = 1,N(K) YV(I) = 2./RAD(I)**2 X1(I) = SIN(2.*ANG(I)) X2(I) = COS(2.*ANG(I)) X3(I) = 1.0 10 CONTINUE NN = N(K) CALL LSS3(NN,YV,X1,X2,X3,CO1,CO2,CON) N(K) = NN C C *** recover the fit parameters ABS to avoid interuption in batch C ANGU(K) = -0.5*ATAN2(CO1,CO2) ANGU(K) = 180.*ANGU(K)/PI + 90.0 CC = SQRT(ABS(CO1*CO1+CO2*CO2)) B(K) = SQRT(ABS(2./ (CON+CC))) A(K) = SQRT(ABS(2./ (CON-CC))) E(K) = A(K)/B(K) EE = B(K)/A(K) EX(K) = SQRT(ABS(1.-EE*EE)) C C *** retrieve found datapoints and ellips in mask frame C NN = N(K) CALL PLMASK(K,RAD,ANG,NN,NP,NPIX,START,STEP,CENTER, + MADRID(PNTRM),CON,CO1,CO2) N(K) = NN 20 CONTINUE 30 CONTINUE C C *** Print table C CALL STTPUT(' ',ISTAT) CALL STTPUT('*** Summary of the Results ***',ISTAT) CALL STTPUT('Isophote #Points A B '// + 'A/B exc. angle',ISTAT) DO 40 J = 1,NGAL WRITE (TEXT,8000) ISOP1(J),N(J),A(J),B(J),E(J),EX(J),ANGU(J) CALL STTPUT(TEXT,ISTAT) 40 CONTINUE CALL STSEPI STOP ' ' C 8000 FORMAT (G9.3,2X,I8,5F10.2) C END C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.IDENTIFICATION: POLAR C.PURPOSE: Creates 2 virtual images containing the polar coordinates C with respect to CENTER using descripters of a MIDAS image C.AUTHOR: E.A.Valentijn ESO/Garching C.INPUT/OUTPUT: NP,NPIX,START,STEP standard descriptor parameters C radii go to image R C angles go to image ANGLE C.VERSION: 821028 EAV tested and works ok C.VERSION: 880707 RHW Converted to standard fortran C -------------------------------------------------------------------- SUBROUTINE POLAR(NP,NPIX,START,STEP,CENTER,R,ANGLE) C IMPLICIT NONE C COMMON /VMR/MADRID(1) C INTEGER MADRID INTEGER NP,NPIX(2),NX,NY,I,IX,IY REAL X,Y,XSTRT,YSTRT,XSTEP,YSTEP,XCT,YCT,CENTER(2) REAL R(NP),ANGLE(NP) DOUBLE PRECISION START(2),STEP(2) C XSTRT = SNGL(START(1)) YSTRT = SNGL(START(2)) XSTEP = SNGL(STEP(1)) YSTEP = SNGL(STEP(2)) XCT = CENTER(1) YCT = CENTER(2) NX = NPIX(1) NY = NPIX(2) I = 0 C DO 20 IY = 1,NY Y = YSTRT + ((IY-1)*YSTEP) DO 10 IX = 1,NX X = XSTRT + ((IX-1)*XSTEP) I = I + 1 R(I) = SQRT((X-XCT)**2+ (Y-YCT)**2) IF (R(I).EQ.0) THEN ! avoid singularity at center R(I) = 0.D-10 ANGLE(I) = 0.0D-10 ELSE ANGLE(I) = ATAN2(X-XCT,Y-YCT) END IF 10 CONTINUE 20 CONTINUE RETURN END C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.IDENTIFICATION: FITEL C.AUTHOR: Edwin A. Valentijn C.PURPOSE: Selects out of image a on int and radius and fills rad and ang C.INPUT: NP total number of pixels C A image intensities input C R polar coord image of raddi input C ANGLE polar coord image of angles input C RAD radii of found points output C ANG angles of found points output C ISOP1 isophotal level in % of background keyword C ICAL isophotal window in pro mille of background keyword C SKY background level to be added before fitting keyword C IRLIM maximum radius used for search keyword C N number of points found output C.VERSION: 880707 Converted to standard fortran C -------------------------------------------------------------------- SUBROUTINE FITEL(NP,A,R,ANGLE,RAD,ANG,ISOP1,ICAL,SKY,IRLIM,N) C IMPLICIT NONE C COMMON /VMR/MADRID(1) C INTEGER MADRID INTEGER NP,N,ICAL,NI,ISTAT REAL A(NP),R(NP),ANGLE(NP) REAL RAD(500),ANG(500),ISOP1 REAL SKY(2),RISOP,LOLIM,UPLIM,CAL1,IRLIM CHARACTER TEXT*80 C RISOP = ISOP1 CAL1 = (ICAL*ISOP1)/1000. LOLIM = SKY(1) + RISOP - CAL1 UPLIM = SKY(1) + RISOP + CAL1 C C *** search int array for filtered values C simultaneously R array for radii C if true put corresponding R and ANGLE in RAD and ANG C N = 0 DO 10 NI = 1,NP IF (R(NI).LT.IRLIM .AND. A(NI).LT.UPLIM .AND. 2 A(NI).GT.LOLIM) THEN N = N + 1 IF (N.LE.500) THEN RAD(N) = R(NI) ANG(N) = ANGLE(NI) END IF END IF 10 CONTINUE WRITE (TEXT,9000) ISOP1,ICAL,LOLIM,UPLIM,IRLIM,N CALL STTPUT(TEXT,ISTAT) C RETURN C 9000 FORMAT ('Isoph=',G9.3,' Window=',I4,' Lim=',F8.2,1X,F8.2, + ' Rad=',F8.2,' #Pnts=',I4) C END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.IDENTIFICATION: LSS3 C.PURPOSE: least square solution of a system of equations with 3 unknown: C A * X + B * Y + C * Z = CNS C.AUTHOR: Edwin A. Valentijn C.INPUT.OUTPUT: N : number of equations C CNS: array containing the constants C AAR: coef. of X C BAR: coef. of Y C CAR: coef. of Z C X : least square solution of first unknown C Y : second C Z : third C.VERSION: 880707 RHW conversion to standard fortran C ----------------------------------------------------------------- SUBROUTINE LSS3(N,CNS,AAR,BAR,CAR,X,Y,Z) C IMPLICIT NONE C COMMON /VMR/MADRID(1) C INTEGER MADRID INTEGER N,I REAL CNS(N),AAR(N),BAR(N),CAR(N),X,Y,Z DOUBLE PRECISION AA,AB,AC,AN,BB,BC,BN,CC,CN,DA,DB,DC,DN DOUBLE PRECISION FBB,FBN,FBC,FCC,FCN C IF (N.LT.3) THEN ! Check if there are more then 3 equations RETURN ENDIF C C *** initiate the sum - variables C AA = 0.0D+0 AB = 0.0D+0 AC = 0.0D+0 AN = 0.0D+0 BB = 0.0D+0 BC = 0.0D+0 BN = 0.0D+0 CC = 0.0D+0 CN = 0.0D+0 C DO 10 I = 1,N ! start the summing DA = DBLE(AAR(I)) DB = DBLE(BAR(I)) DC = DBLE(CAR(I)) DN = DBLE(CNS(I)) AA = AA + DA*DA AB = AB + DA*DB AC = AC + DA*DC AN = AN + DA*DN BB = BB + DB*DB BC = BC + DB*DC BN = BN + DB*DN CC = CC + DC*DC CN = CN + DC*DN 10 CONTINUE C C *** first step in solving C FBB = BB - AB*AB/AA FBC = BC - AB*AC/AA FCC = CC - AC*AC/AA FCN = CN - AC*AN/AA FBN = BN - AB*AN/AA C C *** get the unknown C Z = SNGL((FCN-FBC*FBN/FBB)/ (FCC-FBC*FBC/FBB)) Y = SNGL((FBN-FBC*Z)/FBB) X = SNGL((AN-AB*Y-AC*Z)/AA) C RETURN END C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.IDENTIFICATION: PLMASK C.PURPOSE: Finds rectangular image coordinates and labels that in mask C frame(A) of: C set of (N) points in polar coord given by RAD and ANG C ellips given by least square parameters CON,CO1,CO2 C.AUTHOR: Edwin A. Valentijn C.NOTE: mask does not work : instead use normal image C.VERSION: 880707 RHW conversion to standard fortran C ------------------------------------------------------------------- SUBROUTINE PLMASK(K,RAD,ANG,N,NP,NPIX,START,STEP,CENTER,A,CON, + CO1,CO2) C IMPLICIT NONE C COMMON /VMR/MADRID(1) C INTEGER MADRID INTEGER K,N,NP,NPIX(2),I,II,IX,IY REAL A(NP) REAL RAD(500),ANG(500),CENTER(2) REAL R,TH,CON,CO1,CO2,XX,YY DOUBLE PRECISION START(2),STEP(2) C C *** initialize Mask frame C IF (K.EQ.1) THEN DO 10 I = 1,NP A(I) = 0. 10 CONTINUE END IF C C *** find indices of points of the polar coord array C DO 20 I = 1,N XX = RAD(I)*SIN(ANG(I)) YY = RAD(I)*COS(ANG(I)) IX = IFIX((XX+CENTER(1)+0.5-SNGL(START(1)))/SNGL(STEP(1))) IY = IFIX((YY+CENTER(2)+0.5-SNGL(START(2)))/SNGL(STEP(2))) II = (IY-1)*NPIX(1) + IX ! put a label in mask A(II) = K*10. 20 CONTINUE C C *** now for the ellips; take 90 points duplication does not harm C data of found points from previous step are lost if they coincide C DO 30 I = 1,91 TH = (I-1)*4.0 TH = 3.1415926*TH/180. R = CON + CO1*SIN(2.*TH) + CO2*COS(2.*TH) R = SQRT(ABS(2.0/R)) XX = R*SIN(TH) YY = R*COS(TH) IX = IFIX((XX+CENTER(1)+0.5-SNGL(START(1)))/SNGL(STEP(1))) IY = IFIX((YY+CENTER(2)+0.5-SNGL(START(2)))/SNGL(STEP(2))) II = (IY-1)*NPIX(1) + IX A(II) = K*10. 30 CONTINUE RETURN END