C  This program creates the unformatted file containing the logical array       
C  which gives the excluded zones.                                              
                                                                                
      program write_lmask                                                       
      logical*1 lmask(41167)                                                    
      do i = 1, 41167                                                           
         lmask(i) = .false.                                                     
      enddo                                                                     
      open (2,file='exclude.lis',status='old')                                  
      do i = 1, 41167                                                           
         read (2,*,end=10) j                                                    
         lmask(j) = .true.                                                      
      enddo                                                                     
 10   open( unit=10, file='lmask.lis', form='unformatted', status='unknown')    
      write (10) lmask           !array of lune bins                            
      close (10)                                                                
      stop                                                                      
      end                                                                       
                                                                                
                                                                                
                                                                                
C  This subroutine reads in the excluded region mask.                           
      subroutine init_lmask                                                     
      logical*1 lmask(41167)                                                    
      common /exclude_mask/ lmask                                               
      open( unit=10, file='lmask.lis', form='unformatted', status='old')        
      read (10) lmask           !array of lune bins                             
      close (10)                                                                
      return                                                                    
      end                                                                       
                                                                                
                                                                                
C  This logical function returns .true. if the source is in the excluded        
C  regions, and .false. otherwise. Input ra and dec are real, in hours and      
C  degrees.                                                                     
                                                                                
      logical function exclude(ra,dec)                                          
      double precision dl2, db2, dqr, dqd, delon, delat                         
      logical*1 lmask(41167)                                                    
      common /exclude_mask/ lmask                                               
      call eqgal(ra,dec,tl2,tb2)                                                
      db2 = dble(tb2)                                                           
      dl2 = dble(tl2)                                                           
      call g2qe(dl2, db2, dqr, dqd, delon, delat, 2)                            
      n = 0                                                                     
      call p2bin(delon, delat, n) !lune bin number of source                    
      exclude = .false.                                                         
      if (lmask(n) .or. abs(tb2) .le. 5) exclude = .true.                       
      return                                                                    
      end                                                                       
                                                                                
                                                                                
C  The remaining subroutines and functions are used in the conversions          
C  between various coordinate systems.                                          
	SUBROUTINE EQGAL(RA, DEC, GLG, GLT)                                            
C  Single precision code, RA in hours, dec in degrees.                          
	DOUBLE PRECISION DRA, DDEC, DGL, DGB, PI, DD                                   
	PI = ACOS(-1.d0)                                                               
	DD = 180.d0/PI                                                                 
	DRA = RA*15.d0/DD                                                              
	DDEC = DEC/DD                                                                  
	CALL CEQTGA(DRA, DDEC, DGL, DGB)                                               
	GLG = DGL*DD                                                                   
	GLT = DGB*DD                                                                   
	IF (GLG .LT. 0.) GLG = 360. + GLG                                              
	RETURN                                                                         
	END                                                                            
                                                                                
       SUBROUTINE G2QE(GLON,GLAT,QLON,QLAT,ELON,ELAT,IOPT)                      
C-------------------------------------------------------------                  
C                                                                               
C  THIS DOUBLE PRECISION ROUTINE CONVERTS FROM GALACTIC CO-ORDINATES            
C  TO EQUATORIAL AND ECLIPTIC CO-ORDINATES (IF IOPT=2), OR TO                   
C  JUST EQUATORIAL CO-ORDINATES (IF IOPT=1)                                     
C                                                                               
C  ALL CO-ORDINATES ARE IN DEGREES.                                             
C                                                                               
C--------------------------------------------------------------                 
       DOUBLE PRECISION GLON,GLAT,QLON,QLAT,ELON,ELAT,RA2DEG,                   
     # RADLON,RADLAT,ZLON,ZLAT,ANGPOS,ANGSGN                                    
       INTEGER IOPT                                                             
C                                                                               
       DATA RA2DEG/5.72957795D1/                                                
C                                                                               
       RADLON=GLON/RA2DEG                                                       
       RADLAT=GLAT/RA2DEG                                                       
       CALL CEQFGA(RADLON,RADLAT,ZLON,ZLAT)                                     
       QLON=RA2DEG*ANGPOS(ZLON)                                                 
       QLAT=RA2DEG*ANGSGN(ZLAT)                                                 
       IF(QLON.GE.3.6D2)QLON=0.0D0                                              
C                                                                               
       IF(IOPT.EQ.1)RETURN                                                      
C                                                                               
       CALL Q2E(QLON,QLAT,ELON,ELAT)                                            
       RETURN                                                                   
       END                                                                      
                                                                                
       SUBROUTINE P2BIN(ELON,ELAT,NBIN)                                         
       IMPLICIT DOUBLE PRECISION (A-H,P-Z)                                      
