C @(#)extell.for	19.1 (ES0-DMD) 02/25/03 13:30:45
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
C.IDENTIFICATION:
C  SUBROUTINE  EXTELL            VERSION 1.0             MAY  24, 1985
C  P.GROSBOL                     ESO-GARCHING
C
C.PURPOSE
C  EXTRACE A AZIMUTHAL PROFILE OF A GALAXY ALONG AN ELLIPSE 
C
C.ALGORITHM
C  POINTS A INTERPOLATED BILINEAR ALONG AN ELLIPSE. THE NO. OF POINTS
C  IS A POWER OF 2 IN ORDER TO USE A FFT. 
C
C.INPUT/OUTPUT
C  THE SUBROUTINE IS CALLED AS
C
C     CALL EXTELL( AR, NX, NY, XC, YC, R, PA, AI, NM, PROF, NP)
C
C  WHERE THE PARAMETERS ARE
C     AR(NX,NY) :  ARRAY CONTAINING THE IMAGE                      (I)
C     NX        :  NO. OF X PIXELS                                 (I)
C     NY        :  NO. OF Y PIXELS                                 (I)
C     XC        :  X CENTER OF ELLIPSE IN PIXELS                   (I)
C     YC        :  Y CENTER OF ELLIPSE IN PIXELS                   (I)
C     R         :  MAJOR AXIS OF ELLIPSE IN PIXELS                 (I)
C     PA        :  POSITION ANGLE (MATH) IN RAD.                   (I)
C     AI        :  INCLINATION ANGLE IN RAD.                       (I)
C     NM        :  MAX. NO. OF POINTS TO EXTRACT                   (I)
C     PROF(NM)  :  THE EXTRACTED AZIMUTHAL PROFILE                 (O)
C     NP        :  NO. OF POINTS EXTRACTED                         (O)
C
C----------------------------------------------------------------------
      SUBROUTINE EXTELL(AR,NX,NY,XC,YC,R,PA,AI,NM,PROF,NP)
C
      IMPLICIT  NONE
      INTEGER   NX
      INTEGER   NY
      REAL      AR(NX,NY)
      REAL      XC
      REAL      YC
      REAL      R
      REAL      PA
      REAL      AI
      INTEGER   NM	
      REAL      PROF(NM)
      INTEGER   NP
C
      INTEGER   I, IX, IY
      INTEGER   NP2
      REAL      A, A1, A2, A3, A4
      REAL      DF, CF, SF
      REAL      XS, YS
      REAL      F, FX, FY
      REAL      B
      REAL      V1, V2
C
C     CHECK IF ELLIPSE IS FULLY WITHIN IMAGE - ELSE ERROR : NP=0
C
      NP = 0
      IF (NX.LE.0 .OR. NY.LE.0) RETURN
      IF (R-1.0.GE.MIN(XC,NX-XC) .OR.
     .    R-1.0.GE.MIN(YC,NY-YC)) RETURN
C
C     COMPUTE NO. OF POINTS :  NP = 2**(NP2)
C
      NP  = MIN( NM , MAX( 8 , INT(8.0*R*ATAN(1.0)) ) ) - 1
      NP2 = INT( LOG(FLOAT(NP))/LOG(2.0) ) + 1
      NP  = 2**(NP2)
C
C     SETUP TRANSFORMATION PARAMETERS
C
      A  =  R
      B  =  R * COS( AI )
      A1 =  A * COS( PA )
      A2 = -B * SIN( PA )
      A3 =  A * SIN( PA )
      A4 =  B * COS( PA )
C
C     DO BILINEAR INTERPOLATION ALONG ELLIPSE
C
      F  =  0.0
      DF =  8.0 * ATAN(1.0) / FLOAT(NP)
      DO 100, I = 1,NP
         CF = COS( F )
         SF = SIN( F )
         XS = A1*CF + A2*SF + XC 
         YS = A3*CF + A4*SF + YC 
         IX = INT( XS )
         IY = INT( YS )
         FX = XS - FLOAT(IX)
         FY = YS - FLOAT(IY)
         V1 = AR(IX,IY) + FX * (AR(IX+1,IY) - AR(IX,IY))
         V2 = AR(IX,IY+1) + FX * (AR(IX+1,IY+1) - AR(IX,IY+1))
         PROF(I) = V1 + FY * (V2 - V1)
         F  = F + DF
  100 CONTINUE
      RETURN
      END
