PROGRAM  SEARCH
c
c primitive search program for UJ 1.0
c asks user for field center and box width; then goes to
c catalog and extracts relevant data.  Results are either
c displayed or stored in file.
c written 07-Dec-94 AAH
c
c mods 14-Dec-95 AAH: sign prob on output file, leading
c  zeroes in output, case for .ACC files, RA/Dec formatting,
c  allow all zones.
c modified to work with USNO-A 06-Jun-96 AAH
c mod 01-aug-96 aah fix truncation errors
c mod 01-mar-97 aah usno a1.0
c mod 06-mar-97 aah disk lookup table
c
      IMPLICIT
     $  NONE
      REAL*8
     $  RFACTOR,DFACTOR
      INTEGER*4
     $  HR24,NZONES
      PARAMETER
     $  (RFACTOR = 15.0D00*3600.0D00*100.0D00,
     $  DFACTOR =         3600.0D00*100.0D00,
     $  HR24 = 360*3600*100,
     $  NZONES = 24)
      DOUBLE PRECISION
     $  picon,ra,dec
      INTEGER*4
     $  declo,dechi,ralo,rahi,dbeg(4),dend(4),rbeg(4),rend(4),
     $  icentra,icentdec,zlo(NZONES),zhi(NZONES),zone(NZONES),
     $  zonelo,zonehi,z(4),i,nstars,i1,i2,i3,i4,it
      REAL*4
     $   rah,ram,dm,ras,decs,cosdec,dd,centdecd,dsign,
     $   fullwidth,halfwidth
      LOGICAL
     $   fflag,dflag
      CHARACTER
     $   answer*5, fname*80, path*80, path2*80, dds*1
      COMMON /blka/ nstars,path,path2,fflag,dflag
c zone is the spd zone value for each cat file
c zlo, zhi are declination limits each zone in deg*3600*100
c
      picon = datan(1.D0)/45.0
      do i=1,NZONES
        zone(i) = (i-1)*75
        zlo(i) = (i-1)*750*3600
        zhi(i) = zlo(i) + 750*3600 - 1
      enddo
c
c get user input
c
      print *
      print *
      print *
      print *
      print *
      print *
      print *,' U S N O A    S E A R C H    P R O G R A M'
      print *
      print *
      print *,'         Version 1.3  06-Mar-97'
      print *,'For the following questions, an example of the'
      print *,' proper response format is given in the square brackets.'
      print *
      print *,'Enter file path to catalog files [E:]: '
      read (5,'(a)') path
700   continue
      nstars = 0
      print *,'Enter Right Ascension (J2000) of',
     $        ' field center [hh.mmsss]: '
      read (5,*) ra
      it = idnint(ra*100000)
      i1 = it/100000
      rah = float(i1)
      i2 = (it/1000)-i1*100
      ram = float(i2)
      ras = float(it-i2*1000-100000*i1)/10.
      print *,'Enter Declination (J2000) of field center [+-dd.mmsss]: '
      read (5,*) dec
      dsign = 1.0
      dds = ' '
      if (dec.lt.0) then
        dsign = -1.0
        dec = -dec
        dds = '-'
      endif
      it = idnint(dec*100000)
      i3 = it/100000
      dd = float(i3)
      i4 = (it/1000)-i3*100
      dm = float(i4)
      decs = float(it-i4*1000-i3*100000)/10.
      print *,'Enter width of search box in arcseconds [1800.0]: '
      read (5,*) fullwidth
      print *,'Do you want an output file? [YES/NO]: '
      read (5,'(a)') answer
      if (answer(1:1).eq.'Y'.or.answer(1:1).eq.'y') then
         fflag = .true.
         print *,'Enter output file name: '
         read (5,'(a)') fname
         open (unit=3,file=fname,status='new')
      else
         fflag = .false.
      endif
      print *,'Do you want output on your display? [YES/NO]: '
      read (5,'(a)') answer
      if (answer(1:1).eq.'Y'.or.answer(1:1).eq.'y') then
         dflag = .true.
      else
         dflag = .false.
      endif
c
c convert center coordinates to appropriate units.
c
      icentra = (rah + ram/60. + ras/3600.) * RFACTOR
      centdecd = (dd + dm/60. + decs/3600.) * dsign
      cosdec = cos(centdecd*picon)
      icentdec = (centdecd + 90.) * DFACTOR
c
c calculate boxwidth-related parameters
c
      halfwidth = fullwidth * 50.
      declo = icentdec - halfwidth
      dechi = icentdec + halfwidth
      if (declo.lt.0) declo = 0
      if (dechi.ge.180.0*DFACTOR) dechi = 180.0*DFACTOR - 1
      if (declo.eq.dechi) then
         print *,' Declination out of bounds'
         stop
      endif
      ralo = icentra - int(halfwidth/cosdec)
      rahi = icentra + int(halfwidth/cosdec)
