      PROGRAM IPNPOSITION
C.....TESTS POSITIONS TO SEE HOW MANY IPN LOCALIZATIONS ARE CONSISTENT WITH THEM
      CHARACTER HEADER*80,NAME*10
      DOUBLE PRECISION PI,RADEG,DRAB,
     CDDECB,DRB,DRAIPN1,DDECIPN1,DRIPN1,DDRIPN1,DRAIPN2,
     CDDECIPN2,DRIPN2,DDRIPN2,DWATCHRA,DWATCHDEC,DWATCHR,DCOMPTELRA,
     CDCOMPTELDEC,DCOMPTELR,DEGRETRA,DEGRETDEC,DEGRETR,DPHEBUSRA,
     CDPHEBUSDEC,DPHEBUSR,DSIGMARA,DSIGMADEC,DSIGMAR,DSAXRA,DSAXDEC,
     CDSAXR,DHRAB,DHDECB,DHRB,DKONUSEL1,DKONUSEL2,DPRA,DPDEC,DPR,
     CDORAB,DODECB,DORB,STEP1,STEP2,AREA1,AREA2,FRACTION,AREA1TOT,
     CAREA2TOT,FRACTIONTOT,R,DR,RACIRC,DECCIRC,RCIRC,RATEST,DECTEST,
     CDRA,DDEC,DJD1,DJD2,DJDGRB
      REAL KONUSEL1,KONUSEL2
      COMMON/GEOM/PI,RADEG
      LOGICAL KON,IN,INANN,INANN1,INANN2,INOTHER
      PI = 3.14159265358979323846D0
      RADEG=PI/180.D0
      CALL INIT
      FRACTIONTOT1=0.
      NIN1=0
C.....OPEN INPUT FILE.  FORMAT IS NAME, RA, DEC, DATE RANGE, SEPARATOR:
C     2005nb  
C     183.406708  16.121167 
C     2005 12  1.0  2005 12 18.0
****************************************************************
C     2005mn  
C     57.326833   -.692056 
C     2005 11 15.0  2005 11 24.0
C      OPEN(10,FILE='FILE.input')
1     NCALC=0
      ARTOT=0.
      AR1TOT=0.
      AREA1TOT=0.
      AREA2TOT=0.
      FRACTIONTOT=0.
      NREAD=0
      NTREATED=0
      NBAT=0
      NIPN=0
      NIPN2=0
      NW=0
      NC=0
      NE=0
      NP=0
      NS=0
      N1ANN=0
      N2ANN=0
      N1ANNAND=0
      NEL=0
      NBAT=0
      NSAX=0
      NPHE=0
      NSIG=0
      NWAT=0
      NEGR=0
      NCOM=0
      NHET=0
      NPLAN=0
      NOTHER=0
      NNOINT=0
      NNOINT2=0
      N2ANNAND=0
      NNOINT2AND=0
      N2OTHER=0
      NNOINTN2=0
      NIN=0
      NAME='WEBSITE'
C      READ(10,4,END=5000)NAME
4     FORMAT(A10)
C      PRINT *,NAME
C.....GET RA, DEC POSITION
C      READ(10,*)RATEST,DECTEST
C      PRINT *,'TESTING RA, DEC=',RATEST,DECTEST
C.....READ DATE RANGE 
C      READ(10,*)IYR1,IM1,DY1,IYR2,IM2,DY2
C.....CONVERT TO JULIAN DATES	
C.....READ DATE RANGE AND RA, DEC
      READ *,IYR1,IM1,DY1,IYR2,IM2,DY2,RATEST,DECTEST
      CALL JUL(DY1,IM1,IYR1,DJD1)
      CALL JUL(DY2,IM2,IYR2,DJD2)
      DJD11=DJD1+2415020.
      DJD21=DJD2+2415020.
      PRINT *,'LOOKING FOR LOCALIZATIONS BETWEEN ',IYR1,IM1,DY1,' AND ',
     CIYR2,IM2,DY2
      PRINT *,' (JULIAN DAYS ',DJD11,' TO ',DJD21,' INCLUSIVE)'
      PRINT *,'NOT ALL EVENTS WILL HAVE LOCALIZATIONS'
	PRINT *,'SEND ME AN E-MAIL IF YOU WANT TO HAVE THE DATABASE POPULA
     CTED FOR YOUR SEARCH.  khurley@ssl.berkeley.edu'
	PRINT *,' '
      READ(13,*)NTOTAL
      CLOSE(13)
      READ(1,5)HEADER
5     FORMAT(A80)
C.... READ THE IPN DATA FILE.  THERE ARE 47 VARIABLES.  FORMAT:
C.....DAY,MONTH,YEAR,HOURS,MINUTES,SECONDS,BATSE #,BATSE RA,BATSE DEC,BATSE R,
C.....IPN RA1,IPN DEC1,IPN R1,IPN DR1,IPN RA2,IPN DEC2,IPN R2,IPN DR2,WATCH RA,WATCH DEC,WATCH R,
C.....COMPTEL RA,COMPTEL DEC,COMPTEL R,EGRET RA,EGRET DEC,EGRET R,PHEBUS RA,PHEBUS DEC,PHEBUS R,
C.....SIGMA RA,SIGMA DEC,SIGMA R,SAX RA,SAX DEC,SAX R,HETE RA,HETE DEC,HETE R,KONUS ECLIPTIC1,
C.....KONUS ECLIPTIC2,PLANET BLOCKING RA,PLANET BLOCKING DEC,PLANET  BLOCKING RADIUS,
C.... OTHER LOCALIZATION RA,OTHER LOCALIZATION DEC,OTHER LOCALIZATION RADIUS
10    READ(1,*,END=4000,ERR=30)IDAY,MON,IYR,IHR,IMIN,ISEC,
     CIBNO,RAB,DECB,RB,
     CRAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,DECIPN2,RIPN2,DELRIPN2,
     CWATCHRA,WATCHDEC,WATCHR,
     CCOMPTELRA,COMPTELDEC,COMPTELR,
     CEGRETRA,EGRETDEC,EGRETR,
     CPHEBUSRA,PHEBUSDEC,PHEBUSR,
     CSIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,
     CHRAB,HDECB,HRB,
     CKONUSEL1,KONUSEL2,
     CPRA,PDEC,PR,
     CORAB,ODECB,ORB
      NREAD=NREAD+1
C.....CHECK JULIAN DATE BEFORE PROCEEDING
      DY=IDAY
      CALL JUL(DY,MON,IYR,DJDGRB)
      IF(DJDGRB.LT.DJD1 .OR. DJDGRB.GT.DJD2)GOTO 10
      DRAB=DBLE(RAB)
      DDECB=DBLE(DECB)
      DRB=DBLE(RB)
      DRAIPN1=DBLE(RAIPN1)
      DDECIPN1=DBLE(DECIPN1)
      DRIPN1=DBLE(RIPN1)
      DDRIPN1=DBLE(DELRIPN1)
      DRAIPN2=DBLE(RAIPN2)
      DDECIPN2=DBLE(DECIPN2)
      DRIPN2=DBLE(RIPN2)
      DDRIPN2=DBLE(DELRIPN2)
      DWATCHRA=DBLE(WATCHRA)
      DWATCHDEC=DBLE(WATCHDEC)
      DWATCHR=DBLE(WATCHR)
      DCOMPTELRA=DBLE(COMPTELRA)
      DCOMPTELDEC=DBLE(COMPTELDEC)
      DCOMPTELR=DBLE(COMPTELR)
      DEGRETRA=DBLE(EGRETRA)
      DEGRETDEC=DBLE(EGRETDEC)
      DEGRETR=DBLE(EGRETR)
      DPHEBUSRA=DBLE(PHEBUSRA)
      DPHEBUSDEC=DBLE(PHEBUSDEC)
      DPHEBUSR=DBLE(PHEBUSR)
      DSIGMARA=DBLE(SIGMARA)
      DSIGMADEC=DBLE(SIGMADEC)
      DSIGMAR=DBLE(SIGMAR)
      DSAXRA=DBLE(SAXRA)
      DSAXDEC=DBLE(SAXDEC)
      DSAXR=DBLE(SAXR)
      DHRAB=DBLE(HRAB)
      DHDECB=DBLE(HDECB)
      DHRB=DBLE(HRB)
      DKONUSEL1=DBLE(KONUSEL1)
      DKONUSEL2=DBLE(KONUSEL2)
      DPRA=DBLE(PRA)
      DPDEC=DBLE(PDEC)
      DPR=DBLE(PR)
      DORAB=DBLE(ORAB)
      DODECB=DBLE(ODECB)
      DORB=DBLE(ORB)
      GOTO 50
30    NREAD=NREAD+1
      PRINT *,NREAD,IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,RAB,
     CDECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB
      STOP 'ERROR'
50    CONTINUE
      WRITE(*,51)NREAD,NTOTAL,IDAY,MON,IYR,IHR,IMIN,ISEC
51    FORMAT(' TESTING GRB NUMBER ',I4,' OUT OF ~'I6, ' ON ' 
     C,I2,1X,I2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2)
C.....DO THE SIMPLE CASES FIRST.  THESE ARE THE LOCALIZATIONS THAT CONSIST
C.....OF A SINGLE ERROR CIRCLE OR ANNULUS
C.....CASE 1.  SINGLE IPN ANNULUS
      IF(RAB .EQ. 0. .AND. DECB .EQ. 0. .AND. RB .EQ. 0.  .AND.
     CRAIPN1 .NE. 0. .AND. DECIPN1 .NE. 0. .AND. RIPN1 .NE. 0.  .AND.
     CDELRIPN1 .NE. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .EQ. 0. .AND. WATCHDEC .EQ. 0. .AND. WATCHR .EQ. 0. .AND.
     CCOMPTELRA .EQ. 0. .AND. COMPTELDEC .EQ. 0. .AND. COMPTELR .EQ. 0. 
     C.AND.
     CEGRETRA .EQ. 0. .AND. EGRETDEC .EQ. 0. .AND. EGRETR .EQ. 0. .AND.
     CPHEBUSRA .EQ. 0. .AND. PHEBUSDEC .EQ. 0. .AND. PHEBUSR .EQ. 0. 
     C.AND.
     CSIGMARA .EQ. 0. .AND. SIGMADEC .EQ. 0. .AND. SIGMAR .EQ. 0. .AND.
     CSAXRA .EQ. 0. .AND. SAXDEC .EQ. 0. .AND. SAXR .EQ. 0. .AND.
     CHRAB .EQ. 0. .AND. HDECB .EQ. 0. .AND. HRB .EQ. 0. .AND.
     CKONUSEL1 .EQ. 0. .AND. KONUSEL2 .EQ. 0. .AND.
     CPRA .EQ. 0. .AND. PDEC .EQ. 0. .AND. PR .EQ. 0. .AND.
     CORAB .EQ. 0. .AND. ODECB .EQ. 0. .AND. ORB .EQ. 0.)THEN
      CALL ANNULUS(DRAIPN1,DDECIPN1,DRIPN1,DDRIPN1,AREA1,AREA2,FRACTION,
     CRATEST,DECTEST,IN)
      IF(IN)THEN
      NIN=NIN+1
      WRITE(*,60)IDAY,MON,IYR,IHR,IMIN,ISEC,RAIPN1,DECIPN1,RIPN1,
     CDELRIPN1,AREA1,AREA2,FRACTION
60    FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'RAIPN',
     C4X,'DECIPN',4X,'RIPN',5X,'DRIPN',3X,'STER',5X,'SQ.DEG.'2X,
     C'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,F8.4,1X,
     CF8.4,1X,F8.4,1X,F8.4,1X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS IN THIS ERROR ANNULUS'
      ENDIF
      IF(.NOT. IN)THEN
      WRITE(*,60)IDAY,MON,IYR,IHR,IMIN,ISEC,RAIPN1,DECIPN1,RIPN1,
     CDELRIPN1,AREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATION'
	ENDIF
      PRINT *,'--------------------------------------------------------'
	AREA1TOT=AREA1TOT+AREA1
	AREA2TOT=AREA2TOT+AREA2
	FRACTIONTOT=FRACTIONTOT+FRACTION
	N1ANN=N1ANN+1
	NTREATED=NTREATED+1
	GOTO 10
	ENDIF
C.....CASE 2.  KONUS ECLIPTIC LATITUDE BAND
      IF(RAB .EQ. 0. .AND. DECB .EQ. 0. .AND. RB .EQ. 0.  .AND.
     CRAIPN1 .EQ. 0. .AND. DECIPN1 .EQ. 0. .AND. RIPN1 .EQ. 0.  .AND.
     CDELRIPN1 .EQ. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .EQ. 0. .AND. WATCHDEC .EQ. 0. .AND. WATCHR .EQ. 0. .AND.
     CCOMPTELRA .EQ. 0. .AND. COMPTELDEC .EQ. 0. .AND. COMPTELR .EQ. 0. 
     C.AND.
     CEGRETRA .EQ. 0. .AND. EGRETDEC .EQ. 0. .AND. EGRETR .EQ. 0. .AND.
     CPHEBUSRA .EQ. 0. .AND. PHEBUSDEC .EQ. 0. .AND. PHEBUSR .EQ. 0. 
     C.AND.
     CSIGMARA .EQ. 0. .AND. SIGMADEC .EQ. 0. .AND. SIGMAR .EQ. 0. .AND.
     CSAXRA .EQ. 0. .AND. SAXDEC .EQ. 0. .AND. SAXR .EQ. 0. .AND.
     CHRAB .EQ. 0. .AND. HDECB .EQ. 0. .AND. HRB .EQ. 0. .AND.
     CKONUSEL1 .NE. 0. .AND. KONUSEL2 .NE. 0. .AND.
     CPRA .EQ. 0. .AND. PDEC .EQ. 0. .AND. PR .EQ. 0. .AND.
     CORAB .EQ. 0. .AND. ODECB .EQ. 0. .AND. ORB .EQ. 0.)THEN
	R=90.-DBLE((KONUSEL1+KONUSEL2)/2.)
	IF(R.LT.0.)R=DBLE(-R+90.D0)
	DR=DBLE(ABS(KONUSEL1-KONUSEL2)/2.)
	DRA=270.D0
	DDEC=90.-23.4393
	CALL ANNULUS(DRA,DDEC,R,DR,AREA1,AREA2,FRACTION,RATEST,
     CDECTEST,IN)
	IF(IN)THEN
	NIN=NIN+1
	WRITE(*,70)IDAY,MON,IYR,IHR,IMIN,ISEC,KONUSEL1,KONUSEL2,
     CAREA1,AREA2,FRACTION
70    FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'EL1',6X,'EL2',6X,'STER',5X,'SQ.DEG.'2X,'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,F8.4,1X,
     CF8.4,1X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS IN THIS ECLIPTIC LATITUDE BAND'
	ENDIF
      IF(.NOT. IN)THEN
	WRITE(*,70)IDAY,MON,IYR,IHR,IMIN,ISEC,KONUSEL1,KONUSEL2,
     CAREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATION'
      ENDIF
      PRINT *,'--------------------------------------------------------'
	AREA1TOT=AREA1TOT+AREA1
	AREA2TOT=AREA2TOT+AREA2
	FRACTIONTOT=FRACTIONTOT+FRACTION
	NEL=NEL+1
	NTREATED=NTREATED+1
	GOTO 10
	ENDIF
C.....CASE 3. BATSE ERROR CIRCLE
C.....THERE IS A BATSE BURST WITH RA=0
      IF(RB .NE. 0.  .AND.
     CRAIPN1 .EQ. 0. .AND. DECIPN1 .EQ. 0. .AND. RIPN1 .EQ. 0.  .AND.
     CDELRIPN1 .EQ. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .EQ. 0. .AND. WATCHDEC .EQ. 0. .AND. WATCHR .EQ. 0. .AND.
     CCOMPTELRA .EQ. 0. .AND. COMPTELDEC .EQ. 0. .AND. COMPTELR .EQ. 0. 
     C.AND.
     CEGRETRA .EQ. 0. .AND. EGRETDEC .EQ. 0. .AND. EGRETR .EQ. 0. .AND.
     CPHEBUSRA .EQ. 0. .AND. PHEBUSDEC .EQ. 0. .AND. PHEBUSR .EQ. 0. 
     C.AND.
     CSIGMARA .EQ. 0. .AND. SIGMADEC .EQ. 0. .AND. SIGMAR .EQ. 0. .AND.
     CSAXRA .EQ. 0. .AND. SAXDEC .EQ. 0. .AND. SAXR .EQ. 0. .AND.
     CHRAB .EQ. 0. .AND. HDECB .EQ. 0. .AND. HRB .EQ. 0. .AND.
     CKONUSEL1 .EQ. 0. .AND. KONUSEL2 .EQ. 0. .AND.
     CPRA .EQ. 0. .AND. PDEC .EQ. 0. .AND. PR .EQ. 0. .AND.
     CORAB .EQ. 0. .AND. ODECB .EQ. 0. .AND. ORB .EQ. 0.)THEN
	DRB=DSQRT(DRB**2 + 1.6D0**2)
	CALL CIRCLE(DRAB,DDECB,DRB,AREA1,AREA2,FRACTION,RATEST,DECTEST,IN)
	IF(IN)THEN
	NIN=NIN+1
	WRITE(*,80)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,RAB,DECB,RB,
     CAREA1,AREA2,FRACTION
80    FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'BATSE#',6X,'BATSE RA',6X,'BATSE DEC',6X,'RB',6X,'STER',5X,
     C'SQ.DEG.'2X,'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,I6,6X,F8.4,6X,
     CF8.4,6X,F8.4,1X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS IN THIS ERROR CIRCLE'
	ENDIF
      IF(.NOT. IN)THEN
	WRITE(*,80)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,RAB,DECB,RB,
     CAREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATI
     CON'
      PRINT *,'--------------------------------------------------------'
      ENDIF
	AREA1TOT=AREA1TOT+AREA1
	AREA2TOT=AREA2TOT+AREA2
	FRACTIONTOT=FRACTIONTOT+FRACTION
	NBAT=NBAT+1
	NTREATED=NTREATED+1
	GOTO 10
	ENDIF
C.....CASE 4. BEPPOSAX GRBM ERROR CIRCLE
      IF(RAB .EQ. 0. .AND. DECB .EQ. 0. .AND. RB .EQ. 0.  .AND.
     CRAIPN1 .EQ. 0. .AND. DECIPN1 .EQ. 0. .AND. RIPN1 .EQ. 0.  .AND.
     CDELRIPN1 .EQ. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .EQ. 0. .AND. WATCHDEC .EQ. 0. .AND. WATCHR .EQ. 0. .AND.
     CCOMPTELRA .EQ. 0. .AND. COMPTELDEC .EQ. 0. .AND. COMPTELR .EQ. 0. 
     C.AND.
     CEGRETRA .EQ. 0. .AND. EGRETDEC .EQ. 0. .AND. EGRETR .EQ. 0. .AND.
     CPHEBUSRA .EQ. 0. .AND. PHEBUSDEC .EQ. 0. .AND. PHEBUSR .EQ. 0. 
     C.AND.
     CSIGMARA .EQ. 0. .AND. SIGMADEC .EQ. 0. .AND. SIGMAR .EQ. 0. .AND.
C.....THERE IS A SAX GRBM EVENT WITH DECLINATION ZERO
     CSAXR .NE. 0. .AND.
     CHRAB .EQ. 0. .AND. HDECB .EQ. 0. .AND. HRB .EQ. 0. .AND.
     CKONUSEL1 .EQ. 0. .AND. KONUSEL2 .EQ. 0. .AND.
     CPRA .EQ. 0. .AND. PDEC .EQ. 0. .AND. PR .EQ. 0. .AND.
     CORAB .EQ. 0. .AND. ODECB .EQ. 0. .AND. ORB .EQ. 0.)THEN
	CALL CIRCLE(DSAXRA,DSAXDEC,DSAXR,AREA1,AREA2,FRACTION,RATEST,
     CDECTEST,IN)
	IF(IN)THEN
	NIN=NIN+1
	WRITE(*,90)IDAY,MON,IYR,IHR,IMIN,ISEC,SAXRA,SAXDEC,SAXR,
     CAREA1,AREA2,FRACTION
90    FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'SAXRA',4X,'SAXDEC',4X,'SAXR',4X,'STER',5X,'SQ.DEG.'2X,
     C'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,F8.4,1X,
     CF8.4,1X,F8.4,1X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS IN THIS ERROR CIRCLE'
	ENDIF
      IF(.NOT. IN)THEN
	WRITE(*,90)IDAY,MON,IYR,IHR,IMIN,ISEC,SAXRA,SAXDEC,SAXR,
     CAREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATI
     CON'
      ENDIF
      PRINT *,'--------------------------------------------------------'
	AREA1TOT=AREA1TOT+AREA1
	AREA2TOT=AREA2TOT+AREA2
	FRACTIONTOT=FRACTIONTOT+FRACTION
	NSAX=NSAX+1
	NTREATED=NTREATED+1
	GOTO 10
	ENDIF
