00001 program plot_SeBa
00002
00003 c general plotting program for reduced SeBa output
00004
00005 c parameter(Nmax=36389,Np=30000)
00006 parameter(Nmax=100000,Nr=20,Ni=6)
00007
00008 real all_r(Nmax,Nr),arr1(Nmax),arr2(Nmax)
00009 integer all_i(Nmax,Ni)
00010
00011 real rmin(Nr),rmax(Nr),xmin,xmax,ymin,ymax
00012 integer nbin,istat1,istat2
00013 character*30 label(Nr),xlabel,ylabel,filename,filename2
00014 character*100 title
00015 integer pgopen
00016
00017 data rmin/0.1,1.,1.,0.,0.,0.8,0.5,0.8,0.5,0.1,0.1,
00018 + 1.e-4,0.,0.,0.1,0.05,0.1,0.05,1.,1./
00019 data rmax/0.1,1.e4,1.e3,1.,1.,10.,10.,10.,10.,1.e4,1.e2,
00020 + 1.e2,1.,1.,10.,10.,10.,10.,20.,20./
00021
00022 data label/'T\di\u (Myr)','a\di\u (R\d\(2281)\u)','P\di\u (d)',
00023 + 'q\di\u',
00024 + 'ecc\di\u','M\di\u (M\d\(2281)\u)','R\di\u (R\d\(2281)\u)',
00025 + 'm\di\u (M\d\(2281)\u)','r\di\u (R\d\(2281)\u)',
00026 + 'T\df\u (Myr)','a\df\u (R\d\(2281)\u)','P\df\u (d)',
00027 + 'q\df\u',
00028 + 'ecc\df\u','M\df\u (M\d\(2281)\u)','R\df\u (R\d\(2281)\u)',
00029 + 'm\df\u (M\d\(2281)\u)','r\df\u (R\d\(2281)\u)',
00030 + 'id1','id2'/
00031
00032 title = ' '
00033 nbin=30
00034
00035
00036 write(6,*) 'give filename'
00037 read(5,"(A)") filename
00038
00039 write(6,*) filename
00040 n=0
00041 call read_in(all_r,all_i,filename,n,Nmax,Nr,Ni)
00042
00043 write(6,*) n
00044 do i=1,n
00045 c Pi (in days), qi
00046 all_r(i,3) = 0.116*sqrt(all_r(i,2)**3/(all_r(i,6)+all_r(i,8)))
00047 all_r(i,4) = all_r(i,8)/all_r(i,6)
00048 c Pf (in days), qf
00049 all_r(i,12) = 0.116*sqrt(all_r(i,11)**3
00050 & / (all_r(i,15)+all_r(i,17)))
00051 if (all_r(i,15).eq.0.or.all_r(i,17).eq.0) then
00052 all_r(i,13) = 0
00053 all_r(i,17) = max(all_r(i,15),all_r(i,17))
00054 all_r(i,15) = 0
00055 else
00056 all_r(i,13) = all_r(i,17)/all_r(i,15)
00057 end if
00058 all_r(i,19) = all_i(i,5)
00059 all_r(i,20) = all_i(i,6)
00060 end do
00061
00062 c istat1=pgopen('?')
00063 call pgbeg(1,'?',2,2)
00064 istat1 =1
00065 print*, istat1
00066 c call pgsubp(2,2)
00067 call pgslw(3)
00068 call pgsch(1.5)
00069
00070 10 write(6,*) 'Choose which parameters to plot'
00071 write(6,*) ' 1 Ti 2 ai 3 Pi'
00072 write(6,*) ' 4 qi 5 ecci 6 Mi'
00073 write(6,*) ' 7 Ri 8 mi 9 ri'
00074 write(6,*) '10 Tf 11 af 12 Pf'
00075 write(6,*) '13 qf 14 eccf 15 Mf'
00076 write(6,*) '16 Rf 17 mf 18 rf'
00077 write(6,*) '19 id1 20 id2'
00078
00079 read(5,*) inum1,inum2
00080
00081 if (inum1.eq.inum2) then
00082 do i=1,n
00083 if (inum1.gt.0) then
00084 arr1(i)=all_r(i,inum1)
00085 xmin = rmin(inum1)
00086 xmax = rmax(inum1)
00087 xlabel = label(inum1)
00088 else
00089 arr1(i)=log10(all_r(i,-inum1))
00090 xmin = log10(rmin(-inum1))
00091 xmax = log10(rmax(-inum1))
00092 xlabel = 'log '
00093 end if
00094 end do
00095 20 call pghist(n,arr1,xmin,xmax,nbin,0)
00096 call pglabel(xlabel,'N',title)
00097 write(6,*) 'change limits?'
00098 read(5,"(A)") ch
00099 if (ch.eq.'y'.or.ch.eq.'Y') then
00100 write(6,*) 'give new limits'
00101 read(5,*) xmin,xmax
00102 rmin(-inum1) = xmin
00103 rmax(-inum1) = xmax
00104 goto 20
00105 end if
00106 else
00107 do i=1,n
00108 if (inum1.gt.0) then
00109 arr1(i)=all_r(i,inum1)
00110 xmin = rmin(inum1)
00111 xmax = rmax(inum1)
00112 xlabel = label(inum1)
00113 else
00114 arr1(i)=log10(all_r(i,-inum1))
00115 xmin = log10(rmin(-inum1))
00116 xmax = log10(rmax(-inum1))
00117 xlabel = 'log '
00118 end if
00119 if (inum2.gt.0) then
00120 arr2(i)=all_r(i,inum2)
00121 ymin = rmin(inum2)
00122 ymax = rmax(inum2)
00123 ylabel = label(inum2)
00124 else
00125 arr2(i)=log10(all_r(i,-inum2))
00126 ymin = log10(rmin(-inum2))
00127 ymax = log10(rmax(-inum2))
00128 ylabel = 'log '
00129 end if
00130 end do
00131 30 call spzgray(n,arr1,xmin,xmax,arr2,ymin,
00132 + ymax,nbin,nbin,0.,1.,0)
00133 call pglabel(xlabel,ylabel,title)
00134 write(6,*) 'change limits?'
00135 read(5,"(A)") ch
00136 if (ch.eq.'y'.or.ch.eq.'Y') then
00137 write(6,*) 'give new x-limits'
00138 read(5,*) xmin,xmax
00139 if (inum1.gt.0) then
00140 rmin(inum1) = xmin
00141 rmax(inum1) = xmax
00142 else
00143 rmin(-inum1) = 10**xmin
00144 rmax(-inum1) = 10**xmax
00145 end if
00146 write(6,*) 'give new y-limits'
00147 read(5,*) ymin,ymax
00148 if (inum2.gt.0) then
00149 rmin(inum2) = ymin
00150 rmax(inum2) = ymax
00151 else
00152 rmin(-inum2) = 10**ymin
00153 rmax(-inum2) = 10**ymax
00154 end if
00155 goto 30
00156 end if
00157 end if
00158
00159 write(6,*) 'Stop? (s) or Print (p)'
00160 read(5,"(A)") ch
00161
00162 if (ch.ne.'s'.and.ch.ne.'p') goto 10
00163
00164 if (ch.eq.'p') then
00165 write(6,*) 'give filename'
00166 read(5,*) filename
00167 c istat2 = pgopen(filename
00168 write(6,*) 'give title (optional)'
00169 read(5,"(A)") title
00170
00171 if (inum1.eq.inum2) then
00172 call pghist(n,arr1,xmin,xmax,nbin,0)
00173 call pglabel(xlabel,'N',title)
00174 else
00175 call spzgray(n,arr1,xmin,xmax,arr2,ymin,
00176 + ymax,nbin,nbin,0.,1.,0)
00177 call pglabel(xlabel,ylabel,title)
00178 end if
00179 c call pgclos
00180 print*, istat1, istat2
00181 c call pgslct(istat1)
00182 goto 10
00183 end if
00184
00185 c call pgclos
00186 c call end
00187 end
00188
00189 SUBROUTINE read_in(all_r, all_i, filename, n, Nmax, Nr, Ni)
00190
00191 real all_r(Nmax,Nr)
00192 integer all_i(Nmax,Ni)
00193 character*30 filename
00194
00195 open(10,file=filename,status="old")
00196
00197 n=0
00198 do i=1,Nmax
00199 read(10,*,end=99) all_i(i,1),all_r(i,1),all_r(i,2),
00200 + all_r(i,5),all_i(i,2),all_r(i,6),all_r(i,7),
00201 + all_i(i,3),all_r(i,8),all_r(i,9)
00202 read(10,*,end=99) all_i(i,4),all_r(i,10),all_r(i,11),
00203 + all_r(i,14),all_i(i,5),all_r(i,15),all_r(i,16),
00204 + all_i(i,6),all_r(i,17),all_r(i,18)
00205 n=n+1
00206 end do
00207
00208 99 close(1)
00209 return
00210
00211 end
00212
00213 subroutine spzgray(ndata, xdata, xmin, xmax, ydata, ymin, ymax
00214 1 , nxbin,nybin, iflag)
00215 integer ndata, nxbin, nybin, iflag
00216 real xdata, ydata, xy_data
00217 real xmin, xmax, ymin, ymax
00218 dimension xdata(ndata), ydata(ndata), xy_data(nxbin, nybin)
00219 real dx, dy, tr(6), fg, bg, norm
00220 integer nxmin, nxmax, ix, iy, it
00221 norm=0
00222 nxmin=1
00223 nymin=1
00224 dx = (xmax-xmin)/(1.0*nxbin)
00225 dy = (ymax-ymin)/(1.0*nybin)
00226 bg = 0.0
00227 fg = 1.0
00228 tr(1) = xmin
00229 tr(2) = dx
00230 tr(3) = 0
00231 tr(4) = ymin
00232 tr(5) = 0
00233 tr(6) = dy
00234 print *, ndata, xmin, xmax, ymin, ymax, nxbin,nybin, iflag
00235 print *, dx, dy
00236
00237 do ix=1, nxbin
00238 do iy=1, nybin
00239 xy_data(ix, iy) = 0
00240 enddo
00241 enddo
00242 do i=1, ndata
00243 c print *, i, xdata(i), ydata(i)
00244 ix = 0.5 + (xdata(i)-xmin)/dx
00245 iy = 0.5 + (ydata(i)-ymin)/dy
00246 c print *, ix, iy, xdata(i), ydata(i)
00247 if(ix.ge.1.and.iy.ge.1.and.ix.le.nxbin.and.iy.le.nybin)then
00248 xy_data(ix, iy) = xy_data(ix, iy) + 1
00249 norm = max(norm, xy_data(ix, iy))
00250 c print *, ix, it, xy_data(ix, it), norm
00251 endif
00252 enddo
00253 do ix=1, nxbin
00254 do iy=1, nybin
00255 xy_data(ix, iy) = xy_data(ix, iy)/norm
00256 enddo
00257 enddo
00258
00259 if(iflag.eq.0)then
00260 call pgenv(xmin,xmax, ymin,ymax,0,0)
00261 endif
00262 call pggray(xy_data, nxbin, nybin, nxmin, nxbin, nymin,
00263 1 nybin, fg, bg, tr)
00264 end
00265
00266
00267
