PROGRAM u2access
C
C  pgf77 -o u2access u2access.f u2sub.f sorti.f
C
C  030417 NZ create
C  030418 NZ finish, test,debug
C  030422 NZ extend interactive options
C  030423 NZ star implement pos., err. update by newep
C  030424 NZ finish pos.update
C  030427 NZ back to 23 items, ep0=1975, zmax, open_zfile
C  030429 NZ bfx flag
C  030428 NZ no "READONLY" in OPEN statements 
C  030529 NZ width of box RA, Dec
C  030530 NZ wrap around 0/24h also for box, fix subr.
C  030601-2 NZ finalize output formats, tests, index DA
C
C  - top level program to access UCAC2 data
C  - extract all items for stars
C      in selected areas (RA, Dec box)
C  - formatted (ASCII) output table, various format options
C
C  use Fortran unit 13 to read index and zone files
C  use Fortran unit 20 for output file (or screen dump)

      IMPLICIT NONE
      INTEGER    zmax, dimj, dima
      PARAMETER (zmax= 288       ! maximal zone number
     .          ,dimj= 25 + 4    ! add updated position and errors
     .          ,dima=  30000)   ! max.numb. stars in box

      INTEGER  alls(dima,dimj)   ! all data all selected star
      INTEGER  nx (zmax,240)     ! index = accumul.numb.stars f(zone,RA bin)
      INTEGER  fmt

      INTEGER  i,j,jp,is,zn, nsr,nss,errc, uo, rah,ram,dcd,dcm, un
      REAL*8   ra1,ra2, dc1,dc2, ra0,dc0,ras,dcs, mag1,mag2, newep
     .        ,wra,wdc
      CHARACTER*40 pathu,fnxu, fnout, answer
      CHARACTER*1  a1
      LOGICAL      bf, bfx, was

* default
CC    pathu = '/d1/rl/u2g/'      ! path to orig. UCAC2 zone files
      pathu = '/mnt/cdrom/u2/'   ! CDROM at my Linux box
CC    pathu = 'D:\u2\'           ! sample path for Windows, CD-ROM drive

      fnout = 'u2.tab' 
      un    = 13                 ! Fortran unit number for zone files
      ra1   =  0.0d0
      ra2   = 24.0d0
      dc1   = -90.0d0
      dc2   = -89.9d0
      mag1  =  5.0d0
      mag2  = 18.0d0
      newep = 2000.0d0
      wra   = 2.0d0
      wdc   = 0.1d0
      fmt   = 1
      was   = .FALSE.    ! width in arcsec, else degree

* interactive: basic data, common for entire run of program

      WRITE (*,'(/a)') 'Program to access UCAC2 release data'
      WRITE (*,'( a)') '  hit "enter" for defaults'
      WRITE (*,'( a)') 'option: loop through several areas, but'
      WRITE (*,'( a)') '  output with single format to single file'
      WRITE (*,'( a)') '  per run of this program'
      WRITE (*,'( a)') 'wrap around 24/0h RA allowed'

      WRITE (*,'(/a)') 'path for original UCAC2 files = '
      WRITE (*,'(a,$)') pathu
      READ (*,'(a)') answer
      IF (answer.NE.' ') pathu = answer

      WRITE (*,'(/a)') '  fmt = 1 = all original integer'
      WRITE (*,'( a)') '  fmt = 2 = updated RA, Dec (degree) PM..' 
      WRITE (*,'( a)') '  fmt = 3 = updated RA, Dec (hms) photom..'
      WRITE (*,'( a)') '  fmt = 4 = updated RA, Dec (h,d) short'
      WRITE (*,' (a)') '  fmt = ... user defined in subr. put_stars'
      WRITE (*,'(a,i2,a,$)') 'format for output = ',fmt,'  ' 
      READ (*,'(a)') answer
      IF (answer.NE.' ') READ (answer,*) fmt
      
      WRITE (*,'(/a)') 'file name output table or "s" = screen dump '
      WRITE (*,'(a,$)') fnout
      READ (*,'(a)') answer
      IF (answer.NE.' ') fnout = answer

      WRITE (*,'(/a)')
     . 'select magnitude range (UCAC mags, 5 - 18 = all):'
      WRITE (*,'(2f7.2,a,$)') mag1, mag2,'  '
      READ (*,'(a)') answer
      IF (answer.NE.' ') READ (answer,*) mag1,mag2

      IF (fmt.GT.1) THEN
        WRITE (*,'(/a)') 
     .   'update positions for proper motions to new epoch'
        WRITE (*,'( a)') '2000.0 is default = original UCAC2'
        WRITE (*,'(f7.2,a,$)') newep,'  '
        READ (*,'(a)') answer
        IF (answer.NE.' ') READ (answer,*) newep
      ENDIF

      WRITE (*,'(/a)') 'see readme2.txt for all items'
      WRITE (*,'( a)') '  0 = no sort = default'
      WRITE (*,'( a)') '  1 = RA ' 
      WRITE (*,'( a)') '  2 = Dec' 
      WRITE (*,'( a)') '  3 = UCAC mag' 
      WRITE (*,'( a)') ' ...'
      WRITE (*,'( a)') ' 24 = UCAC ID number'
      WRITE (*,'(a,$)') 'sort by item number: '
      READ (*,'(a)') answer
      IF (answer.NE.' ') READ (answer,*) is