C.....CASE 5. PHEBUS ERROR CIRCLE
      IF(RAB .EQ. 0. .AND. DECB .EQ. 0. .AND. RB .EQ. 0.  .AND.
     CRAIPN1 .EQ. 0. .AND. DECIPN1 .EQ. 0. .AND. RIPN1 .EQ. 0.  .AND.
     CDELRIPN1 .EQ. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .EQ. 0. .AND. WATCHDEC .EQ. 0. .AND. WATCHR .EQ. 0. .AND.
     CCOMPTELRA .EQ. 0. .AND. COMPTELDEC .EQ. 0. .AND. COMPTELR .EQ. 0. 
     C.AND.
     CEGRETRA .EQ. 0. .AND. EGRETDEC .EQ. 0. .AND. EGRETR .EQ. 0. .AND.
     CPHEBUSRA .NE. 0. .AND. PHEBUSDEC .NE. 0. .AND. PHEBUSR .NE. 0. 
     C.AND.
     CSIGMARA .EQ. 0. .AND. SIGMADEC .EQ. 0. .AND. SIGMAR .EQ. 0. .AND.
     CSAXRA .EQ. 0. .AND. SAXDEC .EQ. 0. .AND. SAXR .EQ. 0. .AND.
     CHRAB .EQ. 0. .AND. HDECB .EQ. 0. .AND. HRB .EQ. 0. .AND.
     CKONUSEL1 .EQ. 0. .AND. KONUSEL2 .EQ. 0. .AND.
     CPRA .EQ. 0. .AND. PDEC .EQ. 0. .AND. PR .EQ. 0. .AND.
     CORAB .EQ. 0. .AND. ODECB .EQ. 0. .AND. ORB .EQ. 0.)THEN
	CALL CIRCLE(DPHEBUSRA,DPHEBUSDEC,DPHEBUSR,AREA1,AREA2,FRACTION,
     CRATEST,DECTEST,IN)
      IF(IN)THEN
	NIN=NIN+1
	WRITE(*,100)IDAY,MON,IYR,IHR,IMIN,ISEC,PHEBUSRA,PHEBUSDEC,PHEBUSR,
     CAREA1,AREA2,FRACTION
100   FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'PHEBUSRA',1X,'PHEBUSDEC',1X,'PHEBUSR',1X,'STER',5X,'SQ.DEG.'2X,
     C'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,F8.4,1X,
     CF8.4,1X,F8.4,1X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS IN THIS ERROR CIRCLE'
	ENDIF
      IF(.NOT. IN)THEN
	WRITE(*,100)IDAY,MON,IYR,IHR,IMIN,ISEC,PHEBUSRA,PHEBUSDEC,PHEBUSR,
     CAREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATI
     CON'
      ENDIF
      PRINT *,'--------------------------------------------------------'
	AREA1TOT=AREA1TOT+AREA1
	AREA2TOT=AREA2TOT+AREA2
	FRACTIONTOT=FRACTIONTOT+FRACTION
	NPHE=NPHE+1
	NTREATED=NTREATED+1
	GOTO 10
	ENDIF
C.....CASE 6. SIGMA ERROR CIRCLE
      IF(RAB .EQ. 0. .AND. DECB .EQ. 0. .AND. RB .EQ. 0.  .AND.
     CRAIPN1 .EQ. 0. .AND. DECIPN1 .EQ. 0. .AND. RIPN1 .EQ. 0.  .AND.
     CDELRIPN1 .EQ. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .EQ. 0. .AND. WATCHDEC .EQ. 0. .AND. WATCHR .EQ. 0. .AND.
     CCOMPTELRA .EQ. 0. .AND. COMPTELDEC .EQ. 0. .AND. COMPTELR .EQ. 0. 
     C.AND.
     CEGRETRA .EQ. 0. .AND. EGRETDEC .EQ. 0. .AND. EGRETR .EQ. 0. .AND.
     CPHEBUSRA .EQ. 0. .AND. PHEBUSDEC .EQ. 0. .AND. PHEBUSR .EQ. 0. 
     C.AND.
     CSIGMARA .NE. 0. .AND. SIGMADEC .NE. 0. .AND. SIGMAR .NE. 0. .AND.
     CSAXRA .EQ. 0. .AND. SAXDEC .EQ. 0. .AND. SAXR .EQ. 0. .AND.
     CHRAB .EQ. 0. .AND. HDECB .EQ. 0. .AND. HRB .EQ. 0. .AND.
     CKONUSEL1 .EQ. 0. .AND. KONUSEL2 .EQ. 0. .AND.
     CPRA .EQ. 0. .AND. PDEC .EQ. 0. .AND. PR .EQ. 0. .AND.
     CORAB .EQ. 0. .AND. ODECB .EQ. 0. .AND. ORB .EQ. 0.)THEN
	CALL CIRCLE(DSIGMARA,DSIGMADEC,DSIGMAR,AREA1,AREA2,FRACTION,
     CRATEST,DECTEST,IN)
      IF(IN)THEN
	NIN=NIN+1
	WRITE(*,110)IDAY,MON,IYR,IHR,IMIN,ISEC,SIGMARA,SIGMADEC,SIGMAR,
     CAREA1,AREA2,FRACTION
110   FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'SIGMARA',2X,'SIGMADEC',2X,'SIGMAR',2X,'STER',5X,'SQ.DEG.'2X,
     C'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,F8.4,1X,
     CF8.4,1X,F8.4,1X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS IN THIS ERROR CIRCLE'
	ENDIF
      IF(.NOT. IN)THEN
	WRITE(*,110)IDAY,MON,IYR,IHR,IMIN,ISEC,SIGMARA,SIGMADEC,SIGMAR,
     CAREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATI
     CON'
      ENDIF
      PRINT *,'--------------------------------------------------------'
	AREA1TOT=AREA1TOT+AREA1
	AREA2TOT=AREA2TOT+AREA2
	FRACTIONTOT=FRACTIONTOT+FRACTION
	NSIG=NSIG+1
	NTREATED=NTREATED+1
	GOTO 10
	ENDIF
C.....CASE 7. WATCH ERROR CIRCLE
      IF(RAB .EQ. 0. .AND. DECB .EQ. 0. .AND. RB .EQ. 0.  .AND.
     CRAIPN1 .EQ. 0. .AND. DECIPN1 .EQ. 0. .AND. RIPN1 .EQ. 0.  .AND.
     CDELRIPN1 .EQ. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .NE. 0. .AND. WATCHDEC .NE. 0. .AND. WATCHR .NE. 0. .AND.
     CCOMPTELRA .EQ. 0. .AND. COMPTELDEC .EQ. 0. .AND. COMPTELR .EQ. 0. 
     C.AND.
     CEGRETRA .EQ. 0. .AND. EGRETDEC .EQ. 0. .AND. EGRETR .EQ. 0. .AND.
     CPHEBUSRA .EQ. 0. .AND. PHEBUSDEC .EQ. 0. .AND. PHEBUSR .EQ. 0. 
     C.AND.
     CSIGMARA .EQ. 0. .AND. SIGMADEC .EQ. 0. .AND. SIGMAR .EQ. 0. .AND.
     CSAXRA .EQ. 0. .AND. SAXDEC .EQ. 0. .AND. SAXR .EQ. 0. .AND.
     CHRAB .EQ. 0. .AND. HDECB .EQ. 0. .AND. HRB .EQ. 0. .AND.
     CKONUSEL1 .EQ. 0. .AND. KONUSEL2 .EQ. 0. .AND.
     CPRA .EQ. 0. .AND. PDEC .EQ. 0. .AND. PR .EQ. 0. .AND.
     CORAB .EQ. 0. .AND. ODECB .EQ. 0. .AND. ORB .EQ. 0.)THEN
	CALL CIRCLE(DWATCHRA,DWATCHDEC,DWATCHR,AREA1,AREA2,FRACTION,
     CRATEST,DECTEST,IN)
      IF(IN)THEN
	NIN=NIN+1
	WRITE(*,120)IDAY,MON,IYR,IHR,IMIN,ISEC,WATCHRA,WATCHDEC,WATCHR,
     CAREA1,AREA2,FRACTION
120   FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'WATCHRA',2X,'WATCHDEC',2X,'WATCHR',2X,'STER',5X,'SQ.DEG.'2X,
     C'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,F8.4,1X,
     CF8.4,1X,F8.4,1X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS IN THIS ERROR CIRCLE'
	ENDIF
      IF(.NOT. IN)THEN
	WRITE(*,120)IDAY,MON,IYR,IHR,IMIN,ISEC,WATCHRA,WATCHDEC,WATCHR,
     CAREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATI
     CON'
      ENDIF
      PRINT *,'--------------------------------------------------------'
	AREA1TOT=AREA1TOT+AREA1
	AREA2TOT=AREA2TOT+AREA2
	FRACTIONTOT=FRACTIONTOT+FRACTION
	NWAT=NWAT+1
	NTREATED=NTREATED+1
	GOTO 10
	ENDIF
C.....CASE 8. EGRET ERROR CIRCLE
      IF(RAB .EQ. 0. .AND. DECB .EQ. 0. .AND. RB .EQ. 0.  .AND.
     CRAIPN1 .EQ. 0. .AND. DECIPN1 .EQ. 0. .AND. RIPN1 .EQ. 0.  .AND.
     CDELRIPN1 .EQ. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .EQ. 0. .AND. WATCHDEC .EQ. 0. .AND. WATCHR .EQ. 0. .AND.
     CCOMPTELRA .EQ. 0. .AND. COMPTELDEC .EQ. 0. .AND. COMPTELR .EQ. 0. 
     C.AND.
     CEGRETRA .NE. 0. .AND. EGRETDEC .NE. 0. .AND. EGRETR .NE. 0. .AND.
     CPHEBUSRA .EQ. 0. .AND. PHEBUSDEC .EQ. 0. .AND. PHEBUSR .EQ. 0. 
     C.AND.
     CSIGMARA .EQ. 0. .AND. SIGMADEC .EQ. 0. .AND. SIGMAR .EQ. 0. .AND.
     CSAXRA .EQ. 0. .AND. SAXDEC .EQ. 0. .AND. SAXR .EQ. 0. .AND.
     CHRAB .EQ. 0. .AND. HDECB .EQ. 0. .AND. HRB .EQ. 0. .AND.
     CKONUSEL1 .EQ. 0. .AND. KONUSEL2 .EQ. 0. .AND.
     CPRA .EQ. 0. .AND. PDEC .EQ. 0. .AND. PR .EQ. 0. .AND.
     CORAB .EQ. 0. .AND. ODECB .EQ. 0. .AND. ORB .EQ. 0.)THEN
	CALL CIRCLE(DEGRETRA,DEGRETDEC,DEGRETR,AREA1,AREA2,FRACTION,
     CRATEST,DECTEST,IN)
	IF(IN)THEN
	NIN=NIN+1
	WRITE(*,130)IDAY,MON,IYR,IHR,IMIN,ISEC,EGRETRA,EGRETDEC,EGRETR,
     CAREA1,AREA2,FRACTION
130   FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'EGRETRA',2X,'EGRETDEC',2X,'EGRETR',2X,'STER',5X,'SQ.DEG.'2X,
     C'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,F8.4,1X,
     CF8.4,1X,F8.4,1X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS IN THIS ERROR CIRCLE'
	ENDIF
      IF(.NOT. IN)THEN
	WRITE(*,130)IDAY,MON,IYR,IHR,IMIN,ISEC,EGRETRA,EGRETDEC,EGRETR,
     CAREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATI
     CON'
      ENDIF
      PRINT *,'--------------------------------------------------------'
	AREA1TOT=AREA1TOT+AREA1
	AREA2TOT=AREA2TOT+AREA2
	FRACTIONTOT=FRACTIONTOT+FRACTION
	NEGR=NEGR+1
	NTREATED=NTREATED+1
	GOTO 10
	ENDIF
C.....CASE 9. COMPTEL ERROR CIRCLE
      IF(RAB .EQ. 0. .AND. DECB .EQ. 0. .AND. RB .EQ. 0.  .AND.
     CRAIPN1 .EQ. 0. .AND. DECIPN1 .EQ. 0. .AND. RIPN1 .EQ. 0.  .AND.
     CDELRIPN1 .EQ. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .EQ. 0. .AND. WATCHDEC .EQ. 0. .AND. WATCHR .EQ. 0. .AND.
     CCOMPTELRA .NE. 0. .AND. COMPTELDEC .NE. 0. .AND. COMPTELR .NE. 0. 
     C.AND.
     CEGRETRA .EQ. 0. .AND. EGRETDEC .EQ. 0. .AND. EGRETR .EQ. 0. .AND.
     CPHEBUSRA .EQ. 0. .AND. PHEBUSDEC .EQ. 0. .AND. PHEBUSR .EQ. 0. 
     C.AND.
     CSIGMARA .EQ. 0. .AND. SIGMADEC .EQ. 0. .AND. SIGMAR .EQ. 0. .AND.
     CSAXRA .EQ. 0. .AND. SAXDEC .EQ. 0. .AND. SAXR .EQ. 0. .AND.
     CHRAB .EQ. 0. .AND. HDECB .EQ. 0. .AND. HRB .EQ. 0. .AND.
     CKONUSEL1 .EQ. 0. .AND. KONUSEL2 .EQ. 0. .AND.
     CPRA .EQ. 0. .AND. PDEC .EQ. 0. .AND. PR .EQ. 0. .AND.
     CORAB .EQ. 0. .AND. ODECB .EQ. 0. .AND. ORB .EQ. 0.)THEN
	CALL CIRCLE(DCOMPTELRA,DCOMPTELDEC,DCOMPTELR,AREA1,AREA2,
     CFRACTION,RATEST,DECTEST,IN)
	IF(IN)THEN
	NIN=NIN+1
	WRITE(*,140)IDAY,MON,IYR,IHR,IMIN,ISEC,COMPTELRA,COMPTELDEC,
     CCOMPTELR,AREA1,AREA2,FRACTION
140   FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'COMPTELRA',1X,'COMPTELDEC',1X,'COMPTELR',1X,'STER',5X,'SQ.DEG.'2X,
     C'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,F8.4,1X,
     CF8.4,1X,F8.4,1X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS IN THIS ERROR CIRCLE'
	ENDIF
      IF(.NOT. IN)THEN
	WRITE(*,140)IDAY,MON,IYR,IHR,IMIN,ISEC,COMPTELRA,COMPTELDEC,
     CCOMPTELR,AREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATI
     CON'
      ENDIF
      PRINT *,'--------------------------------------------------------'
      AREA1TOT=AREA1TOT+AREA1
      AREA2TOT=AREA2TOT+AREA2
      FRACTIONTOT=FRACTIONTOT+FRACTION
      NCOM=NCOM+1
	NTREATED=NTREATED+1
	GOTO 10
	ENDIF
C.....CASE 10. HETE ERROR CIRCLE
      IF(RAB .EQ. 0. .AND. DECB .EQ. 0. .AND. RB .EQ. 0.  .AND.
     CRAIPN1 .EQ. 0. .AND. DECIPN1 .EQ. 0. .AND. RIPN1 .EQ. 0.  .AND.
     CDELRIPN1 .EQ. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .EQ. 0. .AND. WATCHDEC .EQ. 0. .AND. WATCHR .EQ. 0. .AND.
     CCOMPTELRA .EQ. 0. .AND. COMPTELDEC .EQ. 0. .AND. COMPTELR .EQ. 0. 
     C.AND.
     CEGRETRA .EQ. 0. .AND. EGRETDEC .EQ. 0. .AND. EGRETR .EQ. 0. .AND.
     CPHEBUSRA .EQ. 0. .AND. PHEBUSDEC .EQ. 0. .AND. PHEBUSR .EQ. 0. 
     C.AND.
     CSIGMARA .EQ. 0. .AND. SIGMADEC .EQ. 0. .AND. SIGMAR .EQ. 0. .AND.
     CSAXRA .EQ. 0. .AND. SAXDEC .EQ. 0. .AND. SAXR .EQ. 0. .AND.
     CHRAB .NE. 0. .AND. HDECB .NE. 0. .AND. HRB .NE. 0. .AND.
     CKONUSEL1 .EQ. 0. .AND. KONUSEL2 .EQ. 0. .AND.
     CPRA .EQ. 0. .AND. PDEC .EQ. 0. .AND. PR .EQ. 0. .AND.
     CORAB .EQ. 0. .AND. ODECB .EQ. 0. .AND. ORB .EQ. 0.)THEN
	CALL CIRCLE(DHRAB,DHDECB,DHRB,AREA1,AREA2,FRACTION,
     CRATEST,DECTEST,IN)
	IF(IN)THEN
	NIN=NIN+1
	WRITE(*,150)IDAY,MON,IYR,IHR,IMIN,ISEC,HRAB,HDECB,
     CHRB,AREA1,AREA2,FRACTION
150   FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'HETE RA',2X,'HETE DEC',1X,'HETE R',3X,'STER',5X,'SQ.DEG.'2X,
     C'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,F8.4,1X,
     CF8.4,1X,F8.4,1X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS IN THIS ERROR CIRCLE'
	ENDIF
      IF(.NOT. IN)THEN
	WRITE(*,150)IDAY,MON,IYR,IHR,IMIN,ISEC,HRAB,HDECB,
     CHRB,AREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATI
     CON'
      ENDIF
      PRINT *,'--------------------------------------------------------'
      AREA1TOT=AREA1TOT+AREA1
      AREA2TOT=AREA2TOT+AREA2
      FRACTIONTOT=FRACTIONTOT+FRACTION
      NHET=NHET+1
      NTREATED=NTREATED+1
      GOTO 10
      ENDIF
C.....CASE 11. PLANET-BLOCKING
      IF(RAB .EQ. 0. .AND. DECB .EQ. 0. .AND. RB .EQ. 0.  .AND.
     CRAIPN1 .EQ. 0. .AND. DECIPN1 .EQ. 0. .AND. RIPN1 .EQ. 0.  .AND.
     CDELRIPN1 .EQ. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .EQ. 0. .AND. WATCHDEC .EQ. 0. .AND. WATCHR .EQ. 0. .AND.
     CCOMPTELRA .EQ. 0. .AND. COMPTELDEC .EQ. 0. .AND. COMPTELR .EQ. 0. 
     C.AND.
     CEGRETRA .EQ. 0. .AND. EGRETDEC .EQ. 0. .AND. EGRETR .EQ. 0. .AND.
     CPHEBUSRA .EQ. 0. .AND. PHEBUSDEC .EQ. 0. .AND. PHEBUSR .EQ. 0. 
     C.AND.
     CSIGMARA .EQ. 0. .AND. SIGMADEC .EQ. 0. .AND. SIGMAR .EQ. 0. .AND.
     CSAXRA .EQ. 0. .AND. SAXDEC .EQ. 0. .AND. SAXR .EQ. 0. .AND.
     CHRAB .EQ. 0. .AND. HDECB .EQ. 0. .AND. HRB .EQ. 0. .AND.
     CKONUSEL1 .EQ. 0. .AND. KONUSEL2 .EQ. 0. .AND.
     CPRA .NE. 0. .AND. PDEC .NE. 0. .AND. PR .NE. 0. .AND.
     CORAB .EQ. 0. .AND. ODECB .EQ. 0. .AND. ORB .EQ. 0.)THEN
      R=180.-PR
      CALL CIRCLE(DPRA,DPDEC,DPR,AREA1,AREA2,FRACTION,
     CRATEST,DECTEST,IN)
      IF(.NOT.IN)THEN
      NIN=NIN+1
      WRITE(*,160)IDAY,MON,IYR,IHR,IMIN,ISEC,PRA,PDEC,
     CPR,AREA1,AREA2,FRACTION
160   FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'PLANET RA',1X,'PLANET DEC',1X,'PLANET R',3X,'STER',5X,'SQ.DEG.'2X,
     C'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,F8.4,2X,
     CF8.4,2X,F8.4,4X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS NOT PLANET-BLOCKED'
      ENDIF
      IF(IN)THEN
      WRITE(*,160)IDAY,MON,IYR,IHR,IMIN,ISEC,PRA,PDEC,
     CPR,AREA1,AREA2,FRACTION
	PRINT *,'TEST POINT IS PLANET-BLOCKED'
      ENDIF
      PRINT *,'--------------------------------------------------------'
      AREA1TOT=AREA1TOT+AREA1
      AREA2TOT=AREA2TOT+AREA2
      FRACTIONTOT=FRACTIONTOT+FRACTION
      NPLAN=NPLAN+1
      NTREATED=NTREATED+1
      GOTO 10
      ENDIF