c
c now divide read up into 4 pieces: possible crossover of
c declination zone, and possible crossover of 0h RA
c
      do i=1,NZONES
        if (declo.ge.zlo(i).and.declo.le.zhi(i)) zonelo = i
        if (dechi.ge.zlo(i).and.dechi.le.zhi(i)) zonehi = i
      enddo
      if (zonelo.eq.zonehi) then
        dbeg(1) = declo
        dend(1) = dechi
        dbeg(2) = 0
        dend(2) = 0
        z(1) = zone(zonelo)
        z(2) = zone(zonelo)
      else
        dbeg(1) = declo
        dend(1) = zhi(zonelo)
        dbeg(2) = zlo(zonehi)
        dend(2) = dechi
        z(1) = zone(zonelo)
        z(2) = zone(zonehi)
      endif
      if (ralo.lt.0) then
        rbeg(1) = HR24 + ralo
        rend(1) = HR24-1
        rbeg(2) = 0
        rend(2) = rahi
      elseif (rahi.gt.HR24) then
        rbeg(1) = ralo
        rend(1) = HR24-1
        rbeg(2) = 0
        rend(2) = rahi - HR24
      else
        rbeg(1) = ralo
        rend(1) = rahi
        rbeg(2) = 0
        rend(2) = 0
      endif
c
c now 4 possible read cases
c
      if (dbeg(2).eq.dend(2).and.rbeg(2).eq.rend(2)) then
        do i=3,4
          rbeg(i) = 0
          rend(i) = 0
          dbeg(i) = 0
          dend(i) = 0
          z(i) = 0
        enddo
      else if (dbeg(2).ne.dend(2).and.rbeg(2).eq.rend(2)) then
        rbeg(2) = rbeg(1)
        rend(2) = rend(1)
        do i=3,4
          rbeg(i) = 0
          rend(i) = 0
          dbeg(i) = 0
          dend(i) = 0
          z(i) = 0
        enddo
      else if (dbeg(2).eq.dend(2).and.rbeg(2).ne.rend(2)) then
        dbeg(2) = dbeg(1)
        dend(2) = dend(1)
        do i=3,4
          rbeg(i) = 0
          rend(i) = 0
          dbeg(i) = 0
          dend(i) = 0
          z(i) = 0
        enddo
      else if (dbeg(2).ne.dend(2).and.rbeg(2).ne.rend(2)) then
        do i=1,2
          rbeg(5-i) = rbeg(i)
          rend(5-i) = rend(i)
          dbeg(i+2) = dbeg(i)
          dend(i+2) = dend(i)
          z(i+2) = z(i)
        enddo
      endif
c
c do the 4 partial reads of the catalog
c
      print *,'Starting catalog search, please be patient...'
      print *
      if(fflag) write(3,903) i1,i2,ras,dds,i3,i4,decs,fullwidth
      if(dflag) write(6,903) i1,i2,ras,dds,i3,i4,decs,fullwidth
903      format(10x,' USNO-A Search'/
     $     ' Field Center(J2000):  RA ',
     $     i2.2,1x,i2.2,f7.3,'  Dec ',a1,i2.2,1x,i2.2,f6.2/
     $     ' Boxwidth (arcsec): ',f7.2/
     $     '    Star        RA(J2000)     DEC(J2000)    Bmag',
     $     '   Rmag   Field GSC? Magerr?')
      do i=1,4
        call readcat (rbeg(i),rend(i),dbeg(i),dend(i),z(i))
      enddo
      print 904,nstars
904   format (i5,' stars were found')
      if (fflag.and.nstars.ne.0) close(3)
      if (fflag.and.nstars.eq.0) close(3,status='delete')
      print *,'Do you want another search? [YES/NO]: '
      read (5,'(a)') answer
      if (answer(1:1).eq.'Y'.or.answer(1:1).eq.'y') goto 700
      stop
      end

      SUBROUTINE READCAT (rbeg,rend,dbeg,dend,zone)
