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
