00001
00002 c
00003 c Copyright (c) 1986,1987,1988,1989,1990,1991,1992,1993,
00004 c by Steve McMillan, Drexel University, Philadelphia, PA.
00005 c
00006 c All rights reserved.
00007 c
00008 c Redistribution and use in source and binary forms are permitted
00009 c provided that the above copyright notice and this paragraph are
00010 c duplicated in all such forms and that any documentation,
00011 c advertising materials, and other materials related to such
00012 c distribution and use acknowledge that the software was developed
00013 c by the author named above.
00014 c
00015 c THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
00016 c IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
00017 c WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
00018 c
00019
00020 c------------------------------------------------------------------------
00021 c
00022 c Histogram-plotting routines: histogram
00023 c histo (used by histogram)
00024 c
00025 c------------------------------------------------------------------------
00026
00027 subroutine histogram(input,istart,nin,x,y,z,nx,ny,nz,
00028 & iherr,ihmode,ihsave,iprompt,*)
00029 save
00030 c
00031 c Draw a histogram from z or y(x).
00032 c
00033 character*(*) input
00034 dimension x(*),y(*),z(*)
00035 c
00036 common/scales/xl,xr,dinchx,ybot,ytop,dinchy,rlen,slen
00037 c
00038 common/draw params/roff,soff,aspect1,xlen,ylen,hs,hn,hp,
00039 & idevset,jbox,iorig
00040 c
00041 common/frame params 1/xmin,xmax,ymin,ymax,modex,modey
00042 character*80 xttl,yttl
00043 common/frame params 2/xttl,yttl
00044 c
00045 parameter (NHISTOMAX=200)
00046 dimension xh(0:NHISTOMAX),yh(0:NHISTOMAX),xp(19)
00047 logical normal
00048 c
00049 common/weights/wt(NHISTOMAX),nw
00050 c
00051 if (input(3:3).eq.'x') then
00052 c
00053 if (nx.le.1) go to 1001
00054 c
00055 c Already have the data in x, y, and the box is drawn.
00056 c
00057 dx2 = .5*(x(2) - x(1))
00058 yh(0) = 0.
00059 xh(0) = max(xl,min(xr,x(1)-1.5*dx2))
00060 do 10180 i=1,nx
00061 xh(i) = max(xl,min(xr,x(i)-dx2))
00062 yh(i) = y(i)
00063 10180 continue
00064 c
00065 if (.not.normal) then
00066 ier = 1
00067 else
00068 ier = nz
00069 end if
00070 c
00071 call histo(xh,yh,nx+1,ihmode,ier*iherr)
00072 return
00073 c
00074 end if
00075 c
00076 if (nz.le.1) go to 1001
00077 c
00078 c Must make the histogram from scratch.
00079 c
00080 call minmax(z,nz,z1,z2)
00081 xmin0 = z1
00082 xmax0 = z2
00083 c
00084 c Determine binning to use:
00085 c
00086 if (input(3:3).ne.'2') then
00087 if (input(3:3).eq.'n'.or.input(4:4).eq.'n') then
00088 normal = .true.
00089 else
00090 normal = .false.
00091 end if
00092 c
00093 ref = -1.e35
00094 io = 1
00095 call readrq(input(istart:nin),2,
00096 & del,ref,dum,dum,*100)
00097 io = 0
00098 c
00099 100 if (io.ne.0.or.ref.lt.-9.e34) then
00100 read(input(istart:nin),*,err=1001,end=1001)del
00101 ref = z1
00102 end if
00103 if (del.le.0.) go to 1001
00104 c
00105 xmin = z1
00106 xmax = z2
00107 c
00108 if (modex.lt.0) then
00109 if (z1.le.0..or.ref.le.0.) go to 1001
00110 z1 = log10(z1)
00111 z2 = log10(z2)
00112 ref = log10(ref)
00113 end if
00114 c
00115 iref = 0.9999*(ref-z1)/del
00116 nh = 0
00117 xh(0) = ref - (iref+1)*del
00118 iref = iref + 2
00119 c
00120 c Notes: xh(i) is the TOP of bin #i.
00121 c ref will be at the BOTTOM of bin #iref
00122 c
00123 180 nh = nh+1
00124 if (nh.gt.NHISTOMAX) then
00125 if (iprompt.eq.1) write(6,'(''Too many bins.'')')
00126 go to 1001
00127 end if
00128 xh(nh) = xh(nh-1) + del
00129 if (xh(nh).le.z2) go to 180
00130 else
00131 if (del.le.0..or.nh.le.0..or.nh.gt.nhmax) go to 1001
00132 end if
00133 c
00134 c Bin the data.
00135 c
00136 do 182 j=0,nh
00137 yh(j) = 0.
00138 182 continue
00139 zbar=0.
00140 var=0.
00141 do 184 i=1,nz
00142 zz = z(i)
00143 if (modex.lt.0) zz = log10(zz)
00144 zbar = zbar + z(i)
00145 var = var + z(i)**2
00146 j = iref + (zz-ref)/del
00147 j = max(1,min(nh,j))
00148 yh(j) = yh(j) + 1.
00149 184 continue
00150 c
00151 if (normal) then
00152 do 185 j = 0,nh
00153 yh(j) = yh(j)/nz
00154 185 continue
00155 end if
00156 c
00157 c Get statistics (percentiles from the binned data)...
00158 c
00159 zbar = zbar/nz
00160 var = var/nz - zbar**2
00161 c
00162 ix=0
00163 sum=0.
00164 do 190 ip=1,19
00165 sp = .05*ip*nz
00166 if (normal) sp = sp/nz
00167 188 ix = ix + 1
00168 if (ix.gt.nh) go to 191
00169 sum = sum + yh(ix)
00170 if (sum.lt.sp) go to 188
00171 xp(ip) = xh(ix) - del*(sum-sp)/yh(ix)
00172 sum = sum - yh(ix)
00173 ix = ix - 1
00174 190 continue
00175 c
00176 c ...and (clip and) draw the histogram.
00177 c
00178 191 call minmax(yh(1),nh,y1,y2)
00179 c
00180 if (modey.gt.0) then
00181 ymin = 0.
00182 else
00183 if (.not.normal) then
00184 ymin = 1.
00185 else
00186 ymin = 1./nz
00187 end if
00188 ymin1 = log10(ymin) - 10.
00189 do 192 j = 0,nh
00190 if (yh(j).lt.ymin) then
00191 yh(j) = ymin1
00192 else
00193 yh(j) = log10(max(ymin,yh(j)))
00194 end if
00195 192 continue
00196 end if
00197 c
00198 if (input(3:3).ne.'1'.and.input(3:3).ne.'2') then
00199 c
00200 ymax = min(10.+y2,1.1*y2)
00201 c
00202 delz = xp(19) - xp(1)
00203 c
00204 if (xp(1)-z1.gt.delz) then
00205 c
00206 c Bottom 5% too spread out:
00207 c
00208 xmin = xp(19) - sqrt(delz*(xp(19)-z1))
00209 ix = abs(xmin - ref)/del
00210 if (xmin.lt.ref) then
00211 xmin = ref - (ix+1)*del
00212 else
00213 xmin = ref + ix * del
00214 end if
00215 if (modex.lt.0) xmin = 10.**z1
00216 end if
00217 c
00218 if (z2-xp(19).gt.delz) then
00219 c
00220 c Top 5% too spread out:
00221 c
00222 xmax = xp(1) + sqrt(delz*(z2-xp(1)))
00223 ix = abs(xmax - ref)/del
00224 if (xmax.lt.ref) then
00225 xmax = ref - ix*del
00226 else
00227 xmax = ref + (ix+1) * del
00228 end if
00229 if (modex.lt.0) xmax = 10.**z2
00230 end if
00231 c
00232 c Expand the range slightly (linear case only):
00233 c
00234 if (modex.gt.0) then
00235 xav = .5*(xmin+xmax)
00236 xmin = xmin - .1*(xav-xmin)
00237 xmax = xmax + .1*(xav-xmin)
00238 end if
00239 c
00240 if (idevset.eq.0) call setup(' ','histogram')
00241 call devon
00242 call clear
00243 c
00244 if (.not.normal) then
00245 yttl = '@H'
00246 else
00247 yttl = ' fraction'
00248 end if
00249 c
00250 call eframe(xmin,xmax,xlen,modex,xttl,
00251 & ymin,ymax,ylen,modey,yttl)
00252 ibox=1
00253 c
00254 c More clipping...
00255 c
00256 do 193 j=0,nh
00257 if (xh(j).lt.xl) xh(j)=xl
00258 if (xh(j).gt.xr) xh(j)=xr
00259 193 continue
00260 end if
00261 c
00262 if (.not.normal) then
00263 ier = 1
00264 else
00265 ier = nz
00266 end if
00267 c
00268 call histo(xh,yh,nh+1,ihmode,ier*iherr)
00269 call devoff
00270 c
00271 if (ihsave.ne.0) then
00272 c
00273 c Save the histogram data (top center of each bar) in y(x).
00274 c
00275 do 200 j=2,nh
00276 x(j-1) = .5*(xh(j-1)+xh(j))
00277 y(j-1) = yh(j)
00278 if (ihsave.eq.2.and.nw.gt.0) z(j-1) = wt(j)
00279 200 continue
00280 nx = nh-1
00281 ny = nx
00282 end if
00283 if (ihsave.eq.2.and.nw.gt.0) nz = nx
00284 c
00285 if (iprompt.eq.1) then
00286 if (modex.lt.0) then
00287 do 194 ip=1,19
00288 xp(ip) = 10.**xp(ip)
00289 194 continue
00290 end if
00291 write(6,195)nz,xmin0,xmax0,zbar,
00292 & sqrt(var),(xp(ip),ip=2,18,2)
00293 195 format(i6,' points: ',
00294 & ' xmin = ',1pe10.3,', xmax = ',e10.3/
00295 & 15x,' mean = ',e10.3,', s.d. = ',e10.3/
00296 & ' %tiles:',5e12.3/8x,4e12.3)
00297 end if
00298 c
00299 return
00300 1001 return 1
00301 c
00302 end
00303
00304
00305 subroutine histo(x,y,n,ihmode,iherr)
00306 save
00307 c
00308 c Draw y(x) as a histogram or as data points with error bars.
00309 c
00310 dimension x(*),y(*)
00311 c
00312 common/frame params 1/xmin,xmax,ymin,ymax,modex,modey
00313 common/scales/xl,xr,dinchx,ybot,ytop,dinchy,rlen,slen
00314 common/draw params/roff,soff,aspect1,xlen,ylen,hs,hn,hp,
00315 & idevset,jbox,iorig
00316 c
00317 parameter (NHISTOMAX=200)
00318 common/weights/wt(NHISTOMAX),nw
00319 c
00320 nw = 0
00321 c
00322 if (iherr.le.0) then
00323 r = (x(1)-xl)*dinchx
00324 so = 0.
00325 do 100 i=2,n
00326 s = max(y(i)-ybot,0.)*dinchy
00327 if (ihmode.eq.0) then
00328 call plot(r,0.,3)
00329 else
00330 call plot(r,min(s,so),3)
00331 end if
00332 call plot(r,max(s,so),2)
00333 call plot(r,s,3)
00334 r = (x(i)-xl)*dinchx
00335 call plot(r,s,2)
00336 so = s
00337 if (r.ge.rlen) go to 101
00338 100 continue
00339 101 if (r.lt.rlen) call plot(r,0.,2)
00340 else
00341 dx = .5*(x(3)-x(2))
00342 if (hp.le.0.) then
00343 nn = 4
00344 hh = .025
00345 else
00346 nn = 10
00347 hh = .5*hp
00348 end if
00349 c
00350 do 250 i = 2,n
00351 if (x(i).ge.xr.or.y(i).lt.ybot) go to 250
00352 r = (x(i)-dx-xl)*dinchx
00353 s = (y(i)-ybot)*dinchy
00354 call ngon(r,s,hh,nn,0.)
00355 c
00356 c Find the actual number of points in the bin, zbin.
00357 c
00358 zbin = y(i)
00359 if (modey.lt.0) zbin = 10.**zbin
00360 zbin = zbin*iherr
00361 c
00362 c Determine the standard error.
00363 c
00364 nw = nw + 1
00365 wt(i-1) = sqrt(zbin)
00366 c
00367 c (Note that wt(1) corresponds to x(2)!)
00368 c
00369 zfac = 1./wt(i-1)
00370 if (zbin.le.1.) zfac = .999
00371 c
00372 if (modey.gt.0) then
00373 su = (y(i)*(1.+zfac) - ybot)*dinchy
00374 208 sl = (y(i)*(1.-zfac) - ybot)*dinchy
00375 else
00376 su = (y(i) + log10(1.+zfac) - ybot)*dinchy
00377 sl = (y(i) + log10(1.-zfac) - ybot)*dinchy
00378 end if
00379 c
00380 c Plot the error bars.
00381 c
00382 call plotin(r-hh,sl,3)
00383 call plotin(r+hh,sl,2)
00384 call plotin(r,sl,3)
00385 call plotin(r,su,2)
00386 call plotin(r-hh,su,3)
00387 call plotin(r+hh,su,2)
00388 c
00389 250 continue
00390 c
00391 end if
00392 c
00393 end
