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