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 ' 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.