* check for byte flip
      CALL chk_byte_flip (pathu,un,bf)

* read index file
      jp = INDEX (pathu,' ') - 1

CC    fnxu = pathu(1:jp) // 'u2index.unf'
CC    OPEN (un,FILE=fnxu,FORM='unformatted',READONLY)
CC    READ (un,ERR=901) nx

* new read of index file, works better accross platforms:
      fnxu = pathu(1:jp) // 'u2index.da' 
      WRITE (*,'(a,a)') 'open index file = ',fnxu
      OPEN (un,FILE=fnxu,ACCESS='direct',RECL=960) ! 240 * 4 byte
      DO zn=1,288
        READ (un,REC=zn,ERR=901) (nx(zn,j),j=1,240)
      ENDDO

      CLOSE(un)
      WRITE (*,'(a)') 'index read successfully'
      DO zn=1,3
        WRITE (*,'(a,6i6)') 'zn, nx...=',zn,(nx(zn,j),j=1,5)
      ENDDO

      CALL nx_byte_flip (nx,zmax,bfx) ! check for byte flip, apply if needed
 
* prepare table output
      IF (fnout.EQ.'s') THEN
        uo = 6
      ELSE
        uo = 20
        OPEN (uo,FILE=fnout)
      ENDIF

* loop boxes 
 101  WRITE (*,'(/a,$)') '== new box: r=range, c=center, q=quit '
      READ (*,'(a)') a1

      IF (a1.EQ.'r'.OR.a1.EQ.'R') THEN
        WRITE (*,'(a,2f8.4,a,$)') 'RA (hour) = ',ra1,ra2,'  '
        READ (*,'(a)') answer
        IF (answer.NE.' ') READ (answer,*) ra1,ra2
        WRITE (*,'(a,2f8.3,a,$)') 'Dec (deg) = ',dc1,dc2,'  '
        READ (*,'(a)') answer
        IF (answer.NE.' ') READ (answer,*) dc1,dc2

      ELSEIF (a1.EQ.'c'.OR.a1.EQ.'C') THEN
        WRITE (*,'(2a)') 'width in degree or enter "a" for arcsec'
     .    ,' after values'
        WRITE (*,'(a,2f7.3,a,$)')
     .     'width of box RA*cosD, Dec =',wra,wdc,'  '
        READ (*,'(a)') answer
        j = INDEX(answer,'a')
        IF (j.GT.0) THEN
          was = .TRUE.
          answer (j:j) = ' '
        ELSE
          was = .FALSE.
        ENDIF

        IF (answer.NE.' ') THEN
          READ (answer,*) wra,wdc
          IF (was) THEN
            wra = wra / 3.6d3    ! arcsec to degree
            wdc = wdc / 3.6d3
          ENDIF
        ENDIF
        IF (wra.LT.0.0d0.OR.wdc.LT.0.0d0) GOTO 101 

        WRITE (*,'(a,$)') 'center RA (hms), Dec (dms) = '
        READ (*,'(a)') answer

        CALL calc_box (answer,wra,wdc,ra1,ra2,dc1,dc2)

        WRITE (*,'(2f10.6,2f10.5)') ra1,ra2,dc1,dc2

      ELSE
        GOTO 99
      ENDIF