C.....CASE 12. OTHER ERROR CIRCLE
      IF(RAB .EQ. 0. .AND. DECB .EQ. 0. .AND. RB .EQ. 0.  .AND.
     CRAIPN1 .EQ. 0. .AND. DECIPN1 .EQ. 0. .AND. RIPN1 .EQ. 0.  .AND.
     CDELRIPN1 .EQ. 0. .AND.
     CRAIPN2 .EQ. 0. .AND. DECIPN2 .EQ. 0. .AND. RIPN2 .EQ. 0. .AND.
     CDELRIPN2 .EQ. 0. .AND.
     CWATCHRA .EQ. 0. .AND. WATCHDEC .EQ. 0. .AND. WATCHR .EQ. 0. .AND.
     CCOMPTELRA .EQ. 0. .AND. COMPTELDEC .EQ. 0. .AND. COMPTELR .EQ. 0. 
     C.AND.
     CEGRETRA .EQ. 0. .AND. EGRETDEC .EQ. 0. .AND. EGRETR .EQ. 0. .AND.
     CPHEBUSRA .EQ. 0. .AND. PHEBUSDEC .EQ. 0. .AND. PHEBUSR .EQ. 0. 
     C.AND.
     CSIGMARA .EQ. 0. .AND. SIGMADEC .EQ. 0. .AND. SIGMAR .EQ. 0. .AND.
     CSAXRA .EQ. 0. .AND. SAXDEC .EQ. 0. .AND. SAXR .EQ. 0. .AND.
     CHRAB .EQ. 0. .AND. HDECB .EQ. 0. .AND. HRB .EQ. 0. .AND.
     CKONUSEL1 .EQ. 0. .AND. KONUSEL2 .EQ. 0. .AND.
     CPRA .EQ. 0. .AND. PDEC .EQ. 0. .AND. PR .EQ. 0. .AND.
     CORAB .NE. 0. .AND. ODECB .NE. 0. .AND. ORB .NE. 0.)THEN
      CALL CIRCLE(DORAB,DODECB,DORB,AREA1,AREA2,FRACTION,
     CRATEST,DECTEST,IN)
      IF(IN)THEN
	NIN=NIN+1
      WRITE(*,170)IDAY,MON,IYR,IHR,IMIN,ISEC,ORAB,ODECB,
     CORB,AREA1,AREA2,FRACTION
170   FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'OTHER RA',1X,'OTHER DEC',1X,'OTHER R',2X,'STER',5X,'SQ.DEG.'2X,
     C'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,F8.4,1X,
     CF8.4,2X,F8.4,1X,1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS IN THE ERROR CIRCLE'
	ENDIF
      IF(.NOT. IN)THEN
      WRITE(*,170)IDAY,MON,IYR,IHR,IMIN,ISEC,ORAB,ODECB,
     CORB,AREA1,AREA2,FRACTION
	PRINT *,'TEST POINT IS NOT IN THE ERROR CIRCLE'
      ENDIF
      PRINT *,'--------------------------------------------------------'
      AREA1TOT=AREA1TOT+AREA1
      AREA2TOT=AREA2TOT+AREA2
      FRACTIONTOT=FRACTIONTOT+FRACTION
      NOTHER=NOTHER+1
      NTREATED=NTREATED+1
      GOTO 10
      ENDIF
C.....ALL THE SIMPLE CASES HAVE NOW BEEN DONE.  START
C.....THE HARDER ONES
C.....CASE 13.  ONE AND ONLY ONE IPN ANNULUS, AND ONE OR MORE OF ANYTHING ELSE 
C.....(BATSE ERROR CIRCLE, KONUS EL, PLANET-BLOCKING, ...)
      IF(RAIPN1.NE.0. .AND. DECIPN1.NE.0. .AND. RIPN1.NE.0. .AND.
     CDELRIPN1.NE.0. .AND.
     CRAIPN2.EQ.0. .AND. DECIPN2.EQ.0. .AND. RIPN2.EQ.0. .AND.
     CDELRIPN2.EQ.0.)THEN
      IF(RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE.0.) GOTO 180 
      IF(WATCHRA.NE.0. AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.) GOTO 180
      IF(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.) 
     CGOTO 180
      IF(EGRETRA.NE.0. AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.) GOTO 180
      IF(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.) 
     CGOTO 180
      IF(SIGMARA.NE.0. AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.) GOTO 180
      IF(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.) GOTO 180
      IF(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.) GOTO 180
      IF(KONUSEL1.NE.0. .OR. KONUSEL2.NE.0.) GOTO 180
      IF(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.) GOTO 180
      IF(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.)GOTO 180
      ENDIF
      GOTO 200
180   CONTINUE
      N1ANNAND=N1ANNAND+1
C.....STEP1 AND STEP2 ARE THE AZIMUTHAL AND RADIAL STEPS THAT
C.....WILL BE USED TO GO AROUND THE ANNULUS AND CALCULATE THE AREA.
C.....THE ANGULAR STEP IS PI*STEP	RADIANS
      STEP1=0.001D0
      STEP2=DDRIPN1/1000.D0
      STEP1=STEP2
      AREA1=0.D0
      CALL GC13(DRAIPN1,DDECIPN1,DRIPN1,DDRIPN1,STEP1,STEP2,DRAB,DDECB,
     CDRB,DWATCHRA,DWATCHDEC,DWATCHR,DCOMPTELRA,DCOMPTELDEC,DCOMPTELR,
     CDEGRETRA,DEGRETDEC,DEGRETR,DPHEBUSRA,DPHEBUSDEC,DPHEBUSR,
     CDSIGMARA,DSIGMADEC,DSIGMAR,DSAXRA,DSAXDEC,DSAXR,DHRAB,DHDECB,DHRB,
     CDKONUSEL1,DKONUSEL2,DPRA,DPDEC,DPR,DORAB,DODECB,DORB,AREA1,AREA2,
     CFRACTION,NAREA,RATEST,DECTEST,INANN,INOTHER)  
      NTREATED=NTREATED+1
      IF(INANN .AND. INOTHER)THEN
      NIN=NIN+1
      WRITE(*,190)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION
190   FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'BATSE#',1X,'BATSE RA',1X,'BATSEDEC',1X,'BATSE R ',1X,'RA IPN1 ',
     C1X,'DEC IPN1',1X,'R IPN1  ',1X,'DR IPN1 ',1X,'RA IPN2 ',1X,
     C'DEC IPN2'
     C,1X,'R IPN2  ',1X,'DR IPN2 ',1X,'WATCH RA',1X,'WATCHDEC',1X,
     C'WATCH R '
     C,1X,'CMPTL RA',1X,'CMPTLDEC',1X,'CMPTL R ',1X,'EGRET RA',1X,
     C'EGRETDEC',1X,'EGRET R ',1X,'PHEBS RA',1X,'PHEBSDEC',1X,
     C'PHEBUS R',1X,'SIGMA RA',1X,'SIGMADEC',1X,'SIGMA R ',1X,
     C'SAX RA  ',1X,'SAX DEC ',1X,'SAX R   ',1X,'HETE RA ',1X,
     C'HETE DEC',1X,
     C'HETE R  ',1X,'KONUSEL1',1X,'KONUSEL2',1X,'PLNET RA',1X,
     C'PLNETDEC',1X,'PLANET R',1X,'OTHER RA',1X,'OTHERDEC',1X,
     C'OTHER R ',1X,'STER',5X,'SQ.DEG.',2X,'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4,1X,I2.2,1H:,I2.2,1H:,I2.2,1X,I6,1X,
     C40(F8.4,1X),1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'TEST POINT IS WITHIN THE ERROR BOX'
      ENDIF
      IF(.NOT. INANN .AND. .NOT. INOTHER)THEN
      WRITE(*,190)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION	
      PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATIONS'
	ENDIF
      IF(INANN .AND. .NOT. INOTHER)THEN
      WRITE(*,190)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION	
      PRINT *,'TEST POINT LIES WITHIN THE ANNULI BUT DOES NOT LIE WITHIN
     C THE OTHER LOCALIZATIONS'
      ENDIF
      IF(.NOT. INANN .AND. INOTHER)THEN
      WRITE(*,190)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE ANNULUS BUT DOES LIE W
     CITHIN THE OTHER LOCALIZATIONS'
	ENDIF
      AREA1TOT=AREA1TOT+AREA1
      AREA2TOT=AREA2TOT+AREA2
      FRACTIONTOT=FRACTIONTOT+FRACTION
      IF(NAREA.EQ.0)THEN
C      WRITE(*,195)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
C     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
C     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
C     CCOMPTELDEC,
C     CCOMPTELR,EGRETRA,EGRETDEC,
C     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
C     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
C     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION
195   FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'BATSE#',1X,'BATSE RA',1X,'BATSEDEC',1X,'BATSE R ',1X,'RA IPN1 ',
     C1X,'DEC IPN1',1X,'R IPN1  ',1X,'DR IPN1 ',1X,'RA IPN2 ',1X,
     C'DEC IPN2'
     C,1X,'R IPN2  ',1X,'DR IPN2 ',1X,'WATCH RA',1X,'WATCHDEC',1X,
     C'WATCH R '
     C,1X,'CMPTL RA',1X,'CMPTLDEC',1X,'CMPTL R ',1X,'EGRET RA',1X,
     C'EGRETDEC',1X,'EGRET R ',1X,'PHEBS RA',1X,'PHEBSDEC',1X,
     C'PHEBUS R',1X,'SIGMA RA',1X,'SIGMADEC',1X,'SIGMA R ',1X,
     C'SAX RA  ',1X,'SAX DEC ',1X,'SAX R   ',1X,'HETE RA ',1X,
     C'HETE DEC',1X,
     C'HETE R  ',1X,'KONUSEL1',1X,'KONUSEL2',1X,'PLNET RA',1X,
     C'PLNETDEC',1X,'PLANET R',1X,'OTHER RA',1X,'OTHERDEC',1X,
     C'OTHER R ',1X,'STER',5X,'SQ.DEG.',2X,'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4.4,1X,I2.2,':',I2.2,':',I2.2,1X,
     CI6,1X,40(F8.4,1X),1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'THERE IS NO INTERSECTION BETWEEN THE IPN ANNULUS AND THE 
     COTHER LOCALIZATIONS; THE AREA OF THIS LOCALIZATION WILL BE TAKEN T
     CO BE THE AREA OF THE IPN ANNULUS'
      NNOINT=NNOINT+1
      IF(INANN)WRITE(*,196)
196   FORMAT(T2,'TEST POINT IS WITHIN THE IPN ANNULUS')
      IF(.NOT. INANN)WRITE(*,197)
197   FORMAT(T2,'TEST POINT IS NOT WITHIN THE IPN ANNULUS')
      IF(INOTHER)WRITE(*,198)
198   FORMAT(T2,'TEST POINT IS WITHIN THE OTHER LOCALIZATION')
      IF(.NOT. INOTHER)WRITE(*,199)
199   FORMAT(T2,'TEST POINT IS NOT WITHIN THE OTHER LOCALIZATION')
      ENDIF
      PRINT *,'--------------------------------------------------------'
      GOTO 10
200   CONTINUE

C.....CASE 14.  THERE ARE TWO IPN ANNULI AND NOTHING ELSE
      IF(RAIPN1.NE.0. .AND. DECIPN1.NE.0. .AND. RIPN1.NE.0. .AND.
     CDELRIPN1.NE.0. .AND.
     CRAIPN2.NE.0. .AND. DECIPN2.NE.0. .AND. RIPN2.NE.0. .AND.
     CDELRIPN2.NE.0. .AND.
     C(RAB.EQ.0. .AND. DECB.EQ.0. .AND. RB.EQ. 0.) .AND.
     C(WATCHRA.EQ.0. AND. WATCHDEC.EQ.0. .AND. WATCHR.EQ.0.).AND.
     C(COMPTELRA.EQ.0. .AND. COMPTELDEC.EQ.0. .AND. COMPTELR.EQ.0.).AND. 
     C(EGRETRA.EQ.0. AND. EGRETDEC.EQ.0. .AND. EGRETR.EQ.0.) .AND.
     C(PHEBUSRA.EQ.0. .AND. PHEBUSDEC.EQ.0. .AND. PHEBUSR.EQ.0.) .AND. 
     C(SIGMARA.EQ.0. AND. SIGMADEC.EQ.0. .AND. SIGMAR.EQ.0.) .AND.
     C(SAXRA.EQ.0. .AND. SAXDEC.EQ.0. .AND. SAXR.EQ.0.) .AND.
     C(HRAB.EQ.0. .AND. HDECB.EQ.0. .AND. HRB.EQ.0.) .AND.
     C(KONUSEL1.EQ.0. .AND. KONUSEL2.EQ.0.) .AND.
     C(PRA.EQ.0. .AND. PDEC.EQ.0. .AND. PR.EQ.0.) .AND.
     C(ORAB.EQ.0. .AND. ODECB.EQ.0. .AND. ORB.EQ.0.))GOTO 280
      GOTO 300
280   CONTINUE
C      PRINT *,NREAD,IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,RAB,
C     CDECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
C     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
C     CCOMPTELDEC,
C     CCOMPTELR,EGRETRA,EGRETDEC,
C     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
C     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
C     CORAB,ODECB,ORB           
      N2ANN=N2ANN+1
C.....STEP1 AND STEP2 ARE THE AZIMUTHAL AND RADIAL STEPS THAT
C.....WILL BE USED TO GO AROUND THE ANNULUS AND CALCULATE THE AREA.
C.....THE ANGULAR STEP IS PI*STEP	RADIANS
      STEP1=DDRIPN1/1000.D0
      STEP2=DDRIPN2/1000.D0
      STEP=DMIN1(STEP1,STEP2)
      STEP1=STEP
      STEP2=STEP
      AREA1=0.D0
      CALL GC14(DRAIPN1,DDECIPN1,DRIPN1,DDRIPN1,DRAIPN2,DDECIPN2,
     CDRIPN2,DDRIPN2,STEP1,STEP2,AREA1,AREA2,FRACTION,NAREA,RATEST,
     CDECTEST,IN)  
      NTREATED=NTREATED+1
 	IF(IN)THEN
	NIN=NIN+1
      WRITE(*,190)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION
	PRINT *,'TEST POINT LIES WITHIN THE TWO IPN ANNULI'
	ENDIF
      IF(.NOT. IN)THEN
      WRITE(*,190)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE TWO IPN ANNULI'
	ENDIF
      PRINT *,'--------------------------------------------------------'
      AREA1TOT=AREA1TOT+AREA1
      AREA2TOT=AREA2TOT+AREA2
      FRACTIONTOT=FRACTIONTOT+FRACTION
      IF(NAREA.EQ.0)THEN
      WRITE(*,195)IDAY,MON,IYR,IHR,IMIN,ISEC
      NNOINT2=NNOINT2+1
      ENDIF
      GOTO 10
300   CONTINUE

C.....CASE 15.  THERE ARE TWO IPN ANNULI AND ANYTHING ELSE
      IF(RAIPN1.NE.0. .AND. DECIPN1.NE.0. .AND. RIPN1.NE.0. .AND.
     CDELRIPN1.NE.0.) THEN
      IF(RAIPN2.NE.0. .AND. DECIPN2.NE.0. .AND. RIPN2.NE.0. .AND.
     CDELRIPN2.NE.0.) THEN
      IF((RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE. 0.) .OR.
     C(WATCHRA.NE.0. AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.).OR.
     C(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.).OR. 
     C(EGRETRA.NE.0. AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.) .OR.
     C(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.) .OR. 
     C(SIGMARA.NE.0. AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.) .OR.
     C(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.) .OR.
     C(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.) .OR.
     C(KONUSEL1.NE.0. .OR. KONUSEL2.NE.0.) .OR.
     C(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.) .OR.
     C(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.))GOTO 380
      ENDIF
      ENDIF
      GOTO 400
380   CONTINUE
C      PRINT *,NREAD,IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,RAB,
C     CDECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
C     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
C     CCOMPTELDEC,
C     CCOMPTELR,EGRETRA,EGRETDEC,
C     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
C     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
C     CORAB,ODECB,ORB           
      N2ANNAND=N2ANNAND+1
C.....STEP1 AND STEP2 ARE THE AZIMUTHAL AND RADIAL STEPS THAT
C.....WILL BE USED TO GO AROUND THE ANNULUS AND CALCULATE THE AREA.
C.....THE ANGULAR STEP IS PI*STEP	RADIANS
      STEP1=DDRIPN1/1000.D0
      STEP2=DDRIPN2/1000.D0
      STEP=DMIN1(STEP1,STEP2)
      STEP1=STEP
      STEP2=STEP
      AREA1=0.D0
      CALL GC15(DRAIPN1,DDECIPN1,DRIPN1,DDRIPN1,DRAIPN2,DDECIPN2,
     CDRIPN2,DDRIPN2,DRAB,DDECB,
     CDRB,DWATCHRA,DWATCHDEC,DWATCHR,DCOMPTELRA,DCOMPTELDEC,DCOMPTELR,
     CDEGRETRA,DEGRETDEC,DEGRETR,DPHEBUSRA,DPHEBUSDEC,DPHEBUSR,
     CDSIGMARA,DSIGMADEC,DSIGMAR,DSAXRA,DSAXDEC,DSAXR,DHRAB,DHDECB,DHRB,
     CDKONUSEL1,DKONUSEL2,DPRA,DPDEC,DPR,DORAB,DODECB,DORB,
     CSTEP1,STEP2,AREA1,AREA2,FRACTION,NAREA,RATEST,DECTEST,INANN1,
     CINANN2,INOTHER)  
      NTREATED=NTREATED+1
      IF(INANN1 .AND. INANN2 .AND. INOTHER)THEN
      NIN=NIN+1
      WRITE(*,190)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION
	PRINT *,'TEST POINT LIES WITHIN THE ANNULI AND THE OTHER LOCALIZAT
     CION'
      ENDIF
	IF(NAREA.NE.0)THEN
      WRITE(*,190)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION
      IF(.NOT.INANN1)PRINT *,'TEST POINT DOES NOT LIE WITHIN THE FIRST A
     CNNULUS'
      IF(.NOT.INANN2)PRINT *,'TEST POINT DOES NOT LIE WITHIN THE SECOND 
     CANNULUS'
      IF(.NOT.INOTHER)PRINT *,'TEST POINT DOES NOT LIE WITHIN THE OTHER 
     CLOCALIZATION'
      ENDIF
      AREA1TOT=AREA1TOT+AREA1
      AREA2TOT=AREA2TOT+AREA2
      FRACTIONTOT=FRACTIONTOT+FRACTION
      IF(NAREA.EQ.0)THEN
C.....GET THE AREA OF THE INTERSECTION OF THE TWO ANNULI
      CALL GC14(DRAIPN1,DDECIPN1,DRIPN1,DDRIPN1,DRAIPN2,DDECIPN2,
     CDRIPN2,DDRIPN2,STEP1,STEP2,AREA1,AREA2,FRACTION,NAREA,RATEST,
     CDECTEST,IN)  
      AREA1TOT=AREA1TOT+AREA1
      AREA2TOT=AREA2TOT+AREA2
      FRACTIONTOT=FRACTIONTOT+FRACTION
      WRITE(*,190)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION
      IF(.NOT.INANN1)PRINT *,'TEST POINT DOES NOT LIE WITHIN THE FIRST A
     CNNULUS'
      IF(.NOT.INANN2)PRINT *,'TEST POINT DOES NOT LIE WITHIN THE SECOND 
     CANNULUS'
      IF(.NOT.INOTHER)PRINT *,'TEST POINT DOES NOT LIE WITHIN THE OTHER 
     CLOCALIZATION'