c
c actual zone file reader.
c
      INTEGER*4
     $   NERA,NMAX,CAT_REC_SIZE,NZONES
      REAL*8
     $   RFACTOR,DFACTOR
      PARAMETER
     $   (NERA = 96,
     $   NZONES = 24,
     $   NMAX = 10*1024*1024,
     $   RFACTOR = 15.0D00*3600.0D00*100.0D00,
     $   DFACTOR =         3600.0D00*100.0D00,
     $   CAT_REC_SIZE = 12)  ! 3 longwords, or 12 bytes/record
      INTEGER*4
     $   rbeg,rend,dbeg,dend,zone,ira,idec,ibuf,
     $   rah,ram,dd,dm,nstars,diskno(NZONES)
     $   i, n, estart(NERA), elength(NERA),
     $   high, med, low, q, ifield
      DOUBLE PRECISION
     $  era(NERA+1), r, d
      REAL*4
     $   ras,ds,bmag,rmag
      CHARACTER
     $   zname*12, catname*12, path*80, dsign*1,
     $   cgsc*1, cmagerr*1, flag*1, path2*80, dup*1
      LOGICAL
     $   fflag,dflag,gscflg,bothmag,gsconly,magerr,bmaghere,
     $   rmaghere
      COMMON /blka/ nstars,path,path2,fflag,dflag
      DATA diskno/1,1,6,5,3,2,1,4,6,5,7,10,
     $     8,7,8,9,9,4,10,3,2,6,2,3/
c
c first, test for the 'unread' condition
c
      IF (rbeg.eq.rend.or.dbeg.eq.dend) RETURN
c
c inform user about disk swap
c
      i = zone/75 + 1
      write (6,9944) diskno(i)
9944  format (' Insert disk ',i2,
     $ ', then press <cr>:')
      read (5,'(a)') flag
c
c read the accelerator file
c this is a kludge, since the file was written on a Unix
c machine with LF between lines, but not CR/LF as DOS prefers.
c
      write (zname,9004) zone
9004  format ('ZONE',i4.4,'.ACC')
        OPEN (
c     *        access='sequential',
     *      access='transparent',
c     *        carriagecontrol='list',
     *        err=210,
c     *        form='formatted',
     +        form='unformatted',
     *        recl=1,
     *        file=path(1:nblank(path))//zname,
     *        status='old',
     *        unit=1
     *       )
c read file length
        n = 0
        DO i=1,NERA
c          READ (1,9002,rec=i) era(i),estart(i),elength(i)
c          READ (1,9002) era(i),estart(i),elength(i)
c9002      FORMAT (f5.2, 2i12)
          CALL accread(i,era(i),estart(i),elength(i))
          n = n + elength(i)
        ENDDO
        CLOSE (1)
c tests on file length.  if zero, no stars, return
        IF (n.le.0) RETURN
        indx = 1
        era(NERA+1) = 24.0
        ralo = rbeg / RFACTOR
        DO i=1,NERA
          if (ralo.ge.era(i).and.ralo.lt.era(i+1)) then
            indx=i
            goto 20
          endif
        ENDDO
20    continue
c
c now do actual search.
c
          write (catname,9005) zone
9005      format('ZONE',i4.4,'.CAT')
          OPEN (
     *          access='transparent',
c     *          carriagecontrol='none',
     *          form='unformatted',
     *          file=path(1:nblank(path))//catname,
c     *          readonly,
     *          err=190,
     *          recl=1,
c     *          recordtype='fixed',
c     *          shared,
     *          status='old',
     *          unit=1
     *         )
c
c the file read is system dependent.  The catalog was written on
c a Silicon Graphics workstation, which uses most-significant-byte-first
c byte ordering.  For IBM PCs, you will have to do a byte swap here.
c Some Fortran compilers will do the byte ordering automatically
c for you if you specify the input file type in the open statement.
c
      do 150 i=estart(indx),n
        call bufread(i,ira,idec,ibuf)
        if (ira.gt.rend) goto 200
        if (ira.lt.rbeg.or.idec.lt.dbeg.or.idec.gt.dend) goto 150
           nstars = nstars+1
           r = ira/RFACTOR
           rah = r
           r = (r - float(rah)) * 60.
           ram = r
           ras = (r - float(ram)) * 60.
           d = idec/DFACTOR - 90.0D00
           dsign = '+'
           if (d.lt.0.0) then
             dsign = '-'
             d = abs(d)
           endif
           dd = d
           d = (d - float(dd)) * 60.
           dm = d
           ds = (d - float(dm)) * 60.
           if (ibuf.lt.0) then
             dup = 'D'
             gscflg = .true.
           else
             dup = ' '
             gscflg = .false.
           endif
c
           q = iabs(ibuf)
           if (q/1000000000 .ne. 0) then
             magerr = .true.
           else
             magerr = .false.
           endif
c
           high = q/1000000000
           q = q-high*1000000000
           high = q/1000000
           q = q-high*1000000
           med = q/1000
           low = q-med*1000
           rmag = 0.1*low
           bmag = 0.1*med
           bmaghere = .true.
           rmaghere = .true.
           if (low.eq.999) rmaghere = .false.
           ifield = high
           if (ifield.eq.0) then
             gsconly = .true.
             rmaghere = .true.
             bmaghere = .false.
           else
             gsconly = .false.
           endif
           if (bmaghere.and.rmaghere) then
             bothmag = .true.
           else
             bothmag = .false.
           endif