* output info
      WRITE (uo,'(/a, f10.3)') 'epoch      = ',newep
      WRITE (uo,'( a, i10  )') 'sort item  = ',is 
      WRITE (uo,'( a, i10  )') 'output fmt = ',fmt
      WRITE (uo,'( a,2f10.2)') 'mag range  = ',mag1,mag2
      WRITE (uo,'( a,f10.6,f10.5)') 'min RA, DE = ',ra1,dc1 
      WRITE (uo,'( a,f10.6,f10.5)') 'max RA, DE = ',ra2,dc2 

* get stars 
      CALL get_stars (bf,ra1,ra2,dc1,dc2,mag1,mag2
     .               ,newep,pathu,nx,zmax,dima,dimj
     .               ,alls,nss,nsr,errc)

      WRITE (*,'(a,3i7)') 
     .   'stars read, slected, errors = ',nsr,nss,errc

* sort
      IF (is.GE.1.AND.is.LE.dimj) THEN
        WRITE (*,'(a,i3)') 'start sort by column = ',is
        CALL sorti (alls,is,dimj,1,nss,dima,dimj)
      ENDIF

* output
      CALL put_stars (alls,dima,dimj,nss,fmt,uo)
      WRITE (*,'(a,i6,a)') 'output complete',nss,' stars'

      GOTO 101         ! more boxes

* the end
 901  WRITE (*,'(a)') 'error in reading index file'
      STOP

  99  WRITE (*,'(a,/)') 'thanks for using <u2access>, have a nice day!'
 
      END  ! main <u2access>

************************************************************************

      SUBROUTINE calc_box (answer,wra,wdc,ra1,ra2,dc1,dc2)
C
C  input : answer = string with RA, Dec in hms format
C          wra    = width of box (degree) along RA (*cosDec)
C          wdc    = width of box (degree) along Dec
C  output: ra1,ra2= range in RA (hour)
C          dc1,dc2= range in Declination (degree)

      IMPLICIT NONE
      CHARACTER*(*) answer
      REAL*8  wra,wdc, ra1,ra2,dc1,dc2
      REAL*8  degrad, cosdec, hwr,hwd, ras,dcs, ra0,dc0
      INTEGER j, rah,ram,dcd,dcm
      LOGICAL negdec

      degrad = DATAN (1.0d0) / (4.5d1)    ! degree to radian
      j = INDEX (answer,'-')

      IF (j.GT.0) THEN
        negdec = .TRUE.
        answer (j:j) = ' '
      ELSE
        negdec = .FALSE.
      ENDIF
      READ (answer,*) rah,ram,ras, dcd,dcm,dcs

      ra0 = rah + ram/60.0d0 + ras/3.6d3
      dc0 = dcd + dcm/60.0d0 + dcs/3.6d3
      IF (negdec) dc0 = -dc0
      cosdec = DCOS (dc0 * degrad)

      hwr = wra / (cosdec * 30.0d0)  ! half width in hour for RA
      hwd = wdc / 2.0d0              ! half width in Dec

      ra1 = ra0 - hwr
      ra2 = ra0 + hwr
      IF (ra1.LT. 0.0d0) ra1 = ra1 + 24.0d0
      IF (ra2.GT.24.0d0) ra2 = ra2 - 24.0d0

      dc1 = dc0 - hwd
      dc2 = dc0 + hwd

      END  ! subr. <calc_box>

************************************************************************

      SUBROUTINE get_stars (bf,ra1,ra2,dc1,dc2,mag1,mag2
     .                     ,newep,pathu,nx,zmax,dima,dimj
     .                     ,alls,nss,nsr,errc)