CCC	Original Module modified to work in degrees instead of radians.             
CCC                                                                             
CCC       DATA R2D/57.2957795/                                                  
       INTEGER MAXBIN(182)                                                      
       DATA MAXBIN                                                              
     #  /0,1,7,19,37,62,93,130,173,223,279,341,409,483,563,                     
     #  650,743,842,947,1058,1175,1298,1427,1561,1701,1847,1999,                
     #  2156,2319,2488,2662,2841,3026,3216,3412,3613,3819,4030,                 
     #  4246,4467,4693,4924,5160,5400,5645,5895,6149,6407,6670,                 
     #  6937,7208,7483,7762,8045,8332,8623,8917,9215,9516,9821,                 
     #  10129,10440,10754,11071,11391,11714,12040,12368,12699,                  
     #  13032,13368,13706,14046,14388,14732,15078,15425,15774,                  
     #  16124,16476,16829,17183,17538,17894,18251,18609,18967,                  
     #  19326,19685,20044,20403,20763,21122,21481,21840,22199,                  
     #  22557,22915,23272,23628,23983,24337,24690,25042,25392,                  
     #  25741,26088,26434,26778,27120,27460,27798,28134,28467,                  
     #  28798,29126,29452,29775,30095,30412,30726,31037,31345,                  
     #  31650,31951,32249,32543,32834,33121,33404,33683,33958,                  
     #  34229,34496,34759,35017,35271,35521,35766,36006,36242,                  
     #  36473,36699,36920,37136,37347,37553,37754,37950,38140,                  
     #  38325,38505,38679,38848,39011,39168,39320,39466,39606,                  
     #  39740,39869,39992,40109,40220,40325,40424,40517,40604,                  
     #  40684,40758,40826,40888,40944,40994,41037,41074,41105,                  
     #  41130,41148,41160,41166,41167/                                          
C                                                                               
CCC       I=90.0-ELAT*R2D+2.5                                                   
       I=90.0-ELAT+2.5                                                          
       N=MAXBIN(I)-MAXBIN(I-1)                                                  
       SIZE=360.0/REAL(N)                                                       
CCC       J=ELON*R2D/SIZE+1                                                     
       J=ELON/SIZE+1                                                            
       NBIN=MAXBIN(I-1)+J                                                       
       RETURN                                                                   
       END                                                                      
                                                                                
      SUBROUTINE CEQTEC(Y,RA,DEC,LNG,LAT)                                       
C       EQUATORIAL COORDINATES AT EPOCH Y TO ECLIPTIC COORDINATES               
        DOUBLE PRECISION Y,RA,DEC,LNG,LAT,E1,E2,E3,E4,PIBY2,                    
     #  DTR,T,E,EQP,PEQ,EQ,ANGSGN                                               
C                                                                               
      DATA E1,E2,E3,E4/2.3452294D1,-1.30125D-2,-1.64D-6,5.03D-7/                
      DATA PIBY2/1.5707963268D0/                                                
                                                                                
      DTR=PIBY2/90.0                                                            
      T=(Y-1900.0)/100.0                                                        
      E=DTR*(E1+T*(E2+T*(E3+T*E4)))                                             
      CALL SPTR1(RA+PIBY2,PIBY2-DEC,E,EQP,PEQ,EQ)                               
      LNG=ANGSGN(PIBY2-PEQ)                                                     
      LAT=PIBY2-EQ                                                              
                                                                                
      RETURN                                                                    
      END                                                                       
                                                                                
      SUBROUTINE CEQTGA(RA,DEC,GLG,GLT)                                         
C       RA AND DEC 1950.0 TO GALACTIC LONGITUDE AND LATITUDE                    
       DOUBLE PRECISION RA,DEC,GLG,GLT,PIBY2,RAG,DECG,PGL,DTR,                  
     # PGQ,GQP,GQ,ANGSGN                                                        
      DATA PIBY2/1.5707963268D0/                                                
      DATA RAG,DECG,PGL /1.9225D2, 2.74D1, 1.23D2/                              
C                                                                               
      DTR=PIBY2/90.0                                                            
      CALL SPTR1(RA-RAG*DTR,PIBY2-DECG*DTR,PIBY2-DEC,PGQ,GQP,GQ)                
      GLG=ANGSGN(PGL*DTR-PGQ)                                                   
      GLT=PIBY2-GQ                                                              
      RETURN                                                                    
      END                                                                       
                                                                                
      SUBROUTINE SPTR1(A,AB,AC,B,C,BC)                                          
C       TO SOLVE A SPHERICAL TRIANGLE GIVEN TWO SIDES AND THE INCLUDED          
C       ANGLE                                                                   
       DOUBLE PRECISION A,AB,AC,B,C,BC,Z,SA,CA,SAB,CAB,SAC,CAC,                 
     # CBC,CS,SS,SBCSQ,SBC,SB,CB,SC,CC,DSIN,DCOS,DSQRT,DABS                     
C                                                                               
      DATA Z/0.0D0/                                                             
C                                                                               
      SA=DSIN(A)                                                                
      CA=DCOS(A)                                                                
      SAB=DSIN(AB)                                                              
      CAB=DCOS(AB)                                                              
      SAC=DSIN(AC)                                                              
      CAC=DCOS(AC)                                                              
      CBC=CAB*CAC+SAB*SAC*CA                                                    
C                                                                               
      CS=DCOS(5.0D-1*A)*DSIN(AB-AC)                                             
      SS=DSIN(5.0D-1*A)*DSIN(AB+AC)                                             
      SBCSQ=CS*CS+SS*SS+SAB*SAB*SAC*SAC*SA*SA                                   
      SBC=DSQRT(SBCSQ)                                                          
