PROGRAM u2rdplot
C
C pgf77 -o u2rdplot u2rdplot.f u2sub.f -L /usr/local/pgplot -lpgplot /usr/X11R6/lib/libX11.so
C
C    - read controll data from u2rdplot.in file, no interactive
C    - read RA and Dec from UCAC2 release zone files 
C        (assume 1 RA range, no cross 24/0 hours)
C    - use 3 colors for 3 different magnitude bins
C    - plot RA (hours) and DC (degree) distribution on sky, many stars
C    - options for different projections (pole, box) -- future: tang.proj.
C    - loop several plots, automatic rename of pgplot.ps files
C
C  030226 NZ taken from ../agk2/radc_plot.f 
C  030427 NZ debug/test
C  030429 NZ automatic title, pgdev 
C  030528 NZ don't use READONLY in OPEN statements
C  030531 NZ fix 24/0 RA subr.call

      IMPLICIT NONE

      INTEGER    dim1, dim2, dim3, zmax
      PARAMETER (dim1=  10000    ! max. number selected stars/zone,color
     .          ,dim2=  50000
     .          ,dim3= 300000 
     .          ,zmax=    288 )  ! largest zone number UCAC2

*  variables for PGPLOT:
      REAL*4     xp1(dim1), yp1(dim1)   ! plot points per zone, color
     .          ,xp2(dim2), yp2(dim2)
     .          ,xp3(dim3), yp3(dim3)

      INTEGER    lwt,lws    ! line width = 1...6
      INTEGER    coltxt     ! color number text, coordinates
     .          ,col1,col2,col3  ! colors for stars (3 mag bins)
      INTEGER    symbol     ! symbol number
      INTEGER    just       ! =1 --> scale equal on x,y
      INTEGER    axis       ! see expl. of PGENV

      REAL*4     st, ss     ! character size
      REAL*4     xmin,xmax, ymin,ymax      ! range of x,y values (box)
      REAL*4     xt, yt
      CHARACTER*40  xtit0,ytit0, xtit1,ytit1, title*50
      CHARACTER*8   pgdev
      INTEGER    PGBEG      ! function

* variables for UCAC2
      INTEGER  nx(zmax,240), idat(25)
      INTEGER  i1(2),i2(2), r1m(2),r2m(2)
      INTEGER  un, i0, j1,j2, zn,z1,z2, recn
      INTEGER  d1m,d2m,nz, nr
      LOGICAL  bf, errflg, only_rd

* other variables
      REAL*8     ra1,ra2,dc1,dc2
      REAL*4     masrad, spd
      INTEGER    mag1,mag2,mag3,mag4, ptyp, nplot 
      INTEGER    i,j, jb, nrd, n0,n1,n2,n3,n4, n0t,n1t,n2t,n3t,n4t 
      CHARACTER*40  fnin, fnindex, pathu, fnu2, fnbase, fnplot
      CHARACTER*60  cmd, a1*1
      INTEGER  ierr, SYSTEM

* read control data
      WRITE (*,'(/a/)') 'plot distribution of stars on sky'
      fnin   = 'u2rdplot.in'
CC    OPEN (11,FILE=fnin, READONLY)
      OPEN (11,FILE=fnin)

      READ (11,*) lwt, lws 
      READ (11,*) st,  ss
      READ (11,*) coltxt
      READ (11,*) symbol
      READ (11,*) mag1,mag2,mag3,mag4
      READ (11,*) col1,col2,col3 
      READ (11,*) fnindex 
      READ (11,*) pathu  
      READ (11,*) fnbase 

* prepare
      masrad= ATAN (1.0) / (45.0 * 3.6e6)   ! mas to radian
      jb = INDEX (fnbase,' ') - 1

      WRITE (title,'(a,4i5,a,3i2)')
     .  'UCAC2, mag=',mag1,mag2,mag3,mag4,'  col=',col1,col2,col3

      pgdev= '/cps'  ! color PS plot 
      axis = 0       ! 0 = normal, else logarithmic ...
      nplot= 0       ! count plots
      xtit0= 'center = South Pole'
      ytit0= 'unit = degree'
      xtit1= 'RA (hour)'
      ytit1= 'Dec (degree)'

      un = 12           ! unit for UCAC2 files
      bf      = .FALSE. ! byte flip
      only_rd = .TRUE.  ! no write to zone files here

* begin log file
      OPEN (9,FILE='u2rdplot.log')
      WRITE (9,'(a)') title

* read index
      WRITE (*,'(a,a)') 'open index file = ',fnindex