C      WRITE(*,395)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
C     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
C     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
C     CCOMPTELDEC,
C     CCOMPTELR,EGRETRA,EGRETDEC,
C     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
C     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
C     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION
395   FORMAT(T2,'DD',1X,'MM',1X,'YYYY',1X,'HH',1X,'MM',1X,'SS',1X,
     C'BATSE#',1X,'BATSE RA',1X,'BATSEDEC',1X,'BATSE R ',1X,'RA IPN1 ',
     C1X,'DEC IPN1',1X,'R IPN1  ',1X,'DR IPN1 ',1X,'RA IPN2 ',1X,
     C'DEC IPN2'
     C,1X,'R IPN2  ',1X,'DR IPN2 ',1X,'WATCH RA',1X,'WATCHDEC',1X,
     C'WATCH R '
     C,1X,'CMPTL RA',1X,'CMPTLDEC',1X,'CMPTL R ',1X,'EGRET RA',1X,
     C'EGRETDEC',1X,'EGRET R ',1X,'PHEBS RA',1X,'PHEBSDEC',1X,
     C'PHEBUS R',1X,'SIGMA RA',1X,'SIGMADEC',1X,'SIGMA R ',1X,
     C'SAX RA  ',1X,'SAX DEC ',1X,'SAX R   ',1X,'HETE RA ',1X,
     C'HETE DEC',1X,
     C'HETE R  ',1X,'KONUSEL1',1X,'KONUSEL2',1X,'PLNET RA',1X,
     C'PLNETDEC',1X,'PLANET R',1X,'OTHER RA',1X,'OTHERDEC',1X,
     C'OTHER R ',2X,'STER',5X,'SQ.DEG.'2X,
     C'FRAC.OF SKY',/,
     CT2,I2.2,1X,I2.2,1X,I4.4,1X,I2.2,':',I2.2,':',I2.2,1X,
     CI6,1X,40(F8.4,1X),1PE8.2,1X,1PE8.2,1X,1PE8.2)
      PRINT *,'THERE IS NO INTERSECTION BETWEEN THE IPN ANNULI AND THE O
     CTHER LOCALIZATIONS; THE AREA OF THIS LOCALIZATION WILL BE TAKEN TO 
     C BE THE AREA OF THE INTERSECTIONS OF THE IPN ANNULI'
      NNOINT2AND=NNOINT2AND+1
      ENDIF
      PRINT *,'--------------------------------------------------------'
      GOTO 10
400   CONTINUE
C.....CASE 16. LAST CASE.  THERE ARE NO IPN ANNULI.  THE LOCALIZATION CONSISTS
C.....OF AT LEAST TWO OF ANYTHING ELSE, LIKE:
C.....KONUS EL + EARTH BLOCKING
C.....OTHER LOCALIZATION + SAX CIRCLE
C.....EARTH BLOCKING + OTHER LOCALIZATION
C.....ETC.
      IF(RAIPN1.EQ.0. .AND. DECIPN1.EQ.0. .AND. RIPN1.EQ.0. .AND.
     CDELRIPN1.EQ.0.) THEN
      IF(RAIPN2.EQ.0. .AND. DECIPN2.EQ.0. .AND. RIPN2.EQ.0. .AND.
     CDELRIPN2.EQ.0.) THEN
C.....FIND THE SMALLEST ERROR CIRCLE IN THE LOCALIZATIONS; IGNORE KONUS
C.....ECLIPTIC LATITUDES FOR NOW
      RMIN=1.E9
      KON=.FALSE.
      NLOC=0
      IF(RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE. 0.) THEN
      NLOC=NLOC+1
      IF(RB.LT.RMIN)THEN
      RACIRC=DRAB
      DECCIRC=DDECB
      RCIRC=DRB
      RMIN=RCIRC
      ENDIF
      ENDIF
      IF(WATCHRA.NE.0. AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.) THEN
      NLOC=NLOC+1
      IF(WATCHR.LT.RMIN)THEN
      RACIRC=DWATCHRA
      DECCIRC=DWATCHDEC
      RCIRC=DWATCHR
      RMIN=RCIRC
      ENDIF
      ENDIF
      IF(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.)	
     CTHEN
      NLOC=NLOC+1
      IF(COMPTELR.LT.RMIN)THEN
      RACIRC=DCOMPTELRA
      DECCIRC=DCOMPTELDEC
      RCIRC=DCOMPTELR
      RMIN=RCIRC
      ENDIF
      ENDIF	 
      IF(EGRETRA.NE.0. AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.)THEN 
      NLOC=NLOC+1
      IF(EGRETR.LT.RMIN)THEN
      RACIRC=DEGRETRA
      DECCIRC=EGRETDEC
      RCIRC=DEGRETR
      RMIN=RCIRC
      ENDIF
      ENDIF
      IF(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.)THEN 
      NLOC=NLOC+1
      IF(PHEBUSR.LT.RMIN)THEN
      RACIRC=DPHEBUSRA
      DECCIRD=DPHEBUSDEC
      RCIRC=DPHEBUSR
      RMIN=RCIRC
      ENDIF	 
      ENDIF
      IF(SIGMARA.NE.0. AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.)THEN 
      NLOC=NLOC+1
      IF(SIGMAR.LT.RMIN)THEN
      RACIRC=DSIGMARA
      DECCIRC=SIGMADEC
      RCIRC=DSIGMAR
      RMIN=RCIRC
      ENDIF
      ENDIF
      IF(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.) THEN
      NLOC=NLOC+1
      IF(SAXR.LT.RMIN)THEN
      RACIRC=DSAXRA
      DECCIRC=DSAXDEC
      RCIRC=DSAXR
      RMIN=RCIRC
      ENDIF
      ENDIF
      IF(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.)THEN 
      NLOC=NLOC+1
      IF(HRB.LT.RMIN)THEN
      RACIRC=DHRAB
      DECCIRC=DHDECB
      RCIRC=DHRB
      RMIN=RCIRC
      ENDIF
      ENDIF
C.....DON'T USE PLANET-BLOCKING TO DEFINE RMIN!
      IF(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.) THEN
      NLOC=NLOC+1
      ENDIF
      IF(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.) THEN
      NLOC=NLOC+1
      IF(ORB.LT.RMIN)THEN
      RACIRC=DORAB
      DECCIRC=DODECB
      RCIRC=DORB
      RMIN=RCIRC
      ENDIF
      ENDIF
C.....IF RMIN IS STILL NOT DEFINED, USE THE KONUS EL TO DEFINE IT
      IF(KONUSEL1.NE.0. .OR. KONUSEL2.NE.0.)THEN
      NLOC=NLOC+1
      IF(RMIN.GT.1.E8)THEN
      RACIRC=270.0
      DECCIRC=66.561
      RMIN=10.
      KON=.TRUE.
      ENDIF
      ENDIF
      IF(NLOC.LT.2)THEN
C      PRINT *,NREAD,IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,RAB,
C     CDECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
C     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
C     CCOMPTELDEC,
C     CCOMPTELR,EGRETRA,EGRETDEC,
C     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
C     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
C     CORAB,ODECB,ORB         
      STOP 'NLOC'
      ENDIF
      GOTO 410
      ENDIF
      ENDIF
C      PRINT *,NREAD,IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,RAB,
C     CDECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
C     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
C     CCOMPTELDEC,
C     CCOMPTELR,EGRETRA,EGRETDEC,
C     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
C     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
C     CORAB,ODECB,ORB         
      STOP 'NO CASE 16'
410	N2OTHER=N2OTHER+1
C.....STEP1 AND STEP2 ARE THE AZIMUTHAL AND RADIAL STEPS THAT
C.....WILL BE USED TO GO AROUND THE SMALLEST LOCALIZATION CIRCLE
C.....AND CALCULATE THE AREA.
C.....THE ANGULAR STEP IS PI*STEP	RADIANS
      STEP1=RMIN/1000.D0
      STEP2=STEP1
      AREA1=0.D0
C	PRINT *,IDAY,MON,IYR,IHR,IMIN,ISEC
      CALL GC16(RACIRC,DECCIRC,RCIRC,DRAB,DDECB,
     CDRB,DWATCHRA,DWATCHDEC,DWATCHR,DCOMPTELRA,DCOMPTELDEC,DCOMPTELR,
     CDEGRETRA,DEGRETDEC,DEGRETR,DPHEBUSRA,DPHEBUSDEC,DPHEBUSR,
     CDSIGMARA,DSIGMADEC,DSIGMAR,DSAXRA,DSAXDEC,DSAXR,DHRAB,DHDECB,DHRB,
     CDKONUSEL1,DKONUSEL2,DPRA,DPDEC,DPR,DORAB,DODECB,DORB,
     CSTEP1,STEP2,AREA1,AREA2,FRACTION,NAREA,KON,RATEST,DECTEST,IN)  
      NTREATED=NTREATED+1
      IF(IN)THEN
      NIN=NIN+1
      WRITE(*,190)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION
	PRINT *,'TEST POINT LIES WITHIN THE LOCALIZATIONS'
      ENDIF
      IF(.NOT. IN)THEN
      WRITE(*,190)IDAY,MON,IYR,IHR,IMIN,ISEC,IBNO,
     CRAB,DECB,RB,RAIPN1,DECIPN1,RIPN1,DELRIPN1,RAIPN2,
     CDECIPN2,RIPN2,DELRIPN2,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,
     CCOMPTELDEC,
     CCOMPTELR,EGRETRA,EGRETDEC,
     CEGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,SIGMARA,SIGMADEC,SIGMAR,
     CSAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,KONUSEL1,KONUSEL2,PRA,PDEC,PR,
     CORAB,ODECB,ORB,AREA1,AREA2,FRACTION
	PRINT *,'TEST POINT DOES NOT LIE WITHIN THE LOCALIZATIONS'
      ENDIF
      AREA1TOT=AREA1TOT+AREA1
      AREA2TOT=AREA2TOT+AREA2
      FRACTIONTOT=FRACTIONTOT+FRACTION
      IF(NAREA.EQ.0)THEN
      WRITE(*,415)IDAY,MON,IYR,IHR,IMIN,ISEC
415   FORMAT(I2.2,'/',I2.2,'/',I4.4,1X,I2.2,':',I2.2,':',I2.2,1X,
     C'THERE IS NO INTERSECTION BETWEEN THE LOCALIZATIONS')
      NNOINTN2=NNOINTN2+1
      ENDIF
      PRINT *,'--------------------------------------------------------'
      GOTO 10
4000  CONTINUE
      BURSTSPERDAY=NTREATED/(DJD21-DJD11+1)
      WRITE(*,4100)NREAD,NTREATED,BURSTSPERDAY,AREA1TOT,AREA2TOT,
     CFRACTIONTOT