C                                                                               
      SB=SAC*SA                                                                 
      CB=CAC*SAB-CAB*SAC*CA                                                     
      SC=SAB*SA                                                                 
      CC=CAB*SAC-CAC*SAB*CA                                                     
C                                                                               
      IF ((SB.NE.Z.OR.CB.NE.Z).AND.(SC.NE.Z.OR.CC.NE.Z)) GO TO 5                
C       SPECIAL CASES                                                           
      SB=SA                                                                     
      CB=CA                                                                     
      SC=0.0                                                                    
      CC=-1.0                                                                   
      IF (DABS(SAB)+DABS(SAC).GT.DABS(SA)) GO TO 5                              
      CB=-CAB*CA                                                                
      CC=CAB                                                                    
C                                                                               
    5 B=DATAN2(SB,CB)                                                           
      C=DATAN2(SC,CC)                                                           
      BC=DATAN2(SBC,CBC)                                                        
      RETURN                                                                    
      END                                                                       
      FUNCTION ANGSGN(X)                                                        
C       TO PUT AN ANGLE INTO THE RANGE -PI TO PI                                
C       IE  -PI LE ANGSGN LT PI                                                 
       DOUBLE PRECISION ANGSGN,X,PI,Y,DMOD                                      
                                                                                
C                                                                               
      DATA PI/3.1415926536D0/                                                   
C                                                                               
      Y=DMOD(X+PI,2.0*PI)                                                       
      IF (Y) 11,10,10                                                           
   11 Y=Y+2.0*PI                                                                
   10 ANGSGN=Y-PI                                                               
      RETURN                                                                    
      END                                                                       
                                                                                
       SUBROUTINE Q2E(QLON,QLAT,ELON,ELAT)                                      
C---------------------------------------------------------------------          
C                                                                               
C  THIS DOUBLE PRECISION ROUTINE CONVERTS FROM EQUATORIAL CO-ORDINATES          
C  TO ECLIPTIC CO-ORDINATES. ALL POSITIONS ARE IN DEGREES.                      
C                                                                               
C---------------------------------------------------------------------          
       DOUBLE PRECISION QLON,QLAT,ELON,ELAT,RA2DEG,RADLON,RADLAT,               
     # ZLON,ZLAT,EPOCH,ANGPOS,ANGSGN                                            
C                                                                               
       DATA RA2DEG/5.72957795D1/,EPOCH/1.95D3/                                  
C                                                                               
       RADLON=QLON/RA2DEG                                                       
       RADLAT=QLAT/RA2DEG                                                       
       CALL CEQTEC(EPOCH,RADLON,RADLAT,ZLON,ZLAT)                               
       ELON=RA2DEG*ANGPOS(ZLON)                                                 
       ELAT=RA2DEG*ANGSGN(ZLAT)                                                 
       IF(ELON.GE.3.6D2)ELON=0.0D0                                              
       RETURN                                                                   
       END                                                                      
                                                                                
      SUBROUTINE CEQFGA(GLG,GLT,RA,DEC)                                         
C       RA AND DEC 1950.0 FROM GALACTIC LONGITUDE AND LATITUDE                  
       DOUBLE PRECISION GLG,GLT,RA,DEC,RAG,DECG,PGL,PIBY2,                      
     # DTR,GQP,QPG,PQ,ANGSGN                                                    
C                                                                               
      DATA RAG,DECG,PGL /1.9225D2, 2.74D1, 1.23D2/                              
      DATA PIBY2/1.5707963268D0/                                                
C                                                                               
      DTR=PIBY2/90.0                                                            
      CALL SPTR1(PGL*DTR-GLG,PIBY2-GLT,PIBY2-DECG*DTR,GQP,QPG,PQ)               
      RA=ANGSGN(QPG+RAG*DTR)                                                    
      DEC=PIBY2-PQ                                                              
                                                                                
      RETURN                                                                    
      END                                                                       
                                                                                
      FUNCTION ANGPOS(X)                                                        
C       TO PUT AN ANGLE INTO THE RANGE ZERO TO TWOPI                            
C       IE  0 LE ANGPOS LT TWOPI                                                
       DOUBLE PRECISION ANGPOS,X,TWOPI,AMOD1                                    
C                                                                               
      DATA TWOPI/6.2831853072D0/                                                
                                                                                
      ANGPOS=AMOD1(X,TWOPI)                                                     
                                                                                
      RETURN                                                                    
      END                                                                       
      FUNCTION AMOD1(X,Y)                                                       
C       AMOD1 HAS THE SAME SIGN AS THE 2ND ARGUMENT                             
C       AMOD HAS THE SAME SIGN AS THE 1ST ARGUMENT                              
       DOUBLE PRECISION AMOD1,X,Y,DMOD                                          
                                                                                
      AMOD1=DMOD(X,Y)                                                           
      IF (AMOD1/Y) 11,10,10                                                     
   11 AMOD1=AMOD1+Y                                                             
C                                                                               
   10 RETURN                                                                    
      END                                                                       
