C @(#)fitell2.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
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C.IDENTIFICATION: ELLIPS.FOR
C.PURPOSE:   Fit ellipses to 2-dim images,retrieve parameters relevant 
C            to isophotal twisti
C.LANGUAGE:  F77+ESOext
C.AUTHOR:    Edwin A. Valentijn
C.KEYWORDS:  Ellipse fitting,isophotal twisting, display fit
C.ALGORITHM: Make polar coord frames, search for points selected on radius and
C            intensity, least square fit to the polar coord of thee points,
C            retrieve ellips parameters, display points + fit
C.VERSION:   821101 EAV  creation
C.VERSION:   871002 RHW  conversion to standard interfaces
C.VERSION:   871015 RHW  conversion to standard fortran
C 
C 021205	last modif
C 
C ---------------------------------------------------------------------
      PROGRAM FITELL2
C
      IMPLICIT  NONE
C
      INTEGER   MADRID
      INTEGER   I,J,K,NN,IAV,ICAL,IPOL,ISTAT
      INTEGER   MAXLEV,NAXIS,NGAL,NP,NX,NY
      INTEGER   STAT,XSTEP,YSTEP
      INTEGER   NA,NF,NM,NR
      INTEGER   NPIX(2),N(10)
      INTEGER   KUN,KNUL
      INTEGER*8   PNTRA,PNTRF,PNTRM,PNTRR
C
      DOUBLE PRECISION  START(2),STEP(2)
      REAL      IRLIM,XSTRT,YSTRT,CENTER(2),SKY
      REAL      RAD(500),ANG(500),YV(500),X1(500),X2(500),X3(500)
      REAL      CUTS(4),ISOP1(10),ISOP
      REAL      CO1,CO2,CON,CC,B(10),A(10),E(10),EX(10),ANGU(10),PI,EE
C
      CHARACTER IDENT*72,CUNIT*48,TEXT*80,CINPUT*72,INSTRM*72
      CHARACTER INFR*60,RADII*8,ANGLI*8,MASK*8,CTYPE
C
      INCLUDE   'MID_INCLUDE:ST_DEF.INC/NOLIST'
      COMMON    /VMR/MADRID(1)
      INCLUDE   'MID_INCLUDE:ST_DAT.INC/NOLIST'
C
      DATA      PI/3.1415926/
      DATA      MAXLEV/10/
      DATA      CTYPE/'R'/
C
C *** start code
C
      CALL STSPRO('FITELL2 ')                        ! set up Midas environment
C
C ** get keywords
C  
      CALL STKRDC('IN_A',1,1,60,IAV,INFR,KUN,KNUL,STAT)       !  get frame name
      CALL STKRDI('INPUTI',1,1,IAV,IPOL,KUN,KNUL,STAT)         !  polar frame ?
      CALL STKRDI('INPUTI',2,1,IAV,ICAL,KUN,KNUL,STAT)      !  isophotal window
      CALL STKRDC('INPUTC',1,1,72,IAV,CINPUT,KUN,KNUL,STAT)         !  contours
      CALL GETINP(CTYPE,CINPUT,MAXLEV,ISOP1,NGAL,STAT)
      IF (NGAL.GT.1) THEN
         CALL SORT1(ISOP1,NGAL)  
      ENDIF
      CALL STKRDR('INPUTR',1,2,IAV,CENTER,KUN,KNUL,STAT)              !  center
      CALL STKRDR('INPUTR',3,1,IAV,IRLIM,KUN,KNUL,STAT) !  radius around center
      CALL STKRDR('INPUTR',4,1,IAV,SKY,KUN,KNUL,STAT)         !  sky background
C
C *** get the image
C
      CALL STIGET(INFR,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXIS,NPIX,START,STEP,IDENT,
     +            CUNIT,PNTRF,NF,STAT)
C
C  initialize coordinates, parameters
C
      NX     = NPIX(1)
      NY     = NPIX(2)
      NP     = NX*NY
      XSTRT  = SNGL(START(1))
      YSTRT  = SNGL(START(2))
      XSTEP  = SNGL(STEP(1))
      YSTEP  = SNGL(STEP(2))
C
C *** polar coordinates
C     create two images with polar coord: radii and angles
C
      RADII  = 'SCRRADII'
      ANGLI  = 'SCRANGLI'
C
C *** but enable run of program with preexisting polar frames
C
      IF (IPOL.EQ.1) THEN
          CALL STIGET(RADII,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +                1,NAXIS,NPIX,START,STEP,
     +                INSTRM,CUNIT,PNTRR,NR,STAT)
          CALL STIGET(ANGLI,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +                1,NAXIS,NPIX,START,STEP,
     +                INSTRM,CUNIT,PNTRA,NA,STAT)
      ELSE
          CALL STIPUT(RADII,10,1,1,NAXIS,NPIX,START,STEP,
     +                'POLAR-RAD',CUNIT,PNTRR,NR,STAT)
          CALL STIPUT(ANGLI,10,1,1,NAXIS,NPIX,START,STEP,
     +                'POL-ANG',CUNIT,PNTRA,NA,STAT)
          CALL POLAR(NP,NPIX,START,STEP,CENTER,
     +               MADRID(PNTRR),MADRID(PNTRA))
      END IF
C
C     create image for retrieved points and ellipses
C
      MASK   = 'ELLIPS  '
      CUTS(1) = 0.
      CUTS(3) = 0.
      CUTS(2) = NGAL*10.
      CUTS(4) = NGAL*10.
      CALL STIPUT(MASK,10,1,1,NAXIS,NPIX,START,STEP,
     +            'ELLIPS',CUNIT,PNTRM,NM,STAT)
      CALL STDWRR(NM,'LHCUTS',CUTS,1,4,KUN,STAT)
      CALL STDWRC(NM,'INSTRUM',1,INSTRM,1,72,KUN,STAT)
C
C *** Enter main loop
C
      DO 30 K = 1,NGAL
C
C *** Search for points wthin radius and isophotal window
C     For found points retrieve only their polar coord in RAD and ANG
C
         NN = N(K)
         ISOP = ISOP1(K)
         CALL FITEL(NP,MADRID(PNTRF),MADRID(PNTRR),
     +              MADRID(PNTRA),RAD,ANG,ISOP,
     +              ICAL,SKY,IRLIM,NN)
         N(K) = NN
         ISOP1(K) = ISOP
C
C *** Check on the actual number of points found
C
         IF (N(K).GT.499) THEN
            CALL STTPUT('*** WARNING: More then 500 points found '//
     +                   ' results unreliable',ISTAT)
            CALL STTPUT('*** ADVICE:  make window smaller',ISTAT)
            N(K) = 500
         END IF
C
C *** Continue the proces when no or few points were found
C
         IF (N(K).LT.6) THEN
            A(K)  = 0.
            B(K)  = 0.
            E(K)  = 0.
            EX(K) = 0.
            ANGU(K) = 0.
            GO TO 20
         END IF
C
C *** set up for least square fit
C
         DO 10 I = 1,N(K)
            YV(I)  = 2./RAD(I)**2
            X1(I)  = SIN(2.*ANG(I))
            X2(I)  = COS(2.*ANG(I))
            X3(I)  = 1.0
   10    CONTINUE
         NN = N(K)
         CALL LSS3(NN,YV,X1,X2,X3,CO1,CO2,CON)
         N(K) = NN
C
C *** recover the fit parameters  ABS to avoid interuption in batch
C
         ANGU(K) = -0.5*ATAN2(CO1,CO2)
         ANGU(K) = 180.*ANGU(K)/PI + 90.0
         CC       = SQRT(ABS(CO1*CO1+CO2*CO2))
         B(K)    = SQRT(ABS(2./ (CON+CC)))
         A(K)    = SQRT(ABS(2./ (CON-CC)))
         E(K)    = A(K)/B(K)
         EE       = B(K)/A(K)
         EX(K)   = SQRT(ABS(1.-EE*EE))
C
C *** retrieve found datapoints and ellips in mask frame
C
         NN = N(K)
         CALL PLMASK(K,RAD,ANG,NN,NP,NPIX,START,STEP,CENTER,
     +               MADRID(PNTRM),CON,CO1,CO2)
         N(K) = NN
   20    CONTINUE
   30 CONTINUE
C
C *** Print table
C
      CALL STTPUT(' ',ISTAT)
      CALL STTPUT('*** Summary of the Results ***',ISTAT)
      CALL STTPUT('Isophote      #Points     A         B         '//
     +            'A/B      exc.     angle',ISTAT)
      DO 40 J = 1,NGAL
         WRITE (TEXT,8000) ISOP1(J),N(J),A(J),B(J),E(J),EX(J),ANGU(J)
         CALL STTPUT(TEXT,ISTAT)
   40 CONTINUE

      CALL STSEPI
      STOP ' '
C
 8000 FORMAT (G9.3,2X,I8,5F10.2)
C
      END


C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENTIFICATION: POLAR
C.PURPOSE:        Creates 2 virtual images containing the polar coordinates 
C                 with respect to CENTER using descripters of a MIDAS image
C.AUTHOR:         E.A.Valentijn  ESO/Garching
C.INPUT/OUTPUT:   NP,NPIX,START,STEP standard descriptor parameters
C                 radii  go to image R
C                 angles go to image ANGLE
C.VERSION:        821028  EAV  tested and works ok
C.VERSION:        880707  RHW  Converted to standard fortran
C --------------------------------------------------------------------
      SUBROUTINE POLAR(NP,NPIX,START,STEP,CENTER,R,ANGLE)
C
      IMPLICIT NONE
C
      COMMON /VMR/MADRID(1)
C
      INTEGER   MADRID
      INTEGER   NP,NPIX(2),NX,NY,I,IX,IY
      REAL      X,Y,XSTRT,YSTRT,XSTEP,YSTEP,XCT,YCT,CENTER(2)
      REAL      R(NP),ANGLE(NP)
      DOUBLE PRECISION START(2),STEP(2)
C
      XSTRT  = SNGL(START(1))
      YSTRT  = SNGL(START(2))
      XSTEP  = SNGL(STEP(1))
      YSTEP  = SNGL(STEP(2))
      XCT    = CENTER(1)
      YCT    = CENTER(2)
      NX     = NPIX(1)
      NY     = NPIX(2)
      I      = 0
C
      DO 20 IY = 1,NY
         Y      = YSTRT + ((IY-1)*YSTEP)
         DO 10 IX = 1,NX
            X      = XSTRT + ((IX-1)*XSTEP)
            I      = I + 1
            R(I)   = SQRT((X-XCT)**2+ (Y-YCT)**2)
            IF (R(I).EQ.0) THEN                !   avoid singularity at center
               R(I)   = 0.D-10
               ANGLE(I) = 0.0D-10
            ELSE
               ANGLE(I) = ATAN2(X-XCT,Y-YCT)
            END IF
   10    CONTINUE
   20 CONTINUE
      RETURN
      END



C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENTIFICATION: FITEL
C.AUTHOR:   Edwin A. Valentijn
C.PURPOSE:  Selects out of image a on int and radius and fills rad and ang
C.INPUT:    NP    total number of pixels
C           A     image intensities                           input
C           R     polar coord image of raddi                  input
C           ANGLE polar coord image of angles                 input
C           RAD   radii of found points                       output
C           ANG   angles of found points                      output
C           ISOP1 isophotal level in % of background          keyword
C           ICAL  isophotal window in pro mille of background keyword
C           SKY   background level to be added before fitting keyword
C           IRLIM maximum radius used for search              keyword
C           N     number of points found                      output
C.VERSION:  880707 Converted to standard fortran
C --------------------------------------------------------------------
      SUBROUTINE FITEL(NP,A,R,ANGLE,RAD,ANG,ISOP1,ICAL,SKY,IRLIM,N)
C
      IMPLICIT NONE
C
      COMMON /VMR/MADRID(1)
C
      INTEGER   MADRID
      INTEGER   NP,N,ICAL,NI,ISTAT
      REAL      A(NP),R(NP),ANGLE(NP)
      REAL      RAD(500),ANG(500),ISOP1
      REAL      SKY(2),RISOP,LOLIM,UPLIM,CAL1,IRLIM
      CHARACTER TEXT*80
C
      RISOP = ISOP1
      CAL1  = (ICAL*ISOP1)/1000.
      LOLIM = SKY(1) + RISOP - CAL1
      UPLIM = SKY(1) + RISOP + CAL1
C
C *** search int array for filtered values
C     simultaneously R array for radii
C     if true put corresponding R and ANGLE in RAD and ANG
C
      N = 0
      DO 10 NI = 1,NP
         IF (R(NI).LT.IRLIM .AND. A(NI).LT.UPLIM .AND. 
     2      A(NI).GT.LOLIM) THEN
            N      = N + 1
            IF (N.LE.500) THEN
               RAD(N) = R(NI)
               ANG(N) = ANGLE(NI)
            END IF
         END IF
   10 CONTINUE

      WRITE (TEXT,9000) ISOP1,ICAL,LOLIM,UPLIM,IRLIM,N
      CALL STTPUT(TEXT,ISTAT)
C
      RETURN
C
 9000 FORMAT ('Isoph=',G9.3,'  Window=',I4,'  Lim=',F8.2,1X,F8.2,
     +        '  Rad=',F8.2,'  #Pnts=',I4)
C
      END


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENTIFICATION: LSS3
C.PURPOSE:        least square solution of a system of equations with 3 unknown:
C                 A * X  +  B * Y  +  C * Z  =  CNS
C.AUTHOR:         Edwin A. Valentijn
C.INPUT.OUTPUT:   N  :   number of equations
C                 CNS:   array containing the constants
C                 AAR:   coef.  of X
C                 BAR:   coef.  of Y
C                 CAR:   coef.  of Z
C                 X  :   least square solution of first unknown
C                 Y  :                            second
C                 Z  :                            third
C.VERSION:        880707 RHW  conversion to standard fortran
C -----------------------------------------------------------------
      SUBROUTINE LSS3(N,CNS,AAR,BAR,CAR,X,Y,Z)
C
      IMPLICIT NONE
C
      COMMON /VMR/MADRID(1)
C 
      INTEGER   MADRID
      INTEGER   N,I
      REAL      CNS(N),AAR(N),BAR(N),CAR(N),X,Y,Z
      DOUBLE PRECISION AA,AB,AC,AN,BB,BC,BN,CC,CN,DA,DB,DC,DN
      DOUBLE PRECISION FBB,FBN,FBC,FCC,FCN
C
      IF (N.LT.3) THEN               !  Check if there are more then 3 equations
         RETURN             
      ENDIF
C
C *** initiate the sum - variables
C
      AA     = 0.0D+0
      AB     = 0.0D+0
      AC     = 0.0D+0
      AN     = 0.0D+0
      BB     = 0.0D+0
      BC     = 0.0D+0
      BN     = 0.0D+0
      CC     = 0.0D+0
      CN     = 0.0D+0
C  
      DO 10 I = 1,N                                      !   start the summing
         DA     = DBLE(AAR(I))
         DB     = DBLE(BAR(I))
         DC     = DBLE(CAR(I))
         DN     = DBLE(CNS(I))
         AA     = AA + DA*DA
         AB     = AB + DA*DB
         AC     = AC + DA*DC
         AN     = AN + DA*DN
         BB     = BB + DB*DB
         BC     = BC + DB*DC
         BN     = BN + DB*DN
         CC     = CC + DC*DC
         CN     = CN + DC*DN
   10 CONTINUE
C
C *** first step in solving
C
      FBB    = BB - AB*AB/AA
      FBC    = BC - AB*AC/AA
      FCC    = CC - AC*AC/AA
      FCN    = CN - AC*AN/AA
      FBN    = BN - AB*AN/AA
C
C *** get the unknown
C
      Z      = SNGL((FCN-FBC*FBN/FBB)/ (FCC-FBC*FBC/FBB))
      Y      = SNGL((FBN-FBC*Z)/FBB)
      X      = SNGL((AN-AB*Y-AC*Z)/AA)
C
      RETURN
      END


C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENTIFICATION: PLMASK
C.PURPOSE:        Finds rectangular image coordinates and labels that in mask
C                 frame(A) of:
C                 set of (N) points in polar coord given by RAD and ANG
C                 ellips given by least square parameters CON,CO1,CO2
C.AUTHOR:         Edwin A. Valentijn
C.NOTE:           mask does not work : instead use normal image
C.VERSION:        880707  RHW  conversion to standard fortran
C -------------------------------------------------------------------
      SUBROUTINE PLMASK(K,RAD,ANG,N,NP,NPIX,START,STEP,CENTER,A,CON,
     +                  CO1,CO2)
C
      IMPLICIT NONE
C
      COMMON   /VMR/MADRID(1)
C
      INTEGER   MADRID
      INTEGER   K,N,NP,NPIX(2),I,II,IX,IY
      REAL      A(NP)
      REAL      RAD(500),ANG(500),CENTER(2)
      REAL      R,TH,CON,CO1,CO2,XX,YY
      DOUBLE PRECISION START(2),STEP(2)
C
C *** initialize Mask frame
C
      IF (K.EQ.1) THEN
         DO 10 I = 1,NP
             A(I)   = 0.
   10    CONTINUE
      END IF
C
C *** find indices of points of the polar coord array
C
      DO 20 I = 1,N
         XX    = RAD(I)*SIN(ANG(I))
         YY    = RAD(I)*COS(ANG(I))
         IX    = IFIX((XX+CENTER(1)+0.5-SNGL(START(1)))/SNGL(STEP(1)))
         IY    = IFIX((YY+CENTER(2)+0.5-SNGL(START(2)))/SNGL(STEP(2)))
         II    = (IY-1)*NPIX(1) + IX                   !    put a label in mask
         A(II) = K*10.
   20 CONTINUE
C
C *** now for the ellips; take 90 points duplication does not harm
C     data of found points from previous step are lost if they coincide
C
      DO 30 I = 1,91
         TH    = (I-1)*4.0
         TH    = 3.1415926*TH/180.
         R     = CON + CO1*SIN(2.*TH) + CO2*COS(2.*TH)
         R     = SQRT(ABS(2.0/R))
         XX    = R*SIN(TH)
         YY    = R*COS(TH)
         IX    = IFIX((XX+CENTER(1)+0.5-SNGL(START(1)))/SNGL(STEP(1)))
         IY    = IFIX((YY+CENTER(2)+0.5-SNGL(START(2)))/SNGL(STEP(2)))
         II    = (IY-1)*NPIX(1) + IX
         A(II) = K*10.
   30 CONTINUE
      
      RETURN
      END