4100  FORMAT(///,T2,'SUMMARY:',/,
     CT2,'THERE WERE A TOTAL OF ',I4,' LOCALIZATIONS IN THE DATABASE',/
     C,I4,' LOCALIZATIONS FELL WITHIN THE SPECIFIED TIME PERIOD (',F5.2
     C,' BURSTS PER DAY)',/,
     C'THEY OCCUPIED ',1PE8.2,' TOTAL STERADIANS',/,
     C'OR ',1PE8.2,' TOTAL SQUARE DEGREES',/,
     C'THIS COVERS A FRACTION OF THE TOTAL SKY EQUAL TO ',1PE8.2,' SO TH
     CIS IS THE NUMBER OF CHANCE ASSOCIATIONS EXPECTED')
      WRITE(*,4150)NIN
4150  FORMAT('THE POSITION TESTED FELL WITHIN ',I3,' ERROR BOXES')
      WRITE(*,4200) N1ANN,NEL,NBAT,NSAX,NPHE,NSIG,NWAT,NEGR,
     CNCOM,NHET,NPLAN,NOTHER
4200  FORMAT(///,T2,'SOME STATISTICS:',/,
     CT2,'NUMBER OF LOCALIZATIONS WHERE THERE IS ONLY:',/,
     C'A SINGLE IPN ANNULUS ',I4,/,
     C'AN ECLIPTIC LATITUDE BAND ',I4,/,
     C'A BATSE/FERMI ERROR CIRCLE ',I4,/,
     C'A BEPPOSAX ERROR CIRCLE ',I4,/,
     C'A PHEBUS ERROR CIRCLE ',I4,/,
     C'A SIGMA ERROR CIRCLE ',I4,/,
     C'A WATCH ERROR CIRCLE ',I4,/,
     C'AN EGRET ERROR CIRCLE ',I4,/,
     C'A COMPTEL ERROR CIRCLE ',I4,/,
     C'A HETE ERROR CIRCLE ',I4,/,
     C'PLANET-BLOCKING ',I4,/,
     C'AN OTHER ERROR CIRCLE ',I4)
      WRITE(*,4300)N1ANNAND
4300  FORMAT(/,T2,'NUMBER OF LOCALIZATIONS WHERE THERE IS ONE AND ONLY
     CONE IPN ANNULUS AND ANYTHING ELSE IS ',I4)
      WRITE(*,4310)NNOINT
4310  FORMAT(T2,'OF THEM, THE NUMBER WHERE THERE IS NO INTERSECTION IS 
     C ',I4)
      WRITE(*,4320)N2ANN
4320  FORMAT(/,T2,'NUMBER OF LOCALIZATIONS WHERE THERE ARE ONLY TWO IPN
     C ANNULI AND NOTHING ELSE IS ',I4)
      WRITE(*,4330)NNOINT2
4330  FORMAT(T2,'OF THEM, THE NUMBER WHERE THERE IS NO INTERSECTION IS
     C ',I4)
      WRITE(*,4340)N2ANNAND
4340  FORMAT(/,T2,'NUMBER OF LOCALIZATIONS WHERE THERE ARE TWO IPN ANNUL
     CI AND OTHER LOCALIZATION INFORMATION IS ',I4)
      WRITE(*,4350)NNOINT2AND
4350  FORMAT(T2,'OF THEM, THE NUMBER WHERE THERE IS NO INTERSECTION IS
     C ',I4)
      WRITE(*,4360)N2OTHER
4360  FORMAT(/,T2,'NUMBER OF LOCALIZATIONS WHERE THERE ARE NO IPN ANNULI
     C AND TWO OR MORE OTHER LOCALIZATIONS IS ',I4)
      WRITE(*,4330)NNOINTN2
      REWIND 1
      WRITE(*,4400)
4400  FORMAT('******************************************************END*
     C***************************************************************')
      WRITE(12,4500)NAME,NTREATED,FRACTIONTOT,NIN
	FRACTIONTOT1=FRACTIONTOT1+FRACTIONTOT
	NIN1=NIN1+NIN
4500  FORMAT(A10,/,I4,' EVENTS TREATED', /,D8.3,' EVENTS EXPECTED',/,I4
     C' EVENTS CONSISTENT WITH GIVEN POSITION')
C.....READ SEPARATOR
C      READ(10,4,END=5000)SEPARATOR
C.....UPDATE RATEST,DECTEST
C     	GOTO 1
5000  WRITE(12,4600)FRACTIONTOT1,NIN1
4600  FORMAT('TOTAL EVENTS EXPECTED',10X,D7.3,1X,I5,' EVENTS FOUND')
      OPEN(13,FILE='TOTALNUMBER',ACTION='WRITE')
      WRITE(13,*)NREAD
      CLOSE(13)
      STOP
      END
      SUBROUTINE INIT
C.....OPEN FILES
C      OPEN(1,FILE='I:/FIL/ACCESS/chronological2.csv',STATUS='OLD')
      OPEN(1,FILE='chronological2.csv',STATUS='OLD')
      OPEN(2,FILE='ipnpositions',ACTION='WRITE',STATUS='REPLACE')
      OPEN(12,FILE='summary',ACTION='WRITE',STATUS='REPLACE')
      OPEN(13,FILE='TOTALNUMBER',ACTION='READ')
      RETURN
      END
      SUBROUTINE ANNULUS(RAIPN,DECIPN,R,DR,AREA1,AREA2,FRACTION,RATEST,
     CDECTEST,IN)
C.....CALCULATES THE AREA OF AN ANNULUS IN STERADIANS, SQUARE DEGREES,
C.....AND FRACTION OF THE SKY.  ANNULUS RADIUS IS R +/- DR
      COMMON/GEOM/PI,RADEG
      DOUBLE PRECISION PI,RADEG,R,DR,AREA1,AREA2,FRACTION,RAIPN,
     CDECIPN,R1,R2,RATEST,DECTEST
	LOGICAL IN
      AREA1=2.D0*PI*(DCOS((R-DR)*RADEG)-DCOS((R+DR)*RADEG))
      AREA2=AREA1/(RADEG**2)
      FRACTION=AREA1/(4.D0*PI)
	R1=R-DR
	R2=R+DR
      CALL ANGLE(RAIPN,DECIPN,RATEST,DECTEST,R1,R2,DIST,IN)
	RETURN
      END
      SUBROUTINE CIRCLE(DRAB,DDECB,R,AREA1,AREA2,FRACTION,RATEST,
     CDECTEST,IN)
C.....CALCULATES THE AREA OF A SMALL CIRCLE IN STERADIANS, SQUARE DEGREES,
C.....AND FRACTION OF THE SKY.  
      COMMON/GEOM/PI,RADEG
      DOUBLE PRECISION PI,RADEG,R,AREA1,AREA2,FRACTION,RATEST,DECTEST,
     CDRAB,DDECB,R1,R2
	LOGICAL IN
      AREA1=2.D0*PI*(1.-DCOS(R*RADEG))
      AREA2=AREA1/(RADEG**2)
      FRACTION=AREA1/(4.D0*PI)
	R1=0.D0
	R2=R
      CALL ANGLE(DRAB,DDECB,RATEST,DECTEST,R1,R2,DIST,IN)
      RETURN
      END

      SUBROUTINE ANGLE(RA1,DEC1,RA2,DEC2,R1,R2,DIST,IN)
C.....FINDS THE ANGULAR DISTANCE DIST BETWEEN POINTS RA1,DEC1, AND RA2,DEC2.
C.....IF R1<DIST<R2, POINT IS WITHIN THE ERROR ANNULUS, AND IN=.TRUE.
      LOGICAL IN
	DOUBLE PRECISION PI,RD,RA1,DEC1,RA2,DEC2,R1,R2
      COMMON/GEOM/PI,RD
	IN=.FALSE.
      DIST=DACOS(DSIN(DEC1*RD)*DSIN(DEC2*RD) +
     CDCOS(DEC1*RD)*DCOS(DEC2*RD)*DCOS((RA1-RA2)*RD))
      DIST=DIST/RD
      IF(DIST.GE.R1 .AND. DIST.LE.R2)IN=.TRUE.
      IF(DIST.GE.R1 .AND. DIST.LE.R2)THEN
C	PRINT *,R1,R2,DIST
	ENDIF
      RETURN
      END
      SUBROUTINE FIND (IJ,I,IK,IL,*,*)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      DOUBLE PRECISION NUM1,NUM2
      COMMON/INTER/RA(10,10,10,10,10),DEC(10,10,10,10,10),THETA(10,10),
     CDTHETA(10,10)
      COMMON/VECT/V(10,3),DV(10,3),VDIFF(10,10,3),DT(10,10),DDT(10,10)
      COMMON/ITOT/ITOT
      RD=3.1415926536/180.
      PI=3.1415926536
      P2=2*PI
C.....FINDS 2 INTERSECTION POINTS OF 2 SMALL CIRCLES (CIRCLES IJ TO I AND IJ
C.....TO IK) ON CELESTIAL SPHERE. FIRST RETURN IS NORMAL PATH, 2ND IS FOR
C.....CASE WHERE INTERSECTION CAN'T BE CALCULATED
C.....REFER TO "SIGNE 2MP COMPUTER PROGRAMS" BINDER FOR DETAILS
      IF(DTHETA(IJ,I).EQ.0. .OR. DTHETA(IK,IL).EQ.0.) GOTO 1
      GOTO 5
1     WRITE(*,2) DTHETA(IJ,I),DTHETA(IK,IL)
2     FORMAT('\r',T2'S/R FIND CALLED FOR 'I3'-'I3','I3'-'I3' ',
     C'INTERSECTION. NO DELTAT GIVEN FOR ONE OR BOTH. CANT CALCULATE') 
C      pause
      RETURN 2
5     CONTINUE
C.....RA AND DEC OF CENTERS OF CIRCLES
      AO=VDIFF(I,IJ,2)*RD
      AOP=VDIFF(IK ,IL,2)*RD
      DO=VDIFF(I,IJ,3)*RD
      DOP=VDIFF(IK ,IL,3)*RD 
C.....IJ AND I REFER TO THE IPN ANNULUS.  DO THE CALCULATION FOR BOTH
C.....ANNULI.
C.....IK AND IL REFER TO THE BATSE ERROR CIRCLE.  DO THE CALCULATION
C.....ONLY FOR THE OUTER ANNULUS  
      DO 200 K=2,3
      GO TO(10,20,30,40,50) K
10    CHI=THETA(I,IJ)
      CHIP=THETA(IK ,IL)
      GO TO 60
20    CHI=THETA(I,IJ)+DTHETA(I,IJ)
      CHIP=THETA(IK ,IL)+DTHETA(IK ,IL)
      GO TO 60
30    CHI=THETA(I,IJ)-DTHETA(I,IJ)
      CHIP=THETA(IK ,IL)+DTHETA(IK ,IL)
      GO TO 60
40    CHI=THETA(I,IJ)+DTHETA(I,IJ)
      CHIP=THETA(IK ,IL)-DTHETA(IK ,IL)
      GO TO 60
50    CHI=THETA(I,IJ)-DTHETA(I,IJ)
      CHIP=THETA(IK ,IL)-DTHETA(IK ,IL)
60    CONTINUE
C.....SINES AND COSINES OF ANGLES
      COSDO=DCOS(DO)
      COSDOP=DCOS(DOP)
      SINDO=DSIN(DO)
      SINDOP=DSIN(DOP)
      COSCHI=DCOS(CHI)
      COSCHIP=DCOS(CHIP)
      COSAOAP=DCOS(AO-AOP)
      SINAOAP=DSIN(AO-AOP)
C.....COEFFICIENTS OF QUADRATIC EQUATION
      A=-(COSDO*COSDOP*SINAOAP)**2
     C+2*SINDO*SINDOP*COSDO*COSDOP*COSAOAP
     C-(SINDOP*COSDO)**2   -   (SINDO*COSDOP)**2
      B=-2*(COSCHI*SINDOP + COSCHIP*SINDO)*COSDO*COSDOP*COSAOAP
     C+2*COSCHIP*SINDOP*COSDO**2 + 2*COSCHI*SINDO*COSDOP**2
      C=(COSDO*COSDOP*SINAOAP)**2+ 2*COSCHI*COSCHIP*COSDO*COSDOP*COSAOAP
     C-(COSCHIP*COSDO)**2 - (COSCHI*COSDOP)**2
      B24AC=B**2-4*A*C    
      IF(B24AC.LT.0. .AND. ITOT.GT.2) GOTO70
      GOTO 90
70    WRITE(*,71) IJ,I,IK,IL
71    FORMAT('\r',T2'MISSING OR GRAZING INTERSECTION BETWEEN ',
     C''I3'-'I3' ANNULUS AND 'I3'-'I3' ANNULUS')
      RETURN 2
C.....IF ITOT=2 AND THERE IS NO INTERSECTION POINT, DEFINE
C.....THE INTERSECTION AS THE CENTER OF THE BATSE CIRCLE
90    IF(B24AC.LT.0. .AND. ITOT.EQ.2) THEN 
C****      R1=VDIFF(3,4,2)
C****      R2=VDIFF(3,4,2)
C****      D1=VDIFF(3,4,3)
C****      D2=VDIFF(3,4,3)
      RETURN 2
C****      GOTO 100
      ENDIF
      NUM1=-B+DSQRT(B24AC)
      NUM2=-B-DSQRT(B24AC)
      SIND1=NUM1/(2*A)
      SIND2=NUM2/(2*A)
      D1=DASIN(SIND1)
      D2=DASIN(SIND2)
      COSD1=DCOS(D1)
      COSD2=DCOS(D2)
      ACOS1=DACOS((COSCHI-SIND1*SINDO)/(COSD1*COSDO))
      ACOS2=DACOS((COSCHI-SIND2*SINDO)/(COSD2*COSDO))
      A1NEG=MOD(AO-ACOS1+P2,P2)
      A1POS=MOD(AO+ACOS1+P2,P2)
      A2NEG=MOD(AO-ACOS2+P2,P2)
      A2POS=MOD(AO+ACOS2+P2,P2)
C.....FOR A CIRCLE, THERE ARE TWO ALPHAS (E.G. A1NEG,A1POS) FOR A GIVEN
C.....DELTA(D1).  THE FOLLOWING CALCULATION ASSURES THAT THE ALPHA TAKEN FOR
C.....THE INTERSECTION POINT LIES ON THE SECOND (PRIMED) CIRCLE
      ACOSP1=DACOS((COSCHIP-SIND1*SINDOP)/(COSD1*COSDOP))
      ACOSP2=DACOS((COSCHIP-SIND2*SINDOP)/(COSD2*COSDOP))
      A1PNEG=MOD(AOP-ACOSP1+P2,P2)
      A1PPOS=MOD(AOP+ACOSP1+P2,P2)
      A2PNEG=MOD(AOP-ACOSP2+P2,P2)
      A2PPOS=MOD(AOP+ACOSP2+P2,P2)
C;....INTERSECTION POINT 1
      IF(ABS(A1NEG-A1PNEG).LT.1.D-5) R1=A1NEG/RD
      IF(ABS(A1NEG-A1PPOS).LT.1.D-5) R1=A1NEG/RD
      IF(ABS(A1POS-A1PNEG).LT.1.D-5)R1=A1POS/RD
      IF(ABS(A1POS-A1PPOS).LT.1.D-5)R1=A1POS/RD
C)....INTERSECTION POINT 2
      IF(ABS(A2NEG-A2PNEG).LT.1.D-5) R2=A2NEG/RD
      IF(ABS(A2NEG-A2PPOS).LT.1.D-5) R2=A2NEG/RD
      IF(ABS(A2POS-A2PNEG).LT.1.D-5)R2=A2POS/RD
      IF(ABS(A2POS-A2PPOS).LT.1.D-5)R2=A2POS/RD
      D1=D1/RD
      D2=D2/RD
100   K1=2*K-1
      K2=2*K
C****      CALL ABER(R1,D1,R1,D1)
C****      CALL ABER(R2,D2,R2,D2)
      RA(K1,I,IJ,IK,IL)=R1
      DEC(K1,I,IJ,IK,IL)=D1
      RA(K1,IJ,I,IK,IL)=R1
      DEC(K1,IJ,I,IK,IL)=D1
      RA(K1,I,IJ,IL,IK)=R1
      DEC(K1,I,IJ,IL,IK)=D1
      RA(K1,IJ,I,IL,IK)=R1
      DEC(K1,IJ,I,IL,IK)=D1
      RA(K1,IK,IL,I,IJ)=R1
      DEC(K1,IK,IL,I,IJ)=D1
      RA(K1,IL,IK,I,IJ)=R1
      DEC(K1,IL,IK,I,IJ)=D1
      RA(K1,IL,IK,IJ,I)=R1
      DEC(K1,IL,IK,IJ,I)=D1
      RA(K1,IK,IL,IJ,I)=R1
      DEC(K1,IK,IL,IJ,I)=D1
      RA(K2,I,IJ,IK,IL)=R2
      DEC(K2,I,IJ,IK,IL)=D2
      RA(K2,IJ,I,IK,IL)=R2
      DEC(K2,IJ,I,IK,IL)=D2
      RA(K2,I,IJ,IL,IK)=R2
      DEC(K2,I,IJ,IL,IK)=D2
      RA(K2,IJ,I,IL,IK)=R2
      DEC(K2,IJ,I,IL,IK)=D2
      RA(K2,IK,IL,I,IJ)=R2
      DEC(K2,IK,IL,I,IJ)=D2
      RA(K2,IL,IK,I,IJ)=R2
      DEC(K2,IL,IK,I,IJ)=D2
      RA(K2,IL,IK,IJ,I)=R2
      DEC(K2,IL,IK,IJ,I)=D2
      RA(K2,IK,IL,IJ,I)=R2
      DEC(K2,IK,IL,IJ,I)=D2
200   CONTINUE
      RETURN 1
      END
      SUBROUTINE AREA(RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE)
C.....CALCULATES THE AREA OF A 4 SIDED POLYGON IN SPHERICAL TRIGONOMETRY BY
C.....BREAKING IT DOWN INTO TWO TRIANGLES	WITH EQUAL AREAS. RA,DEC DEFINE
C.....THE CORNERS, IN DEGREES. ARE IS IN SQUARE ARCMINUTES.  THIS IS
C.....FOR A POLYGON DEFINED BY SEGMENTS OF GREAT CIRCLES
      DOUBLE PRECISION ANG1,ANG2,ANG3,RA1,RA2,RA3,DEC1,DEC2,
     CDEC3,PI,RD,ARE
      COMMON/GEOM/PI,RD
C.....SOLVE FOR THE INNER ANGLES OF THE TRIANGLE
      CALL SPHANGLE(RA1,DEC1,RA2,DEC2,RA3,DEC3,ANG1)
	CALL SPHANGLE(RA2,DEC2,RA3,DEC3,RA1,DEC1,ANG2)
	CALL SPHANGLE(RA3,DEC3,RA1,DEC1,RA2,DEC2,ANG3)
	    PRINT *,'ANG1,ANG2,ANG3=',ANG1,ANG2,ANG3
      ARE=2*3600.D0*3282.8D0*((ANG1+ANG2+ANG3)*RD-PI)
      RETURN
      END
      SUBROUTINE SPHANGLE(A,B,C,D,E,F,ANG)
C.....SOLVES THE SPHERICAL TRIANGLE GIVEN BY THE 3 COORDINATES
C.....(A,B), (C,D), (E,F) IN DEGREES FOR THE ANGLE (A,B)/(C,D)/(E,F),
C.....ALSO IN DEGREES
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      COMMON/GEOM/PI,RD
C.....FIRST CALCULATE THE ANGLES BETWEEN (A,B) AND (C,D), (C,D) AND
C.....(E,F), AND (E,F) AND (A,B)
      CALL ANGLE1(A,B,C,D,A1)
      CALL ANGLE1(C,D,E,F,A2)
      CALL ANGLE1(E,F,A,B,A3)
	     PRINT *,'A1,A2,A3= ',A1,A2,A3
C      COSC=(DCOS(A2*RD)-DCOS(A1*RD)*DCOS(A3*RD))/
C     C(DSIN(A1*RD)*DSIN(A3*RD))
      COSC=(DCOS(A3*RD)-DCOS(A1*RD)*DCOS(A2*RD))/
     C(DSIN(A1*RD)*DSIN(A2*RD))
      ANG=DACOS(COSC)/RD
      RETURN
      END
      SUBROUTINE ANGLE1(A1,D1,A2,D2,ANG)
C.....CALCULATES THE ANGLE ANG (IN DEGREES) BETWEEN TWO POINTS ON THE CELESTIAL SPHERE
C.....A1,D1 AND A2,D2 (IN DEGREES)
      DOUBLE PRECISION ANG,A1,D1,A2,D2,PI,RD,ARG
      COMMON/GEOM/PI,RD
      ARG=DSIN(D1*RD)*DSIN(D2*RD)+DCOS(D1*RD)*DCOS(D2*RD)
     C*DCOS(A1*RD-A2*RD)
	IF(ARG.GT.1.D0 .OR. ARG.LT.-1.D0)THEN
	ANG=0.D0
	RETURN
	ENDIF
      ANG=DACOS(DSIN(D1*RD)*DSIN(D2*RD)+DCOS(D1*RD)*DCOS(D2*RD)
     C*DCOS(A1*RD-A2*RD))
      ANG=ANG/RD
      RETURN
      END
      SUBROUTINE JUL(DY,IMON,IYR,DJD)
C.....DJD IS THE NUMBER OF JULIAN DAYS SINCE 1900 JAN 0.5	(2415020)
C.....S/R IS FROM "ASTRONOMY WITH YOUR PERSONAL COMPUTER" BY P.
C.....DUFFETT-SMITH
      REAL MN1
      DOUBLE PRECISION DJD,C
      INTEGER YR1
      FNITG(W)=SIGN(1.,W)*INT(ABS(W))
1110  MN1=IMON
      YR1=IYR
      B=0
1115  IF (YR1.LT.0) YR1=YR1+1
1120  IF (IMON.LT.3) THEN
      MN1=IMON+12
      YR1=YR1-1
      ENDIF
1125  IF (IYR.LT.1582) GOTO 1145
1130  IF (IYR.EQ.1582 .AND. IMON.LT.10) GOTO 1145
1135  IF (IYR.EQ.1582 .AND. IMON.EQ.10 .AND. DY.LT.15.D0) GOTO 1145
1140  A=INT(YR1/100)
      B=2-A+INT(A/4)
1145  IF (YR1.LT.0) GOTO 1155
1150  C=INT(365.25*YR1)-694025
      GOTO 1160
1155  C=FNITG((365.25*YR1)-0.75)-694025
1160  D=INT(30.6001*(MN1+1))
1165  DJD=B+C+D+DY-0.5
1170  RETURN
      END
      SUBROUTINE GC13(RA,DEC,R,DR,STEP1,STEP2,RAB,DECB,RB,
     CWATCHRA,WATCHDEC,WATCHR,COMPTELRA,COMPTELDEC,COMPTELR,
     CEGRETRA,EGRETDEC,EGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,
     CSIGMARA,SIGMADEC,SIGMAR,SAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,
     CKONUSEL1,KONUSEL2,PRA,PDEC,PR,ORAB,ODECB,ORB,AREA1,AREA2,FRACTION,
     CNAREA,RATEST,DECTEST,INANN,INOTHER)
C.....ADAPTED FROM SUBROUTINE GC IN BMAP AND GC13 IN HETEIPN.  
C.....THIS SUBROUTINE FINDS THE RA, DEC
C.....COORDINATES ALONG AN IPN ANNULUS, GCA, GCD, GCAA, GCDD.  THE ANNULUS IS
C.....CENTERED AT RA, DEC, AND HAS RADIUS R +/- DR. 
C.....IN THIS VERSION, THE PARAMETER T IS ADJUSTED SO THAT POINTS ALONG THE
C.....CIRCLE (IN THE AZIMUTHAL DIRECTION) ARE FOUND FOR A TOTAL OF 1/STEP1 EQUAL ANGULAR STEPS.
C.....POINTS ALONG THE RADIAL DIRECTION ARE INCREMENTED EVERY 1/STEP2 DEGREES.  
C.....EACH POINT IS TESTED TO SEE WHETHER IT IS INSIDE THE OTHER LOCALIZATIONS
C.....IF IT IS, THE AREA IS INCREMENTED.
C.....THIS SUBROUTINE IS FOR ONE AND ONLY ONE IPN ANNULUS, PLUS ANY OTHER
C.....LOCALIZATION
C.....THIS VERSION TESTS THE POINT RATEST,DECTEST TO SEE WHETHER IT IS WITHIN
C.....ALL THE LOCALIZATIONS OF THE BURST
      DOUBLE PRECISION PI,RD,X,RA1,RA2,RA3,DEC1,DEC2,DEC3,DIST,
     CGCA,GCD,GCAA,GCDD,RA,DEC,R,DR,STEP1,STEP2,RAB,DECB,RB,
     CWATCHRA,WATCHDEC,WATCHR,COMPTELRA,COMPTELDEC,COMPTELR,
     CEGRETRA,EGRETDEC,EGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,
     CSIGMARA,SIGMADEC,SIGMAR,SAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,
     CKONUSEL1,KONUSEL2,PRA,PDEC,PR,ORAB,ODECB,ORB,RAEL,DECEL,
     CAREA1,AREA2,FRACTION,ARE,CHI,D,DTHETA,DTHETA2,THETA,THETA2,T,
     CSCHI,CCHI,SD,CD,ALPH,DELT,T2,CHI1,SCHI1,CCHI1,ANG,RATEST,DECTEST,
     CA1,A2,FRAC
      COMMON/GEOM/PI,RD
      LOGICAL INANN,INOTHER,IN
C.....DEFINE IN TO BE TRUE, AND TEST EVERY LOCALIZATION; CHANGE IT TO 
C.....FALSE IF IT FAILS ANY TEST
      INANN=.TRUE.
      INOTHER=.TRUE.
      AREA1=0.
      NAREA=0
      IF(R.EQ.0.)THEN
      PRINT *,'RA,DEC,R=',RA,DEC,R
      STOP 'R=0 IN GC13'
      ENDIF
C.....CHECK WHETHER THE TEST POINT IS WITHIN THE ANNULUS
      CALL ANNULUS(RA,DEC,R,DR,AREA1,AREA2,FRACTION,RATEST,
     CDECTEST,INANN)
20    CONTINUE
C.....THETA IS THE AZIMUTHAL ANGLE (AROUND THE CIRCLE). DTHETA IS THE ANGULAR STEP
C.....IN THE AZIMUTHAL DIRECTION
      NTIMES=1
C.....THIS CALCULATES THE RA, DEC VALUES WHICH ARE WITHIN THE LIMITS OF THE ERROR BOX
C.....THIS FIRST LOOP TRACES THE RIGHT HALF OF A CIRCLE FROM SMALLER DEC
C.....VALUES TO LARGER.  T GOES FROM +1 TO -1
C.....FIRST STEP: OUTER RADIUS OF THE ANNULUS
C.....CHI IS THE HALF ANGLE OF THE ANNULUS (RADIANS),
C.....AND THE NORTH POLE IS AT RA=X DEGREES,DEC=D RADIANS
      CHI=(R-DR)*RD
      X=RA
      D=DEC*RD
      DTHETA=PI*STEP1
      DTHETA2=PI*STEP2
5     THETA=0.D0
      T=1.D0
      IMAX=INT(1./STEP1)+1
C.....THIS CALCULATES THE RA, DEC VALUES WHICH ARE WITHIN THE LIMITS OF THE ERROR BOX
C.....THIS FIRST LOOP TRACES THE RIGHT HALF OF A CIRCLE FROM SMALLER DEC
C.....VALUES TO LARGER.  T GOES FROM +1 TO -1
C.....FIRST STEP: OUTER RADIUS OF THE ANNULUS
10    SCHI=DSIN(CHI)
      CCHI=DCOS(CHI)
      SD=DSIN(D)
      CD=DCOS(D)
      DO 100 I=1,IMAX
      THETA=THETA+DTHETA
50    IF(T.EQ.-1.D0) GO TO 105
      T=DCOS(THETA)
      IF(T.LT.-1.D0) T=-1.D0
      GCA=ALPH(SCHI,CCHI,SD,CD,T) + X
      IF(GCA.GT.360.D0)GCA=GCA-360.D0
      GCD=DELT(SCHI,CCHI,SD,CD,T) 
C.....CHECK TO SEE IF THIS POINT IS WITHIN THE BATSE, EGRET, COMPTEL, ...
C.....ERROR CIRCLES, WITHIN THE KONUS ECLIPTIC LATITUDE BAND, OUTSIDE
C.....THE PLANET-BLOCKING CIRCLE, ETC.  IF THE POINT IS INSIDE, CONTINUE
C.....TESTING; IF OUTSIDE, GET ANOTHER POINT
      IF(RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,RAB,DECB,DIST)
      IF(DIST.GT.SQRT(RB**2+2.56))GOTO 100
      ENDIF
      IF(WATCHRA.NE.0. .AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,WATCHRA,WATCHDEC,DIST)
      IF(DIST.GT.WATCHR)GOTO 100
      ENDIF
      IF(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.)
     CTHEN
      CALL ANGLE1(GCA,GCD,COMPTELRA,COMPTELDEC,DIST)
      IF(DIST .GT. COMPTELR)GOTO 100
      ENDIF
      IF(EGRETRA.NE.0. .AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,EGRETRA,EGRETDEC,DIST)
      IF(DIST.GT.EGRETR)GOTO 100
      ENDIF
      IF(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,PHEBUSRA,PHEBUSDEC,DIST)
      IF(DIST.GT.PHEBUSR)GOTO 100
      ENDIF
      IF(SIGMARA.NE.0. .AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,SIGMARA,SIGMADEC,DIST)
      IF(DIST.GT.SIGMAR)GOTO 100
      ENDIF
      IF(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,SAXRA,SAXDEC,DIST)
      IF(DIST.GT.SAXR)GOTO 100
      ENDIF
      IF(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,HRAB,HDECB,DIST)
      IF(DIST.GT.HRB)GOTO 100
      ENDIF
      IF(KONUSEL1 .NE.0. .OR. KONUSEL2.NE.0.)THEN
      RAEL=270.D0
      DECEL=66.55D0
      REL=90.D0-(KONUSEL1+KONUSEL2)/2.
      DREL=ABS((90.-KONUSEL1)-REL)
      CALL ANGLE1(GCA,GCD,RAEL,DECEL,DIST)
      IF(DIST.GT. (REL+DREL) .OR. DIST .LT. (REL-DREL) )GOTO 100
      ENDIF
      IF(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,PRA,PDEC,DIST)
C.....THIS ANGLE HAS TO BE GREATER THAN PR IF THE POINT IS OUTSIDE THE PLANET
      IF(DIST .LT. PR)GOTO 100
      ENDIF
      IF(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,ORAB,ODECB,DIST)
      IF(DIST.GT.ORB)GOTO 100
      ENDIF
C.....THE POINT HAS PASSED ALL THE TESTS.  INCREMENT THE AREA.
C.....THE ELEMENTAL AREA IS DEFINED BY 3 RA, DEC POINTS.  THE
C.....FIRST TWO ARE ON THE CIRCLE.  THE THIRD IS OBTAINED BY
C.....INCREMENTING THE RADIAL DISTANCE
      NAREA=NAREA+1
      RA1=GCA
      DEC1=GCD
      THETA2=THETA+DTHETA		     
      T2=DCOS(THETA2)
      IF(T2.LT.-1.) T2=-1.
      RA2=ALPH(SCHI,CCHI,SD,CD,T2) + X
      IF(RA2.GT.360.D0)RA2=RA2-360.D0
      DEC2=DELT(SCHI,CCHI,SD,CD,T2) 
      CHI1=CHI+DTHETA2
      SCHI1=DSIN(CHI1)
      CCHI1=DCOS(CHI1)
      RA3=ALPH(SCHI1,CCHI1,SD,CD,T2) + X
      IF(RA3.GT.360.D0)RA3=RA3-360.D0
      DEC3=DELT(SCHI1,CCHI1,SD,CD,T2) 
C      CALL AREA(RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE)
C	   PRINT *,RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE
C.....APPROXIMATE THE ELEMENTAL AREA USING A PLANE TRIANGLE
      CALL ANGLE1(RA1,DEC1,RA2,DEC2,ANG)
      ARE=ANG*DTHETA2/RD
C	   PRINT *,RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE
      AREA1=AREA1+ARE*3600. 
100   CONTINUE
105   THETA=PI
      T=COS(THETA+DTHETA)
C.....THIS SECOND LOOP TRACES THE LEFT HALF OF A CIRCLE GOING FROM LARGER
C.....DEC VALUES TO SMALLER.  T GOES FROM -1 TO +1
      T=-1.D0
      DO 200 J=1,IMAX
      THETA=THETA+DTHETA
      IF(T.EQ.1.D0) GO TO 202
      T=DCOS(THETA)
      IF(T.GT.1.D0) T=1.D0
      GCAA=360.D0-ALPH(SCHI,CCHI,SD,CD,T)+X
      IF(GCAA.GT.360.D0)GCAA=GCAA-360.D0
      GCDD=DELT(SCHI,CCHI,SD,CD,T)
C.....CHECK TO SEE IF THIS POINT IS WITHIN THE BATSE, EGRET, COMPTEL, ...
C.....ERROR CIRCLES, WITHIN THE KONUS ECLIPTIC LATITUDE BAND, OUTSIDE
C.....THE PLANET-BLOCKING CIRCLE, ETC.  IF THE POINT IS INSIDE, CONTINUE
C.....TESTING; IF OUTSIDE, GET ANOTHER POINT
      IF(RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,RAB,DECB,DIST)
      IF(DIST.GT.SQRT(RB**2+2.56))GOTO 200
      ENDIF
      IF(WATCHRA.NE.0. .AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,WATCHRA,WATCHDEC,DIST)
      IF(DIST.GT.WATCHR)GOTO 200
      ENDIF
      IF(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.)
     CTHEN
      CALL ANGLE1(GCAA,GCDD,COMPTELRA,COMPTELDEC,DIST)
      IF(DIST .GT. COMPTELR)GOTO 200
      ENDIF
      IF(EGRETRA.NE.0. .AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,EGRETRA,EGRETDEC,DIST)
      IF(DIST.GT.EGRETR)GOTO 200
      ENDIF
      IF(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,PHEBUSRA,PHEBUSDEC,DIST)
      IF(DIST.GT.PHEBUSR)GOTO 200
      ENDIF
      IF(SIGMARA.NE.0. .AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,SIGMARA,SIGMADEC,DIST)
      IF(DIST.GT.SIGMAR)GOTO 200
      ENDIF
      IF(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,SAXRA,SAXDEC,DIST)
      IF(DIST.GT.SAXR)GOTO 200
      ENDIF
      IF(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,HRAB,HDECB,DIST)
      IF(DIST.GT.HRB)GOTO 200
      ENDIF
      IF(KONUSEL1 .NE.0. .OR. KONUSEL2.NE.0.)THEN
      RAEL=270.D0
      DECEL=66.55D0
      REL=90.-(KONUSEL1+KONUSEL2)/2.
      DREL=ABS((90.-KONUSEL1)-REL)
      CALL ANGLE1(GCAA,GCDD,RAEL,DECEL,DIST)
      IF(DIST.GT. (REL+DREL) .OR. DIST .LT. (REL-DREL) )GOTO 200
      ENDIF
      IF(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,PRA,PDEC,DIST)
C.....THIS ANGLE HAS TO BE GREATER THAN PR IF THE POINT IS OUTSIDE THE PLANET
      IF(DIST .LT. PR)GOTO 200
      ENDIF
      IF(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,ORAB,ODECB,DIST)
      IF(DIST.GT.ORB)GOTO 200
      ENDIF
C.....THE POINT HAS PASSED ALL THE TESTS.  INCREMENT THE AREA.
C.....THE ELEMENTAL AREA IS DEFINED BY 3 RA, DEC POINTS.
C.....RA1, DEC1 IS THE POINT JUST TREATED
C.....RA2, DEC2 IS THE POINT ALONG THE CIRCLE
C.....RA3, DEC3 IS THE POINT IN THE RADIAL DIRECTION
      NAREA=NAREA+1
      RA1=GCAA
      DEC1=GCDD
      THETA2=THETA+DTHETA
      T2=DCOS(THETA2)
      IF(T2.LT.-1.D0) T2=-1.D0
      RA2=360.D0-ALPH(SCHI,CCHI,SD,CD,T2) + X
      IF(RA2.GT.360.D0)RA2=RA2-360.D0
      DEC2=DELT(SCHI,CCHI,SD,CD,T2) 
      CHI1=CHI+DTHETA2
      SCHI1=DSIN(CHI1)
      CCHI1=DCOS(CHI1)
      RA3=360.D0-ALPH(SCHI1,CCHI1,SD,CD,T2) + X
      IF(RA3.GT.360.D0)RA3=RA3-360.D0
      DEC3=DELT(SCHI1,CCHI1,SD,CD,T2) 
C      RA3=RA2
C      DEC3=DEC2+DTHETA2
C      CALL AREA(RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE)
C.....APPROXIMATE THE ELEMENTAL AREA USING A PLANE TRIANGLE
      CALL ANGLE1(RA1,DEC1,RA2,DEC2,ANG)
      ARE=ANG*DTHETA2/RD
C	   PRINT *,RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE
      AREA1=AREA1+ARE*3600. 
200   CONTINUE
C.....NOW INCREMENT THE RADIAL ANGLE AND GO AROUND THE CIRCLE AGAIN
202   CHI=CHI+DTHETA2
      J=J+1
      IF(CHI.GE.(R+DR)*RD)GOTO 205
      X=RA
      D=DEC*RD
      GOTO 5
C.....AREA1 IS IN SQUARE ARCMINUTES.  CONVERT IT TO STERADIANS.
C.....AREA2 IS IN SQUARE DEGREES
205   AREA2=AREA1/3600.D0
      AREA1=AREA1*RD**2/3600.D0
      FRACTION=AREA1/(4*PI)
C.....AREA CALCULATION HAS BEEN DONE.  NOW CHECK TO SEE IF THE TEST
C.....POINT IS WITHIN ALL THE LOCALIZATIONS.  IF IT IS NOT WITHIN ANY
C.....ONE OF THEM, RETURN.
C.....IF THERE IS NO INTERSECTION BETWEEN THE IPN ANNULUS AND THE OTHER
C.....LOCALIZATION, TAKE THE AREA TO BE THE ANNULUS AREA
C.....CHECK WHETHER THE TEST POINT IS WITHIN THE ANNULUS
      CALL ANNULUS(RA,DEC,R,DR,A1,A2,FRAC,RATEST,
     CDECTEST,INANN)
      IF(NAREA.EQ.0)THEN
      AREA1=A1
      AREA2=A2
      FRACTION=FRAC
      ENDIF
      IF(RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,RAB,DECB,DIST)
      IF(DIST.GT.SQRT(RB**2+2.56))INOTHER=.FALSE.
      IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(WATCHRA.NE.0. .AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,WATCHRA,WATCHDEC,DIST)
      IF(DIST.GT.WATCHR)INOTHER=.FALSE.
      IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.)
     CTHEN
      CALL ANGLE1(RATEST,DECTEST,COMPTELRA,COMPTELDEC,DIST)
      IF(DIST .GT. COMPTELR)INOTHER=.FALSE.
      IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(EGRETRA.NE.0. .AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,EGRETRA,EGRETDEC,DIST)
      IF(DIST.GT.EGRETR)INOTHER=.FALSE.
      IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,PHEBUSRA,PHEBUSDEC,DIST)
      IF(DIST.GT.PHEBUSR)INOTHER=.FALSE.
      IF(.NOT.IN)RETURN
      ENDIF
      IF(SIGMARA.NE.0. .AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,SIGMARA,SIGMADEC,DIST)
      IF(DIST.GT.SIGMAR)INOTHER=.FALSE.
      IF(.NOT.IN)RETURN
      ENDIF
      IF(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,SAXRA,SAXDEC,DIST)
      IF(DIST.GT.SAXR)INOTHER=.FALSE.
      IF(.NOT.IN)RETURN
      ENDIF
      IF(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,HRAB,HDECB,DIST)
      IF(DIST.GT.HRB)INOTHER=.FALSE.
      IF(.NOT.IN)RETURN
      ENDIF
      IF(KONUSEL1 .NE.0. .OR. KONUSEL2.NE.0.)THEN
      RAEL=270.D0
      DECEL=66.55D0
      REL=90.-(KONUSEL1+KONUSEL2)/2.
      DREL=ABS((90.-KONUSEL1)-REL)
      CALL ANGLE1(RATEST,DECTEST,RAEL,DECEL,DIST)
      IF(DIST.GT. (REL+DREL) .OR. DIST .LT. (REL-DREL) )
     CINOTHER=.FALSE.
      IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,PRA,PDEC,DIST)
C.....THIS ANGLE HAS TO BE GREATER THAN PR IF THE POINT IS OUTSIDE THE PLANET
      IF(DIST .LT. PR)INOTHER=.FALSE.
      IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,ORAB,ODECB,DIST)
      IF(DIST.GT.ORB)INOTHER=.FALSE.
      IF(.NOT.INOTHER)RETURN
      ENDIF
      RETURN
      END
      SUBROUTINE GC14(RAIPN1,DECIPN1,RIPN1,DRIPN1,RAIPN2,DECIPN2,
     CRIPN2,DRIPN2,STEP1,STEP2,AREA1,AREA2,FRACTION,NAREA,
     CRATEST,DECTEST,IN)
C.....ADAPTED FROM SUBROUTINE GC IN BMAP AND GC13 IN HETEIPN.  
C.....THIS SUBROUTINE FINDS THE RA, DEC
C.....COORDINATES ALONG AN IPN ANNULUS, GCA, GCD, GCAA, GCDD.  THE ANNULUS IS
C.....CENTERED AT RA, DEC, AND HAS RADIUS R +/- DR. 
C.....IN THIS VERSION, THE PARAMETER T IS ADJUSTED SO THAT POINTS ALONG THE
C.....CIRCLE (IN THE AZIMUTHAL DIRECTION) ARE FOUND FOR A TOTAL OF 1/STEP1 EQUAL ANGULAR STEPS.
C.....POINTS ALONG THE RADIAL DIRECTION ARE INCREMENTED EVERY 1/STEP2 DEGREES.  
C.....EACH POINT IS TESTED TO SEE WHETHER IT IS INSIDE THE OTHER LOCALIZATIONS
C.....IF IT IS, THE AREA IS INCREMENTED.
C.....THIS SUBROUTINE IS FOR TWO IPN ANNULI AND NO OTHER LOCALIZATION INFORMATION
C.....IN THIS VERSION, A CHECK IS DONE TO SEE IF THE TEST POINT RATEST, DECTEST,
C.....IS WITHIN BOTH ANNULI
      DOUBLE PRECISION PI,RD,X,RA1,RA2,RA3,DEC1,DEC2,DEC3,DIST,
     CGCA,GCD,GCAA,GCDD,RAIPN1,DECIPN1,RIPN1,DRIPN1,RAIPN2,DECIPN2,
     CRIPN2,DRIPN2,STEP1,STEP2,RA4,DEC4,R4,DR4,
     CAREA1,AREA2,FRACTION,ARE,CHI,D,DTHETA,DTHETA2,THETA,THETA2,T,
     CSCHI,CCHI,SD,CD,ALPH,DELT,T2,CHI1,SCHI1,CCHI1,ANG,RATEST,DECTEST,
     CA1,A2,FRAC
      COMMON/GEOM/PI,RD
	LOGICAL IN
      AREA1=0.
	NAREA=0.
C.....DEFINE IN TO BE TRUE
      IN=.TRUE.
      IF(RIPN1.EQ.0. .OR. RIPN2.EQ.0.)THEN
      PRINT *,'IPN ANNULUS HAS ZERO WIDTH IN GC14'
      STOP
      ENDIF
20    CONTINUE
C.....THETA IS THE AZIMUTHAL ANGLE (AROUND THE CIRCLE). DTHETA IS THE ANGULAR STEP
C.....IN THE AZIMUTHAL DIRECTION
      NTIMES=1
C.....THIS CALCULATES THE RA, DEC VALUES WHICH ARE WITHIN THE LIMITS OF THE ERROR BOX
C.....THIS FIRST LOOP TRACES THE RIGHT HALF OF A CIRCLE FROM SMALLER DEC
C.....VALUES TO LARGER.  T GOES FROM +1 TO -1
C.....FIRST STEP: OUTER RADIUS OF THE ANNULUS
C.....CHI IS THE HALF ANGLE OF THE ANNULUS (RADIANS),
C.....AND THE NORTH POLE IS AT RA=X DEGREES,DEC=D RADIANS
C.....USE THE NARROWER OF THE TWO IPN ANNULI
	RA=RAIPN1
	DEC=DECIPN1
	R=RIPN1
	DR=DRIPN1
	RA4=RAIPN2
	DEC4=DECIPN2
	R4=RIPN2
	DR4=DRIPN2
	IF(DRIPN2.LT.DRIPN1)THEN
	RA=RAIPN2
	DEC=DECIPN2
	R=RIPN2
	DR=DRIPN2
	RA4=RAIPN1
	DEC4=DECIPN1
	R4=RIPN1
	DR4=DRIPN1
	ENDIF
C	      PRINT *,RA,DEC,R,DR,RA4,DEC4,R4,DR4
      CHI=(R-DR)*RD
      X=RA
      D=DEC*RD
      DTHETA=PI*STEP1
      DTHETA2=PI*STEP2
5     THETA=0.D0
      T=1.D0
      IMAX=INT(1./STEP1)+1
C.....THIS CALCULATES THE RA, DEC VALUES WHICH ARE WITHIN THE LIMITS OF THE ERROR BOX
C.....THIS FIRST LOOP TRACES THE RIGHT HALF OF A CIRCLE FROM SMALLER DEC
C.....VALUES TO LARGER.  T GOES FROM +1 TO -1
C.....FIRST STEP: OUTER RADIUS OF THE ANNULUS
10    SCHI=DSIN(CHI)
      CCHI=DCOS(CHI)
      SD=DSIN(D)
      CD=DCOS(D)
      DO 100 I=1,IMAX
      THETA=THETA+DTHETA
50    IF(T.EQ.-1.D0) GO TO 105
      T=DCOS(THETA)
      IF(T.LT.-1.D0) T=-1.D0
      GCA=ALPH(SCHI,CCHI,SD,CD,T) + X
      IF(GCA.GT.360.D0)GCA=GCA-360.D0
      GCD=DELT(SCHI,CCHI,SD,CD,T) 
C.....CHECK TO SEE IF THIS POINT IS WITHIN THE OTHER IPN ANNULUS.
C.....IF THE POINT IS INSIDE, CONTINUE
C.....TESTING; IF OUTSIDE, GET ANOTHER POINT
      CALL ANGLE1(GCA,GCD,RA4,DEC4,DIST)
      IF(DIST.GT.R4+DR4 .OR. DIST.LT.R4-DR4)GOTO 100 
C.....THE POINT HAS PASSED ALL THE TESTS.  INCREMENT THE AREA.
C.....THE ELEMENTAL AREA IS DEFINED BY 3 RA, DEC POINTS.  THE
C.....FIRST TWO ARE ON THE CIRCLE.  THE THIRD IS OBTAINED BY
C.....INCREMENTING THE RADIAL DISTANCE
      NAREA=NAREA+1
      RA1=GCA
      DEC1=GCD
      THETA2=THETA+DTHETA		     
      T2=DCOS(THETA2)
      IF(T2.LT.-1.) T2=-1.
      RA2=ALPH(SCHI,CCHI,SD,CD,T2) + X
      IF(RA2.GT.360.D0)RA2=RA2-360.D0
      DEC2=DELT(SCHI,CCHI,SD,CD,T2) 
      CHI1=CHI+DTHETA2
      SCHI1=DSIN(CHI1)
      CCHI1=DCOS(CHI1)
      RA3=ALPH(SCHI1,CCHI1,SD,CD,T2) + X
      IF(RA3.GT.360.D0)RA3=RA3-360.D0
      DEC3=DELT(SCHI1,CCHI1,SD,CD,T2) 
C      CALL AREA(RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE)
C	   PRINT *,RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE
C.....APPROXIMATE THE ELEMENTAL AREA USING A PLANE TRIANGLE
      CALL ANGLE1(RA1,DEC1,RA2,DEC2,ANG)
	ARE=ANG*DTHETA2/RD
      AREA1=AREA1+ARE*3600. 
100   CONTINUE
105   THETA=PI
      T=COS(THETA+DTHETA)
C.....THIS SECOND LOOP TRACES THE LEFT HALF OF A CIRCLE GOING FROM LARGER
C.....DEC VALUES TO SMALLER.  T GOES FROM -1 TO +1
      T=-1.D0
      DO 200 J=1,IMAX
      THETA=THETA+DTHETA
      IF(T.EQ.1.D0) GO TO 202
      T=DCOS(THETA)
      IF(T.GT.1.D0) T=1.D0
      GCAA=360.D0-ALPH(SCHI,CCHI,SD,CD,T)+X
      IF(GCAA.GT.360.D0)GCAA=GCAA-360.D0
      GCDD=DELT(SCHI,CCHI,SD,CD,T)
C.....CHECK TO SEE IF THIS POINT IS WITHIN THE OTHER IPN ANNULUS
C.....IF THE POINT IS INSIDE, CONTINUE
C.....TESTING; IF OUTSIDE, GET ANOTHER POINT
      CALL ANGLE1(GCAA,GCDD,RA4,DEC4,DIST)
      IF(DIST.GT.R4+DR4 .OR. DIST.LT.R4-DR4)GOTO 200 
C.....THE POINT HAS PASSED ALL THE TESTS.  INCREMENT THE AREA.
C.....THE ELEMENTAL AREA IS DEFINED BY 3 RA, DEC POINTS.
C.....RA1, DEC1 IS THE POINT JUST TREATED
C.....RA2, DEC2 IS THE POINT ALONG THE CIRCLE
C.....RA3, DEC3 IS THE POINT IN THE RADIAL DIRECTION
      NAREA=NAREA+1
      RA1=GCAA
      DEC1=GCDD
      THETA2=THETA+DTHETA
      T2=DCOS(THETA2)
      IF(T2.LT.-1.D0) T2=-1.D0
      RA2=360.D0-ALPH(SCHI,CCHI,SD,CD,T2) + X
      IF(RA2.GT.360.D0)RA2=RA2-360.D0
      DEC2=DELT(SCHI,CCHI,SD,CD,T2) 
      CHI1=CHI+DTHETA2
      SCHI1=DSIN(CHI1)
      CCHI1=DCOS(CHI1)
      RA3=360.D0-ALPH(SCHI1,CCHI1,SD,CD,T2) + X
      IF(RA3.GT.360.D0)RA3=RA3-360.D0
      DEC3=DELT(SCHI1,CCHI1,SD,CD,T2) 
      CALL ANGLE1(RA1,DEC1,RA2,DEC2,ANG)
	ARE=ANG*DTHETA2/RD
      AREA1=AREA1+ARE*3600. 
200   CONTINUE
C.....NOW INCREMENT THE RADIAL ANGLE AND GO AROUND THE CIRCLE AGAIN
202   CHI=CHI+DTHETA2
      J=J+1
      IF(CHI.GE.(R+DR)*RD)GOTO 205
      X=RA
      D=DEC*RD
      GOTO 5
C.....AREA1 IS IN SQUARE ARCMINUTES.  CONVERT IT TO STERADIANS.
C.....AREA2 IS IN SQUARE DEGREES
205   AREA2=AREA1/3600.D0
      AREA1=AREA1*RD**2/3600.D0
      FRACTION=AREA1/(4*PI)
C.....CHECK WHETHER RATEST, DECTEST IS WITHIN BOTH ANNULI; CHANGE THE
C.....VARIABLES AREA1,AREA2,AND FRACTION SO THEY DON'T OVERWRITE THE
C.....ONES JUST FOUND
      CALL ANNULUS(RAIPN1,DECIPN1,RIPN1,DRIPN1,A1,A2,FRAC,RATEST,
     CDECTEST,IN)
	IF(.NOT.IN)RETURN
      CALL ANNULUS(RAIPN2,DECIPN2,RIPN2,DRIPN2,A1,A2,FRAC,RATEST,
     CDECTEST,IN)
      RETURN
      END
      SUBROUTINE GC15(RAIPN1,DECIPN1,RIPN1,DRIPN1,RAIPN2,DECIPN2,
     CRIPN2,DRIPN2,RAB,DECB,
     CRB,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,COMPTELDEC,COMPTELR,
     CEGRETRA,EGRETDEC,EGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,
     CSIGMARA,SIGMADEC,SIGMAR,SAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,
     CKONUSEL1,KONUSEL2,PRA,PDEC,PR,ORAB,ODECB,ORB,
     CSTEP1,STEP2,AREA1,AREA2,FRACTION,NAREA,RATEST,DECTEST,INANN1,
     CINANN2,INOTHER)
C.....ADAPTED FROM SUBROUTINE GC IN BMAP AND GC13 IN HETEIPN.  
C.....THIS SUBROUTINE FINDS THE RA, DEC
C.....COORDINATES ALONG AN IPN ANNULUS, GCA, GCD, GCAA, GCDD.  THE ANNULUS IS
C.....CENTERED AT RA, DEC, AND HAS RADIUS R +/- DR. 
C.....IN THIS VERSION, THE PARAMETER T IS ADJUSTED SO THAT POINTS ALONG THE
C.....CIRCLE (IN THE AZIMUTHAL DIRECTION) ARE FOUND FOR A TOTAL OF 1/STEP1 EQUAL ANGULAR STEPS.
C.....POINTS ALONG THE RADIAL DIRECTION ARE INCREMENTED EVERY 1/STEP2 DEGREES.  
C.....EACH POINT IS TESTED TO SEE WHETHER IT IS INSIDE THE OTHER LOCALIZATIONS
C.....IF IT IS, THE AREA IS INCREMENTED.
C.....THIS SUBROUTINE IS FOR TWO IPN ANNULI AND ANY OTHER LOCALIZATION INFORMATION
C.....IT ALSO CHECKS THE POSITION RATEST, DECTEST TO SEE WHETHER IT IS IN THE ERROR
C.....BOX
      DOUBLE PRECISION PI,RD,X,RA1,RA2,RA3,DEC1,DEC2,DEC3,DIST,
     CGCA,GCD,GCAA,GCDD,RAIPN1,DECIPN1,RIPN1,DRIPN1,RAIPN2,DECIPN2,
     CRIPN2,DRIPN2,STEP1,STEP2,RA4,DEC4,R4,DR4,
     CRAB,DECB,RB,
     CWATCHRA,WATCHDEC,WATCHR,COMPTELRA,COMPTELDEC,COMPTELR,
     CEGRETRA,EGRETDEC,EGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,
     CSIGMARA,SIGMADEC,SIGMAR,SAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,
     CKONUSEL1,KONUSEL2,PRA,PDEC,PR,ORAB,ODECB,ORB,RAEL,DECEL,
     CAREA1,AREA2,FRACTION,ARE,CHI,D,DTHETA,DTHETA2,THETA,THETA2,T,
     CSCHI,CCHI,SD,CD,ALPH,DELT,T2,CHI1,SCHI1,CCHI1,ANG,RATEST,DECTEST,
     CA1,A2,FRAC
      COMMON/GEOM/PI,RD
      LOGICAL IN,INANN1,INANN2,INOTHER
      AREA1=0.
	NAREA=0
C.....DEFINE IN TO BE TRUE
      IN=.TRUE.
      IF(RIPN1.EQ.0. .OR. RIPN2.EQ.0.)THEN
      PRINT *,'IPN ANNULUS HAS ZERO WIDTH IN GC15'
      STOP
      ENDIF
20    CONTINUE
C.....THETA IS THE AZIMUTHAL ANGLE (AROUND THE CIRCLE). DTHETA IS THE ANGULAR STEP
C.....IN THE AZIMUTHAL DIRECTION
      NTIMES=1
C.....THIS CALCULATES THE RA, DEC VALUES WHICH ARE WITHIN THE LIMITS OF THE ERROR BOX
C.....THIS FIRST LOOP TRACES THE RIGHT HALF OF A CIRCLE FROM SMALLER DEC
C.....VALUES TO LARGER.  T GOES FROM +1 TO -1
C.....FIRST STEP: OUTER RADIUS OF THE ANNULUS
C.....CHI IS THE HALF ANGLE OF THE ANNULUS (RADIANS),
C.....AND THE NORTH POLE IS AT RA=X DEGREES,DEC=D RADIANS
C.....USE THE NARROWER OF THE TWO IPN ANNULI
      RA=RAIPN1
      DEC=DECIPN1
      R=RIPN1
      DR=DRIPN1
      RA4=RAIPN2
      DEC4=DECIPN2
      R4=RIPN2
      DR4=DRIPN2
      IF(DRIPN2.LT.DRIPN1)THEN
      RA=RAIPN2
      DEC=DECIPN2
      R=RIPN2
      DR=DRIPN2
      RA4=RAIPN1
      DEC4=DECIPN1
      R4=RIPN1
      DR4=DRIPN1
      ENDIF
C	      PRINT *,RA,DEC,R,DR,RA4,DEC4,R4,DR4
      CHI=(R-DR)*RD
      X=RA
      D=DEC*RD
      DTHETA=PI*STEP1
      DTHETA2=PI*STEP2
5     THETA=0.D0
      T=1.D0
      IMAX=INT(1./STEP1)+1
C.....THIS CALCULATES THE RA, DEC VALUES WHICH ARE WITHIN THE LIMITS OF THE ERROR BOX
C.....THIS FIRST LOOP TRACES THE RIGHT HALF OF A CIRCLE FROM SMALLER DEC
C.....VALUES TO LARGER.  T GOES FROM +1 TO -1
C.....FIRST STEP: OUTER RADIUS OF THE ANNULUS
10    SCHI=DSIN(CHI)
      CCHI=DCOS(CHI)
      SD=DSIN(D)
      CD=DCOS(D)
      DO 100 I=1,IMAX
      THETA=THETA+DTHETA
50    IF(T.EQ.-1.D0) GO TO 105
      T=DCOS(THETA)
      IF(T.LT.-1.D0) T=-1.D0
      GCA=ALPH(SCHI,CCHI,SD,CD,T) + X
      IF(GCA.GT.360.D0)GCA=GCA-360.D0
      GCD=DELT(SCHI,CCHI,SD,CD,T) 
C.....CHECK TO SEE IF THIS POINT IS WITHIN THE OTHER IPN ANNULUS.
C.....IF THE POINT IS INSIDE, CONTINUE
C.....TESTING; IF OUTSIDE, GET ANOTHER POINT
      CALL ANGLE1(GCA,GCD,RA4,DEC4,DIST)
      IF(DIST.GT.R4+DR4 .OR. DIST.LT.R4-DR4)GOTO 100 
C.....CHECK TO SEE IF THIS POINT IS WITHIN THE BATSE, EGRET, COMPTEL, ...
C.....ERROR CIRCLES, WITHIN THE KONUS ECLIPTIC LATITUDE BAND, OUTSIDE
C.....THE PLANET-BLOCKING CIRCLE, ETC.  IF THE POINT IS INSIDE, CONTINUE
C.....TESTING; IF OUTSIDE, GET ANOTHER POINT
      IF(RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,RAB,DECB,DIST)
      IF(DIST.GT.SQRT(RB**2+2.56))GOTO 100
      ENDIF
      IF(WATCHRA.NE.0. .AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,WATCHRA,WATCHDEC,DIST)
      IF(DIST.GT.WATCHR)GOTO 100
      ENDIF
      IF(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.)
     CTHEN
      CALL ANGLE1(GCA,GCD,COMPTELRA,COMPTELDEC,DIST)
      IF(DIST .GT. COMPTELR)GOTO 100
      ENDIF
      IF(EGRETRA.NE.0. .AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,EGRETRA,EGRETDEC,DIST)
      IF(DIST.GT.EGRETR)GOTO 100
      ENDIF
      IF(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,PHEBUSRA,PHEBUSDEC,DIST)
      IF(DIST.GT.PHEBUSR)GOTO 100
      ENDIF
      IF(SIGMARA.NE.0. .AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,SIGMARA,SIGMADEC,DIST)
      IF(DIST.GT.SIGMAR)GOTO 100
      ENDIF
      IF(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,SAXRA,SAXDEC,DIST)
      IF(DIST.GT.SAXR)GOTO 100
      ENDIF
      IF(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,HRAB,HDECB,DIST)
      IF(DIST.GT.HRB)GOTO 100
      ENDIF
      IF(KONUSEL1 .NE.0. .OR. KONUSEL2.NE.0.)THEN
      RAEL=270.D0
      DECEL=66.55D0
      REL=90.D0-(KONUSEL1+KONUSEL2)/2.
      DREL=ABS((90.-KONUSEL1)-REL)
      CALL ANGLE1(GCA,GCD,RAEL,DECEL,DIST)
      IF(DIST.GT. (REL+DREL) .OR. DIST .LT. (REL-DREL) )GOTO 100
      ENDIF
      IF(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,PRA,PDEC,DIST)
C.....THIS ANGLE HAS TO BE GREATER THAN PR IF THE POINT IS OUTSIDE THE PLANET
      IF(DIST .LT. PR)GOTO 100
      ENDIF
      IF(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,ORAB,ODECB,DIST)
      IF(DIST.GT.ORB)GOTO 100
      ENDIF
C.....THE POINT HAS PASSED ALL THE TESTS.  INCREMENT THE AREA.
C.....THE ELEMENTAL AREA IS DEFINED BY 3 RA, DEC POINTS.  THE
C.....FIRST TWO ARE ON THE CIRCLE.  THE THIRD IS OBTAINED BY
C.....INCREMENTING THE RADIAL DISTANCE
      NAREA=NAREA+1
      RA1=GCA
      DEC1=GCD
      THETA2=THETA+DTHETA		     
      T2=DCOS(THETA2)
      IF(T2.LT.-1.) T2=-1.
      RA2=ALPH(SCHI,CCHI,SD,CD,T2) + X
      IF(RA2.GT.360.D0)RA2=RA2-360.D0
      DEC2=DELT(SCHI,CCHI,SD,CD,T2) 
      CHI1=CHI+DTHETA2
      SCHI1=DSIN(CHI1)
      CCHI1=DCOS(CHI1)
      RA3=ALPH(SCHI1,CCHI1,SD,CD,T2) + X
      IF(RA3.GT.360.D0)RA3=RA3-360.D0
      DEC3=DELT(SCHI1,CCHI1,SD,CD,T2) 
C      CALL AREA(RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE)
C	   PRINT *,RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE
C.....APPROXIMATE THE ELEMENTAL AREA USING A PLANE TRIANGLE
      CALL ANGLE1(RA1,DEC1,RA2,DEC2,ANG)
	ARE=ANG*DTHETA2/RD
      AREA1=AREA1+ARE*3600. 
100   CONTINUE
105   THETA=PI
      T=COS(THETA+DTHETA)
C.....THIS SECOND LOOP TRACES THE LEFT HALF OF A CIRCLE GOING FROM LARGER
C.....DEC VALUES TO SMALLER.  T GOES FROM -1 TO +1
      T=-1.D0
      DO 200 J=1,IMAX
      THETA=THETA+DTHETA
      IF(T.EQ.1.D0) GO TO 202
      T=DCOS(THETA)
      IF(T.GT.1.D0) T=1.D0
      GCAA=360.D0-ALPH(SCHI,CCHI,SD,CD,T)+X
      IF(GCAA.GT.360.D0)GCAA=GCAA-360.D0
      GCDD=DELT(SCHI,CCHI,SD,CD,T)
C.....CHECK TO SEE IF THIS POINT IS WITHIN THE OTHER IPN ANNULUS
C.....IF THE POINT IS INSIDE, CONTINUE
C.....TESTING; IF OUTSIDE, GET ANOTHER POINT
      CALL ANGLE1(GCAA,GCDD,RA4,DEC4,DIST)
      IF(DIST.GT.R4+DR4 .OR. DIST.LT.R4-DR4)GOTO 200 
C.....CHECK TO SEE IF THIS POINT IS WITHIN THE BATSE, EGRET, COMPTEL, ...
C.....ERROR CIRCLES, WITHIN THE KONUS ECLIPTIC LATITUDE BAND, OUTSIDE
C.....THE PLANET-BLOCKING CIRCLE, ETC.  IF THE POINT IS INSIDE, CONTINUE
C.....TESTING; IF OUTSIDE, GET ANOTHER POINT
      IF(RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,RAB,DECB,DIST)
      IF(DIST.GT.SQRT(RB**2+2.56))GOTO 200
      ENDIF
      IF(WATCHRA.NE.0. .AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,WATCHRA,WATCHDEC,DIST)
      IF(DIST.GT.WATCHR)GOTO 200
      ENDIF
      IF(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.)
     CTHEN
      CALL ANGLE1(GCAA,GCDD,COMPTELRA,COMPTELDEC,DIST)
      IF(DIST .GT. COMPTELR)GOTO 200
      ENDIF
      IF(EGRETRA.NE.0. .AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,EGRETRA,EGRETDEC,DIST)
      IF(DIST.GT.EGRETR)GOTO 200
      ENDIF
      IF(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,PHEBUSRA,PHEBUSDEC,DIST)
      IF(DIST.GT.PHEBUSR)GOTO 200
      ENDIF
      IF(SIGMARA.NE.0. .AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,SIGMARA,SIGMADEC,DIST)
      IF(DIST.GT.SIGMAR)GOTO 200
      ENDIF
      IF(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,SAXRA,SAXDEC,DIST)
      IF(DIST.GT.SAXR)GOTO 200
      ENDIF
      IF(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,HRAB,HDECB,DIST)
      IF(DIST.GT.HRB)GOTO 200
      ENDIF
      IF(KONUSEL1 .NE.0. .OR. KONUSEL2.NE.0.)THEN
      RAEL=270.D0
      DECEL=66.55D0
      REL=90.-(KONUSEL1+KONUSEL2)/2.
      DREL=ABS((90.-KONUSEL1)-REL)
      CALL ANGLE1(GCAA,GCDD,RAEL,DECEL,DIST)
      IF(DIST.GT. (REL+DREL) .OR. DIST .LT. (REL-DREL) )GOTO 200
      ENDIF
      IF(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,PRA,PDEC,DIST)
C.....THIS ANGLE HAS TO BE GREATER THAN PR IF THE POINT IS OUTSIDE THE PLANET
      IF(DIST .LT. PR)GOTO 200
      ENDIF
      IF(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,ORAB,ODECB,DIST)
      IF(DIST.GT.ORB)GOTO 200
      ENDIF
C.....THE POINT HAS PASSED ALL THE TESTS.  INCREMENT THE AREA.
C.....THE ELEMENTAL AREA IS DEFINED BY 3 RA, DEC POINTS.
C.....RA1, DEC1 IS THE POINT JUST TREATED
C.....RA2, DEC2 IS THE POINT ALONG THE CIRCLE
C.....RA3, DEC3 IS THE POINT IN THE RADIAL DIRECTION
      NAREA=NAREA+1
      RA1=GCAA
      DEC1=GCDD
      THETA2=THETA+DTHETA
      T2=DCOS(THETA2)
      IF(T2.LT.-1.D0) T2=-1.D0
      RA2=360.D0-ALPH(SCHI,CCHI,SD,CD,T2) + X
      IF(RA2.GT.360.D0)RA2=RA2-360.D0
      DEC2=DELT(SCHI,CCHI,SD,CD,T2) 
      CHI1=CHI+DTHETA2
      SCHI1=DSIN(CHI1)
      CCHI1=DCOS(CHI1)
      RA3=360.D0-ALPH(SCHI1,CCHI1,SD,CD,T2) + X
      IF(RA3.GT.360.D0)RA3=RA3-360.D0
      DEC3=DELT(SCHI1,CCHI1,SD,CD,T2) 
      CALL ANGLE1(RA1,DEC1,RA2,DEC2,ANG)
      ARE=ANG*DTHETA2/RD
      AREA1=AREA1+ARE*3600. 
200   CONTINUE
C.....NOW INCREMENT THE RADIAL ANGLE AND GO AROUND THE CIRCLE AGAIN
202   CHI=CHI+DTHETA2
      J=J+1
      IF(CHI.GE.(R+DR)*RD)GOTO 205
      X=RA
      D=DEC*RD
      GOTO 5
C.....AREA1 IS IN SQUARE ARCMINUTES.  CONVERT IT TO STERADIANS.
C.....AREA2 IS IN SQUARE DEGREES
205   AREA2=AREA1/3600.D0
      AREA1=AREA1*RD**2/3600.D0
      FRACTION=AREA1/(4*PI)
C.....NOW CHECK TO SEE IF RATEST, DECTEST IS WITHIN ALL THE LOCALIZATIONS
C.....CHECK WHETHER RATEST, DECTEST IS WITHIN BOTH ANNULI; CHANGE THE
C.....VARIABLES AREA1,AREA2,AND FRACTION SO THEY DON'T OVERWRITE THE
C.....ONES JUST FOUND
      CALL ANNULUS(RAIPN1,DECIPN1,RIPN1,DRIPN1,A1,A2,FRAC,RATEST,
     CDECTEST,INANN1)
      IF(.NOT.INANN1)RETURN
      CALL ANNULUS(RAIPN2,DECIPN2,RIPN2,DRIPN2,A1,A2,FRAC,RATEST,
     CDECTEST,INANN2)
      IF(.NOT.INANN2)RETURN
      IF(RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,RAB,DECB,DIST)
      IF(DIST.GT.SQRT(RB**2+2.56))INOTHER=.FALSE.
      IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(WATCHRA.NE.0. .AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,WATCHRA,WATCHDEC,DIST)
      IF(DIST.GT.WATCHR)INOTHER=.FALSE.
	IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.)
     CTHEN
      CALL ANGLE1(RATEST,DECTEST,COMPTELRA,COMPTELDEC,DIST)
      IF(DIST .GT. COMPTELR)INOTHER=.FALSE.
	IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(EGRETRA.NE.0. .AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,EGRETRA,EGRETDEC,DIST)
      IF(DIST.GT.EGRETR)INOTHER=.FALSE.
	IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,PHEBUSRA,PHEBUSDEC,DIST)
      IF(DIST.GT.PHEBUSR)INOTHER=.FALSE.
	IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(SIGMARA.NE.0. .AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,SIGMARA,SIGMADEC,DIST)
      IF(DIST.GT.SIGMAR)INOTHER=.FALSE.
	IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,SAXRA,SAXDEC,DIST)
      IF(DIST.GT.SAXR)INOTHER=.FALSE.
	IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,HRAB,HDECB,DIST)
      IF(DIST.GT.HRB)INOTHER=.FALSE.
	IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(KONUSEL1 .NE.0. .OR. KONUSEL2.NE.0.)THEN
      RAEL=270.D0
      DECEL=66.55D0
      REL=90.-(KONUSEL1+KONUSEL2)/2.
      DREL=ABS((90.-KONUSEL1)-REL)
      CALL ANGLE1(RATEST,DECTEST,RAEL,DECEL,DIST)
      IF(DIST.GT. (REL+DREL) .OR. DIST .LT. (REL-DREL) )INOTHER=.FALSE.
	IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,PRA,PDEC,DIST)
