C @(#)intell.for	19.1 (ESO-DMD) 02/25/03 13:31:23
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