CC    OPEN (13,FILE=fnindex,FORM='unformatted',READONLY)
      OPEN (13,FILE=fnindex,FORM='unformatted')
      READ (13,ERR=99) nx
      CLOSE(13)
      DO zn=1,80,10
        WRITE (*,'(a,i3,2i9)') 
     .   'zn, nx(1, 240 = ',zn, nx(zn,1), nx(zn,240)
      ENDDO
      WRITE (*,'(a,$)') 'index read successfully  '
      READ (*,'(a)') a1

* loop plots

 101  READ (11,*,END=99) ra1,ra2,dc1,dc2, ptyp
      nplot = nplot + 1
      nrd = 0           ! all records read
      n0t = 0 
      n1t = 0 
      n2t = 0 
      n3t = 0
      n4t = 0

      WRITE (*,'(/a,i2)') 'start plot = ',nplot
      WRITE (*,'(a,4f7.2,i2,a)') 'ra1,2, dc1,2, ptyp = '
     .   ,ra1,ra2,dc1,dc2, ptyp,'  '
      WRITE (9,'(/a,4f7.2,i2,a)') 'ra1,2, dc1,2, ptyp = '
     .   ,ra1,ra2,dc1,dc2, ptyp,'  '

      IF (ptyp.EQ.0) THEN
        WRITE (*,'(a,4f7.2)') '-- pole projection',ra1,ra2,dc1,dc2
        just = 1    ! 1 = equal scale on both axis
      ELSEIF (ptyp.EQ.1) THEN
        WRITE (*,'(a,4f7.2)') '-- box projection',ra1,ra2,dc1,dc2
        just = 0    ! 1 = equal scale on both axis
      ELSE
        WRITE (*,'(a,i3)') '-- invalid projection type',ptyp 
        STOP
      ENDIF

      CALL get_zone_range (dc1,dc2,zmax, d1m,d2m,z1,z2,nz)
      CALL get_ra_range (ra1,ra2, r1m,r2m,i1,i2,nr)

      WRITE (9,'(a,4i5)') 'z1,z2, i1,2 = ',z1,z2,i1(1),i2(1)
      WRITE (*,'(a,4i5,a,$)') 
     .  'z1,z2, i1,2 = ',z1,z2,i1(1),i2(1),'  '
      READ (*,'(a)') a1 

* prepare plot
      IF (PGBEG(0,pgdev,1,1).NE.1) THEN
        PRINT*,'-- PGBEG failed'
        STOP
      ENDIF

      CALL PGBBUF              ! begin buffer
      CALL PGSLW (lwt)         ! set line width
      CALL PGSCH (st)          ! set character size
      CALL PGSCI (coltxt)      ! set color

      IF (ptyp.EQ.0) THEN      ! pole plot
        spd = 90.0 + dc2
        xmin= -spd
        xmax=  spd
        ymin= -spd
        ymax=  spd
        CALL PGENV  (xmin,xmax, ymin,ymax, just, axis)  ! draw box ...
        CALL PGLAB  (xtit0, ytit0, title)               ! text output
        xt = 0.8 * spd
        yt = 0.9 * spd
        CALL PGPTXT ( xt, yt, 0.0, 0.0, ' 3h')
        CALL PGPTXT ( xt,-yt, 0.0, 0.0, '21h')
        xt = 0.9 * spd
        yt = 0.9 * spd
        CALL PGPTXT (-xt, yt, 0.0, 0.0, ' 9h')
        CALL PGPTXT (-xt,-yt, 0.0, 0.0, '15h')

      ELSE                     ! box
        xmin= ra1 
        xmax= ra2 
        ymin= dc1 
        ymax= dc2 
        CALL PGENV  (xmin,xmax, ymin,ymax, just, axis)  ! draw box ...
        CALL PGLAB  (xtit1, ytit1, title)               ! text output
      ENDIF

* loop zones
      DO zn=z1,z2
        n0 = 0 
        n1 = 0 
        n2 = 0 
        n3 = 0
        n4 = 0

        CALL open_zfile (pathu,un,zn,only_rd)

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

        IF (i1(1).EQ.1) THEN  ! assume 1 RA range
          j1 = 1              ! start with first star on file
        ELSE
          j1 = nx (zn,i1(1)-1) - i0 + 1   ! i1-1 = previous bin
        ENDIF

        j2 = nx (zn,i2(1)) - i0

        WRITE (9,'(a,i3,2i4,3i9)') 
     .    'zn, i1,2, i0,j1,j2 = ', zn,i1(1),i2(1),i0,j1,j2

