C @(#)intell.for 17.1.1.1 (ESO-DMD) 01/25/02 17:19:21 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 C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C.IDENTIFICATION: INTELL C.PURPOSE: Integrate pixel intensities within ellipse in a 2-dim image C.LANGUAGE: F77+ESOext C.AUTHOR: Andris Lauberts C.KEYWORDS: Integration, ellipse C.ALGORITHM: Search for pixel centers selected within ellipse, C compare two cases with slightly different ellipses, C derive flux by linear interpolation in differential flux C.VERSION: 890717 AL creation C C 010201 last modif C C --------------------------------------------------------------------- PROGRAM INTELL C IMPLICIT NONE C INTEGER MADRID(1) INTEGER IAV INTEGER IMF INTEGER NAXIS INTEGER ISTAT INTEGER KUN,KNUL INTEGER NPIX(2),FLAG INTEGER*8 IP C REAL ELLPAR(6),FLUX REAL STA(2), STE(2) DOUBLE PRECISION STEP(2),START(2) C CHARACTER IDENT*72,CUNIT*64,TEXT*80 CHARACTER INFR*60 C INCLUDE 'MID_INCLUDE:ST_DEF.INC/NOLIST' COMMON /VMR/MADRID INCLUDE 'MID_INCLUDE:ST_DAT.INC/NOLIST' C C *** start code CALL STSPRO('INTELL') ! set up Midas environment C C ** get keywords CALL STKRDC('IN_A',1,1,60,IAV,INFR,KUN,KNUL,ISTAT) ! get inframe name CALL STKRDR('INPUTR',1,5,IAV,ELLPAR,KUN,KNUL,ISTAT) ! ellipse params CALL STKRDI('INPUTI',1,1,IAV,FLAG,KUN,KNUL,ISTAT) ! flag C C *** get the image CALL STIGET(INFR,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, 2 2,NAXIS,NPIX,START,STEP,IDENT,CUNIT,IP,IMF,ISTAT) C C ** get descriptor IF (FLAG.EQ.1) THEN CALL STDRDR(IMF,'ELLPAR',1,5,IAV,ELLPAR,KUN,KNUL,ISTAT) ENDIF C STA(1) = START(1) STA(2) = START(2) STE(1) = STEP(1) STE(2) = STEP(2) C C *** Search for pixel centers (points) within ellipse C Simultaneous search for points within slightly larger ellipse C Study differential increase in pixel area, ellipse area and flux C Linear interpolation to get flux from pixels covering same total C area (approx same region) as the nominal ellipse. C CALL INTEGR(MADRID(IP),NPIX,STA,STE,ELLPAR,FLUX) C C *** Print table CALL STTPUT(' ',ISTAT) CALL STTPUT('*** Summary of the Results ***',ISTAT) CALL STTPUT('Center Semi-diameters '// +'Posangle Flux',ISTAT) ELLPAR(6)=FLUX WRITE (TEXT,9000) ELLPAR 9000 FORMAT (5F10.2,G15.5) CALL STTPUT(TEXT,ISTAT) C C ** update ELLPAR with FLUX in 6th position CALL STDWRR(IMF,'ELLPAR',ELLPAR,1,6,KUN,ISTAT) C 999 CALL STSEPI END C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.IDENTIFICATION: INTEGR C.AUTHOR: A. Lauberts C.PURPOSE: Integrates image pixel intensities within limiting ellipse C.INPUT PIX inframe image intensities C NPIX number of pixels/axis C START start X, Y coords C STEP step size in X, Y C ELLPAR ellipse parameters: C center X,Y, semi-diams, posangle C FLUX integrated intensit output C--------------------------------------------------------------------- SUBROUTINE INTEGR (PIX,NPIX,START,STEP,ELLPAR,FLUX) IMPLICIT NONE C INTEGER I, J, ISTAT INTEGER NPIX(2),MF,MG,N REAL PIX(1) REAL START(2),STEP(2),ELLPAR(6),FLUX REAL A,A1,A2,ABF,ABG,AF,AG REAL B,B1,B2,C1,C2,C3,D1,D2,D3 REAL F,FI,G,GI REAL P,P2,PI,PI2,PIH,PN REAL R2,S2,X,Y,Y2 REAL STC1,STC2,STP1,STP2 CHARACTER TEXT*80 C C ** define two ellipse equation functions C F(X,Y) = C1*X*X + C2*X*Y + C3*Y*Y G(X,Y) = D1*X*X + D2*X*Y + D3*Y*Y C C Let ellipse equation in quadratic form be C C C1*X*X + C2*X*Y + C3*Y*Y = 1 C C Let ellipse equation in normal form be C C r*U*U + s*V*V = 1 C C with 1/r = A*A and 1/s = B*B. C C (U,V) and (X,Y) coordinate systems are transformed by relative rotation C C U = X cos P + Y sin P C V =-X sin P + Y cos P C C using this transformation and identifying the two ellipse equation forms C gives C's expressed in A, B and P: C C 2*C1 = r+s + (r-s) cos 2P C C2 = (r-s) sin 2P C 2*C3 = r+s - (r-s) cos 2P C C ** initialize parameters C PI=3.1415926 PIH=PI/2. PI2=PI*2. C A = ELLPAR(3) B = ELLPAR(4) P = ELLPAR(5) P2 = P*PI/90. C C also try slightly larger concentric ellipse G C with semi-major diameter increased 0.5 pixel C A1 = A+STEP(1)/2. B1 = B*A1/A C C areas of ellipses F and G C ABF = A*B*PI ABG = A1*B1*PI C C constants C's of ellipse F in quadratic form C A2 = 1./A**2 B2 = 1./B**2 R2 = A2+B2 S2 = A2-B2 C1 = (R2 + S2*COS(P2))/2. C2 = S2*SIN(P2) C3 = (R2 - S2*COS(P2))/2. C C constants D's of ellipse G C A2 = 1./A1**2 B2 = 1./B1**2 R2 = A2+B2 S2 = A2-B2 D1 = (R2 + S2*COS(P2))/2. D2 = S2*SIN(P2) D3 = (R2 - S2*COS(P2))/2. C C *** search pixels within specified ellipses F and G C MF = 0 MG = 0 N = 0 FI= 0. GI= 0. STC1=START(1)-ELLPAR(1) STC2=START(2)-ELLPAR(2) STP1=STEP(1) STP2=STEP(2) DO J = 1,NPIX(2) Y = STC2+(J-1)*STP2 Y2= Y**2 DO I = 1,NPIX(1) X = STC1+(I-1)*STP1 N = N+1 IF (G(X,Y).LT.1.) THEN MG = MG + 1 PN = PIX(N) GI = GI + PN IF (F(X,Y).LT.1.) THEN MF = MF + 1 FI = FI + PN ENDIF ENDIF ENDDO ENDDO AF = MF*STEP(1)*STEP(2) AG = MG*STEP(1)*STEP(2) C C interpolate flux at approx pixel area = nominal ellipse area C FLUX = FI + (GI-FI)*((ABF-AF)/(AG-AF)) C WRITE (TEXT,9000) AF,ABF,FI CALL STTPUT(TEXT,ISTAT) WRITE (TEXT,9000) AG,ABG,GI CALL STTPUT(TEXT,ISTAT) WRITE (TEXT,9001) FLUX CALL STTPUT(TEXT,ISTAT) 9000 FORMAT ('Pix area=',G12.5,', Ell area=',G12.5,', Flux=',G12.5) 9001 FORMAT ('Interpolated flux=',G12.5) END