C @(#)fitell3.for	19.1 (ESO-DMD) 02/25/03 13:31:22
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
        PROGRAM ELLIPS
C*****************************************************************************
C*
C*... Author             :  M.Schloetelburg based on formulas of
C*                          R.Bender & C.Moellenhof published in
C*                          Astron. Astrophys. 177,71 (1987)
C*
C*... Version            :  11-December-1990
C*
C*... Purpose            :  fit ellipses to isophotes
C*
C 
C 010201		last modif
C 
C******************************************************************************
 
        IMPLICIT         NONE
C
        INTEGER          NPIX(2),MADRID(1)
        INTEGER          KUN, KNUL, NX, NY
        INTEGER          ISTAT, IAV
        INTEGER          N, NID, NP, NL, NO, NAXIS
        INTEGER*8        IO,IP

        REAL             CSTEP
        REAL             CENTER(2),CUTS(2),RANGE(2)
        REAL             A(500),B(500),C(500),X(500),Y(500)

        DOUBLE PRECISION START(2),STEP(2)

        LOGICAL          OPT

        CHARACTER        INIMAG*60,RESULT*60
        CHARACTER        IDENT*72,CUNIT*64,OPTION*1
        CHARACTER*80     TEXT

        INCLUDE          'MID_INCLUDE:ST_DEF.INC/NOLIST'
        COMMON           /PARA/A,B,C,X,Y
        COMMON           /VMR/MADRID
        INCLUDE          'MID_INCLUDE:ST_DAT.INC/NOLIST'
C
  777   FORMAT('*** FATAL: invalid cuts ...')
C--
C--     initialize MIDAS environment
        CALL STSPRO('ELLIPSES')
C--
C--     get keys
        CALL STKRDC('P1',1,1,60,IAV,INIMAG,KUN,KNUL,ISTAT)
        CALL STKRDC('P2',1,1,60,IAV,RESULT,KUN,KNUL,ISTAT)
        CALL STKRDR('INPUTR',1,1,IAV,CSTEP,KUN,KNUL,ISTAT)
        CALL STKRDR('INPUTR',2,2,IAV,CUTS,KUN,KNUL,ISTAT)
        CALL STKRDR('INPUTR',4,2,IAV,RANGE,KUN,KNUL,ISTAT)
        CALL STKRDR('INPUTR',6,2,IAV,CENTER,KUN,KNUL,ISTAT)
        CALL STKRDC('P7',1,1,1,IAV,OPTION,KUN,KNUL,ISTAT)
C--
C--     map input frame
        CALL STIGET(INIMAG,10,0,1,2,NAXIS,NPIX,START,STEP,
     1              IDENT,CUNIT,IP,NP,ISTAT)
        NX = NPIX(1)
        NY = NPIX(2)
C--
C--     get secton to exclude in fit
        IF (RANGE(1).LT.0..OR.RANGE(2).LT.0.) THEN
          RANGE(1) = RANGE(1) + 360.
          RANGE(2) = RANGE(2) + 360.
        ENDIF
C--
C--     get lowest and highest isophote levels
        IF (CUTS(2).GT.CUTS(1)) GOTO 888
          CALL STDRDR(NP,'LHCUTS',1,2,IAV,CUTS,KUN,KNUL,ISTAT)
        IF (CUTS(2).GT.CUTS(1)) GOTO 888
          CALL STDRDR(NP,'LHCUTS',3,2,IAV,CUTS,KUN,KNUL,ISTAT)
        IF (CUTS(2).GT.CUTS(1)) GOTO 888
          WRITE (TEXT,777)
          CALL  STTPUT(TEXT,ISTAT)
C--
C--     convert center coordinates from world to pixel frame
  888   CENTER(1) = (CENTER(1) - START(1))/STEP(1) + 1.
        CENTER(2) = (CENTER(2) - START(2))/STEP(2) + 1.
C--
C--     create output table 'RESULT'.TBL
        CALL TBTINI(RESULT,1,F_O_MODE,7,500,NID,ISTAT)

        CALL TBCINI(NID,10,1,'G15.7','counts','Z',N,ISTAT)
        CALL TBCINI(NID,10,1,'G15.7','pixels','A',N,ISTAT)
        CALL TBCINI(NID,10,1,'G15.7','pixels','B',N,ISTAT)
        CALL TBCINI(NID,10,1,'G15.7','degree','PA',N,ISTAT)
        CALL TBCINI(NID,10,1,'G15.7','pixels','X',N,ISTAT)
        CALL TBCINI(NID,10,1,'G15.7','pixels','Y',N,ISTAT)
        CALL TBCINI(NID,10,1,'G15.7','number','NOISE',N,ISTAT)
C--
C--     fit contours and store parameters in 'RESULT'.TBL
        OPT = .false.
        IF (OPTION(1:1).EQ.'Y'.OR.OPTION(1:1).EQ.'y') OPT = .true.
C--
C--     fit contours and store parameters in 'RESULT'.TBL
        CALL CFIT(MADRID(IP),NX,NY,CUTS,CENTER,CSTEP,
     2            RANGE,NL,NID,OPT)
C--
C--     map output frame
        CALL STIPUT(RESULT,10,1,1,2,NPIX,START,STEP,
     1              IDENT,CUNIT,IO,NO,ISTAT)
C--
C--     copy descriptors and close input frame
        CALL STDCOP(NP,NO,1,' ',ISTAT)
        CALL STFCLO(NP,ISTAT)
C--
C--     reconstruct model image from fitted parameters
        WRITE (TEXT,100)
 100    FORMAT ('*** INFO: reconstructing output image now ...')
        CALL STTPUT(TEXT,ISTAT)
        CALL FILLBDF(MADRID(IO),NX,NY,NL,CUTS,CSTEP)
        CALL STFCLO(NO,ISTAT)
        CALL STSEPI
        END



        SUBROUTINE CFIT (M,NPX,NPY,EZ,CENTER,CSTEP,RANGE,NL,NID,OPT)
C*****************************************************************************
C*
C*                            fit isophotes
C*
C*****************************************************************************
        IMPLICIT      NONE
C
        INTEGER       NPX, NPY
        REAL          M(NPX,NPY)
        REAL          EZ(2)
        INTEGER       HG, VG
        REAL          CENTER(2)
        REAL          CSTEP
        REAL          RANGE(2)
        INTEGER       NL
        INTEGER       NID
        LOGICAL       OPT
C
        INTEGER       I, L, ITURN, IMAX, IW, ISTAT
        REAL          HABS(71),VABS(71),Z(71),G(2)
        REAL          A(500),B(500),C(500),X(500),Y(500)
        REAL          WZ, VV, EV
        REAL          AFIT0, BFIT0 
        REAL          AFIT, BFIT, HFIT, VFIT, WFIT
        REAL          HSIN, HCOS, VSIN, VCOS
        REAL          CEP, SEP, EPS
        REAL          T, P, Q, O, U, R
        REAL          H, V, HB, VB
        REAL          BW, BW0, BX
        REAL          USIGN, PI
        CHARACTER     TEXT*80
C
        COMMON        /PARA/A,B,C,X,Y
C
        DATA IW/71/               ! points per isophote
        DATA G/0.1,0.1/           ! quality limits of fit
        DATA IMAX/25/             ! maximal # iterations
        DATA PI/3.1415926535897932/
C
  112   FORMAT('*** WARNING: Fit parameters outside quality limits !')
  113   FORMAT('*** WARNING: maximum 500 isophotes !')
  400   FORMAT(17X,'parameters of each contour')
  401   FORMAT(17X,26('-'))
  500   FORMAT(1X,'No',1X,'level ',5X,'1.axis',5X,'2.axis',
     +         6X,'angle',7X,'x-cen',6X,'y-cen',6x,'sigma')
  600   FORMAT(1X,I3,3(1X,G10.4),1X,G8.2,3(1X,G10.4))
C
        EZ(2) = EZ(2) - EZ(1)     ! subtract background from maximum
        WZ    = REAL(IW)
        EPS   = 0.0
        VV    = 0.

        DO 1 L = 1,IW             ! reset array
           Z(L) = 0.
   1    CONTINUE
  
        WRITE(TEXT,400)
        CALL STTPUT(TEXT,ISTAT)
        WRITE(TEXT,401)
        CALL STTPUT(TEXT,ISTAT)
        WRITE(TEXT,500)
        CALL STTPUT(TEXT,ISTAT)

        AFIT = 1.                 ! start with circle
        BFIT = 1.                 ! approximately
 
        DO 10 L = 1,500           ! isophote loop (not more then 500)
          ITURN = 0
          EV    = EZ(2) * 10.**(-CSTEP*REAL(L))

          IF (EV.LT.EZ(1)) THEN
             NL = L - 1
             CALL TBTCLO (NID,ISTAT)
             RETURN
          ENDIF

 100      AFIT0 = AFIT
          BFIT0 = BFIT
          HFIT  = 0.
          VFIT  = 0.
          HSIN  = 0.
          HCOS  = 0.
          VSIN  = 0.
          VCOS  = 0.

          DO 20 I = 1,IW          ! angle loop
             T = 360.*REAL(I)/WZ - 0.1      ! polar angle of reference points
             P = T - EPS
             Q = EPS
 
C            IF (ABS(COSD(P)).GT.0.000001) THEN      ! to avoid tan(0.)
C               Q = EPS + ATAN2D((BFIT0*TAND(P)),AFIT0)
C            ENDIF
             IF (ABS(COS(P)/180*PI).GT.0.000001) THEN      ! to avoid tan(0.)
                Q = EPS + ATAN2((BFIT0*TAN(P/180*PI)),AFIT0)*180/PI
             ENDIF

             IF (P.LT.-90.) Q = Q - 180.
             IF (P.GT.+90.) Q = Q + 180.
             IF (P.GT.270.) Q = Q + 180.
             IF (Q.LT.0.)   Q = Q + 360.
             IF (Q.GT.360.) Q = Q - 360.

             O = Q
C            P  = COSD(O)
C            Q  = SIND(O)
             P  = COS(O/180*PI)
             Q  = SIN(O/180*PI)
             U  = AMAX1(Z(I),0.) 
             BW = EZ(2)
 
 200         U     = U + 1.
             BW0   = BW
             USIGN = 1.

             IF (O.GT.RANGE(1).AND.O.LT.RANGE(2)) THEN
                USIGN = -USIGN
             ENDIF

             H  = CENTER(1) + USIGN*U*P
             V  = CENTER(2) + USIGN*U*Q
             HG = INT(H)
             VG = INT(V)

             IF (HG.GE.NPX .OR. HG.LT.1) THEN
                USIGN = -USIGN
             ENDIF

             H  = CENTER(1) + USIGN*U*P
             V  = CENTER(2) + USIGN*U*Q
             HG = INT(H)
             VG = INT(V)

             IF (VG.GE.NPY .OR. VG.LT.1) THEN
                USIGN = -USIGN
             ENDIF

             H  = CENTER(1) + USIGN*U*P
             V  = CENTER(2) + USIGN*U*Q
             HG = INT(H)
             VG = INT(V)

             IF (HG.GE.NPX .OR. HG.LT.1) THEN
                NL = L - 1
                CALL TBTCLO(NID,ISTAT)
                RETURN
             ENDIF
             IF (VG.GE.NPY .OR. VG.LT.1) THEN
                NL = L - 1
                CALL TBTCLO(NID,ISTAT)
                RETURN
             ENDIF

             HB = H - REAL(HG)
             VB = V - REAL(VG)
             BW = (1.-VB)*((1.-HB)*M(HG,VG) + HB*M(HG+1,VG)) +
     +            VB *((1.-HB)*M(HG,VG+1) + HB*M(HG+1,VG+1))
             BW = BW - EZ(1)

             IF (OPT) THEN
                USIGN = -USIGN
                H     = CENTER(1) + USIGN*U*P
                V     = CENTER(2) + USIGN*U*Q
                HG    = INT(H)
                VG    = INT(V)

                IF (HG.GE.NPX .OR. HG.LT.1) THEN
                   USIGN = -USIGN
                ENDIF

                H = CENTER(1) + USIGN*U*P
                V = CENTER(2) + USIGN*U*Q
                HG = INT(H)
                VG = INT(V)
 
                IF (VG.GE.NPY .OR. VG.LT.1) THEN
                   USIGN = -USIGN
                ENDIF

                H  = CENTER(1) + USIGN*U*P
                V  = CENTER(2) + USIGN*U*Q
                HG = INT(H)
                VG = INT(V)

                IF (HG.GE.NPX .OR. HG.LT.1) THEN      
                   NL = L - 1
                   CALL TBTCLO(NID,ISTAT)
                   RETURN
                ENDIF
                IF (VG.GE.NPY .OR. VG.LT.1) THEN      
                   NL = L - 1
                   CALL TBTCLO(NID,ISTAT)
                  RETURN
                ENDIF

                HB = H - REAL(HG)
                VB = V - REAL(VG)
                BX = (1.-VB)*((1.-HB)*M(HG,VG) + HB*M(HG+1,VG)) +
     +               VB *((1.-HB)*M(HG,VG+1) + HB*M(HG+1,VG+1))
                BX = BX - EZ(1)
                BW = AMIN1 (BW,BX)
             ENDIF

             IF (BW.GT.EV) GO TO 200

             R       = U - (EV - BW)/(BW0 - BW)
             HABS(I) = R*P + CENTER(1)
             VABS(I) = R*Q + CENTER(2)
C            HSIN    = HSIN + HABS(I)*SIND(T)
C            HCOS    = HCOS + HABS(I)*COSD(T)
C            VSIN    = VSIN + VABS(I)*SIND(T)
C            VCOS    = VCOS + VABS(I)*COSD(T)
             HSIN    = HSIN + HABS(I)*SIN(T/180*PI)
             HCOS    = HCOS + HABS(I)*COS(T/180*PI)
             VSIN    = VSIN + VABS(I)*SIN(T/180*PI)
             VCOS    = VCOS + VABS(I)*COS(T/180*PI)
             HFIT    = HFIT + HABS(I)/WZ
             VFIT    = VFIT + VABS(I)/WZ
             Z(I)    = U - 2. - 4.*VV
   20     CONTINUE     
          
C         EPS  = 0.5*MOD(ATAN2D((HSIN+VCOS),(HCOS-VSIN))+180.,180.)
C         CEP  = COSD(EPS)
C         SEP  = SIND(EPS)
          EPS  = 0.5*180.0/PI*
     2               MOD(ATAN2((HSIN+VCOS),(HCOS-VSIN))+PI,PI)
          CEP  = COS(EPS/180*PI)
          SEP  = SIN(EPS/180*PI)
          P    = HCOS*CEP*CEP + (HSIN+VCOS)*CEP*SEP + VSIN*SEP*SEP
          Q    = HCOS*SEP*SEP - (VCOS+HSIN)*SEP*CEP + VSIN*CEP*CEP   
          P    = 2.*P/WZ
          Q    = 2.*Q/WZ
          AFIT = AMAX1(P,Q)
          BFIT = AMIN1(P,Q)

          IF (P.LT.Q) EPS = EPS + 90.

          P     = AFIT / BFIT
          Q     = AFIT0/BFIT0
          ITURN = ITURN + 1

          IF (ITURN.GT.IMAX) GOTO 111

          IF (ABS(Q/P - 1.).GT.G(1))  GO TO 100

          P         = ABS(HFIT - CENTER(1))
          Q         = ABS(VFIT - CENTER(2))
          CENTER(1) = HFIT
          CENTER(2) = VFIT

          IF (P.GT.G(2) .OR. Q.GT.G(2)) GO TO 100

  111     IF (ITURN.GT.IMAX) THEN
             WRITE(TEXT,112)
             CALL STTPUT(TEXT,ISTAT)
          ENDIF
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C+
C+                    compute standart deviation
C+
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
           V = 0.
C          CEP = COSD(EPS)
C          SEP = SIND(EPS)
           CEP = COS(EPS/180*PI)
           SEP = SIN(EPS/180*PI)
 
           DO 30 I = 1,IW
              H = 360.*REAL(I)/WZ - EPS
C             P = CEP*COSD(H)*AFIT - SEP*SIND(H)*BFIT
C             Q = SEP*COSD(H)*AFIT + CEP*SIND(H)*BFIT
              P = CEP*COS(H/180*PI)*AFIT - SEP*SIN(H/180*PI)*BFIT
              Q = SEP*COS(H/180*PI)*AFIT + CEP*SIN(H/180*PI)*BFIT
              R = HABS(I) - HFIT 
              T = VABS(I) - VFIT
              P = (P - R)*(P - R)
              Q = (Q - T)*(Q - T)
              V = V + (P + Q)/(WZ - 1.)
   30      CONTINUE

           VV   = SQRT(V/(AFIT*BFIT))
           WFIT = EPS - 90.
           A(L) = AFIT
           B(L) = BFIT
           C(L) = WFIT
           X(L) = HFIT
           Y(L) = VFIT

           CALL TBEWRR(NID,L,1,EV,ISTAT)
           CALL TBEWRR(NID,L,2,AFIT,ISTAT)
           CALL TBEWRR(NID,L,3,BFIT,ISTAT)
           CALL TBEWRR(NID,L,4,WFIT,ISTAT)
           CALL TBEWRR(NID,L,5,HFIT,ISTAT)
           CALL TBEWRR(NID,L,6,VFIT,ISTAT)
           CALL TBEWRR(NID,L,7,VV,ISTAT)
           WRITE(TEXT,600) L,EV,AFIT,BFIT,WFIT,HFIT,VFIT,VV
           CALL STTPUT(TEXT,ISTAT)
   10   CONTINUE
C
        WRITE(TEXT,113)
        CALL STTPUT(TEXT,ISTAT)
        NL = L - 1
        CALL TBTCLO(NID,ISTAT)
        RETURN
        END



        SUBROUTINE FILLBDF (K,NPX,NPY,NL,EZ,CSTEP)
C*****************************************************************************
C*
C*                  fill output image
C*
C*****************************************************************************
        IMPLICIT      NONE
        INTEGER       NPX,NPY
        REAL          K(NPX,NPY)
        INTEGER       NL
        REAL          EZ(2)
        REAL          CSTEP
        REAL          A(500),B(500),C(500),X(500),Y(500)
C
        INTEGER       IMIN, IMAX, JMIN, JMAX
        INTEGER       J, I, L
        REAL          U, S, T, V, O, P, Q, D, F, H, W
        REAL          PI
C
        COMMON        /PARA/A,B,C,X,Y
C
        DATA PI/3.1415926535897932/

        IMIN = NINT(X(1) - A(NL) - 1.)      ! there is clearly no
        IMAX = NINT(X(1) + A(NL) + 1.)      ! reconstruction
        JMIN = NINT(Y(1) - A(NL) - 1.)      ! outside of the
        JMAX = NINT(Y(1) + A(NL) + 1.)      ! largest [A(NL)] isophote
        IMIN = MAX(IMIN,1)
        IMAX = MIN(IMAX,NPX)
        JMIN = MAX(JMIN,1)
        JMAX = MIN(JMAX,NPY)

        DO J = 1,NPY
           DO I = 1,NPX
              K(I,J) = 0.           ! initialize output image
           ENDDO
        ENDDO

        L = NL
        DO J = JMIN,JMAX
          DO I = IMIN,IMAX

            U = Y(L) - REAL(J)              ! shift to
            S = X(L) - REAL(I)              ! ellipses frame
C           T = S*COSD(C(L)) + U*SIND(C(L)) ! rotate to
C           V = U*COSD(C(L)) - S*SIND(C(L)) ! ellipses frame
            T = S*COS(C(L)/180*PI) + U*SIN(C(L)/180*PI) ! rotate to
            V = U*COS(C(L)/180*PI) - S*SIN(C(L)/180*PI) ! ellipses frame
            O = V/A(L)                      ! scale to
            P = T/B(L)                      ! ellipses frame

  900       Q = SQRT(O*O + P*P)             ! 'radius' of ellipses
            D = REAL(L)                     ! find the next
            L = L + NINT(SIGN(1.,(Q-1.)))   ! neigbouring
            L = MAX(L,1)                    ! ellipses, if
            L = MIN(L,NL)                   ! there is
            F = REAL(L)                     ! one.

            IF (ABS(D-F).GT.0.5) THEN       ! switch if R is outside

              U = Y(L) - REAL(J)              ! shift to
              S = X(L) - REAL(I)              ! ellipses frame
C             T = S*COSD(C(L)) + U*SIND(C(L)) ! rotate to
C             V = U*COSD(C(L)) - S*SIND(C(L)) ! ellipses frame
              T = S*COS(C(L)/180*PI) + U*SIN(C(L)/180*PI) ! rotate to
              V = U*COS(C(L)/180*PI) - S*SIN(C(L)/180*PI) ! ellipses frame
              O = V/A(L)                      ! scale to
              P = T/B(L)                      ! ellipses frame
              H = SQRT(O*O + P*P)             ! 'radius' of ellipses
              W = (1.-H)/(Q-H)                ! interpolate

              IF (W.LT.0. .OR. W.GT.1.) GO TO 900 ! not the right intervall

              F = F + (D-F)*W

              K(I,J) = EZ(2)*10.**(-CSTEP*F)+EZ(1)

            ENDIF
          ENDDO
        ENDDO

        RETURN
        END