c
           if (dflag) then
             if (gscflg.or.gsconly) then
               cgsc = 'x'
             else
               cgsc = ' '
             endif
             if (magerr) then
               cmagerr = 'x'
             else
               cmagerr = ' '
             endif
             if (gsconly) then
               write (6,900) zone,i,rah,ram,ras,dsign,dd,dm,ds,
     $           rmag,cgsc,cmagerr
             else if (bothmag) then
               write (6,901) zone,i,rah,ram,ras,dsign,dd,dm,ds,
     $           bmag,rmag,ifield,cgsc,cmagerr
             else
               write (6,902) zone,i,rah,ram,ras,dsign,dd,dm,ds,
     $           bmag,ifield,cgsc,cmagerr
             endif
           endif
900       format (1x,i4.4,'-',i8.8,2x,i2.2,1x,i2.2,f7.3,2x,a1,i2.2,1x,
     $        i2.2,1x,f5.2,f7.2,15x,a1,6x,a1)
901       format (1x,i4.4,'-',i8.8,2x,i2.2,1x,i2.2,f7.3,2x,a1,i2.2,1x,
     $        i2.2,1x,f5.2,f7.2,f7.2,3x,i4.4,1x,a1,6x,a1)
902       format (1x,i4.4,'-',i8.8,2x,i2.2,1x,i2.2,f7.3,2x,a1,i2.2,1x,
     $        i2.2,1x,f5.2,f7.2,10x,i4.4,1x,a1,6x,a1)
           if (fflag) then
             if (gscflg.or.gsconly) then
               cgsc = 'x'
             else
               cgsc = ' '
             endif
             if (magerr) then
               cmagerr = 'x'
             else
               cmagerr = ' '
             endif
             if (gsconly) then
               write (3,900) zone,i,rah,ram,ras,dsign,dd,dm,ds,
     $           rmag,cgsc,cmagerr
             else if (bothmag) then
               write (3,901) zone,i,rah,ram,ras,dsign,dd,dm,ds,
     $           bmag,rmag,ifield,cgsc,cmagerr
             else
               write (3,902) zone,i,rah,ram,ras,dsign,dd,dm,ds,
     $           bmag,ifield,cgsc,cmagerr
             endif
           endif
150     continue
200   continue
      close(1)
      return
190   continue
      print *,'.cat file not found:',path(1:nblank(path))
     $   //catname
      stop
210   continue
      print *,'.acc file not found:',path(1:nblank(path))
     $   //zname
      stop
      END

      SUBROUTINE ACCREAD (indx,era,estart,elength)
c
c unix file reader
c not transportable!
c
      DOUBLE PRECISION era
      INTEGER*4 estart,elength
      INTEGER  indx,i,irec,j
      CHARACTER buf*29
c
      irec = 30*(indx-1) + 1
      do i=irec,irec+28
        j = i-irec+1
        read (1,rec=i) buf(j:j)
      enddo
      read (buf,900) era,estart,elength
900   format(f5.2,2i12)
      RETURN
      END

      SUBROUTINE BUFREAD (indx,ira,idec,ibuf)
c
c IBM PC direct file reader
c NOT TRANSPORTABLE!!!!
c skips over the 12-byte direct file header
c reads record, does byte swap
c
      INTEGER*4
     $    indx,ira,idec,ibuf,i,irec,d1,d2,d3
      CHARACTER*1 buf(12)
      DATA d1,d2,d3/256,65536,16777216/
c
      irec = 12*(indx-1)+1
      do i=irec,irec+11
        read(1,rec=i) buf(i-irec+1)
      enddo
      ira = ichar(buf(4))+ichar(buf(3))*d1+
     $      ichar(buf(2))*d2+ichar(buf(1))*d3
      idec = ichar(buf(8))+ichar(buf(7))*d1+
     $       ichar(buf(6))*d2+ichar(buf(5))*d3
      ibuf = ichar(buf(12))+ichar(buf(11))*d1+
     $       ichar(buf(10))*d2+ichar(buf(9))*d3
      return
      end

      FUNCTION NBLANK (string)
c
c find first blank character
c
      INTEGER
     *       nblank, i
      CHARACTER*1
     *       string(*)
c
      nblank = 0
      DO i=1,256
        nblank = nblank + 1
        IF (string(nblank).eq.' ') GOTO 100
      ENDDO
100   continue
      nblank = nblank - 1
      RETURN
      END