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