C @(#)sgfpi.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 SGFPI               VERSION 1.4                FEB   9, 1988
C  P.GROSBOL                   ESO - GARCHING
C
C.PURPOSE:
C  PROGRAM FITS THE POSITION AND INCLINATION ANGLE OF A GALAXY SO
C  THAT THE PHASES OF 2ND AND 4TH HARMONIC HAVE THE MINIMUM SIGMA.
C
C.INPUT/OUTPUT:
C  THE PROGRAM RUNS UNDER MIDAS AND USES THE FOLLOWING KEYWORDS:
C      IN_A      :  NAME OF IMAGE FILE                               (I)
C      POSI(4)   :  REAL INPUT PARAMETERS DEFINING GALAXY POSITION  (I/O)
C                   1:X CENTER POSITION, 2:Y CENTER POSITION,
C                   3:POSITON ANGLE (DEG), 4:INCLINATION ANGLE (DEG)
C      RANGE(3)  :  REAL INPUT PARAMETERS DEFINING RADIAL RANGE      (I)
C                   1:INNER RADIUS, 2:OUTER RADIUS, 3:RADIAL STEP
C      REGION(2) :  THE REGION OF PA AND I TO SEARCH FOR A MINIMUM   (I)
C                   1:PA REGION (DEG),  2:I REGION (DEG)
C
C.VERSION
C 
C 021031	last modif
C 
C---------------------------------------------------------------------
      PROGRAM SGFPI
C
      IMPLICIT   NONE
      INTEGER    MR
      INTEGER    MS
      INTEGER    MFC
      PARAMETER  (MR=100)
      PARAMETER  (MS=7)
      PARAMETER  (MFC=5)
C
      INTEGER    NPIX(3),IMNO,KUNIT(4),INULL
      INTEGER    MADRID(1)
      INTEGER    IERR, IAV, IP, II, I
      INTEGER    IX, IDF
      INTEGER    MC, NR, NA
      INTEGER*8  PNTR
C
      REAL       AMP(MFC),PHA(MFC),PH(MFC),
     ,           RANGE(3),POSI(4),REG(2),
     ,           RL(MR),P2(MR),P4(MR),W(MR),
     ,           A(MS,MS),V(MS,MS),S(MS,MS),DA(MS,MS),DB(MS,MS),
     ,           SA(MS,MS),SB(MS,MS),RR(MS,MS),DIA(MS)
      REAL       PIQ, PI, PAS, PA, PAR, PI2, PIH
      REAL       DTR, DR, DPA, DAI, DIF
      REAL       RI, RO, RM, R, RP, RD
      REAL       XC, YC
      REAL       AIS, AIR
      REAL       A2, A4
      REAL       AN, AI, AD
      REAL       B2, B4
      REAL       R2, R4
      REAL       BD
      REAL       SM
      REAL       SA2, SB2, SA4, SB4
      REAL       SAD, SBD, SS2, SS4
      REAL       XM, WM
C
      DOUBLE PRECISION START(3),STEP(3)
C
      LOGICAL    LFPHA
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('SGFPI')
C
C     READ PARAMETERS
C
      CALL STKRDC('IN_A',1,1,60,IAV,INAME,KUNIT,INULL,IERR)
      CALL STKRDR('INPUTR',5,3,IAV,RANGE,KUNIT,INULL,IERR)
      CALL STKRDR('INPUTR',1,4,IAV,POSI,KUNIT,INULL,IERR)
      CALL STKRDR('INPUTR',8,2,IAV,REG,KUNIT,INULL,IERR)
C
C     OPEN IMAGE FRAME
C
      CALL STIGET(INAME,D_R4_FORMAT,F_I_MODE,1,3,NA,NPIX,START,
     2            STEP,IDEN,UNIT,PNTR,IMNO,IERR)
      IF (NA.NE.2 .OR. IERR.NE.0 .OR.
     .    ABS(SNGL(STEP(1))/SNGL(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
      LFPHA = .TRUE.
      PIQ = ATAN(1.0)
      PIH = 2.0 * PIQ
      PI  = 4.0 * PIQ
      PI2 = 8.0 * PIQ
      DTR = PI / 180.0
      RI  = MAX( 0.0,ABS(RANGE(1)) )
      RO  = MAX( RI,ABS(RANGE(2)) )
      RM  = MIN( ABS(POSI(1)-SNGL(START(1))), 
     ,           ABS(POSI(2)-SNGL(START(2))) )
      RM  = MIN( RM,ABS(POSI(1)-SNGL(START(1))-NPIX(1)*SNGL(STEP(1))) )
      RM  = MIN( RM,ABS(POSI(2)-SNGL(START(2))-NPIX(2)*SNGL(STEP(2))) )
      RO  = MIN( RM-2.0*ABS(SNGL(STEP(1))),RO )
      DR  = ABS( RANGE(3) )
      XC  = (POSI(1)-SNGL(START(1)))/SNGL(STEP(1)) + 1.0
      YC  = (POSI(2)-SNGL(START(2)))/SNGL(STEP(2)) + 1.0
      NR  = MIN( MR,INT((RO-RI)/DR) )
      DPA = REG(1) / (MS-1)
      DAI = REG(2) / (MS-1)
      PAS = POSI(3) - 0.5*REG(1)
      AIS = POSI(4) - 0.5*REG(2)
      DO 60, I = 1,MS
         DIA(I) = AIS + (I-1)*DAI
   60 CONTINUE
C
C     GO THROUGH ALL POINTS IN THE MESH
C
      PA  = PAS
      DO 1000, IP = 1,MS
         AI  = AIS
         DO 1100, II = 1,MS
C
C     GET FOURIER COEFFICIENTS FROM EACH RADIUS
C
            XM = 0.0
            SM = 0.0
            WM = 0.0
            R  = RI
            PAR = DTR * PA
            AIR = DTR * AI
            DO 1200, I = 1,NR
               RP = R / SNGL(STEP(1))
               RL(I) = LOG(R)
               P2(I) = 0.0
               P4(I) = 0.0
               W(I)  = 0.0
               CALL APFFTC(MADRID(PNTR),NPIX(1),NPIX(2),XC,YC,
     ,                     RP,PAR,AIR,MFC,MC,AMP,PHA,AN)
               IF (MC.GE.MFC) THEN
                  IF (LFPHA) THEN
                     LFPHA = .FALSE.
                     PH(1) = 0.0
                     DO 1250, IX = 2, MFC
                        DIF = PI2 / FLOAT(IX-1)
                        IDF = INT((PHA(IX)-PH(IX-1))/DIF+100.5)-100
                        PH(IX) = PHA(IX) - FLOAT(IDF)*DIF
 1250                CONTINUE
                  ENDIF
                  CALL PHAMOD(PHA,PH,MFC)
                  P2(I)  = PHA(3)
                  P4(I)  = PHA(5)
                  W(I)   = 1.0
                  DIF = ATAN(TAN(2.0*(P2(I)-P4(I)))) / DTR
                  XM  = XM + DIF
                  SM  = SM + DIF*DIF
                  WM  = WM + 1.0
               ENDIF
               R  = R + DR
 1200       CONTINUE
C
C     COMPUTE MEAN AND RMS 
C
            CALL LFIT(RL,P2,W,NR,1,A2,SA2,B2,SB2,R2)
            CALL LFIT(RL,P4,W,NR,1,A4,SA4,B4,SB4,R4)
            CALL LFIT(P2,P4,W,NR,1,AD,SAD,BD,SBD,RD)
            XM = XM / WM
            SM = SQRT( MAX(0.0,SM/WM - XM*XM) )
            A(IP,II) = XM
            S(IP,II) = SM
            DA(IP,II) = AD
            DB(IP,II) = BD
            V(IP,II) = RD
            SS2 = 0.0
            SS4 = 0.0
            DO 1300, I = 1,NR
               SS2 = SS2 + (P2(I)-A2-B2*RL(I))**2
               SS4 = SS4 + (P4(I)-A4-B4*RL(I))**2
 1300       CONTINUE
            SA(IP,II) = SQRT(SS2/NR)
            SB(IP,II) = SQRT(SS4/NR)
            RR(IP,II) = R2
            AI = AI + DAI
 1100    CONTINUE
         PA = PA + DPA
 1000 CONTINUE
C
C     ESTIMATE POSITION OF MINIMUM
C
      WRITE(MSG,620) DIA
  620 FORMAT(' INCL. : ',5X,7F8.2)
      CALL STTPUT(MSG,IERR)
      PA = PAS
      CALL STTPUT('REGR. COEF. P2-P4',IERR)
      DO 2000, IP = 1,MS
         WRITE(MSG,640) PA,(V(IP,II),II=1,MS)
  640    FORMAT(' PA : ',F7.1,3X,7F8.4)
  642    FORMAT(' PA : ',F7.1,3X,7F8.2)
         CALL STTPUT(MSG,IERR)
         PA = PA + DPA
 2000 CONTINUE
      CALL STTPUT('SIGMA(DIF.)',IERR)
      PA = PAS
      DO 2100, IP = 1,MS
         WRITE(MSG,642) PA,(S(IP,II),II=1,MS)
         CALL STTPUT(MSG,IERR)
         PA = PA + DPA
 2100 CONTINUE
      CALL STTPUT('MEAN(DIF.)',IERR)
      PA = PAS
      DO 2200, IP = 1,MS
         WRITE(MSG,642) PA,(A(IP,II),II=1,MS)
         CALL STTPUT(MSG,IERR)
         PA = PA + DPA
 2200 CONTINUE
      CALL STTPUT('B(P2-P4)',IERR)
      PA = PAS
      DO 2300, IP = 1,MS
         WRITE(MSG,640) PA,(DB(IP,II),II=1,MS)
         CALL STTPUT(MSG,IERR)
         PA = PA + DPA
 2300 CONTINUE
      CALL STTPUT('SIGMA P2',IERR)
      PA = PAS
      DO 2400, IP = 1,MS
         WRITE(MSG,640) PA,(SA(IP,II),II=1,MS)
         CALL STTPUT(MSG,IERR)
         PA = PA + DPA
 2400 CONTINUE
      CALL STTPUT('SIGMA P4',IERR)
      PA = PAS
      DO 2500, IP = 1,MS
         WRITE(MSG,640) PA,(SB(IP,II),II=1,MS)
         CALL STTPUT(MSG,IERR)
         PA = PA + DPA
 2500 CONTINUE
      CALL STTPUT('REGR. COEF. R2',IERR)
      PA = PAS
      DO 2600, IP = 1,MS
         WRITE(MSG,640) PA,(RR(IP,II),II=1,MS)
         CALL STTPUT(MSG,IERR)
         PA = PA + DPA
 2600 CONTINUE
C
C     FINISHED - EXIT
C
      CALL STSEPI
      END