* read data
        DO recn= j1,j2                 ! assume 1 RA range
          CALL read_u2line (un,recn,bf,idat,errflg)

          IF (errflg) THEN
            WRITE (*,'(a,i7,i4)') 
     .       'read error, recn, zn = ',recn,zn
            GOTO 71
          ENDIF

	  nrd = nrd + 1

          IF (nrd.LE.5)               ! test printout
     .      WRITE (9,'(a,i8,2i12,i5)') 'nrd, ra,dc,m: ',
     .           nrd, idat(1),idat(2),idat(3) 

          IF (idat(3).GT.mag4) THEN   ! faintest limit
            n4 = n4 + 1               ! out of range, do not plot

          ELSEIF (idat(3).GT.mag3) THEN
            n3 = n3 + 1               ! within mag bin 3
            IF (n3.GT.dim3) THEN
              WRITE (*,'(a,i8)') 'dim3 too small, n3 = ',n3
              STOP
            ENDIF
            CALL calc_xy (ptyp,idat(1),idat(2),masrad,xp3(n3),yp3(n3)) 
            
          ELSEIF (idat(3).GT.mag2) THEN
            n2 = n2 + 1               ! within mag bin 2
            IF (n2.GT.dim2) THEN
              WRITE (*,'(a,i8)') 'dim2 too small, n2 = ',n2
              STOP
            ENDIF
            CALL calc_xy (ptyp,idat(1),idat(2),masrad,xp2(n2),yp2(n2)) 
            
          ELSEIF (idat(3).GT.mag1) THEN
            n1 = n1 + 1               ! within mag bin 1
            IF (n1.GT.dim1) THEN
              WRITE (*,'(a,i8)') 'dim1 too small, n1 = ',n1
              STOP
            ENDIF
            CALL calc_xy (ptyp,idat(1),idat(2),masrad,xp1(n1),yp1(n1)) 
            
          ELSE
            n0 = n0 + 1               ! too bright, do not plot

          ENDIF   ! mag bins
        ENDDO     ! RA range

  71    CONTINUE  ! after read error

* plot each zone

        CALL PGSLW (lws)         ! set line width
        CALL PGSCH (ss)          ! set character size

        CALL PGSCI (col1)  
        CALL PGPT  (n1, xp1,yp1, symbol) 
        CALL PGSCI (col2)  
        CALL PGPT  (n2, xp2,yp2, symbol) 
        CALL PGSCI (col3)  
        CALL PGPT  (n3, xp3,yp3, symbol) 

* sum counters
        n0t = n0t + n0
        n1t = n1t + n1
        n2t = n2t + n2
        n3t = n3t + n3
        n4t = n4t + n4
      ENDDO       ! Dec zones

* end plot
      CALL PGIDEN               ! user id, date, time
      CALL PGEBUF               ! end buffer
      CALL PGEND                ! end PGPLOT

      WRITE (*,'(/a)')     'number of star entire plot: ' 
      WRITE (*,'(a, i9)')  ' - read in :', nrd 
      WRITE (*,'(a,2i9)')  ' - out l/h :', n0t,n4t
      WRITE (*,'(a,3i9/)') ' - for plot:', n1t,n2t,n3t

      WRITE (9,'(/a)')     'number of star entire plot: ' 
      WRITE (9,'(a, i9)')  ' - read in :', nrd 
      WRITE (9,'(a,2i9)')  ' - out l/h :', n0t,n4t
      WRITE (9,'(a,3i9/)') ' - for plot:', n1t,n2t,n3t

      WRITE (cmd,'(a,a,i2.2,a)') 
     .  'mv pgplot.ps ',fnbase(1:jb),nplot,'.ps'
      WRITE (*,'(a)') cmd

      ierr = SYSTEM (cmd)

      GOTO 101                  ! next plot

 99   WRITE (*,'(/a/)') 'end of program <u2rdplot>'

      END

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

      SUBROUTINE calc_xy (protyp,ra,dc,masrad, x,y)
C
      IMPLICIT NONE
      INTEGER  protyp, ra,dc
      REAL*4   masrad, x,y, rarad, spd

      IF (protyp.EQ.0) THEN        ! projection = pole plot
        rarad = ra * masrad        ! mas -> radian
        spd   = 90.0 + dc / 3.6e6  ! mas -> degree -> south pole distance
        x = spd * COS (rarad)
        y = spd * SIN (rarad)
 
      ELSE
        x = ra / 5.4e7             ! mas -> hour
        y = dc / 3.6e6 
      ENDIF

      END   ! subr. <calc_xy>