C
C  get all items from all stars within specified RA,Dec box
C
C  input : bf       = .TRUE. if flip of bytes is required
C          ra1,ra2  = RA range (hour, decimal)
C          dc1,dc2  = declination range (degree, decimal)
C          mag1,mag2= magnitude rane (UCAC 1/100 mag)
C          newep    = requested epoch for positions (year)
C          pathu    = path name for UCAC2 zone files (max. 40 char.)
C          nx       = index array = accumulated numb.stars f(zone, RA bin)
C          zmax     = largest zone number (0.5 deg steps from South Pole)
C          dima     = max. number of stars to retrieve
C          dimj     = number of items per star incl. updated pos.+error
C
C  output: alls     = all items for all retrieved stars
C          nss      = number of stars selected = within specified box
C          nsr      = number of star records read
C          errc     = number of errors while reading zone files

      IMPLICIT NONE
      LOGICAL  bf
      REAL*8   ra1,ra2,dc1,dc2, mag1,mag2, newep
      CHARACTER*(*) pathu
      INTEGER  zmax,dima,dimj, nb, nsr,nss,errc
      INTEGER  nx (zmax,240)     ! accumul.numb.stars = f(zone,RA bin)
      INTEGER  alls(dima,dimj)   ! all integer data each star

      INTEGER  idat(25)          ! data items individual star

      INTEGER  i,j,ju,zn,rr,uz, d1m,d2m,z1,z2,nz, mi1,mi2
     .        ,r1m(2),r2m(2),i1(2),i2(2),nr, j1,j2, recn,n0
     .        ,ran,dcn, sxn,syn
 
      LOGICAL  errflg, only_rd

* prepare
      only_rd = .TRUE.           ! no write to zone files
      mi1 = IDNINT (mag1 * 1.0d2)
      mi2 = IDNINT (mag2 * 1.0d2)

* calculate range 
      CALL get_zone_range (dc1,dc2,zmax,d1m,d2m,z1,z2,nz)

      IF (nz.LT.1) THEN
        nss = 0
        RETURN                   ! exit, declination invalid
      ENDIF

      CALL get_ra_range (ra1,ra2,r1m,r2m,i1,i2,nr)

* initialize
      uz = 13    ! unit number for zone files
      nsr = 0
      nss = 0 
      errc= 0
      ju = INDEX (pathu,' ') - 1

* loop all zones
      DO zn = z1,z2
        WRITE (*,'(a,i3)') 'open file for zone = ',zn 
        CALL open_zfile (pathu,uz,zn,only_rd)

        IF (zn.EQ.1) THEN
          n0 = 0              ! no stars in "previous" zone
        ELSE
          n0 = nx (zn-1,240)  ! = numb.of stars up to end of previous zone
        ENDIF

* rr    = 1 or 2     = numb. of RA ranges, normal=1, else crossover 24/0 hour
* i1,i2 = 1,...,240  = bin number (0.1 hour) along RA
* j1,j2 = 1,... many = record number on zone file (1 record = 1 star)

        DO rr = 1,nr          ! 1 or 2 RA ranges possible
          IF (i1(rr).EQ.1) THEN
            j1 = 1            ! start with first star on file
          ELSE
            j1 = nx (zn,i1(rr)-1) - n0 + 1   ! i1-1 = previous bin
          ENDIF

          j2 = nx (zn,i2(rr)) - n0 

CC        WRITE (*,'(a,i3,i9,3i5,2i9)') 
CC   .     'zn, n0, rr,i1,i2, j1,j2 = ',zn,n0,rr,i1(rr),i2(rr),j1,j2

          DO recn= j1,j2
            CALL read_u2line (uz,recn,bf,idat,errflg) 

            IF (errflg) THEN
              errc = errc + 1  ! count read errors

            ELSE
              nsr = nsr + 1    ! count number of stars read
              IF (idat(3).GE.mi1.AND.idat(3).LE.mi2) THEN   ! mag range
              IF (idat(1).GE.r1m(rr).AND.idat(1).LE.r2m(rr).AND.
     .            idat(2).GE.d1m.AND.idat(2).LE.d2m) THEN
                nss = nss + 1
                IF (nss.GT.dima) THEN
                  WRITE (*,'(/a,i8)') 
     .             'WARNING: too many stars, take first dima =',dima
                  RETURN
                ENDIF

                DO j=1,25     ! current version: only 23 items on file
                  alls (nss,j) = idat(j)    ! still keep data structure
                ENDDO
* here generate on the fly, instead read from file:
                alls (nss,24) = n0 + recn   ! official UCAC2 star ID

                CALL apply_pm (idat,newep, ran,dcn,sxn,syn)
                alls (nss,26) = ran
                alls (nss,27) = dcn
                alls (nss,28) = sxn
                alls (nss,29) = syn

              ENDIF    ! within box
              ENDIF    ! within mag range
            ENDIF      ! case read o.k. or error

          ENDDO        ! all stars within a zone
        ENDDO          ! 1 or 2 RA ranges 
      ENDDO            ! all zones

      END   ! subr. <get_stars>

