C @(#)sgepi.for	19.1 (ESO-DMD) 02/25/03 13:31:24
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
C.IDENTIFICATION:
C  PROGRAM SGEPI               VERSION 1.0                SEP   8, 1988
C  P.GROSBOL                   ESO - GARCHING
C
C.PURPOSE:
C  PROGRAM ESTIMATES THE POSITION AND INCLINATION ANGLE OF A GALAXY.
C  THIS IS DONE BY USE THE VARIATION OF THE 2ND FOURIER HARMONIC.
C
C.INPUT/OUTPUT:
C  THE PROGRAM RUNS UNDER MIDAS AND USES THE FOLLOWING KEYWORDS:
C    IN_A       :  NAME OF IMAGE FILE                               (I)
C    INPUTR(1-5):  REAL INPUT PARAMETERS DEFINING RADIAL RANGE      (I)
C                  1:INNER RADIUS, 2:OUTER RADIUS, 3:RADIAL STEP
C                  4:X CENTER POSITION, 5:Y CENTER POSITION,
C    OUTPUTR(1-4): REAL OUTPUT PARAMETERS GIVING                    (O)
C                  1:X CENTER POSITION, 2:Y CENTER POSITION,
C                  3:POSITON ANGLE (DEG), 4:INCLINATION ANGLE (DEG)
C
C.VERSION:
C  Included in the ESO-MIDAS version                          910208 RHW
C 
C 010201		last modif
C 
C---------------------------------------------------------------------
      PROGRAM SGEPI
C
      IMPLICIT   NONE
      INTEGER    MR
      INTEGER    MF
      INTEGER    MFD
      INTEGER    MIS
      PARAMETER  (MR=100)
      PARAMETER  (MF=6)
      PARAMETER  (MFD=MF+1)
      PARAMETER  (MIS=15)
C
      INTEGER    MADRID(1)
      INTEGER    IERR, IAV
      INTEGER    NA
      INTEGER    MC
      INTEGER    I
      INTEGER    NR
      INTEGER    NPIX(3),IMNO,KUNIT(4),INULL
      INTEGER*8    PNTR
      
      REAL       PI, PA, PM
      REAL       R, RI, RO, RM, RP
      REAL       DR, DTR, D24, DI
      REAL       XC, YC, XM
      REAL       AN
      REAL       SM, SA, SB
      REAL       WM
      REAL       AIR
      REAL       A, AI
      REAL       B
      REAL       XS, YS, SX, SY
C
      REAL       AMP(MFD),PHA(MFD),RANGE(3),POSI(4),
     2           AMP2(MR),PHA2(MR),AMP4(MR),PHA4(MR),
     3           W(MR),XV(MR),YV2(MR),YV4(MR)
      DOUBLE PRECISION START(3), STEP(3)
C
      CHARACTER*60 INAME
      CHARACTER*64 UNIT
      CHARACTER*72 IDEN
      CHARACTER*80 MSG
C
      INCLUDE  'MID_INCLUDE:ST_DEF.INC'
      COMMON   /VMR/MADRID
      INCLUDE  'MID_INCLUDE:ST_DAT.INC'
C
C     INITIATE MIDAS
C
      CALL STSPRO('SGEPI')
C
C     READ PARAMETERS
C
      CALL STKRDC('IN_A',1,1,60,IAV,INAME,KUNIT,INULL,IERR)
      CALL STKRDR('INPUTR',1,3,IAV,RANGE,KUNIT,INULL,IERR)
      CALL STKRDR('INPUTR',4,2,IAV,POSI,KUNIT,INULL,IERR)
C
C     OPEN IMAGE FRAME
C
      CALL STIGET(INAME,D_R4_FORMAT,F_I_MODE,1,3,NA,NPIX,
     ,                START,STEP,IDEN,UNIT,PNTR,IMNO,IERR)
      IF (NA.NE.2 .OR. IERR.NE.0 .OR.
     .    ABS(STEP(1)/STEP(2)-1.0).GT.1E-5) THEN
         CALL STTPUT('*** FATAL: WRONG IMAGE FORMAT',IERR)
         CALL STSEPI
      ENDIF
C
C     CHECK PARAMETERS AND CONVERT TO PIXEL VALUES
C
      PI  = 4.0 * ATAN(1.0)
      DTR = PI / 180.0
      XS  = START(1)
      YS  = START(2)
      SX  = STEP(1)
      SY  = STEP(2)
      RI  = MAX(0.0, ABS(RANGE(1)))
      RO  = MAX(RI, ABS(RANGE(2)))
      RM  = MIN(ABS(POSI(1)-XS), ABS(POSI(2)-YS))
      RM  = MIN(RM, ABS(POSI(1)-XS-NPIX(1)*SX))
      RM  = MIN(RM, ABS(POSI(2)-YS-NPIX(2)*SY))
      RO  = MIN(RM-2.0*ABS(SX), RO)
      DR  = ABS(RANGE(3))
      XC  = (POSI(1)-XS)/SX + 1.0
      YC  = (POSI(2)-YS)/SY + 1.0
      NR = MIN(MR, INT((RO-RI)/DR))
C
C     GET FOURIER COEFFICIENTS FROM EACH RADIUS
C
      R  = RI
      DO 100, I = 1,NR
         RP = R / STEP(1)
         CALL APFFTC(MADRID(PNTR),NPIX(1),NPIX(2),XC,YC,RP,0.0,0.0,
     ,               MFD,MC,AMP,PHA,AN)
         AMP2(I) = 0.0
         PHA2(I) = 0.0
         AMP4(I) = 0.0
         PHA4(I) = 0.0
         W(I)    = 0.0
         IF (MC.GT.0) THEN
            AMP2(I) = AMP(3) / AMP(1)
            PHA2(I) = PHA(3) / (2.0*DTR)
            AMP4(I) = AMP(5) / AMP(1)
            PHA4(I) = PHA(5) / (4.0*DTR)
            D24     = COS( 4.0*DTR*(PHA2(I)-PHA4(I)) )
            W(I)    = MAX( 0.0 , D24*MIN( 10.0,AMP2(3)/AN) )
         ENDIF
         WRITE(MSG,605) R,AMP2(I),PHA2(I),PHA4(I),D24,AN,W(I)
  605    FORMAT(' R,A,P(2,4,6),W :',F6.1,F7.3,2F6.1,2F7.3,F7.3)
         CALL STTPUT(MSG,IERR)
         R  = R + DR
  100 CONTINUE
C
C    COMPUTE ESTIMATE OF POSITION ANGLE
C
      PM = 0.0
      RM = 0.0
      SM = 0.0
      WM = 0.0
      DO 200, I = 1,NR
         RM = RM + W(I)*I
         XM = XM + W(I)*PHA2(I)
         WM = WM + W(I)
  200 CONTINUE
      PA = XM / WM
      RM = (RM/WM - 1.0)*DR + RI
      WRITE(MSG,620) RM,PA
  620 FORMAT('RADIUS, P.A. : ',2F8.1)
      CALL STTPUT(MSG,IERR)
C
C     ESTIMATE INCLINATION ANGLE WITH THIS POSITION ANGLE
C
      RP = RM / STEP(1)
      AI = 0.0
      DI = 5.0 
      PA = DTR * PA
      DO 300, I = 1,MIS
         AIR = DTR * AI
         CALL APFFTC(MADRID(PNTR),NPIX(1),NPIX(2),XC,YC,RP,PA,AIR,
     ,               MFD,MC,AMP,PHA,AN)
         W(I)    = 0.0
         IF (MC.GT.0) THEN
            XV(I) = COS( AIR )
            YV2(I) = AMP(3) * COS( PHA(3) ) / AMP(1)
            YV4(I) = AMP(5) * COS( PHA(5) ) / AMP(1)
            W(I)  = 2.0*COS(ATAN(4.0*YV2(I)))
         ENDIF
         WRITE(MSG,630) AI,XV(I),YV2(I),YV4(I),W(I)
  630    FORMAT(' R,A2,P2,W :',F6.1,F8.3,2F8.3,F8.3)
         CALL STTPUT(MSG,IERR)
         AI = AI + DI
  300 CONTINUE
C
C     FIT LINE TO GET ESTIMATE FOR INCLINATION ANGLE
C
      CALL LFIT(XV,YV2,W,MIS,1,A,SA,B,SB,R)
      AI = ACOS( -A/B ) / DTR
      PA = PA / DTR
      WRITE(MSG,640) R,PA,AI
  640 FORMAT('RADIUS, PA, I : ',F8.1,2F7.1)
      CALL STTPUT(MSG,IERR)
      CALL STKWRR('OUTPUTR',POSI,1,2,KUNIT,IERR)
      CALL STKWRR('OUTPUTR',PA,3,1,KUNIT,IERR)
      CALL STKWRR('OUTPUTR',AI,4,1,KUNIT,IERR)
C
C     FINISHED - EXIT
C
      CALL STSEPI
      END
