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