************************************************************************

      SUBROUTINE put_stars (alls,dima,dimj,nss,fmt,uo)
C
C  input : alls  = array with selected stars
C          dima  = max.numb. of stars
C          dimj  = number of items per star
C          nss   = number of selected stars
C          fmt   = format modus for output (= 1 to 4 here)
C          uo    = Fortran unit number for output file
C  output : items to file with unit number uo (can be screen)

      IMPLICIT NONE
      INTEGER  dima,dimj,nss, fmt,uo
      INTEGER  alls(dima,dimj)

      INTEGER  i,j
      REAL*8   mag, epr,epd, pmr,pmrc,pmd, degrad, cosdec
     .        ,ran,dcn, magj, magh, magk, epmr,epmd, qra,qde
      CHARACTER*13 cra,cdc

      degrad = DATAN (1.0d0) / (4.5d1)   ! degree to radian

* header line
      IF (fmt.EQ.1) THEN
        WRITE (uo,'(/4a)') '        RA         DE'
     .   ,' Umag eRc eDc no eUp nc cf epocR epocD'
     .   ,'    pmRA  pmDE epR epD qRA qDE   2MASS ID'
     .   ,'  2m_J  2m_H 2m_Ks phf ccf UCAC2 ID'
      ELSEIF (fmt.EQ.2) THEN
        WRITE (uo,'(/3a)') '   RA (deg)    DE (deg)'
     .   ,' U2mag eRA eDE     epRA     epDC    pm RA'
     .   ,' pmRAcD   pmDE epmR epmD UCAC2 ID'
      ELSEIF (fmt.EQ.3) THEN
        WRITE (uo,'(/3a)') '           RA            DE'
     .   ,' eRA eDE no Upe nc   qRA   qDE U2mag'
     .   ,'   Jmag   Hmag   Kmag phf ccf UCAC2 ID'
      ELSEIF (fmt.EQ.4) THEN
        WRITE (uo,'(/2a)') '  RA (hour)    DE (deg)'
     .   ,' U2mag eRA eDE UCAC2 ID'
      ELSE
        WRITE (*,'(a,i3)') '** invalid format  fmt = ',fmt
        RETURN
      ENDIF

* table with data
      DO i=1,nss
        IF (fmt.EQ.1) THEN               ! all original integer 
          WRITE (uo,'(i10,i11,i5,2i4,i3,i4,2i3,2i6,i8,i6
     .               ,2i4,2i4,i11,3i6,2i4.3,i9.8)')
     .      (alls(i,j),j=1,24)

        ELSEIF (fmt.EQ.2) THEN
          ran = alls(i,26) / 3.6d6       ! updated mas -> degree
          dcn = alls(i,27) / 3.6d6       ! updated mas -> degree
          mag = alls(i, 3) * 1.0d-2
          epr = alls(i,10) * 1.0d-3 + 1975.0d0
          epd = alls(i,11) * 1.0d-3 + 1975.0d0
          cosdec = DCOS (dcn * degrad)
          pmr = alls(i,12) * 0.1d0       ! 0.1 mas/yr -> mas/yr
          pmrc= pmr * cosdec
          pmd = alls(i,13) * 0.1d0
          epmr= alls(i,14) * 0.1d0
          epmd= alls(i,15) * 0.1d0

          WRITE (uo,'(f11.7,f12.7,f6.2,2i4,2f9.3
     .               ,f9.1,2f7.1,2f5.1,i9.8)') 
     .      ran,dcn,mag, alls(i,28),alls(i,29), epr,epd
     .     ,pmr,pmrc,pmd, epmr,epmd, alls(i,24)

        ELSEIF (fmt.EQ.3) THEN
          ran = alls(i,26) * 1.0d-3   ! updated pos. mas -> arcsec 
          dcn = alls(i,27) * 1.0d-3   ! updated pos. mas -> arcsec
          mag = alls(i, 3) * 1.0d-2
          magj= alls(i,19) * 1.0d-3   ! from 2MASS
          magh= alls(i,20) * 1.0d-3
          magk= alls(i,21) * 1.0d-3
          qra = alls(i,16) / 20.0d0   ! quality PM RA, scatter/model
          qde = alls(i,17) / 20.0d0

          CALL as2hms (ran,dcn,cra,cdc) 

          WRITE (uo,'(a,1x,a,2i4,i3,i4,i3,3f6.2,3f7.3,2i4.3,i9.8)')
     .     cra,cdc, alls(i,28),alls(i,29),alls(i,6),alls(i,7),alls(i,8)
     .    ,qra,qde,mag,magj,magh,magk, alls(i,22),alls(i,23),alls(i,24)

        ELSEIF (fmt.EQ.4) THEN 
          ran = alls(i,26) / 5.4d7       ! updated mas -> hour
          dcn = alls(i,27) / 3.6d6       ! updated mas -> degree
          mag = alls(i, 3) * 1.0d-2
          WRITE (uo,'(f11.8,f12.7,f6.2,2i4,i9.8)')
     .      ran,dcn,mag, alls(i,28),alls(i,29), alls(i,24)
        ENDIF

      ENDDO   ! all stars
    
      END   ! subr. <put_stars>