C.....THIS ANGLE HAS TO BE GREATER THAN PR IF THE POINT IS OUTSIDE THE PLANET
      IF(DIST .LT. PR)INOTHER=.FALSE.
	IF(.NOT.INOTHER)RETURN
      ENDIF
      IF(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,ORAB,ODECB,DIST)
      IF(DIST.GT.ORB)INOTHER=.FALSE.
	IF(.NOT.INOTHER)RETURN
      ENDIF
      RETURN
      END
      SUBROUTINE GC16(RACIRC,DECCIRC,RCIRC,RAB,DECB,
     CRB,WATCHRA,WATCHDEC,WATCHR,COMPTELRA,COMPTELDEC,COMPTELR,
     CEGRETRA,EGRETDEC,EGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,
     CSIGMARA,SIGMADEC,SIGMAR,SAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,
     CKONUSEL1,KONUSEL2,PRA,PDEC,PR,ORAB,ODECB,ORB,
     CSTEP1,STEP2,AREA1,AREA2,FRACTION,NAREA,KON,RATEST,DECTEST,IN)
C.....ADAPTED FROM SUBROUTINE GC13  
C.....THIS SUBROUTINE FINDS THE RA, DEC
C.....COORDINATES ALONG INSIDE A CIRCLE, GCA, GCD, GCAA, GCDD.  THE CIRCLE IS
C.....CENTERED AT RA, DEC, AND HAS RADIUS R +/- DR. 
C.....IN THIS VERSION, THE PARAMETER T IS ADJUSTED SO THAT POINTS ALONG THE
C.....CIRCLE (IN THE AZIMUTHAL DIRECTION) ARE FOUND FOR A TOTAL OF 1/STEP1 EQUAL ANGULAR STEPS.
C.....POINTS ALONG THE RADIAL DIRECTION ARE INCREMENTED EVERY 1/STEP2 DEGREES.  
C.....EACH POINT IS TESTED TO SEE WHETHER IT IS INSIDE THE OTHER LOCALIZATIONS
C.....IF IT IS, THE AREA IS INCREMENTED.
C.....THIS SUBROUTINE IS FOR THE CASE WHERE THERE ARE NO IPN ANNULI.  THE LOCALIZATION CONSISTS
C.....OF AT LEAST TWO OF ANYTHING ELSE, LIKE:
C.....KONUS EL + EARTH BLOCKING
C.....OTHER LOCALIZATION + SAX CIRCLE
C.....EARTH BLOCKING + OTHER LOCALIZATION
C.....ETC.
C.....IT ALSO TESTS THE POSITION RATEST,DECTEST TO SEE WHETHER IT IS IN
C.....THE LOCALIZATIONS
      DOUBLE PRECISION PI,RD,X,RA1,RA2,RA3,DEC1,DEC2,DEC3,DIST,
     CGCA,GCD,GCAA,GCDD,STEP1,STEP2,RAB,DECB,RB,
     CWATCHRA,WATCHDEC,WATCHR,COMPTELRA,COMPTELDEC,COMPTELR,
     CEGRETRA,EGRETDEC,EGRETR,PHEBUSRA,PHEBUSDEC,PHEBUSR,
     CSIGMARA,SIGMADEC,SIGMAR,SAXRA,SAXDEC,SAXR,HRAB,HDECB,HRB,
     CKONUSEL1,KONUSEL2,PRA,PDEC,PR,ORAB,ODECB,ORB,RAEL,DECEL,
     CAREA1,AREA2,FRACTION,ARE,CHI,D,DTHETA,DTHETA2,THETA,THETA2,T,
     CSCHI,CCHI,SD,CD,ALPH,DELT,T2,CHI1,SCHI1,CCHI1,ANG,RACIRC,
     CDECCIRC,RCIRC,ELMIN,RATEST,DECTEST
      COMMON/GEOM/PI,RD
      LOGICAL KON,IN
      AREA1=0.
      NAREA=0
