ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c match a list of ra,dec with the spm4 catalog c this program works with the rsNN.asc version of the spm4 c c Note: modify "aroot" to reflect the catalog directory path c c 2009-11-05 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc implicit real*8 (a-h,o-z) parameter (nmax=9e6) ! max no. of positions in input list dimension ras(nmax),decs(nmax) character*20 label,labels(nmax) integer*4 id,ira2,idec2 integer*2 ifld,ibmag,ivmag,ira1,idec1,ierr,ipoch,irest character*90 fnin,fnout,aword character*8 afile,afind dimension lines(0:360) character*144 aline integer labs(2),lrs(2),lds(2) integer lrhs(2),lrms(2),lrss(2),ldds(2),ldms(2),ldss(2) c character*2 aroot data aroot/'./'/ c c----- parse command line c pi=4.0d0*atan(1.0d0) raddeg=pi/180.0d0 c iarg=iargc() if(iarg.lt.3) then print*,'usage: match_spm4 input.file tolsec output.file' stop endif c call getarg(1,fnin) c call getarg(2,aword) read(aword,*)tolsec tol=tolsec/3600.0d0 ! matching tolerance in degrees c fudge=1.20d0 del=fudge*tol delsec=fudge*tolsec c call getarg(3,fnout) c c----- get ra,dec format character string for input list c open(10,file=fnin,status='old') c print*,'Enter aaaaa, rrr, ddd (in degrees)' print*,'or AAAAAA, HH MM SSS, +DD MM SSS' c do i=1,5 read(10,'(a)',end=7)aword kk=length(aword) if(aword(kk:kk).eq.achar(13)) kk=kk-1 write(*,'(a)')aword(1:kk) enddo 7 continue c rewind(10) c read(*,'(a90)')aword c call parse(aword,ibab,labs,lrs,lds, + lrhs,lrms,lrss,lsig,ldds,ldms,ldss,lmax) c c----- store ra,dec from input list file c i=0 decmin=+95.0d0 decmax=-95.0d0 decnorth=-20.0d0+del 8 continue read(10,'(a)',end=9)aword(1:lmax) call decrypt(aword,ibab,labs,lrs,lds, + lrhs,lrms,lrss,lsig,ldds,ldms,ldss,label,ra,dec) if(ra.lt.0.0d0) goto 8 if(dec.gt.decnorth) goto 8 ! star is above northern border i=i+1 ras(i)=ra decs(i)=dec if(dec.lt.decmin) decmin=dec if(dec.gt.decmax) decmax=dec labels(i)=label if(i.gt.nmax) stop 'nmax not big enough!' goto 8 9 continue c close(10) imax=i print*,imax,' targets read in and stored' c decmin=decmin-tol decmax=decmax+tol decval=abs(decmin) ibeg=int(decval) decval=abs(decmax) iend=int(decval) c c----- search spm4 dec strips in order, extract and write matches c open(22,file=fnout,status='unknown') c c nout=0 fudge=1.20d0 del=fudge*tol delsec=fudge*tolsec c do kfile=ibeg,iend,-1 c write(afile,100)kfile,'.asc' 100 format('rs',i2.2,a4) print*,'processing dec band ',afile,' ...' open(12,file=aroot//afile,status='old',err=20 + ,access='direct',recl=144) write(afind,100)kfile,'.ind' open(13,file=aroot//afind,status='old') do k=1,360 read(13,*)ii,line,ralim lines(ii)=line enddo close(13) lines(0)=1 declo=-1.0d0*(kfile+1) dechi=-1.0d0*kfile c do i=1,imax tra=ras(i) tdec=decs(i) beta=tra*raddeg gamma=tdec*raddeg if(tdec.ge.declo-del.and.tdec.le.dechi+del) then xx=-tolsec yy=0.0d0 call stdrad(beta,gamma,xx,yy,rar,decr) ra=rar/raddeg west=ra xx=tolsec yy=0.0d0 call stdrad(beta,gamma,xx,yy,rar,decr) ra=rar/raddeg east=ra isplit=0 if(west.gt.east) isplit=1 if(kfile.eq.89) then isplit=0 west=0.0d0 east=360.0d0 else xx=-tolsec yy=0.0d0 call stdrad(beta,gamma,xx,yy,rar,decr) ra=rar/raddeg west=ra xx=tolsec yy=0.0d0 call stdrad(beta,gamma,xx,yy,rar,decr) ra=rar/raddeg east=ra isplit=0 if(west.gt.east) isplit=1 endif decval=abs(decmin) kbeg=int(decval) decval=abs(decmax) kend=int(decval) if(isplit.eq.1) then ir1=1 ideg=int(east+1.0d0) if(ideg.gt.lines(360)) ideg=360 ir2=lines(ideg) ideg=int(west) ir3=lines(ideg) ir4=lines(360) else ir1=1 ir2=0 ideg=int(west) if(ideg.lt.0) ideg=0 ir3=lines(ideg) ideg=int(east+1.0d0) if(ideg.gt.360) ideg=360 ir4=lines(ideg) endif do ir=ir1,ir4 if(ir.le.ir2.or.ir.ge.ir3) then read(12,rec=ir)aline read(aline,201)ra,dec,ex,ey,pma,pmd,epmx,epmy,b,v,ibv, + ta,t1,t2,mp,np,nc,ifld,imas,igc,fj,fh,fk 201 format(2f12.7,2f6.1,2f9.2,2f7.2,2f6.2,i3.2, + 3f6.2,i3,i2,i2,1x,i3.3,i7.7,i3.2,3f7.3, + f7.3,1x,a20) rar=ra*raddeg decr=dec*raddeg call radstd(rar,decr,beta,gamma,xi,eta) sep=sqrt(xi*xi+eta*eta) if(sep.le.tolsec) then nout=nout+1 write(22,201)ra,dec,ex,ey,pma,pmd,epmx,epmy,b,v,ibv, + ta,t1,t2,mp,np,nc,ifld,imas,igc,fj,fh,fk, + sep,labels(i) endif endif enddo endif enddo c close(12) print*,' current matches total = ',nout 20 continue enddo c close(11) c print*,nout,' grand total matches extracted' c c----- all done c stop end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine parse(aword,ibab,labs,lrs,lds, + lrhs,lrms,lrss,lsig,ldds,ldms,ldss,lmax) c implicit real*8 (a-h,o-z) character*90 aword integer labs(2),lrs(2),lds(2) integer lrhs(2),lrms(2),lrss(2),ldds(2),ldms(2),ldss(2) c c----- c ibab=0 kk=index(aword,'a') if(kk.eq.0) ibab=1 c if(ibab.eq.0) then do i=1,90 if(aword(i:i).eq.'a') labs(2)=i if(aword(i:i).eq.'r') lrs(2)=i if(aword(i:i).eq.'d') lds(2)=i enddo do i=90,1,-1 if(aword(i:i).eq.'a') labs(1)=i if(aword(i:i).eq.'r') lrs(1)=i if(aword(i:i).eq.'d') lds(1)=i enddo endif c if(ibab.eq.1) then lsig=index(aword,'+') if(lsig.eq.0) stop 'no + sign!' do i=1,90 if(aword(i:i).eq.'A') labs(2)=i enddo do i=90,1,-1 if(aword(i:i).eq.'A') labs(1)=i enddo do i=1,lsig if(aword(i:i).eq.'H') lrhs(2)=i if(aword(i:i).eq.'M') lrms(2)=i if(aword(i:i).eq.'S') lrss(2)=i enddo do i=lsig,1,-1 if(aword(i:i).eq.'H') lrhs(1)=i if(aword(i:i).eq.'M') lrms(1)=i if(aword(i:i).eq.'S') lrss(1)=i enddo do i=lsig,90 if(aword(i:i).eq.'D') ldds(2)=i if(aword(i:i).eq.'M') ldms(2)=i if(aword(i:i).eq.'S') ldss(2)=i enddo do i=90,lsig,-1 if(aword(i:i).eq.'D') ldds(1)=i if(aword(i:i).eq.'M') ldms(1)=i if(aword(i:i).eq.'S') ldss(1)=i enddo endif c c----- c lmax=0 lmax=max(lmax,labs(2)) lmax=max(lmax,ldss(2)) lmax=max(lmax,lds(2)) c c----- c return end ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc subroutine decrypt(aword,ibab,labs,lrs,lds, + lrhs,lrms,lrss,lsig,ldds,ldms,ldss,label,ra,dec) c implicit real*8 (a-h,o-z) character*90 aword character*20 label,blank character*1 asig integer labs(2),lrs(2),lds(2) integer lrhs(2),lrms(2),lrss(2),ldds(2),ldms(2),ldss(2) data blank/' '/ c c----- c ra=-1.0d0 c label=blank i1=labs(1) i2=labs(2) label(1:i2-i1+1)=aword(i1:i2) c if(ibab.eq.0) then i1=lrs(1) i2=lrs(2) read(aword(i1:i2),*,err=90)ra i1=lds(1) i2=lds(2) read(aword(i1:i2),*)dec endif c if(ibab.eq.1) then i1=lrhs(1) i2=lrhs(2) read(aword(i1:i2),*,err=90)irh i1=lrms(1) i2=lrms(2) read(aword(i1:i2),*)irm i1=lrss(1) i2=lrss(2) read(aword(i1:i2),*)rsec asig=aword(lsig:lsig) i1=ldds(1) i2=ldds(2) read(aword(i1:i2),*)idd i1=ldms(1) i2=ldms(2) read(aword(i1:i2),*)idm i1=ldss(1) i2=ldss(2) read(aword(i1:i2),*)dsec ra=15.0d0*(irh+irm/60.0d0+rsec/3600.0d0) dec=1.0d0*(idd+idm/60.0d0+dsec/3600.0d0) if(asig.eq.'-') dec=-dec endif c c----- c 90 return end CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C RADSTD C C CONVERTS RA AND DEC TO XI AND ETA VALUES C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE RADSTD(RA,DEC,BETA,GAMMA,XI,ETA) C IMPLICIT real*8 (A-H,O-Z) C FF=206264.8062470964D0 C DIV=dSIN(DEC)*dSIN(GAMMA)+dCOS(DEC)*dCOS(GAMMA)*dCOS(RA-BETA) DIV=FF/DIV XI=dCOS(DEC)*dSIN(RA-BETA)*DIV ETA=(dSIN(DEC)*dCOS(GAMMA)-dCOS(DEC)*dSIN(GAMMA) #*dCOS(RA-BETA))*DIV C RETURN END CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C STDRAD C C CONVERTS XI AND ETA TO POSITIONS IN RADIANS C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE STDRAD(BETA,GAMMA,FFX,FFY,ALPHA,DELTA) C IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PI=4D0*DATAN(1D0) FF=206264.8062470964D0 FX=FFX/FF FY=FFY/FF DSG=SIN(GAMMA) DCG=COS(GAMMA) aux=DCG-FY*DSG ALPHA=BETA+DATAN(FX/AUX) if(aux.lt.0d0) alpha=alpha+pi IF(ALPHA.LT.0D0) ALPHA=ALPHA+2D0*PI IF(ALPHA.GT.2D0*PI) ALPHA=ALPHA-2d0*PI SINIS=(DSG+FY*DCG)/SQRT(FX*FX+FY*FY+1.0D0) COSIS=SQRT(1.0D0-SINIS*SINIS) DELTA=ATAN(SINIS/COSIS) C RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc function length(fn) c integer length character*(*) fn c kk=len(fn) c length=0 do i=kk,1,-1 if(fn(i:i).ne.' ') then length=i return endif enddo c return end