************************************************************************

      SUBROUTINE pos_up (ra,dc,pmra,pmdc,newep, ran,dcn)
C
C  input : ra, dc    = RA, Dec (mas) at J2000 epoch
C          pmra,pmdc = proper motion RA, Dec (0.1 mas/yr)
C                      no cos(Dec) for pmra, i.e. large values near pole
C          newep     = new epoch (year)
C  output: ran,dcn   = RA, Dec (mas) at new epoch

      IMPLICIT NONE
      INTEGER  ra,dc,pmra,pmdc,ran,dcn
      REAL*8   newep, dt, dra, ddc

      dt  = newep - 2000.0d0
      dra = DFLOAT(pmra) * 0.1d0 * dt
      ddc = DFLOAT(pmdc) * 0.1d0 * dt
 
      ran = ra + IDNINT (dra)
      IF (ran.GT.1296000000) ran = ran - 1296000000   ! 24 hour in mas
      IF (ran.LT.         0) ran = ran + 1296000000   ! restrict to 0..24 range

      dcn = dc + IDNINT (ddc)     ! assume no jump over pole

CC    WRITE (20,'(a,f9.3,2f9.1)') 'pos_up: dt,dra,ddc = ',dt,dra,ddc

      END    ! subr. <pos_up>

************************************************************************

      SUBROUTINE pos_error (sigx,sigy,cepx,cepy,spmx,spmy,newep
     .                     ,sxn,syn)
C
C     x = RA * cosDec,  y = Dec  component
C
C input : sigx, sigy = position error (mas) at central epoch
C         cepx, cepy = central epoch (1/1000 yr from 1975)
C         spmx, spmy = error of proper motion component (0.1 mas/yr)
C         newep      = requested epoch for position error (year)
C output: sxn, syn   = position error at new epoch (mas)

      IMPLICIT NONE
      INTEGER  sigx,sigy, spmx,spmy, cepx,cepy, sxn,syn
      REAL*8   newep, dtx,dty

      dtx = newep - DFLOAT(cepx) * 1.0d-3 - 1975.0d0
      dty = newep - DFLOAT(cepy) * 1.0d-3 - 1975.0d0

      sxn = IDNINT (DSQRT (sigx**2 + (spmx * 0.1d0 * dtx)**2))
      syn = IDNINT (DSQRT (sigy**2 + (spmy * 0.1d0 * dty)**2))

CC    WRITE (20,'(a,2f9.3)') 'pos_error: dtx,dty = ',dtx,dty

      END   ! subr. <pos_error>

************************************************************************

      SUBROUTINE apply_pm (idat,newep, ran,dcn,sxn,syn)
C
C input : idat = integer vector with original UCAC2 data for 1 star
C         newep= requested new epoch (year)
C output: ran, dcn = new position, updated for proper motion (mas)
C         sxn, syn = error of new position (at new epoch) (mas)

      IMPLICIT NONE
      INTEGER  idat(25), ran,dcn,sxn,syn
      REAL*8   newep

      CALL pos_up (idat(1),idat(2),idat(12),idat(13),newep, ran,dcn)

      CALL pos_error (idat( 4),idat( 5), idat(10),idat(11)
     .               ,idat(14),idat(15), newep, sxn,syn)

      END   ! subr. <apply_pm>

************************************************************************