C.....DEFINE IN TO BE TRUE
      IN=.TRUE.
	IF(KON)GOTO 20
      IF(RCIRC.EQ.0.)THEN
      PRINT *,'RACIRC,DECCIRC,RCIRC=',RACIRC,DECCIRC,RCIRC
      STOP 'RCIRC=0 IN GC16'
      ENDIF
20    CONTINUE
C.....THETA IS THE AZIMUTHAL ANGLE (AROUND THE CIRCLE). DTHETA IS THE ANGULAR STEP
C.....IN THE AZIMUTHAL DIRECTION
      NTIMES=1
C.....THIS CALCULATES THE RA, DEC VALUES WHICH ARE WITHIN THE LIMITS OF THE ERROR CIRCLE
C.....THIS FIRST LOOP TRACES THE RIGHT HALF OF A CIRCLE FROM SMALLER DEC
C.....VALUES TO LARGER.  T GOES FROM +1 TO -1
C.....FIRST STEP: OUTER RADIUS OF THE ANNULUS
C.....CHI IS THE HALF ANGLE OF THE ANNULUS (RADIANS),
C.....AND THE NORTH POLE IS AT RA=X DEGREES,DEC=D RADIANS
C.....WORK FROM CENTER OF CIRCLE OUTWARDS
      DTHETA2=PI*STEP2
      CHI=DTHETA2
      X=RACIRC
