
CD-Section: Integration and Soft Limit Routines
C A set of routines is available to measure integrated intensities for
C a Laue image. They were written by Hao Quan to determine integrated
C intensities of reflections for use in improving the Laue soft limits
C 'lambda-min' and 'dmin' and routines are included for determining these
C improved values. They make use of the Laue Data Module (LDM)
C routines and work with an image read into memory. 
CD-end

C
C LAST UPDATE:
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C**************************
C**     LFN_SAINT        **
C**************************
C
CD-Routine: Standalone Integration Routine - LFN_SAINT
C This routine is used to measure integrated intensities from a Laue image
C stored in memory. It may be used in a one pass mode to carry out box
C integration or, more usually, in a two pass mode to carry out profile fitted
C integration. The Laue Reflection List (LRL) must be generated for the
C plate being integrated prior to calling LFN_SAINT. Note also that
C when profile fitted integration is done for a plate, the spot size
C parameters used will be those used in the profile determination pass
C which might not necessarily be the same as those defined for the plate
C in question. A special option is available to carry out a profile
C fitted integration for an individual spot for diagnostic purposes.
C The LRL must be sorted on 'xfd' or 'yfd' before calling the routine.
C The best order is that which will progress through the image in the
C order that the image data is stored in memory.
CD-end
CD-Fortran:
      SUBROUTINE LFN_SAINT (IPACK, IPLATE, IPASS, IMG, ITYPE, IORDER,
     +                     NF_OFF, NAX1_RASTS, NAX2_RASTS, PROMIN,
     +                     NTHE, PROF_ROT, NOVLD, PROGRESS, CANCEL,
     +                     IERR, ERRSTR)
CD-end
      INCLUDE 'implicit_none.finc'
C
C PURPOSE: Stand-alone integration program (full image version)
C 
C AUTHOR: Q. Hao (Modified by J.W. Campbell)
C
C ARGUMENTS:
C
      INTEGER IPACK
      INTEGER IPLATE 
      INTEGER IPASS
      INTEGER IMG(*)
      INTEGER ITYPE
      INTEGER IORDER
      INTEGER NF_OFF
      INTEGER NAX1_RASTS
      INTEGER NAX2_RASTS
      REAL PROMIN
      INTEGER NTHE
      LOGICAL PROF_ROT
      INTEGER NOVLD
      EXTERNAL PROGRESS
      LOGICAL CANCEL
      EXTERNAL CANCEL
      INTEGER IERR
      CHARACTER*(*) ERRSTR
