SUBROUTINE EVPOLAR( STATUS )
*+
* Name:
* EVPOLAR
* Purpose:
* Bins 2 lists into a 1 or 2 dimensional object
* Language:
* Starlink Fortran
* Type of Module:
* ASTERIX task
* Invocation:
* CALL EVPOLAR( STATUS )
* Arguments:
* STATUS = INTEGER (Given and Returned)
* The global status.
* Description:
* A binned polar dataset is produced from 2 lists from the input dataset.
*
* The azimuthal axis is always regularly spaced, and is specified by
* the number of bins. The radial axis may be regularly or irregularly
* spaced. If regular, bin width and number of bins are specified : if
* irregular, the bin bounds.
*
* All bins are INCLUSIVE of their lower bound and EXCLUSIVE of their
* upper bound.
*
* The way QUALITY lists are handled is controlled by 2 parameters QVAL,
* and QKEEP. Values present in the quality list > QVAL are treated as bad
* quality values. If QKEEP is true, then bad events are written to the
* output DATA_ARRAY, and the corresponding element of the output QUALITY
* array is set to bad (i.e.1). If QKEEP is false, then all bad events are
* simply ignored, and no output QUALITY array is produced.
* Usage:
* evpolar {parameter_usage}
* Environment Parameters:
* INP = CHAR (read)
* Name of input EVDS
* LISTS = CHAR (read)
* Index No(s) of input list(s)
* QVAL = INTEGER (read)
* Event quality > QVAL = bad
* QKEEP = LOGICAL (read)
* Produce output QUALITY array?
* REG = LOGICAL (read)
* Regular or irregular radial bins
* ABINSIZE = REAL (read)
* Azimuthal bin size
* RBINSIZE = REAL (read)
* Regular radial bin width
* RRANGE = CHAR (read)
* Irregular bin bounds
* NAZ, NRAD = INTEGER (read)
* Numbers of bins
* NORMALISE = LOGICAL (read)
* Normalise output array
* OUT = CHAR (read)
* Name of output dataset
* Examples:
* {routine_example_text}
* {routine_example_description}
* Pitfalls:
* {pitfall_description}...
* Notes:
* {routine_notes}...
* Prior Requirements:
* {routine_prior_requirements}...
* Side Effects:
* {routine_side_effects}...
* Algorithm:
* {algorithm_description}...
* Accuracy:
* {routine_accuracy}
* Timing:
* {routine_timing}
* Implementation Status:
* {routine_implementation_status}
* External Routines Used:
* {name_of_facility_or_package}:
* {routine_used}...
* Implementation Deficiencies:
* {routine_deficiencies}...
* References:
* {task_references}...
* Keywords:
* evpolar, usage:public
* Copyright:
* Copyright (C) University of Birmingham, 1995
* Authors:
* DJA: David J. Allan (Jet-X, University of Birmingham)
* {enter_new_authors_here}
* History:
* 21 Feb 1991 V1.4-0 (DJA):
* Original version.
* 11 Sep 1993 V1.7-0 (DJA):
* EVPOLAR_CAREA renamed GEO_CAREA
* 25 Feb 1994 V1.7-1 (DJA):
* Use BIT_ routines to do bit manipulations
* 24 Nov 1994 V1.8-0 (DJA):
* Now use USI for user interface
* 21 Apr 1995 V1.8-1 (DJA):
* Updated data interface
* 18 Aug 1995 V2.0-0 (DJA):
* ADI port of event handling
* 21 Feb 1996 V2.0-1 (DJA):
* Removed ATAN2D for Linux port
* {enter_changes_here}
* Bugs:
* {note_any_bugs_here}
*-
* Type Definitions:
IMPLICIT NONE ! No implicit typing
* Global Constants:
INCLUDE 'SAE_PAR' ! Standard SAE constants
INCLUDE 'ADI_PAR'
INCLUDE 'QUAL_PAR'
* Status:
INTEGER STATUS ! Global status
* External References:
EXTERNAL CHR_LEN
INTEGER CHR_LEN
EXTERNAL CHR_SIMLR
LOGICAL CHR_SIMLR
* Local Constants:
INTEGER MXRANGE ! 2 x max No of irregular bins
PARAMETER (MXRANGE = 1000)
INTEGER MXTEXT
PARAMETER (MXTEXT = 7)
CHARACTER*30 VERSION
PARAMETER ( VERSION = 'EVPOLAR Version V2.0-1' )
* Local Variables:
CHARACTER*20 BNAME(2) ! List names
CHARACTER*40 BUNIT(2) ! List units
CHARACTER*80 HTEXT(MXTEXT) ! History text
CHARACTER*40 ORUNIT ! Output radial unit
CHARACTER*80 TEXT ! Temporary text string
REAL ABINSZ ! Azimuthal bin size
REAL LOX,HIX,LOY,HIY ! Field min and max values
REAL MAXR ! Maximum radial distance
REAL RANGE(MXRANGE) ! Irregular bin ranges.
REAL RBINSZ ! Radial bin size
REAL SPARR(2) ! Spaced array data
REAL X0, Y0 ! Centre of polar distribution
INTEGER BADQUAL ! Exclude events with quality
! > this value.
INTEGER BLID(2) ! List identifiers
INTEGER BPTR(2) ! Binning list pointers
INTEGER I ! Loop counter
INTEGER IFID ! Input dataset id
INTEGER INDEX(2) ! Selected lists
INTEGER INLIST ! Number of lists in input EVDS
INTEGER IQPTR ! Input quality data
INTEGER LEN ! Length of TEMP
INTEGER NAZ, NRAD ! # azimuthal/radial bins
INTEGER NDIM ! Output dimensionality
INTEGER NEV ! Number of events in input
INTEGER NINDEX ! Number of lists selected.
INTEGER NLINE ! # history lines used
INTEGER NREJ ! # events not binned
INTEGER OAPTR, OWPTR ! Output irregular axis data
INTEGER ODIMS(2) ! Output dataset dimensions
INTEGER ODPTR ! Output data pointer
INTEGER OFID ! Output dataset id
INTEGER ONELM ! # output bins
INTEGER OQPTR ! Output quality pointer
INTEGER WPTR ! Pointer to workspace array
LOGICAL NORMALISE ! Normalise output bins?
LOGICAL QKEEP ! If true produce output
! quality array.
LOGICAL QUALITY ! Is QUALITY list present?
LOGICAL REG ! Allow irregular binning?
*.
* Check inherited global status.
IF ( STATUS .NE. SAI__OK ) RETURN
* Version id
CALL MSG_PRNT( VERSION )
* Initialise ASTERIX
CALL AST_INIT()
* Get event dataset
CALL USI_ASSOC( 'INP', 'EventDS', 'READ', IFID, STATUS )
IF ( STATUS .NE. SAI__OK ) GOTO 99
* Find all valid lists in INP. Display them to user.
CALL MSG_BLNK()
CALL MSG_PRNT( 'The available lists are :' )
CALL MSG_BLNK()
CALL EDI_DISP( IFID, STATUS )
IF ( STATUS .NE. SAI__OK ) GOTO 99
* Get number of events and lists
CALL EDI_GETNS( IFID, NEV, INLIST, STATUS )
* Select the lists
CALL MSG_BLNK()
CALL MSG_PRNT( 'Select two lists to bin using index'/
: /' numbers above, eg. 1 2. The direction' )
CALL MSG_PRNT( 'of zero azimuth is that of increasing list'/
: /' 1 value.')
CALL MSG_BLNK()
CALL EDI_SELCT( 'LISTS', INLIST, 2, 2, INDEX, NINDEX, STATUS )
IF ( STATUS .NE. SAI__OK ) GOTO 99
* Is there a quality list?
CALL EDI_CHK( IFID, 'QUALITY', QUALITY, STATUS )
IF ( QUALITY ) THEN
* Map it
CALL EDI_MAPI( IFID, 'QUALITY', 'READ', 0, 0, IQPTR, STATUS )
* Get quality processing mode
CALL USI_GET0I( 'QVAL', BADQUAL, STATUS )
CALL USI_GET0L( 'QKEEP', QKEEP, STATUS )
END IF
IF ( STATUS .NE. SAI__OK ) GOTO 99
* Get info on selected lists
DO I = 1, 2
* Make sure user didn't select quality
C IF ( QUALITY ) THEN
C CALL HDX_SAME( LLOC(INDEX(I)), QLOC, SAME, STATUS )
C IF ( SAME ) THEN
C STATUS = SAI__ERROR
C CALL ERR_REP( ' ', 'ERROR : Cannot bin using QUALITY list',
C : STATUS )
C GOTO 99
C END IF
C END IF
* Locate the list
CALL EDI_IDX( IFID, INDEX(I), BLID(I), STATUS )
* Get its name
CALL ADI_CGET0C( BLID(I), 'Name', BNAME(I), STATUS )
* Map it
CALL EDI_MAPR( IFID, BNAME(I), 'READ', 0, 0, BPTR(I), STATUS )
* Get units
CALL ADI_CGET0C( BLID(I), 'Units', BUNIT(I), STATUS )
IF ( STATUS .NE. SAI__OK ) THEN
CALL ERR_ANNUL( STATUS )
BUNIT(I) = ' '
END IF
END DO
* Get field extrema
CALL ADI_CGET0R( BLID(1), 'Min', LOX, STATUS )
CALL ADI_CGET0R( BLID(1), 'Max', HIX, STATUS )
CALL ADI_CGET0R( BLID(2), 'Min', LOY, STATUS )
CALL ADI_CGET0R( BLID(2), 'Max', HIY, STATUS )
* Decide on radial units
ORUNIT = BUNIT(1)
IF ( .NOT. CHR_SIMLR( BUNIT(1), BUNIT(2) ) ) THEN
CALL MSG_PRNT( 'WARNING : List units are different, '/
: /'output may be rubbish' )
END IF
* Announce axis ranges
CALL MSG_SETC( 'AX', BNAME(1) )
CALL MSG_SETR( 'LO', LOX )
CALL MSG_SETR( 'HI', HIX )
CALL MSG_PRNT( 'Axis ^AX ranges from ^LO to ^HI' )
CALL MSG_SETC( 'AX', BNAME(2) )
CALL MSG_SETR( 'LO', LOY )
CALL MSG_SETR( 'HI', HIY )
CALL MSG_PRNT( 'Axis ^AX ranges from ^LO to ^HI' )
* Get polar centre
CALL USI_GET0R( 'X0', X0, STATUS )
CALL USI_GET0R( 'Y0', Y0, STATUS )
IF ( STATUS .NE. SAI__OK ) GOTO 99
* Map workspace
CALL DYN_MAPR( 1, 2*NEV, WPTR, STATUS )
IF ( STATUS .NE. SAI__OK ) GOTO 99
* Find azimuth and radius for each bin
CALL EVPOLAR_RECPOL( NEV, X0, Y0, %VAL(BPTR(1)), %VAL(BPTR(2)),
: %VAL(WPTR), MAXR, STATUS )
* Get azimuthal bin size
CALL USI_GET0R( 'ABINSIZE', ABINSZ, STATUS )
IF ( STATUS .NE. SAI__OK ) GOTO 99
NAZ = 360.0 / ABINSZ
CALL MSG_SETI( 'NAZ', NAZ )
CALL MSG_PRNT( 'There will be ^NAZ azimuthal bins in the output' )
* Inform user of maximum radius, and get radial bins
CALL USI_GET0L( 'REG', REG, STATUS )
IF ( STATUS .NE. SAI__OK ) GOTO 99
CALL MSG_SETR( 'RMAX', MAXR )
CALL MSG_SETC( 'RUNIT', ORUNIT )
CALL MSG_PRNT( 'Radial distance varies from 0 to ^RMAX ^RUNIT' )
IF ( .NOT. REG ) THEN
CALL MSG_PRNT( 'You must specify INCREASING ranges e.g.'/
: /' 0:10:20' )
CALL PRS_GETRANGES( 'RRANGE', MXRANGE, 1, 0.0, MAXR,
: RANGE, NRAD, STATUS )
ELSE
* Get binsize
CALL USI_GET0R( 'RBINSIZE', RBINSZ, STATUS )
IF ( STATUS .NE. SAI__OK ) GOTO 99
IF ( RBINSZ .LE. 0.0 ) THEN
CALL MSG_PRNT( 'ERROR : Zero or negative radial bin width' )
STATUS = SAI__ERROR
END IF
* Get number of radial bins
CALL USI_DEF0I( 'NRAD', MAX(1,NINT(MAXR/RBINSZ)), STATUS )
CALL USI_GET0I( 'NRAD', NRAD, STATUS )
END IF
IF ( STATUS .NE. SAI__OK ) GOTO 99
* Decide output dimensionality
IF ( ( NAZ .EQ. 1 ) .OR. ( NRAD .EQ. 1 ) ) THEN
NDIM = 1
ODIMS(1) = NRAD
ELSE
NDIM = 2
ODIMS(1) = NRAD
ODIMS(2) = NAZ
END IF
* Create output dataset
CALL USI_CREAT( 'OUT', ADI__NULLID, OFID, STATUS )
CALL BDI_LINK( 'BinDS', NDIM, ODIMS, 'REAL', OFID, STATUS )
CALL BDI_MAPR( OFID, 'Data', 'WRITE/ZERO', ODPTR, STATUS )
* Create AXIS structure
IF ( REG ) THEN
SPARR(1) = 0.0
SPARR(2) = RBINSZ
CALL BDI_AXPUT1R( OFID, 1, 'SpacedData', 2, SPARR, STATUS )
CALL BDI_AXPUT0R( OFID, 1, 'ScalarWidth', RBINSZ, STATUS )
ELSE
CALL BDI_AXMAPR( OFID, 1, 'Data', 'WRITE', OAPTR, STATUS )
CALL BDI_AXMAPR( OFID, 1, 'Width', 'WRITE', OWPTR, STATUS )
CALL ARR_BND2CWR( NRAD, RANGE, %VAL(OAPTR), %VAL(OWPTR),
: STATUS )
END IF
CALL BDI_AXPUT0C( OFID, 1, 'Label', 'Radius', STATUS )
IF ( ORUNIT .GT. ' ' ) THEN
CALL BDI_AXPUT0C( OFID, 1, 'Units', ORUNIT, STATUS )
END IF
IF ( NDIM .EQ. 2 ) THEN
SPARR(1) = 0.5*ABINSZ
SPARR(2) = ABINSZ
CALL BDI_AXPUT1R( OFID, 2, 'SpacedData', 2, SPARR, STATUS )
CALL BDI_AXPUT0R( OFID, 2, 'ScalarWidth', ABINSZ, STATUS )
CALL BDI_AXPUT0C( OFID, 2, 'Label', 'Azimuth', STATUS )
CALL BDI_AXPUT0C( OFID, 2, 'Units', 'degree', STATUS )
END IF
IF ( STATUS .NE. SAI__OK ) GOTO 99
* Initialise quality to data missing
CALL ARR_SUMDIM( NDIM, ODIMS, ONELM )
IF ( QUALITY .AND. QKEEP ) THEN
CALL BDI_MAPUB( OFID, 'Quality', 'WRITE/QMISSING', OQPTR,
: STATUS )
CALL BDI_PUT0UB( OFID, 'QualityMask', QUAL__MASK, STATUS )
END IF
* Bin the events
IF ( REG ) THEN
CALL EVPOLAR_REGBIN( NAZ, NRAD, RBINSZ, NEV, %VAL(WPTR),
: QUALITY, QKEEP, %VAL(IQPTR), BADQUAL,
: %VAL(ODPTR), %VAL(OQPTR), NREJ, STATUS )
ELSE
CALL EVPOLAR_IRREGBIN( NAZ, NRAD, RANGE, NEV,
: %VAL(WPTR), QUALITY, QKEEP, %VAL(IQPTR),
: BADQUAL, %VAL(ODPTR), %VAL(OQPTR), NREJ,
: STATUS )
END IF
* Inform of events binned
CALL MSG_SETI( 'BIN', NEV - NREJ )
CALL MSG_SETI( 'NIN', NEV )
CALL MSG_PRNT( '^BIN events were binned out of ^NIN input.' )
* Normalise
NORMALISE = .FALSE.
IF ( NINDEX .EQ. 2 ) THEN
CALL USI_GET0L( 'NORM', NORMALISE, STATUS )
END IF
IF ( NORMALISE .AND. ( STATUS .EQ. SAI__OK ) ) THEN
IF ( REG ) THEN
CALL EVPOLAR_REGNORM( NAZ, NRAD, RBINSZ, X0, Y0, LOX, HIX,
: LOY, HIY, %VAL(ODPTR), STATUS )
ELSE
CALL EVPOLAR_IRREGNORM( NAZ, NRAD, RANGE, X0, Y0, LOX, HIX,
: LOY, HIY, %VAL(ODPTR), STATUS )
END IF
END IF
CALL BDI_AXPUT0L( OFID, 1, 'Normalised', NORMALISE, STATUS )
IF ( NDIM .EQ. 2 ) THEN
CALL BDI_AXPUT0L( OFID, 2, 'Normalised', NORMALISE, STATUS )
END IF
* Copy ancillary stuff
CALL BDI_COPY( IFID, 'Title,Label,Units', OFID, ' ', STATUS )
CALL UDI_COPANC( IFID, 'grf', OFID, STATUS )
IF ( STATUS .NE. SAI__OK ) GOTO 99
* Re-write data units
IF ( NORMALISE ) THEN
CALL MSG_SETC( 'UN', ORUNIT )
CALL MSG_MAKE( 'count/^UN**2', TEXT, LEN )
CALL BDI_PUT0C( OFID, 'Units', TEXT(:LEN), STATUS )
ELSE
CALL BDI_PUT0C( OFID, 'Units', 'count', STATUS )
END IF
* Update history
CALL HSI_COPY( IFID, OFID, STATUS )
CALL HSI_ADD( OFID, VERSION, STATUS )
NLINE = MXTEXT
HTEXT(1) = 'Input {INP}'
HTEXT(2) = 'Binned dataset created from the LISTs : '//BNAME(1)
HTEXT(3) = ' '//BNAME(2)
WRITE (HTEXT(4), '(12X,A,I,A)') 'Containing data on ',
: NEV, ' events.'
IF ( QUALITY ) THEN
IF ( QKEEP ) THEN
HTEXT(5) = 'QUALITY array created'
ELSE
HTEXT(5) = 'Bad quality events excluded.'
END IF
CALL USI_TEXT( 5, HTEXT, NLINE, STATUS )
ELSE
CALL USI_TEXT( 4, HTEXT, NLINE, STATUS )
END IF
* Write expanded text
CALL HSI_PTXT( OFID, NLINE, HTEXT, STATUS )
* Free the selected lists
DO I = 1, 2
CALL ADI_ERASE( BLID(I), STATUS )
CALL EDI_UNMAP( IFID, BNAME(I), STATUS )
END DO
* Tidy up
99 CALL AST_CLOSE
CALL AST_ERR( STATUS )
END
*+ EVPOLAR_RECPOL - Perform rectangular to polar conversion
SUBROUTINE EVPOLAR_RECPOL( NEV, X0, Y0, XV, YV, POL,MAXR, STATUS )
:
* Description :
*
* Method :
*
* Authors :
*
* David J. Allan (BHVAD::DJA)
*
* History :
*
* 21 Feb 91 : Original
*
* Type definitions :
*
IMPLICIT NONE
*
* Global constants :
*
INCLUDE 'SAE_PAR'
INCLUDE 'MATH_PAR'
*
* Import :
*
INTEGER NEV ! Number of events
REAL X0, Y0 ! Position of pole
REAL XV(*), YV(*) ! Rectangular coords
*
* Export :
*
REAL POL(2,NEV) ! Polar coordinates
REAL MAXR ! Maximum radius
*
* Status :
*
INTEGER STATUS
*
* Local variables :
*
REAL DX, DY ! Offsets from pole
INTEGER I ! Loop over events
*-
IF ( STATUS .EQ. SAI__OK ) THEN
MAXR = -1.0
DO I = 1, NEV
DX = XV(I) - X0
DY = YV(I) - Y0
POL(1,I) = SQRT(DX*DX+DY*DY)
POL(2,I) = ATAN2(DX,DY)*MATH__RTOD
IF ( POL(2,I) .LT. 0 ) POL(2,I) = POL(2,I) + 360.0
IF ( POL(1,I) .GT. MAXR ) MAXR = POL(1,I)
END DO
END IF
END
*+ EVPOLAR_REGBIN - Regular bin polar binning
SUBROUTINE EVPOLAR_REGBIN( NAZ, NRAD, RBINSZ, NEV, POL,
: GOTQUAL, QKEEP, IQUAL, BADQUAL, DATA, QUAL,
: NREJ, STATUS )
*
* Description :
*
* Method :
*
* FOR each event
* Get radial bin
* IF valid
* if 2D get azimuthal bin
* IF input_quality present and good
* data() = data() + 1
* if qual() not already bad, set good
* ELSE if create output quality
* qual() = bad
* ELSE
* nreject = nreject + 1
*
* Authors :
*
* David J. Allan (BHVAD::DJA)
*
* History :
*
* 21 Feb 91 : Original
*
* Type definitions :
*
IMPLICIT NONE
*
* Global constants :
*
INCLUDE 'SAE_PAR'
INCLUDE 'QUAL_PAR'
*
* Import :
*
INTEGER NAZ, NRAD ! # azimuthal,radial bins
REAL RBINSZ ! Radial bin size
INTEGER NEV ! Number of events
REAL POL(2,NEV) ! Polar coordinates
LOGICAL GOTQUAL ! Got input quality list?
LOGICAL QKEEP ! Create output quality?
INTEGER IQUAL(NEV) ! Input quality values
INTEGER BADQUAL ! Quality threshold
*
* Export :
*
REAL DATA(NRAD,NAZ) ! Output data
BYTE QUAL(NRAD,NAZ) ! Output data
INTEGER NREJ ! Events not binned
*
* Status :
*
INTEGER STATUS
*
* Local variables :
*
REAL ABINSZ ! Azimuthal binsize
INTEGER ABIN ! Azimuthal bin
INTEGER I ! Loop over events
INTEGER RBIN ! Radial bin
*-
IF ( STATUS .EQ. SAI__OK ) THEN
NREJ = 0
IF ( NAZ .EQ. 1 ) THEN
IF ( GOTQUAL ) THEN
DO I = 1, NEV
RBIN = INT( POL(1,I) / RBINSZ ) + 1
IF ( RBIN .GT. NRAD ) THEN
NREJ = NREJ + 1
ELSE
IF ( IQUAL(I) .LE. BADQUAL ) THEN
DATA(RBIN,1) = DATA(RBIN,1) + 1.0
IF ( QUAL(RBIN,1) .NE. QUAL__BAD )
: QUAL(RBIN,1) = QUAL__GOOD
ELSE IF ( QKEEP ) THEN
DATA(RBIN,1) = DATA(RBIN,1) + 1.0
QUAL(RBIN,1) = QUAL__BAD
END IF
END IF
END DO
ELSE
DO I = 1, NEV
RBIN = INT( POL(1,I) / RBINSZ ) + 1
IF ( RBIN .GT. NRAD ) THEN
NREJ = NREJ + 1
ELSE
DATA(RBIN,1) = DATA(RBIN,1) + 1.0
END IF
END DO
END IF
ELSE
* Need radial and azimuthal bins. Azimuth bin value is always valid
ABINSZ = 360.0 / REAL(NAZ)
IF ( GOTQUAL ) THEN
DO I = 1, NEV
RBIN = INT( POL(1,I) / RBINSZ ) + 1
IF ( RBIN .GT. NRAD ) THEN
NREJ = NREJ + 1
ELSE
ABIN = INT( POL(2,I)/ ABINSZ ) + 1
IF ( IQUAL(I) .LE. BADQUAL ) THEN
DATA(RBIN,ABIN) = DATA(RBIN,ABIN) + 1.0
IF ( QUAL(RBIN,ABIN) .NE. QUAL__BAD )
: QUAL(RBIN,ABIN) = QUAL__GOOD
ELSE IF ( QKEEP ) THEN
DATA(RBIN,ABIN) = DATA(RBIN,ABIN) + 1.0
QUAL(RBIN,ABIN) = QUAL__BAD
END IF
END IF
END DO
ELSE
DO I = 1, NEV
RBIN = INT( POL(1,I) / RBINSZ ) + 1
IF ( RBIN .GT. NRAD ) THEN
NREJ = NREJ + 1
ELSE
ABIN = INT( POL(2,I)/ ABINSZ ) + 1
DATA(RBIN,ABIN) = DATA(RBIN,ABIN) + 1.0
END IF
END DO
END IF
END IF
END IF
END
*+ EVPOLAR_IRREGBIN - Irregular polar binning
SUBROUTINE EVPOLAR_IRREGBIN( NAZ, NRAD, RANGE, NEV, POL,
: GOTQUAL, QKEEP, IQUAL, BADQUAL, DATA, QUAL,
: NREJ, STATUS )
*
* Description :
*
* Method :
*
* FOR each event
* Get radial bin
* IF valid
* if 2D get azimuthal bin
* IF input_quality present and good
* data() = data() + 1
* if qual() not already bad, set good
* ELSE if create output quality
* qual() = bad
* ELSE
* nreject = nreject + 1
*
* Authors :
*
* David J. Allan (BHVAD::DJA)
*
* History :
*
* 21 Feb 91 : Original
*
* Type definitions :
*
IMPLICIT NONE
*
* Global constants :
*
INCLUDE 'SAE_PAR'
INCLUDE 'QUAL_PAR'
*
* Import :
*
INTEGER NAZ, NRAD ! # azimuthal,radial bins
REAL RANGE(*) ! Bin boundaries
INTEGER NEV ! Number of events
REAL POL(2,NEV) ! Polar coordinates
LOGICAL GOTQUAL ! Got input quality list?
LOGICAL QKEEP ! Create output quality?
INTEGER IQUAL(NEV) ! Input quality values
INTEGER BADQUAL ! Quality threshold
*
* Export :
*
REAL DATA(NRAD,NAZ) ! Output data
BYTE QUAL(NRAD,NAZ) ! Output data
INTEGER NREJ ! Events not binned
*
* Status :
*
INTEGER STATUS
*
* Local variables :
*
REAL ABINSZ ! Azimuthal binsize
INTEGER ABIN ! Azimuthal bin
INTEGER I ! Loop over events
INTEGER RBIN ! Radial bin
LOGICAL OK ! Valid radial bin?
*-
IF ( STATUS .EQ. SAI__OK ) THEN
NREJ = 0
IF ( NAZ .EQ. 1 ) THEN
IF ( GOTQUAL ) THEN
DO I = 1, NEV
CALL AXIS_RNGIDX( NRAD, RANGE,.FALSE.,POL(1,I), RBIN, OK )
IF ( OK ) THEN
IF ( IQUAL(I) .LE. BADQUAL ) THEN
DATA(RBIN,1) = DATA(RBIN,1) + 1.0
IF ( QUAL(RBIN,1) .NE. QUAL__BAD )
: QUAL(RBIN,1) = QUAL__GOOD
ELSE IF ( QKEEP ) THEN
DATA(RBIN,1) = DATA(RBIN,1) + 1.0
QUAL(RBIN,1) = QUAL__BAD
END IF
ELSE
NREJ = NREJ + 1
END IF
END DO
ELSE
DO I = 1, NEV
CALL AXIS_RNGIDX( NRAD, RANGE,.FALSE.,POL(1,I), RBIN, OK )
IF ( OK ) THEN
DATA(RBIN,1) = DATA(RBIN,1) + 1.0
ELSE
NREJ = NREJ + 1
END IF
END DO
END IF
ELSE
* Need radial and azimuthal bins. Azimuth bin value is always valid
ABINSZ = 360.0 / REAL(NAZ)
IF ( GOTQUAL ) THEN
DO I = 1, NEV
CALL AXIS_RNGIDX( NRAD, RANGE,.FALSE.,POL(1,I), RBIN, OK )
IF ( OK ) THEN
ABIN = INT( POL(2,I)/ ABINSZ ) + 1
IF ( IQUAL(I) .LE. BADQUAL ) THEN
DATA(RBIN,ABIN) = DATA(RBIN,ABIN) + 1.0
IF ( QUAL(RBIN,ABIN) .NE. QUAL__BAD )
: QUAL(RBIN,ABIN) = QUAL__GOOD
ELSE IF ( QKEEP ) THEN
DATA(RBIN,ABIN) = DATA(RBIN,ABIN) + 1.0
QUAL(RBIN,ABIN) = QUAL__BAD
END IF
ELSE
NREJ = NREJ + 1
END IF
END DO
ELSE
DO I = 1, NEV
CALL AXIS_RNGIDX( NRAD, RANGE,.FALSE.,POL(1,I), RBIN, OK )
IF ( OK ) THEN
ABIN = INT( POL(2,I)/ ABINSZ ) + 1
DATA(RBIN,ABIN) = DATA(RBIN,ABIN) + 1.0
ELSE
NREJ = NREJ + 1
END IF
END DO
END IF
END IF
END IF
END
*+ EVPOLAR_REGNORM - Normalise data array
SUBROUTINE EVPOLAR_REGNORM( NAZ, NRAD, RBINSZ, CX, CY, LOX,
: HIX, LOY, HIY, DATA, STATUS )
*
* Description :
*
* Method :
*
* Authors :
*
* David J. Allan (BHVAD::DJA)
*
* History :
*
* 21 Feb 91 : Original
*
* Type definitions :
*
IMPLICIT NONE
*
* Global constants :
*
INCLUDE 'SAE_PAR'
*
* Import :
*
INTEGER NAZ, NRAD ! # azimuthal,radial bins
REAL RBINSZ ! Radial bin size
REAL CX,CY ! Centre of circle
REAL LOX, HIX ! Range of rectangle in X
REAL LOY, HIY ! Range of rectangle in Y
*
* Import / export :
*
REAL DATA(NRAD,NAZ) ! Output data
*
* Status :
*
INTEGER STATUS
*
* Functions :
*
REAL GEO_CAREA
*
* Local variables :
*
REAL FACTOR ! Normalisation factor
REAL LAST_AREA ! Last circle area
REAL RAD ! Current radius
REAL THIS_AREA ! Current circle area
INTEGER I, J ! Loop over bins
*-
IF ( STATUS .EQ. SAI__OK ) THEN
RAD = RBINSZ
LAST_AREA = 0.0
DO I = 1, NRAD
THIS_AREA = GEO_CAREA(RAD,CX,CY,LOX,HIX,LOY,HIY)
FACTOR = ( THIS_AREA - LAST_AREA ) / REAL(NAZ)
DO J = 1, NAZ
DATA(I,J) = DATA(I,J) / FACTOR
END DO
RAD = RAD + RBINSZ
LAST_AREA = THIS_AREA
END DO
END IF
END
*+ EVPOLAR_IRREGNORM - Normalise irregularly data array
SUBROUTINE EVPOLAR_IRREGNORM( NAZ, NRAD, RANGE, CX, CY, LOX,
: HIX, LOY, HIY, DATA, STATUS )
*
* Description :
*
* Method :
*
* Authors :
*
* David J. Allan (BHVAD::DJA)
*
* History :
*
* 21 Feb 91 : Original
*
* Type definitions :
*
IMPLICIT NONE
*
* Global constants :
*
INCLUDE 'SAE_PAR'
*
* Import :
*
INTEGER NAZ, NRAD ! # azimuthal,radial bins
REAL RANGE(2*NRAD) ! Radial bin size
REAL CX,CY ! Centre of circle
REAL LOX, HIX ! Range of rectangle in X
REAL LOY, HIY ! Range of rectangle in Y
*
* Import / export :
*
REAL DATA(NRAD,NAZ) ! Output data
*
* Status :
*
INTEGER STATUS
*
* Functions :
*
REAL GEO_CAREA
*
* Local variables :
*
REAL FACTOR ! Normalisation factor
INTEGER I, J ! Loop over bins
INTEGER IR ! Range index counter
*-
IF ( STATUS .EQ. SAI__OK ) THEN
DO I = 1, NRAD
IR = (I-1)*2+1
FACTOR = ( GEO_CAREA(RANGE(IR+1),CX,CY,LOX,HIX,LOY,HIY)
: - GEO_CAREA(RANGE(IR),CX,CY,LOX,HIX,LOY,HIY) )
: / REAL(NAZ)
DO J = 1, NAZ
DATA(I,J) = DATA(I,J) / FACTOR
END DO
END DO
END IF
END