C.....FOR THE SPECIAL CASE WHERE KONUS EL DEFINES THE REGION,
C.....USE THE SMALLER VALUE OF 90-EL TO START, AND GO TO THE
C.....LARGER ONE
      IF(KON)THEN
	ELMIN=DMIN1(90.D0-KONUSEL1,90.D0-KONUSEL2)
	CHI=ELMIN*RD+DTHETA2
	RCIRC=DMAX1(90.D0-KONUSEL1,90.D0-KONUSEL2)
	ENDIF
      D=DECCIRC*RD
      DTHETA=PI*STEP1
5     THETA=0.D0
      T=1.D0
      IMAX=INT(1./STEP1)+1
10    SCHI=DSIN(CHI)
      CCHI=DCOS(CHI)
      SD=DSIN(D)
      CD=DCOS(D)
      DO 100 I=1,IMAX
      THETA=THETA+DTHETA
50    IF(T.EQ.-1.D0) GO TO 105
      T=DCOS(THETA)
      IF(T.LT.-1.D0) T=-1.D0
      GCA=ALPH(SCHI,CCHI,SD,CD,T) + X
      IF(GCA.GT.360.D0)GCA=GCA-360.D0
      GCD=DELT(SCHI,CCHI,SD,CD,T) 
C.....CHECK TO SEE IF THIS POINT IS WITHIN THE BATSE, EGRET, COMPTEL, ...
C.....ERROR CIRCLES, WITHIN THE KONUS ECLIPTIC LATITUDE BAND, OUTSIDE
C.....THE PLANET-BLOCKING CIRCLE, ETC.  IF THE POINT IS INSIDE, CONTINUE
C.....TESTING; IF OUTSIDE, GET ANOTHER POINT
      IF(RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,RAB,DECB,DIST)
      IF(DIST.GT.SQRT(RB**2+2.56))GOTO 100
      ENDIF
      IF(WATCHRA.NE.0. .AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,WATCHRA,WATCHDEC,DIST)
      IF(DIST.GT.WATCHR)GOTO 100
      ENDIF
      IF(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.)
     CTHEN
      CALL ANGLE1(GCA,GCD,COMPTELRA,COMPTELDEC,DIST)
      IF(DIST .GT. COMPTELR)GOTO 100
      ENDIF
      IF(EGRETRA.NE.0. .AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,EGRETRA,EGRETDEC,DIST)
      IF(DIST.GT.EGRETR)GOTO 100
      ENDIF
      IF(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,PHEBUSRA,PHEBUSDEC,DIST)
      IF(DIST.GT.PHEBUSR)GOTO 100
      ENDIF
      IF(SIGMARA.NE.0. .AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,SIGMARA,SIGMADEC,DIST)
      IF(DIST.GT.SIGMAR)GOTO 100
      ENDIF
      IF(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,SAXRA,SAXDEC,DIST)
      IF(DIST.GT.SAXR)GOTO 100
      ENDIF
      IF(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,HRAB,HDECB,DIST)
      IF(DIST.GT.HRB)GOTO 100
      ENDIF
C.....SKIP THIS TEST IF WE ARE GOING AROUND THE KONUS ANNULUS
      IF(KON)GOTO 55
      IF(KONUSEL1 .NE.0. .OR. KONUSEL2.NE.0.)THEN
      RAEL=270.D0
      DECEL=66.55D0
      REL=90.D0-(KONUSEL1+KONUSEL2)/2.
      DREL=ABS((90.-KONUSEL1)-REL)
      CALL ANGLE1(GCA,GCD,RAEL,DECEL,DIST)
      IF(DIST.GT. (REL+DREL) .OR. DIST .LT. (REL-DREL) )GOTO 100
      ENDIF
55    CONTINUE
      IF(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,PRA,PDEC,DIST)
C.....THIS ANGLE HAS TO BE GREATER THAN PR IF THE POINT IS OUTSIDE THE PLANET
      IF(DIST .LT. PR)GOTO 100
      ENDIF
      IF(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.)THEN
      CALL ANGLE1(GCA,GCD,ORAB,ODECB,DIST)
      IF(DIST.GT.ORB)GOTO 100
      ENDIF
C.....THE POINT HAS PASSED ALL THE TESTS.  INCREMENT THE AREA.
C.....THE ELEMENTAL AREA IS DEFINED BY 3 RA, DEC POINTS.  THE
C.....FIRST TWO ARE ON THE CIRCLE.  THE THIRD IS OBTAINED BY
C.....INCREMENTING THE RADIAL DISTANCE
      NAREA=NAREA+1
      RA1=GCA
      DEC1=GCD
      THETA2=THETA+DTHETA		     
      T2=DCOS(THETA2)
      IF(T2.LT.-1.) T2=-1.
      RA2=ALPH(SCHI,CCHI,SD,CD,T2) + X
      IF(RA2.GT.360.D0)RA2=RA2-360.D0
      DEC2=DELT(SCHI,CCHI,SD,CD,T2) 
      CHI1=CHI+DTHETA2
      SCHI1=DSIN(CHI1)
      CCHI1=DCOS(CHI1)
      RA3=ALPH(SCHI1,CCHI1,SD,CD,T2) + X
      IF(RA3.GT.360.D0)RA3=RA3-360.D0
      DEC3=DELT(SCHI1,CCHI1,SD,CD,T2) 
C      CALL AREA(RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE)
C	   PRINT *,RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE
C.....APPROXIMATE THE ELEMENTAL AREA USING A PLANE TRIANGLE
      CALL ANGLE1(RA1,DEC1,RA2,DEC2,ANG)
      ARE=ANG*DTHETA2/RD
C	   PRINT *,RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE
      AREA1=AREA1+ARE*3600.
100   CONTINUE
105   THETA=PI
      T=COS(THETA+DTHETA)
C.....THIS SECOND LOOP TRACES THE LEFT HALF OF A CIRCLE GOING FROM LARGER
C.....DEC VALUES TO SMALLER.  T GOES FROM -1 TO +1
      T=-1.D0
      DO 200 J=1,IMAX
      THETA=THETA+DTHETA
      IF(T.EQ.1.D0) GO TO 202
      T=DCOS(THETA)
      IF(T.GT.1.D0) T=1.D0
      GCAA=360.D0-ALPH(SCHI,CCHI,SD,CD,T)+X
      IF(GCAA.GT.360.D0)GCAA=GCAA-360.D0
      GCDD=DELT(SCHI,CCHI,SD,CD,T)
C.....CHECK TO SEE IF THIS POINT IS WITHIN THE BATSE, EGRET, COMPTEL, ...
C.....ERROR CIRCLES, WITHIN THE KONUS ECLIPTIC LATITUDE BAND, OUTSIDE
C.....THE PLANET-BLOCKING CIRCLE, ETC.  IF THE POINT IS INSIDE, CONTINUE
C.....TESTING; IF OUTSIDE, GET ANOTHER POINT
      IF(RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,RAB,DECB,DIST)
      IF(DIST.GT.SQRT(RB**2+2.56))GOTO 200
      ENDIF
      IF(WATCHRA.NE.0. .AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,WATCHRA,WATCHDEC,DIST)
      IF(DIST.GT.WATCHR)GOTO 200
      ENDIF
      IF(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.)
     CTHEN
      CALL ANGLE1(GCAA,GCDD,COMPTELRA,COMPTELDEC,DIST)
      IF(DIST .GT. COMPTELR)GOTO 200
      ENDIF
      IF(EGRETRA.NE.0. .AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,EGRETRA,EGRETDEC,DIST)
      IF(DIST.GT.EGRETR)GOTO 200
      ENDIF
      IF(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,PHEBUSRA,PHEBUSDEC,DIST)
      IF(DIST.GT.PHEBUSR)GOTO 200
      ENDIF
      IF(SIGMARA.NE.0. .AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,SIGMARA,SIGMADEC,DIST)
      IF(DIST.GT.SIGMAR)GOTO 200
      ENDIF
      IF(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,SAXRA,SAXDEC,DIST)
      IF(DIST.GT.SAXR)GOTO 200
      ENDIF
      IF(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,HRAB,HDECB,DIST)
      IF(DIST.GT.HRB)GOTO 200
      ENDIF
C.....SKIP THIS TEST IF WE ARE GOING AROUND THE KONUS ANNULUS
      IF(KON)GOTO 155
      IF(KONUSEL1 .NE.0. .OR. KONUSEL2.NE.0.)THEN
      RAEL=270.D0
      DECEL=66.55D0
      REL=90.-(KONUSEL1+KONUSEL2)/2.
      DREL=ABS((90.-KONUSEL1)-REL)
      CALL ANGLE1(GCAA,GCDD,RAEL,DECEL,DIST)
      IF(DIST.GT. (REL+DREL) .OR. DIST .LT. (REL-DREL) )GOTO 200
      ENDIF
155   CONTINUE
      IF(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,PRA,PDEC,DIST)
C.....THIS ANGLE HAS TO BE GREATER THAN PR IF THE POINT IS OUTSIDE THE PLANET
      IF(DIST .LT. PR)GOTO 200
      ENDIF
      IF(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.)THEN
      CALL ANGLE1(GCAA,GCDD,ORAB,ODECB,DIST)
      IF(DIST.GT.ORB)GOTO 200
      ENDIF
C.....THE POINT HAS PASSED ALL THE TESTS.  INCREMENT THE AREA.
C.....THE ELEMENTAL AREA IS DEFINED BY 3 RA, DEC POINTS.
C.....RA1, DEC1 IS THE POINT JUST TREATED
C.....RA2, DEC2 IS THE POINT ALONG THE CIRCLE
C.....RA3, DEC3 IS THE POINT IN THE RADIAL DIRECTION
      NAREA=NAREA+1
      RA1=GCAA
      DEC1=GCDD
      THETA2=THETA+DTHETA
      T2=DCOS(THETA2)
      IF(T2.LT.-1.D0) T2=-1.D0
      RA2=360.D0-ALPH(SCHI,CCHI,SD,CD,T2) + X
      IF(RA2.GT.360.D0)RA2=RA2-360.D0
      DEC2=DELT(SCHI,CCHI,SD,CD,T2) 
      CHI1=CHI+DTHETA2
      SCHI1=DSIN(CHI1)
      CCHI1=DCOS(CHI1)
      RA3=360.D0-ALPH(SCHI1,CCHI1,SD,CD,T2) + X
      IF(RA3.GT.360.D0)RA3=RA3-360.D0
      DEC3=DELT(SCHI1,CCHI1,SD,CD,T2) 
C      RA3=RA2
C      DEC3=DEC2+DTHETA2
C      CALL AREA(RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE)
C.....APPROXIMATE THE ELEMENTAL AREA USING A PLANE TRIANGLE
      CALL ANGLE1(RA1,DEC1,RA2,DEC2,ANG)
      ARE=ANG*DTHETA2/RD
C	   PRINT *,RA1,DEC1,RA2,DEC2,RA3,DEC3,ARE
      AREA1=AREA1+ARE*3600. 
200   CONTINUE
C.....NOW INCREMENT THE RADIAL ANGLE AND GO AROUND THE CIRCLE AGAIN
202   CHI=CHI+DTHETA2
      J=J+1
      IF(CHI.GE.RCIRC*RD)GOTO 205
      X=RACIRC
      D=DECCIRC*RD
      GOTO 5
C.....AREA1 IS IN SQUARE ARCMINUTES.  CONVERT IT TO STERADIANS.
C.....AREA2 IS IN SQUARE DEGREES
205   AREA2=AREA1/3600.D0
      AREA1=AREA1*RD**2/3600.D0
      FRACTION=AREA1/(4*PI)
C.....NOW CHECK THE POINT RATEST,DECTEST
      IF(RAB.NE.0. .AND. DECB.NE.0. .AND. RB.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,RAB,DECB,DIST)
      IF(DIST.GT.SQRT(RB**2+2.56))IN=.FALSE.
	IF(.NOT.IN)RETURN
      ENDIF
      IF(WATCHRA.NE.0. .AND. WATCHDEC.NE.0. .AND. WATCHR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,WATCHRA,WATCHDEC,DIST)
      IF(DIST.GT.WATCHR)IN=.FALSE.
	IF(.NOT.IN)RETURN
      ENDIF
      IF(COMPTELRA.NE.0. .AND. COMPTELDEC.NE.0. .AND. COMPTELR.NE.0.)
     CTHEN
      CALL ANGLE1(RATEST,DECTEST,COMPTELRA,COMPTELDEC,DIST)
      IF(DIST .GT. COMPTELR)IN=.FALSE.
	IF(.NOT.IN)RETURN
      ENDIF
      IF(EGRETRA.NE.0. .AND. EGRETDEC.NE.0. .AND. EGRETR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,EGRETRA,EGRETDEC,DIST)
      IF(DIST.GT.EGRETR)IN=.FALSE.
	IF(.NOT.IN)RETURN
      ENDIF
      IF(PHEBUSRA.NE.0. .AND. PHEBUSDEC.NE.0. .AND. PHEBUSR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,PHEBUSRA,PHEBUSDEC,DIST)
      IF(DIST.GT.PHEBUSR)IN=.FALSE.
	IF(.NOT.IN)RETURN
      ENDIF
      IF(SIGMARA.NE.0. .AND. SIGMADEC.NE.0. .AND. SIGMAR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,SIGMARA,SIGMADEC,DIST)
      IF(DIST.GT.SIGMAR)IN=.FALSE.
	IF(.NOT.IN)RETURN
      ENDIF
      IF(SAXRA.NE.0. .AND. SAXDEC.NE.0. .AND. SAXR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,SAXRA,SAXDEC,DIST)
      IF(DIST.GT.SAXR)IN=.FALSE.
	IF(.NOT.IN)RETURN
      ENDIF
      IF(HRAB.NE.0. .AND. HDECB.NE.0. .AND. HRB.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,HRAB,HDECB,DIST)
      IF(DIST.GT.HRB)IN=.FALSE.
	IF(.NOT.IN)RETURN
      ENDIF
      IF(KONUSEL1 .NE.0. .OR. KONUSEL2.NE.0.)THEN
      RAEL=270.D0
      DECEL=66.55D0
      REL=90.-(KONUSEL1+KONUSEL2)/2.
      DREL=ABS((90.-KONUSEL1)-REL)
      CALL ANGLE1(RATEST,DECTEST,RAEL,DECEL,DIST)
      IF(DIST.GT. (REL+DREL) .OR. DIST .LT. (REL-DREL) )IN=.FALSE.
	IF(.NOT.IN)RETURN
      ENDIF
      IF(PRA.NE.0. .AND. PDEC.NE.0. .AND. PR.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,PRA,PDEC,DIST)
C.....THIS ANGLE HAS TO BE GREATER THAN PR IF THE POINT IS OUTSIDE THE PLANET
      IF(DIST .LT. PR)IN=.FALSE.
	IF(.NOT.IN)RETURN
      ENDIF
      IF(ORAB.NE.0. .AND. ODECB.NE.0. .AND. ORB.NE.0.)THEN
      CALL ANGLE1(RATEST,DECTEST,ORAB,ODECB,DIST)
      IF(DIST.GT.ORB)IN=.FALSE.
	IF(.NOT.IN)RETURN
      ENDIF
      RETURN
      END
      DOUBLE PRECISION FUNCTION ALPH(SCHI,CCHI,SD,CD,T)
C.....PARAMETRIC EQUATION FOR RIGHT ASCENSION VALUES ALONG A CIRCLE
C.....T IS THE PARAMETER - IT VARIES BETWEEN -1 AND +1.  ALPH IS IN
C.....DEGREES
C     IMPLICIT DOUBLE PRECISION (A-Z)
      DOUBLE PRECISION SCHI,CCHI,SD,CD,T,NUM,DEN,A,B,PI,RD
      COMMON/GEOM/PI,RD
C      REAL NUM
      A=T*SD*SCHI
	B=CD*CCHI
      NUM=T*SD*SCHI + CD*CCHI
      DEN=DSQRT(SCHI**2 + (CD*CCHI)**2-(T*CD*SCHI)**2
     C+2*T*CCHI*SCHI *CD*SD)
      IF(DABS(NUM-DEN).LT.1.D-10) GO TO 10
      ALPH=DACOS(DSIGN(DMIN1(1.D0,DABS(NUM/DEN)),NUM))/RD
      RETURN
10    IF(NUM.GT.0.D0) ALPH=0.D0
      IF(NUM.LT.0.D0)ALPH=180.D0
      RETURN
      END
      DOUBLE PRECISION FUNCTION DELT(SCHI,CCHI,SD,CD,T)
C     PARAMETRIC EQUATION FOR DELTA VALUES OF CIRCLE
C     WHOSE CENTER IS AT DELTAO.
C.....DELT IS IN DEGREES
      COMMON/GEOM/PI,RD
C     IMPLICIT DOUBLE PRECISION (A-Z)
      DOUBLE PRECISION SCHI,CCHI,SD,CD,T,PI,RD
      DELT=DASIN(CCHI*SD - T*CD*SCHI)/RD
      RETURN
      END