CD-Parameters:
C   IPACK      i (R)  Pack number
C   IPLATE     i (R)  Plate number
C   IPASS      i (R)  =0  Box integration
C                     =1  forming profile
C                     =2  profile fitting
C                     >10 profile fit individual spot no. (IPASS-10) in
C                         spot list (for diagnostic purposes - intensity
C                         will not be stored in LRL)
C                         Data may be returned using LFN_GET_INTSPOT
C   IMG (*)    i (R)  Image data array (may hold signed integer data or
C                     unsigned byte data packed into an integer array or
C                     unsigned two byte data packed into an integer array    
C   ITYPE      i (R)  Image data type =1 unsigned byte
C                                     =2 unsigned two byte
C                                     =3 signed (full) integer
C                                     =4 'squashed i2'
C   IORDER     i (R)  Order of data in image wrt. the two local axes
C                     =1 +ax1 +ax2
C                     =2 +ax1 -ax2
C                     =3 -ax1 +ax2
C                     =4 -ax1 -ax2
C                     =5 +ax2 +ax1
C                     =6 +ax2 -ax1
C                     =7 -ax2 +ax1
C                     =8 -ax2 -ax1
C   NF_OFF     i (R)  Fast offset between (slow) rasters of image 
C                     (e.g. first dimension of the image data if the data 
C                     are held in a two dimensional array). This value is 
C                     the number of elements for the image data type 
C                     (bytes, 2-bytes or integers as set by ITYPE)
C   NAX1_RASTS i (R)  No. of rasters in the image along the first 
C                     local axis 
C   NAX2_RASTS i (R)  No. of rasters in the image along the second 
C                     local axis
C   PROMIN     r (R)  Minimum intensity to include in profile determination    
C                     e.g. 100 for film, 600 for image plate
C   NTHE       i (R)  Number of theta bins for profile fitting (e.g. 8,
C                     max=17 unless program limit was reset)
C   PROF_ROT   l (R)  Flag to rotate/do-not-rotate profiles data to theta
C                     = 0.0 around the beam direction when forming profiles 
C                     from elliptical spots. (.TRUE./.FALSE.)
C                     No rotation is done when spots are circular. 
C                     Allowing rotation means that account is taken of the 
C                     lengths and theta angles around the beam direction
C                     of the individual spots used to form the profiles
C                     but at the expense of accuracy in the mapping of the
C                     pixel positions done during the rotation. When
C                     using the profile for integrating a spot, the 
C                     reverse rotation is applied to the bin profile. If 
C                     no rotation is done, the profiles are gatherd by 
C                     superimposing the spots used to determine a bin 
C                     profile in the orientations in which they occur and 
C                     the profile is used as gathered when profile fitting 
C                     a spot for integration. Rotation may well be 
C                     beneficial if the spots cover many pixels but should 
C                     be used with caution if the spots only cover a few 
C                     pixels. This flag is only used in a profile forming 
C                     pass. When doing profile_fitted integration the
C                     value of the flag used when forming the profiles
C                     applies.
C   NOVLD      i (R)  Number of overload pixels allowed in a spot for
C                     it still to be integrated as a 'good' spot which
C                     will not be flagged as having a overload intensity.
C   PROGRESS     ( )  Subroutine which will be called at intervals during
C                     the integration which may be used to display its 
C                     progress if required (See below for its arguments).
C                     If not required, the dummy routine LFN_SOFT_PROG
C                     may be passed. Note that the routine must be
C                     declared as EXTERNAL in the calling program.
C   CANCEL     f (R)  Logical function used to interrupt procedure
C                     It has one dummy argument and must return a
C                     logical value of .TRUE. to cancel the procedure
C                     or .FALSE. to let it proceed.
C                     The actual function passed must be declared as
C                     both LOGICAL and EXTERNAL in the calling program.
C                     (The dummy function LFN_DUMMY_CNCL may be used)
C   IERR       i (W)  Error flag =  0, OK
C                                = -1, Invalid pack number
C                                = -2, Invalid plate number
C                                = -3, Invalid image type or order flag
C                                = -4, Byte handling not available
C                                =  1, Box size required based on spot size
C                                      exceeds current program limits
C                                =  2, Too many bins requested for profile 
C                                      fitting
C                                =  3, Invalid pass flag
C                                =  4, LRL not generated for current plate
C                                =  5, No profiles calculated
C                                =  6, Spot not present for individual
C                                      spot integration request
C                                =  7, LRL not sorted on 'xfd' or 'yfd'
C                                = 10, Integration of 'too close to
C                                      integrate' single spot requested
C                                = 11, Cannot integrate - <3 non-overlapped
C                                      background pixels
C                                = 12, Cannot integrate - no non-overlapped/
C                                      non-overloaded peak pixels (may also 
C                                      be <3 non-overlapped background pixels)
C                                = 20, Out of limits box - cannot integrate
C                                      single spot
C                                = 100 Interrupted via 'cancel' routine
C   ERRSTR     c (W)  Error message string (max 80 chars) set if IERR 
C                     is non-zero
C
C
C                    Subroutine PROGRESS (IPASS, NSPOT, ISPOT)
C
C                    This routine called by the integration routine
C                    with the following parameters:-
C
C                        IPASS i (R) Integration pass = 0 box integration,
C                                                     = 1 get profiles,
C                                                     = 2 profile fitted
C                                                         integration
C                        NSPOT i (R) No. of spots to be integrated
C                        ISPOT i (W) Last spot integrated
CD-end

CD-doc:
C The routine sets the integration error code (see LRL_GET_INT, LRL_SET_INT)
C for bad spots as follows (the Bad spot flag will also be set):
C
C bit 0 no non-overlapped/non-overloded peak pixels.
C bit 1 < 20% of peak pixels are non-overlapped/non-overloaded.
C bit 2 <3 non-overlapped background pixels.
C bit 4 < 20% of background pixels are non-overlapped.
C
C The 20% (approx) values are calculated from the expression (NPIX-3)/5 
C where NPIX is the number of pixels is the peak or background as relevant.
CD-end
C
C====== Parameter statements
C 
      INCLUDE 'saint_limits.finc'
C
C====== GLOBALS:
C
      INCLUDE 'lfn_saint.finc'
C
C====== FUNCTION TYPES:
C
      EXTERNAL LFN_SAINT_BDAT
C
C====== LOCALS:
C
      INTEGER I
      INTEGER J
      INTEGER NSPOT
      INTEGER OVPIX
      INTEGER IDATA((2*MAX_HPR+1)*(2*MAX_HPR+1))
      INTEGER IWORK(2*MAX_HPR+1)
      INTEGER ND
      INTEGER IFLAG
      INTEGER JERR 
      INTEGER IH
      INTEGER IK
      INTEGER IL
      INTEGER MULT
      INTEGER NAX
      INTEGER N_HALF
      INTEGER N
      INTEGER NBIN
      INTEGER NTEMP
      INTEGER IYD
      INTEGER IXD
      INTEGER IX
      INTEGER IY
      INTEGER IMTYP
      INTEGER N1OD
      INTEGER NTH
      INTEGER IPK
      INTEGER IPL
      INTEGER KPK
      INTEGER KPL
      INTEGER ISPOT_MIN
      INTEGER ISPOT_MAX
      INTEGER KPASS
      INTEGER KSPOT
      INTEGER NUM_OV
      INTEGER IOVTYP
      INTEGER ISORT
      INTEGER INTCOD
      INTEGER NUM_NOT_BG
      INTEGER NUM_IN_BG
      INTEGER NUM_IN_PK
      INTEGER NPK_MIN
      INTEGER NBG_MIN
      INTEGER NPK_OK
      INTEGER NBG_OK
      REAL BASEOD
      REAL NONL
      REAL CFAC
      REAL BASE
      REAL TRUOD
      REAL TWOPI
      REAL X_CEN
      REAL Y_CEN
      REAL SPOTL
      REAL SPOTW
      REAL SPOTF
      REAL BORDER
      REAL SP_L_CUR
      REAL SP_W_CUR
      REAL SP_F_CUR
      REAL SP_B_CUR
      REAL SP_D_CUR
      REAL SCALE
      REAL MM_RAST_X
      REAL MM_RAST_Y
      REAL CTOF
      REAL RMAX
      REAL X
      REAL Y
      REAL ALMBDA
      REAL DMINTHR2
      REAL OD((2*MAX_HPR+1)*(2*MAX_HPR+1))
      REAL XM
      REAL YM
      REAL RADIUS
      REAL AMAJOR
      REAL AMAJOR_1
      REAL AMAJOR_B
      REAL BMINOR
      REAL BMINOR_1
      REAL BMINOR_B
      REAL THETA
      REAL STHETA
      REAL CTHETA
      REAL XPRIME
      REAL YPRIME
      REAL ELLIP
      REAL ELLIP_1
      REAL ELLIP_B
      REAL SIZ
      REAL AVP
      REAL SUM2
      REAL SUMI2(MAX_PRBIN)
      REAL SIGK
      REAL DELTAM
      REAL X_OV(MAX_OVLPS)
      REAL Y_OV(MAX_OVLPS)
      REAL XPN
      REAL YPN
      CHARACTER*10 TYPSTR
      LOGICAL ERR
      LOGICAL SPAT
      LOGICAL CLOS
      LOGICAL NODFLG
      LOGICAL IOVER
      LOGICAL ROT
      LOGICAL DISTOR
      LOGICAL INTEG
C
C     INTEGER I         !Temp/loop variable
C     INTEGER J         !Temp/loop variable
C     INTEGER NSPOT	!Number of spots
C     INTEGER OVPIX	!Value for the overload pixel
C     INTEGER IDATA()   !INTEGER ARRAY TO RETURN THE DENSITY VALUES
C     INTEGER IWORK()   !Integer work array
C     INTEGER ND	!Number of decimal places for print
C     INTEGER IFLAG 	!Flag = -1 undefined value
C                            =  0 default value 
C                            =  1 value set explicitly 
C                            =  2 value set globally 
C     INTEGER JERR      !Error flag 
C     INTEGER IH	!Reflection 'h' index 
C     INTEGER IK	!Reflection 'k' index
C     INTEGER IL	!Reflection 'l' index
C     INTEGER MULT	!Multiplicity
C     INTEGER NAX       !Size of the section in pixels along both axes
C     INTEGER N_HALF	!NAX/2
C     INTEGER N		!Loop counter over spot number
C     INTEGER NBIN	!Bin number
C     INTEGER NTEMP     !Temporary number
C     INTEGER IXD	!X pixel position of spot centre
C     INTEGER IYD	!Y pixel position of spot centre
C     INTEGER IX	!Loop counter over X pixel
C     INTEGER IY	!Loop counter over Y pixel
C     INTEGER IMTYP	!Code number 1=film, 2=ip
C     INTEGER N1OD	!Pixel value for 1 OD
C     INTEGER NTH       !No. of theta bins being used
C     INTEGER IPK       !Pack no. to use for LDM parameters
C     INTEGER IPL       !Plate no. to use for LDM parameters
C     INTEGER KPK       !Last generated pack no. from LRL
C     INTEGER KPL       !Last generated plate no. from LRL
C     INTEGER ISPOT_MIN !Minimum spot no. for loop
C     INTEGER ISPOT_MAX !Maximum spot no. for loop
C     INTEGER KPASS     !Pass flag 0-2. box int., gather profs., profile int.
C     INTEGER KSPOT     !Spot no. for individual integration (0=all spots)
C     INTEGER IOVTYP    !Overlaps option
C     INTEGER ISORT     !LRL sort order code
C     INTEGER INTCOD    !Integration type code
C     INTEGER NUM_NOT_BG!No. non-background pixels in spot ellipse
C     INTEGER NUM_IN_BG !No. background pixels in spot ellipse
C     INTEGER NUM_IN_PK !No. pixels in peak
C     INTEGER NPK_MIN   !Minimum no. non-overlapped pixels in spot to  
C                       !proceed with spot integration (Declare BAD if fewer
C                       !found. (~20% of peak pixels)
C     INTEGER NBG_MIN   !Minimum no. non-overlapped pixels in border to  
C                       !proceed with spot integration (Declare BAD if fewer
C                       !found. (~20% of bg pixels)
C     INTEGER NPK_OK    !No. non-overlapped peak pixels
C     INTEGER NBG_OK    !No. non-overlapped background pixels
C     REAL BASEOD	!Optical density of the film base (OD unit)
C     REAL NONL		!Non-linearity parameter for film 
C     REAL CFAC		!Temp. parameter
C     REAL BASE		!Optical density of the film base (pixel unit)
C     REAL TRUOD	!True pixel value
C     REAL TWOPI	!2.0*PI
C     REAL X_CEN   !Corrected 'x' centre position in rasters (wrt image origin)
C     REAL Y_CEN   !Corrected 'y' centre position in rasters (wrt image origin)	
C     REAL SPOTL	!Spot length in mm - profs plate
C     REAL SPOTW	!Spot width in mm - profs plate
C     REAL SPOTF	!Spot length factor - profs plate
C     REAL BORDER	!Spot background border width in pixels - profs plate
C     REAL SP_L_CUR     !Spot length in mm - current plate
C     REAL SP_W_CUR     !Spot width in mm - current plate
C     REAL SP_F_CUR     !Spot length factor- current plate 
C     REAL SP_B_CUR     !Spot background border width - current plate
C     REAL SP_D_CUR	!Spot delta mm - current plate
C     REAL SCALE        !Spot expansion factor
C     REAL MM_RAST_X	!Conversion factor for mm to x-rasters
C     REAL MM_RAST_Y	!Conversion factor for mm to y-rasters
C     REAL CTOF		!Crystal to film distance
C     REAL RMAX         !LDM RMAX parameter
C     REAL X		!'x' position on detector in raster unit
C     REAL Y		!'y' position on detector in raster unit
C     REAL ALMBDA	!Reflection lambda value
C     REAL DMINTHR2	!Reflection dmin_thresold2
C     REAL OD()         !Optical density in strip
C     REAL XM 	   	!X offset of ellipse centre from film centre in rasters
C     REAL YM      	!Y offset of ellipse centre from film centre in rasters
C     REAL RADIUS	!Distance of ellipse centre from film centre (rasters)
C     REAL AMAJOR	!Ellipse major axis (corrected) in rasters
C     REAL AMAJOR_1	!Ellipse major axis (including omit pixel) in rasters 
C     REAL AMAJOR_B	!Ellipse major axis (including background) in rasters 
C     REAL BMINOR	!Ellipse minor axis in rasters
C     REAL BMINOR_1	!Ellipse minor axis (including omit pixel) in rasters
C     REAL BMINOR_B	!Ellipse minor axis (including background) in rasters
C     REAL THETA	!angle of ellipse major axis wrt. xf
C     REAL STHETA	!SIN (angle of ellipse major axis wrt. xf)
C     REAL CTHETA	!COS (Angle of ellipse major axis wrt. xf)
C     REAL XPRIME       !x' coordinate of current point wrt ellipse axes 
C                       (rasters)
C     REAL YPRIME       !y' coordinate of current point wrt ellipse axes 
C                       (rasters)
C     REAL ELLIP	!Value of ellipse function x**2/a**+y**2/b**2
C     REAL ELLIP_1	!Value of ellipse function x**2/a**+y**2/b**2 (inc 
C                        omit pixel)
C     REAL ELLIP_B	!Value of ellipse function x**2/a**+y**2/b**2 (inc bg)
C     REAL SIZ          !Maximum (expanded) spot dimension
C     REAL AVP		!Average profile
C     REAL SUM2		!Average squared profile
C     REAL SUMI2()      !Sum of integrated intensities for each bin
C     REAL SIGK		!Standard deviation of Gaussian distribution
C     REAL DELTAM       !Distance beyond which no overlap possible
C     REAL X_OV(MAX_OVLPS) !'x' coordinate of an overlapping spot
C     REAL Y_OV(MAX_OVLPS) !'y' coordinate of an overlapping spot
C     REAL XPN          !'x' pixel position in overlapping spot
C     REAL YPN          !'y' pixel position in overlapping spot
C     CHARACTER*10 TYPSTR     !Parameter as a string 'film' or 'i2'
C     LOGICAL ERR	!ERROR FLAG: .TRUE. ERROR, .FALSE. NO ERROR
C     LOGICAL SPAT	!Flag = .true. spatially overlapped spot, = .false. not
C     LOGICAL CLOS      !Flag = .true. too close to integrate, = .false. not
C     LOGICAL NODFLG	!Flag = .true.  non-spatially overlapped nodal spot 
C                            = .false. not 
C     LOGICAL IOVER	!Flag =.TRUE. overload, = .FALSE. not
C     LOGICAL ROT	!Flag =.TRUE. rotate profile, = .FALSE. not
C     LOGICAL DISTOR    !Distortion corrected coordinates flag
C     LOGICAL INTEG     !Data integrated flag
C
C-------------------------------------------------------------------------------
C
      IERR = 0
      ERRSTR = ' '
      KPASS = IPASS
      IF (IPASS.GT.10) KPASS = 2
      KSPOT = 0
      IF (IPASS.GT.10) KSPOT = IPASS - 10
      ROT = .FALSE.
      KSPOT_LAST = 0
      I_H = 0
      I_K = 0
      I_L = 0
      IF (NTHE.LT.1) NTHE = 1
      IF (NTHE.GT.MAX_PRBIN) THEN
         IERR = 2
         ERRSTR = 'Too many profile bins requested, max ='
         CALL STR_ADDI(ERRSTR,MAX_PRBIN,1)
         RETURN
      END IF
      IF (KPASS.LT.0.OR.KPASS.GT.2) THEN
         IERR = 3
         ERRSTR = 'Invalid pass number flag'
         RETURN
      END IF
C
C====== Check that LRL generated
C
      CALL LRL_GET_PLATE (KPK, KPL)
      IF (IPACK.NE.KPK.OR.IPLATE.NE.KPL) THEN
         IERR = 4
         ERRSTR = 'LRL not generated for requested pack and plate'
         RETURN
      END IF
         
C
C====== Check that profiles calculated if needed
C
      IF (KPASS.GT.1) THEN
         IF (IPROF_PK.EQ.0.OR.IPROF_PL.EQ.0) THEN
            IERR = 5
            ERRSTR = 'No profiles calculated'
            RETURN
         END IF
      END IF
C
C====== Get no. of pack/plate to use; if getting profile fitted intensities
C====== must use same spot parameters as on pack & plate from which profiles
C====== were determined
C
      IF (KPASS.EQ.0) THEN
         IPK = IPACK
         IPL = IPLATE
      ELSE IF (KPASS.EQ.1) THEN
         IPK = IPACK
         IPL = IPLATE
         IPROF_PK = IPACK
         IPROF_PL = IPLATE
         NPROF_BIN = NTHE
         NTH = NTHE
      ELSE IF (KPASS.EQ.2) THEN
         IPK = IPROF_PK
         IPL = IPROF_PL
         NTH = NPROF_BIN
      END IF
C
C====== Determine max box size to use (only need to check first LDM_GET
C====== call with pack/plate parameters
C
      CALL LDM_GET_SPOTL (IPK, IPL, SPOTL, IFLAG, ND, IERR, ERRSTR)
      IF (IERR.NE.0) RETURN
      CALL LDM_GET_SPOTW (IPK, IPL, SPOTW, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_SPOTF (IPK, IPL, SPOTF, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_SPOTB (IPK, IPL, BORDER, IFLAG, ND, IERR,ERRSTR)
      IF (BORDER.LT.1.0) BORDER = 1.0
      CALL LDM_MMRAST (IPK, IPL, MM_RAST_X, MM_RAST_Y)
      CALL LDM_GET_CTOF (IPK, IPL, CTOF, IFLAG, ND, IERR, ERRSTR)
      CALL LDM_GET_RMAX (RMAX, IFLAG, ERRSTR)
      SIZ = SPOTL
      IF (SPOTF.GT.0.0) SIZ = SPOTL*(1.0+SPOTF*RMAX/CTOF)
      IF (SPOTW.GT.SIZ) SIZ = SPOTW
      NAX = NINT(SIZ*MM_RAST_X + 2.0 + 2.0*BORDER + 2.0)
      IF (MOD(NAX,2).EQ.0) NAX = NAX + 1
      N_HALF = NAX/2
      IF (N_HALF.GT.MAX_HPR) THEN
         IERR = 1
         ERRSTR = 'Required box size exceeds program limit of'
         CALL STR_ADDI (ERRSTR,2*MAX_HPR+1,1)
         RETURN
      END IF
      ROT = PROF_ROT
      IF (SPOTW.EQ.0.0) THEN 
         SPOTW = SPOTL
         ROT = .FALSE.
      ENDIF
      IF (KPASS.EQ.1) THEN
         PROF_ROTATE = ROT
      ELSE IF (KPASS.EQ.2) THEN
         ROT = PROF_ROTATE
      END IF 
      IF (KPASS.EQ.1) NPR_HALF = N_HALF
C
C====== Get spot size parameters for current plate
C
      CALL LDM_GET_SPOTL (IPACK, IPLATE, SP_L_CUR, IFLAG, ND, 
     +                    IERR, ERRSTR)
      CALL LDM_GET_SPOTW (IPACK, IPLATE, SP_W_CUR, IFLAG, ND, 
     +                    IERR, ERRSTR)
      CALL LDM_GET_SPOTF (IPACK, IPLATE, SP_F_CUR, IFLAG, ND, 
     +                    IERR, ERRSTR)
      CALL LDM_GET_SPOTB (IPACK, IPLATE, SP_B_CUR, IFLAG, ND, 
     +                    IERR,ERRSTR)
      CALL LDM_GET_SPOTD (IPACK, IPLATE, SP_D_CUR, IFLAG, ND, 
     +                    IERR,ERRSTR)
C
C====== Get overlaps calculation type
C
      CALL LRL_STATUS (NSPOT, ISORT, IOVTYP, DISTOR, INTEG, INTCOD)
C
C====== Determine profile parameters
C
      IF (KPASS.EQ.1) THEN
         CALL LFN_DETPRO (IPK, IPL, N_HALF, IERR, ERRSTR)
         IF (IERR.NE.0) RETURN
      END IF
C
C====== Get some other LDM parameter values
C
      CALL LDM_GET_IMTYP (TYPSTR, IMTYP, IFLAG) 
      CALL LRL_NUMSPOTS (NSPOT) 
      CALL LRL_GET_XYCEN (X_CEN, Y_CEN)
      CALL LDM_GET_OVPIX (OVPIX, IFLAG)
C
C====== Set spot loop limits
C
      IF (KSPOT.EQ.0) THEN
         ISPOT_MIN = 1
         ISPOT_MAX = NSPOT
      ELSE
         IF (KSPOT.LT.1.OR.KSPOT.GT.NSPOT) THEN
            IERR = 6
            ERRSTR = 'Requested spot no. not present in LRL'
            RETURN
         END IF
         ISPOT_MIN = KSPOT
         ISPOT_MAX = KSPOT
         KSPOT_LAST = KSPOT
      END IF
C
C====== (Only need to check error flag for 1'st LDM_GET.. call as
C       check is just on pack/plate number and is the same for all)
C
      IF (IMTYP.EQ.1.AND.ITYPE.NE.1) THEN
         CALL LDM_GET_BASE (BASEOD, IFLAG, ND)
         CALL LDM_GET_NONL (NONL, IFLAG, ND)
         CALL LDM_GET_N1OD (N1OD, IFLAG)
         CFAC=NONL/N1OD
         BASE=BASEOD*N1OD
      END IF
      TWOPI = ATAN2(1.0,1.0)*8.0
      IF (KPASS.EQ.1) THEN
         NUM_IN_SPOT = NUM_IN_PROF
         DO 200 J = 1, NTH
             SUMI2 (J) = 0.0
            DO 100 I = 1, NUM_IN_PROF
               PROF(I, J) = 0.0
               IXP(I) = IXPROF(I)
               IYP(I) = IYPROF(I)
100         CONTINUE
200      CONTINUE
      ENDIF
      IF (IMTYP.EQ.1.AND.ITYPE.EQ.1) THEN
C
C====== Generate look-up table for optical density
C
         CALL LFN_OD_TABLE
         IF (OVPIX.LT.0) OVPIX = 0
         IF (OVPIX.GT.255) OVPIX = 255
         OVPIX = NINT(TABLE(OVPIX))
      ENDIF
C
C====== Read reflection list
C 
      DO 2500 N = ISPOT_MIN, ISPOT_MAX
      CALL LRL_GET (N, IH, IK, IL, X, Y, ALMBDA, DMINTHR2,
     +     MULT, NODFLG, SPAT, JERR)
      IF (SPAT.EQV..TRUE.) THEN
         CALL LRL_GET_OVLP (N, SPAT, CLOS, JERR)
         IF (CLOS) THEN
            IF (KSPOT.GT.0) THEN
               IERR = 10
               ERRSTR = 'Overlapped spot too close to integrate'
               KSPOT_LAST = 0
               RETURN
            ELSE
               CALL LRL_CLR_INT (N)
            END IF
            GO TO 2000
         END IF
      END IF
      ICODE = 0
      BAD = .FALSE.
      IF (KPASS.EQ.1.AND.SPAT.EQV..TRUE.) GOTO 2000
      IXD = NINT(X)
      IYD = NINT(Y)
      IF (IPASS.GT.10) THEN
         I_H = IH
         I_K = IK
         I_L = IL
         XPOS = X
         YPOS = Y
      END IF
C
      CALL LFN_IMG_RECT (IMG, ITYPE, IORDER, NF_OFF, IXD-N_HALF,
     +                   IYD-N_HALF, NAX, NAX, NAX1_RASTS, NAX2_RASTS,
     +                   0, IDATA, IWORK, JERR)
      IF (JERR.EQ.-1) THEN
         IERR = -3
         ERRSTR = 'Invalid image type or order flag'
         KSPOT_LAST = 0
         RETURN
      END IF
      IF (JERR.EQ.-2) THEN
         IERR = -4
         ERRSTR = 'Byte handling not available'
         KSPOT_LAST = 0
         RETURN
      END IF
      IF (JERR.GT.0) THEN
         IF (KSPOT.GT.0) THEN
            IERR = 20
            ERRSTR = 'Integration box out of limits'
            KSPOT_LAST = 0
            RETURN
         ELSE
            GO TO 2000
         END IF
      END IF
      IF (IMTYP.EQ.1) THEN
         IF (ITYPE.EQ.1) THEN
            DO 300 I = 1, NAX*NAX
               OD(I) = TABLE(IDATA(I))
300         CONTINUE
         ELSE
            DO 350 I = 1, NAX*NAX
               
               TRUOD = IDATA(I) - BASE 
               TRUOD = TRUOD*(1.0+CFAC*TRUOD)
               TRUOD = MAX(0.0, TRUOD)
               OD(I) = TRUOD
350         CONTINUE
         END IF
      ELSE
         DO 400 I = 1, NAX*NAX
            OD(I) = IDATA(I)
400      CONTINUE
      ENDIF
C
      XM = X - X_CEN
      YM = Y - Y_CEN
      RADIUS = SQRT(XM**2+YM**2)
      STHETA = YM/RADIUS
      CTHETA = XM/RADIUS
      THETA = AMOD(ATAN2(STHETA,CTHETA)+TWOPI,TWOPI)
      IF (ROT.EQV..FALSE.) THEN
         STHETA = 0.0
         CTHETA = 1.0
      ENDIF
      NBIN = MIN(INT(THETA/(TWOPI/NTH))+1,NTH)
      IF (ROT) THEN
         SCALE = 1.0+(SPOTF*RADIUS)/(CTOF*MM_RAST_X)
      ELSE
         SCALE = 1.0
      END IF  
C
C====== Collect od for profile determination
C
      IF (KPASS.EQ.1) THEN
         DO 700 I = 1, NUM_IN_PROF
            IX = NINT(IXPROF(I)*CTHETA*SCALE - IYPROF(I)*STHETA)
            IY = NINT(IYPROF(I)*CTHETA + IXPROF(I)*STHETA*SCALE)
            PIXV(I) = OD(NAX*(N_HALF + IX) + IY + N_HALF + 1)
            IF (OVPIX.GT.0.AND.NINT(PIXV(I)).GE.OVPIX) GOTO 2000
700      CONTINUE
C
C====== Get spot shape
C
      ELSE
         AMAJOR = SPOTL*MM_RAST_X/2.0*(1.+SPOTF*RADIUS/(CTOF*MM_RAST_X)) 
         AMAJOR_1 = AMAJOR + 1.0
         AMAJOR_B = AMAJOR_1 + BORDER
         BMINOR = SPOTW*MM_RAST_X/2.0
         BMINOR_1 = BMINOR + 1.0
         BMINOR_B = BMINOR_1 + BORDER 
         IOVER = .FALSE.
         NUM_IN_SPOT = 0
         NUM_NOT_BG = 0
         NUM_IN_PK = 0 
         NOVPIX = 0 

         IF(SPAT.EQV..TRUE.) THEN
            CALL LRL_OVALL (IOVTYP, N, NSPOT, 
     +                      SP_L_CUR, SP_W_CUR, SP_F_CUR, 
     +                      SP_B_CUR, SP_D_CUR, RMAX, CTOF, 
     +                      MM_RAST_X, MM_RAST_Y,
     +                      MAX_OVLPS, NUM_OV, X_OV, Y_OV, 
     +                      JERR, ERRSTR)
         END IF

         DO 800 IX = -N_HALF, N_HALF
            DO 750 IY = -N_HALF, N_HALF
               XPRIME = IX*CTHETA + IY*STHETA
               YPRIME = IY*CTHETA - IX*STHETA
               ELLIP = (XPRIME**2/AMAJOR**2) + (YPRIME**2/BMINOR**2)
               ELLIP_1 =(XPRIME**2/AMAJOR_1**2)+(YPRIME**2/BMINOR_1**2)
               ELLIP_B =(XPRIME**2/AMAJOR_B**2)+(YPRIME**2/BMINOR_B**2)
               IF (ELLIP_B.LE.1.0) THEN
                  NUM_IN_SPOT = NUM_IN_SPOT + 1
C
C====== Set 'SHAPE' mask defining background
C
                  SHAPE(NUM_IN_SPOT) = 0
                  PIXV(NUM_IN_SPOT) = OD(NAX*(N_HALF+IX)+IY+N_HALF+1)
                  IXP(NUM_IN_SPOT) = IX 
                  IYP(NUM_IN_SPOT) = IY
C
C====== Get profile for this spot
C
                  IF (KPASS.EQ.2) THEN
       		     NTEMP = INDXPR(NINT(XPRIME/SCALE),NINT(YPRIME))
                     IF (NTEMP.EQ.0) THEN
                        PROF_IND(NUM_IN_SPOT) = 0.0
                     ELSE
                        PROF_IND(NUM_IN_SPOT) = PROF(NTEMP, NBIN)
                     ENDIF
                  ENDIF
C
C====== Set 'SHAPE' mask defining omitted 1pixel border
C
                  IF (ELLIP_1.LE.1.0) THEN
                      SHAPE(NUM_IN_SPOT) = -1
                      NUM_NOT_BG = NUM_NOT_BG + 1
C
C====== Set 'SHAPE' mask defining peak
C
                     IF (ELLIP.LE.1.0) THEN
                        SHAPE(NUM_IN_SPOT) = 10
                        IF (OVPIX.GT.0.AND.NINT(PIXV(NUM_IN_SPOT))
     +                     .GE.OVPIX) THEN
                           NOVPIX = NOVPIX + 1
C
C====== Set 'SHAPE' mask defining saturated pixels
C
                           SHAPE(NUM_IN_SPOT) = -3
                        ELSE
                           NUM_IN_PK = NUM_IN_PK + 1
                        ENDIF
                     ENDIF
                  ENDIF
C
C====== Set 'SHAPE' mask defining omitted overlapping pixels
C
                  IF (SPAT.EQV..TRUE..AND.NUM_OV.GT.0) THEN
                     DO 720 I = 1, NUM_OV
                        XPN = FLOAT(IX) + X - X_OV(I)    
                        YPN = FLOAT(IY) + Y - Y_OV(I)
                        XPRIME = XPN*CTHETA + YPN*STHETA
                        YPRIME = YPN*CTHETA - XPN*STHETA
                        ELLIP = (XPRIME**2/AMAJOR**2) + 
     +                          (YPRIME**2/BMINOR**2)
                        IF (ELLIP.LE.1.0) SHAPE(NUM_IN_SPOT) = -2
720                  CONTINUE
                  END IF

               ENDIF
750         CONTINUE
800      CONTINUE
         IF (NOVPIX.GT.NOVLD) IOVER =.TRUE.
         NUM_IN_BG = NUM_IN_SPOT - NUM_NOT_BG
         NPK_MIN = (NUM_IN_PK-3)/5
         NBG_MIN = (NUM_IN_BG-3)/5 
         NPK_OK = 0
         NBG_OK = 0
         DO 850 I = 1, NUM_IN_SPOT
            IF (SHAPE(I).GT.0) NPK_OK = NPK_OK + 1
            IF (SHAPE(I).EQ.0) NBG_OK = NBG_OK + 1
 850     CONTINUE
         IF (NPK_OK.LT.1) ICODE = ICODE + 1
         IF (NBG_OK.LT.3) ICODE = ICODE + 4
         IF (ICODE.NE.0) THEN
            BAD = .TRUE.
            IF (KSPOT.GT.0) THEN
               IF (ICODE.EQ.1) THEN
                  IERR = 11
                  ERRSTR='Less than 3 non-overlapped background pixels'      
               ELSE 
                  IERR = 12
                  ERRSTR = 
     +               'No non-overlapped/non-overloaded peak pixels'
               END IF
               KSPOT_LAST = 0
               RETURN
            ELSE
               INTEN = 0.0
               SIGI = 0.0  
               IF (KSPOT.EQ.0) CALL LRL_SET_INT (N, INTEN, SIGI, 
     +                              .TRUE., .FALSE., ICODE, JERR)
               IF (KSPOT.EQ.0) CALL LRL_SET_OVPIX (N, NOVPIX, JERR)
            END IF
         END IF
         IF (NPK_OK.LT.NPK_MIN) THEN
            ICODE = ICODE + 2
            BAD = .TRUE.
         END IF
         IF (NBG_OK.LT.NBG_MIN) THEN
            ICODE = ICODE + 8
            BAD = .TRUE.
         END IF
      ENDIF
C
      IF (KPASS.EQ.0) THEN
         AVP = 0.0
         SUM2 = 0.0
         DO 1200 I = 1, NUM_IN_SPOT
            AVP = AVP + PIXV(I)
            SUM2 = SUM2 + PIXV(I)**2
1200     CONTINUE
C
C====== Compute Gaussian standard deviation
C
         AVP = AVP / FLOAT(NUM_IN_SPOT)
         SUM2 = SUM2 / FLOAT(NUM_IN_SPOT)
         SIGK = SQRT(SUM2 - AVP**2)
C
         DO 1300 I = 1, NUM_IN_SPOT
            IF ((PIXV(I)-AVP).GT.SIGK) SHAPE(I) = SHAPE(I) + 1
1300     CONTINUE
      ENDIF
C
C====== Determine background plane: Pass 1
C
      CALL LFN_PLFIT(1)
C
C====== Subtract background and eliminate outlying pixels by correcting LS sums
C
      CALL LFN_PLFSUB(1)
C
C====== Determine background: Pass 2 (outlying pixels have been rejected)
C
      CALL LFN_PLFIT(2)
C
C====== Subtract background and return integrated intensity and standard 
C       deviation
C
      IF (KPASS.EQ.2) THEN
         CALL LFN_PLFSUB(3)
      ELSE
         CALL LFN_PLFSUB(2)
      ENDIF
      IF (KPASS.EQ.1) THEN
         IF (INTEN.LT.PROMIN) GOTO 2000
         DO 1800 I = 1, NUM_IN_PROF
            PROF(I,NBIN) = PROF(I,NBIN) + INTEN*SUBPR(I)
1800     CONTINUE
         SUMI2(NBIN) = SUMI2(NBIN) + INTEN*INTEN
      ELSE
         IF (KSPOT.EQ.0.AND.SIGI.GT.0.0) 
     +      CALL LRL_SET_INT (N, INTEN, SIGI, BAD, IOVER, ICODE, JERR)
         IF (KSPOT.EQ.0.AND.SIGI.GT.0.0) 
     +      CALL LRL_SET_OVPIX (N, NOVPIX, JERR)
      ENDIF
2000  CONTINUE
      CALL PROGRESS (KPASS, NSPOT, N)
      IF (CANCEL(0)) THEN
         IERR = 100
         ERRSTR = 'Integration interrupted'
         KSPOT_LAST = 0
         RETURN
      END IF
2500  CONTINUE
      IF (KPASS.EQ.1) THEN
         DO 2200 J = 1, NTH
            DO 2100 I = 1, NUM_IN_PROF
               PROF(I,J) = PROF(I,J)/SUMI2(J)
               IF (PROF(I,J).LT.0.0000) PROF(I,J) = 0.0000
2100        CONTINUE
2200     CONTINUE
      ELSE
         IF (KSPOT.EQ.0) CALL LRL_SET_IDONE (1) 
      ENDIF
      RETURN
      END
