C @(#)spflx.for	19.1 (ES0-DMD) 02/25/03 13:19:40
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
      SUBROUTINE SPFLXA(XIMAG,NPIX,START,STEP,RMAX,PASR,FONDMX,
     2                  XNORM,XYZI,MASK,XYZ)
C ++++
C.INPUT/OUTPUT:  RMAX   [I] radius max of star analysis
C                PASR   [I] annulus of analysis width
C                XNORM  [I] level of magnitude normalization
C                FONDMX [I] level max for background
C.VERSION:       910115 RHW IMPLICIT NONE included
C ---
      IMPLICIT      NONE
C
      REAL          XIMAG(1)                                                
      INTEGER       NPIX(2)                                                  
      REAL          START(2)
      REAL          STEP(2)
      REAL          RMAX
      REAL          PASR
      REAL          FONDMX                                   
      REAL          XNORM
      REAL          XYZI(2)                                            
      INTEGER       MASK                                                   
      REAL          XYZ(8)                                             
C
      INTEGER       MADRID(1)
      INTEGER       KCLAST
      INTEGER       NO
      INTEGER       IERR
      INTEGER       NPL, NPL1, NPL2, NPL3
      INTEGER       NL, NLC, NL1, NL2
      INTEGER       IOPO, NBM, NPLC, IPL, ITEF
      INTEGER       I, J, IY, NPP
      INTEGER       IPV, KCL
C
      REAL          D1(256),DFLUX(256)
      REAL          AIRE(256)
      REAL          RMOY(256)      
      REAL          XFR, RP, RPOI, RPOJ
      REAL          FLUX, FOND, FLUXV, XMAGS, ROBJ
      REAL          DIAMV, RADIUS, XLOGD, DFCAL
      REAL          SURFAC
      REAL          XPW, YPW, XC0, YC0
      REAL          RMAX2, DXP
      REAL          YFR
C
      SAVE          NO
      COMMON        /VMR/MADRID
      DATA          NO/0/
C
      NO     = NO + 1
      NPL    = NPIX(1)
      NL     = NPIX(2)
      NPL3   = 1      
      IOPO   = 0      
      NBM    = NINT(RMAX) + 1
      PASR   = AMAX1(.7,PASR)
      KCLAST = RMAX/PASR
      SURFAC = ABS(STEP(1)*STEP(2))
C
C *** is star into the frame?
      IERR   = 0                                                             
      XPW    = XYZI(1)                                                    
      YPW    = XYZI(2)                                                    
      XC0    = (XPW-START(1))/STEP(1) + 1.                                    
      YC0    = (YPW-START(2))/STEP(2) + 1.                                    
      NPLC   = NINT(XC0)                                                     
      NPL1   = MAX0(NPLC-NBM,1)                                              
      NPL2   = MIN0(NPLC+NBM,NPL)                                            
      IF (NPLC.LT.1 .OR. NPLC.GT.NPIX(1)) IERR = -4                        
      NLC    = NINT(YC0)                                                     
      NL1    = MAX0(NLC-NBM,1)                                               
      NL2    = MIN0(NLC+NBM,NL)                                              
      IF (NLC.LT.1 .OR. NLC.GT.NPIX(2)) IERR = IERR - 5                    
      IF (IERR.LT.0) GO TO 100                                               
C                                
C *** yes it is
      NPP    = NPL2 - NPL1 + 1                                               
      DXP    = ABS(STEP(1))                                                  
C
C *** boucle de calcul du flux et des parametres y afferant
      DO 10 J = 1,256                                                        
         DFLUX(J) = 0.                                                       
         RMOY(J)  = 0.                                                        
         AIRE(J)  = 0.                                                        
   10 CONTINUE                                                               
C
      RMAX2  = RMAX*RMAX                                                     
      DO 30 IY = NL1,NL2                                                     
         YFR    = FLOAT(IY)                                                  
         RPOJ   = (YFR-YC0)**2                                               
         CALL LIRE(IY,NPL,NPL1,NPL2,NPL3,XIMAG,D1)                           
         DO 20 IPL = 1,NPP                                                   
            IPV    = NPL1 + IPL - 1                                          
            XFR    = FLOAT(IPV)                                              
            RPOI   = (XFR-XC0)**2 + RPOJ                                     
            IF (RPOI.LE.RMAX2) THEN                                          
               RP     = SQRT(RPOI)                                           
               KCL    = IFIX(RP/PASR) + 1                                    
               DFLUX(KCL) = DFLUX(KCL) + D1(IPL)                             
               RMOY(KCL) = RMOY(KCL) + RP                                    
               AIRE(KCL) = AIRE(KCL) + 1.                                    
            END IF                                                           
   20    CONTINUE                                                            
   30 CONTINUE                                                               
C
      FLUX   = 0.                                                            
      DO 40 I = 1,KCLAST                                                     
         IF (AIRE(I).GT.0.) THEN                                             
             RMOY(I) = RMOY(I)/AIRE(I)                                       
             FLUX   = FLUX + DFLUX(I)                                        
             DFLUX(I) = DFLUX(I)/AIRE(I)                                     
         END IF                                                              
   40 CONTINUE                                                               
C
C *** calcul du fond 
      CALL CALFON(DFLUX,KCLAST,FONDMX,FOND,ROBJ,IERR)                        
      IF (IERR.LT.0) GO TO 90                                                
C
C *** calcul du flux corrige par le fond calcule                               
      DO 50 I = 1,KCLAST - 1                                                 
         IF (ROBJ.GE.RMOY(I) .AND. ROBJ.LE.RMOY(I+1)) THEN                   
            ITEF   = I + 1                                                   
         ENDIF                                                               
   50 CONTINUE                                                               
C
      FLUXV  = 0.                                                            
      XMAGS  = 0.                                                            
      DO 60 J = 1,ITEF                                                       
         DFCAL  = DFLUX(J) - FOND                                            
         FLUXV  = FLUXV + DFCAL*AIRE(J)                                      
   60 CONTINUE                                                               
      DO 70 J = ITEF,1,-1                                                    
         DFCAL  = DFLUX(J) - FOND                                            
         IF (ABS(DFCAL).GT..001*FLUXV) GO TO 80                              
         FLUXV  = FLUXV - DFCAL*AIRE(J)                                      
   70 CONTINUE                                                               
      IERR   = -3                                                            
      RADIUS = FLOAT(KCLAST)*DXP                                             
      RETURN
                                                                                
   80 CONTINUE                                                               
      RADIUS = FLOAT(J)*DXP                                                  
      DIAMV  = 2.*RADIUS                                                     
      XLOGD  = ALOG10(DIAMV)                                                 
      FLUXV  = FLUXV*SURFAC                                                  
      IF (FLUXV.GT.0.) THEN                                                  
         XMAGS  = -2.5*ALOG10(FLUXV) + XNORM                                 
      ENDIF                                                                  
C 
      MASK     = 1                                                        
      XYZ(1)   = XPW                                           !  x center
      XYZ(2)   = YPW                                           !  y center
      XYZ(3)   = FLUXV                                   !  integrate flux
      XYZ(4)   = XMAGS                                        !  magnitude
      XYZ(5)   = FOND                                  !  local background
      XYZ(6)   = DIAMV                                         !  diameter
      XYZ(7)   = XLOGD                                 !  log. of diameter
      XYZ(8)   = NO                                    !  reference number
      GO TO 100                                                              
C                                                                            
   90 CONTINUE                                                               
      XYZ(8)   = NO                                    !  reference number
      MASK = IERR
  100 CONTINUE
      RETURN 
      END  


      SUBROUTINE SPFLXI(XIMAG,NPIX,START,STEP,RMAX,PASR,FONDMX,
     2                  XNORM,XYZI,MASK,XYZ)
C ++++ 
C.INPUT/OUTPUT:  RMAX   [I] radius max of star analysis
C                PASR   [I] annulus of analysis width
C                XNORM  [I] level of magnitude normalization
C                FONDMX [I] level max for background 
C ---
      IMPLICIT      NONE
C
      REAL          XIMAG(1)
      INTEGER       NPIX(2)
      REAL          START(2)
      REAL          STEP(2)
      REAL          RMAX
      REAL          PASR
      REAL          FONDMX
      REAL          XNORM
      REAL          XYZI(2)
      INTEGER       MASK
      REAL          XYZ(8)
C
      INTEGER       MADRID(1)
      INTEGER       NL1, NL2
      INTEGER       NPLC, NPL, NL
      INTEGER       NPL3, NPL1, NPL2
      INTEGER       IY
      INTEGER       NPP
      INTEGER       NLC, IPL, IPV
      INTEGER       ISTAT, IST
      INTEGER       IRAD
      INTEGER       KCL
      INTEGER       ITEF 
      INTEGER       I, K, J
      INTEGER       IOPO, NBM, KCLAST
      INTEGER       NO, IERR 
      INTEGER       ACCESS, PLMODE
      INTEGER       KA
C
      REAL          D1(256)
      REAL          XCUR,YCUR
      REAL          FRAME(8)
      REAL          XPW, YPW, XC0, YC0
      REAL          SURFAC
      REAL          DXP, RMAX2
      REAL          ROBJ, RPOI, RP
      REAL          XFR, YFR
      REAL          RPOJ, DFCAL
      REAL          FLUXV, XMAGS
      REAL          DIAMV, RADIUS, XLOGD
      REAL          FOND, FLUX
      REAL          SCALE(3)
      REAL          OFFSET(2)
      REAL          STAR(256,3)
      REAL          ZMAX, ZMIN
      REAL          XOUT, YOUT
C
      CHARACTER*80  TEXT
      CHARACTER*16  LABEL1,LABEL2
C      
      LOGICAL       PRINT

      INCLUDE       'MID_INCLUDE:PLTDEC.INC/NOLIST'
      COMMON        /VMR/MADRID
      SAVE          NO
C
      DATA          NO/0/
      DATA          SCALE/0.0,0.0,0.0/
      DATA          OFFSET/-999.,-999./
      DATA          LABEL1/'Radius'/
      DATA          LABEL2/'Pixel Value'/
C
  905 FORMAT(' User input: ', 41X,2(G12.6))
C
C *** START CODE
      NO     = NO + 1
      NPL    = NPIX(1)
      NL     = NPIX(2)
      NPL3   = 1
      IOPO   = 0
      NBM    = NINT(RMAX) + 1
      PASR   = AMAX1(.7,PASR) 
      KCLAST = RMAX/PASR 
      SURFAC = ABS(STEP(1)*STEP(2)) 
C
C *** is star into the frame? 
      IERR   = 0
      XPW    = XYZI(1)
      YPW    = XYZI(2)
      XC0    = (XPW-START(1))/STEP(1) + 1.
      YC0    = (YPW-START(2))/STEP(2) + 1.
      NPLC   = NINT(XC0)
      NPL1   = MAX0(NPLC-NBM,1)
      NPL2   = MIN0(NPLC+NBM,NPL)
      IF (NPLC.LT.1 .OR. NPLC.GT.NPIX(1)) IERR = -4
      NLC    = NINT(YC0)
      NL1    = MAX0(NLC-NBM,1)
      NL2    = MIN0(NLC+NBM,NL)
      IF (NLC.LT.1 .OR. NLC.GT.NPIX(2)) IERR = IERR - 5
      IF (IERR.LT.0) GO TO 100
C
C *** yes it is
      NPP    = NPL2 - NPL1 + 1                                               
      DXP    = ABS(STEP(1))                                                  
C    
C *** boucle de calcul du flux et des parametres y afferant
      DO 10 K = 1,3
         DO 11 J = 1,256
            STAR(J,K) = 0.
   11    CONTINUE
   10 CONTINUE
C
      RMAX2  = RMAX*RMAX                                                     
      DO 20 IY = NL1,NL2                                                     
         YFR    = FLOAT(IY)                                                  
         RPOJ   = (YFR-YC0)**2                                               
         CALL LIRE(IY,NPL,NPL1,NPL2,NPL3,XIMAG,D1)                           
         DO 21 IPL = 1,NPP                                                   
            IPV    = NPL1 + IPL - 1                                          
            XFR    = FLOAT(IPV)                                              
            RPOI   = (XFR-XC0)**2 + RPOJ                                     
            IF (RPOI.LE.RMAX2) THEN                                          
               RP     = SQRT(RPOI)                                           
               KCL    = IFIX(RP/PASR) + 1                                    
               STAR(KCL,2) = STAR(KCL,2) + D1(IPL)
               STAR(KCL,1) = STAR(KCL,1) + RP
               STAR(KCL,3) = STAR(KCL,3) + 1.
            END IF                                                           
   21    CONTINUE                                                            
   20 CONTINUE                                                               
C                                         
      FLUX = 0.                                                            
      DO 30 I = 1,KCLAST                                                     
         IF (STAR(I,3).GT.0.) THEN
            STAR(I,1) = STAR(I,1)/STAR(I,3)
            FLUX      = FLUX + STAR(I,2)
            STAR(I,2) = STAR(I,2)/STAR(I,3)
         END IF                                                              
   30 CONTINUE                                                               
C
C *** calcul du fond
      CALL CALFON(STAR(1,2),KCLAST,FONDMX,FOND,ROBJ,IERR)
      IF (IERR.LT.0) THEN
         RADIUS = FLOAT(KCLAST)*DXP
         GO TO 66
      ENDIF
C
C *** calcul du flux corrige par le fond calcule                               
      DO 40 I = 1,KCLAST - 1                                                 
         IF (ROBJ.GE.STAR(I,1) .AND. ROBJ.LE.STAR(I+1,1)) THEN
            ITEF   = I + 1                                                   
         ENDIF                                                               
   40 CONTINUE                                                               
C                                                                            
      FLUXV  = 0.                                                            
      XMAGS  = 0.                                                            
      DO 50 J = 1,ITEF                                                       
         DFCAL  = STAR(J,2) - FOND
         FLUXV  = FLUXV + DFCAL*STAR(J,3)
   50 CONTINUE                                                               
C
      DO 60 J = ITEF,1,-1                                                    
         DFCAL  = STAR(J,2) - FOND
         IF (ABS(DFCAL).GT..001*FLUXV) GO TO 65
         FLUXV  = FLUXV - DFCAL*STAR(J,3)
   60 CONTINUE                                                               
      IERR   = -3                                                            
      RADIUS = FLOAT(KCLAST)*DXP 
      GO TO 66
C
   65 CONTINUE                                                               
      RADIUS = FLOAT(J)*DXP
      DIAMV  = 2.*RADIUS                                                     
      XLOGD  = ALOG10(DIAMV)                                                 
      FLUXV  = FLUXV*SURFAC                                                  
      IF (FLUXV.GT.0.) THEN                                                  
         XMAGS  = -2.5*ALOG10(FLUXV) + XNORM                                 
      ENDIF                                                                  
      PRINT = .TRUE.
C
   66 CONTINUE	
      PRINT = .FALSE.
      ZMAX = STAR(1,2)
      DO 70 I = 2,KCLAST/3
         IF (ZMAX.LT.STAR(I,2)) THEN
            ZMAX = STAR(I,2)
         ENDIF
   70 CONTINUE
C
      ZMIN      = FOND
      FRAME(1)  = 0.
      FRAME(2)  = FLOAT(KCLAST)*DXP
      FRAME(5)  = ZMIN*1.1-0.1*ZMAX
      FRAME(6)  = ZMAX*1.1-0.1*ZMIN
      CALL GETAXS('AUTO',FRAME(1))
      CALL GETAXS('AUTO',FRAME(5))
      CALL PTKWRR('XWNDL',4,FRAME(1))
      CALL PTKWRR('YWNDL',4,FRAME(5))
      CALL PTKWRR('SCALE',3,SCALE)
      CALL PTKWRR('OFFSET',2,OFFSET)
C
      PLMODE = 1
      ACCESS = 0
      CALL PTOPEN(' ',' ',ACCESS,PLMODE)
      CALL PLFLUX(NO,STAR(1,1),STAR(1,2),KCLAST,
     2            RADIUS,FOND,XOUT,YOUT,IST)
      CALL PTFRAM(FRAME(1),FRAME(5),LABEL1,LABEL2)
C
C *** loop to set the background and radius
   80 CONTINUE

C *** get cursor
      CALL PTGCUR(XCUR,YCUR,KA,IST)
      IF (KA.EQ.32) THEN
         CALL PTCLOS()
         GO TO 90

      ELSE 
         RADIUS = XCUR
         FOND   = YCUR
         WRITE (TEXT,905) FOND, 2*RADIUS
         CALL STTPUT(TEXT,ISTAT)
         GO TO 80
      ENDIF
C
   90 CONTINUE
      FLUXV  = 0.0
      XMAGS  = 0.0
      IRAD   = NINT(RADIUS)
      DO 95 J = 1,IRAD
         DFCAL = STAR(J,2)-FOND
         FLUXV = FLUXV+DFCAL*STAR(J,3)
   95 CONTINUE
      DIAMV = 2.*RADIUS
      XLOGD = ALOG10(DIAMV)
      FLUXV = FLUXV*SURFAC
C
      IF (FLUXV.GT.0.) THEN
         XMAGS = -2.5*ALOG10(FLUXV)+XNORM
      ENDIF
C
      IF (IST.LE.0 .AND. IERR.LT.0) THEN
         GO TO 999
      ENDIF
C
      MASK   = 1
      XYZ(1) = XPW                               ! x center
      XYZ(2) = YPW                               ! y center
      XYZ(3) = FLUXV	                         ! integrate flux
      XYZ(4) = XMAGS	                         ! magnitude
      XYZ(5) = FOND	                         ! local background
      XYZ(6) = DIAMV	                         ! diameter
      XYZ(7) = XLOGD	                         ! log. of diameter
      XYZ(8) = NO                                ! reference number
      GO TO 998

  999 CONTINUE
      XYZ(8) = NO                                ! reference number
      MASK = IERR
  998 CONTINUE

  100 CONTINUE

      RETURN
      END





