ccc
c  SUBROUTINE PLOT3D
c
c  This routine plots three dimensional surfaces on a rectangular grid with
c  hidden line removal.  The array f(nx,ny) contains the values of the surface
c  above a given grid point.  The arrays x(nx,ny) and y(nx,ny) are internal 
c  work arrays containing the x and y coordinate of the projected surface.
c  The variables nx and ny are the dimensions of the arrays.  The variables
c  xw and yw are the size of the projected plot in inches.  The variables
c  xview, yview and zview are the position at which the surface is to be
c  viewed in which the grid lies in a rectangle bounded by the diagonal 
c  vertices (0,0) and (xw,yw).  The variable line determines which directions
c  the grid lines are plotted.  For line = 1, only those grid lines parallel
c  to the x-axis are plotted.  For line = 2, only those grid lines parallel
c  to the y-axis are plotted.  For line = 0 all grid lines are plotted.
c  The variable scale determines the scaling factor for the surface to be 
c  plotted.  If scale = 0, the routine will choose a suitable scaling factor.
c
c
	subroutine plot3d(f,x,y,nx,ny,xw,yw,xview,yview,zview,line,scale)
c
	dimension f(nx,ny),x(nx,ny),y(nx,ny)
c
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, pdx, pdy, dxb, dyb, ixcur, iycur, xmido, ymido, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
	common/pscale/xv,yv,dx,dy,ca,sa,ly,lx,xscale,xn,yscale,yn
	common/color/isubch
c
	data small,pt5/1.e-06,.5/
c
	if(line.le.2)goto 10
	line=line-3
	goto 305
c
c  normalize constants
c
10	xb=1.
	yb=yw/xw
	xv=xview/xw
	yv=yview/xw
	h=zview/xw
	ir=1
	nxx=nx+1
	nyy=ny+1
c
c  determine projection plane
c
	a=atan2(.5*xb-xv,.5*yb-yv)
	sa=sin(a)
	ca=cos(a)
c
c  y component of the vertices of the projected grid
c
	y2=xb*sa
	y3=yb*ca
	y4=y2+y3
	yc=xv*sa+yv*ca
	yn=amin1(0.,y2,y3,y4)
	if(yc.ge.yn)yc=yn-.01*abs(yb)
	yc=-yc
	xc=xv*ca-yv*sa
	dx=xb/(nx-1)
	dxx=1./dx
	vx=dxx*xv+1.00001
	lx=vx
	dy=yb/(ny-1)
	dyy=1./dy
	vy=dyy*yv+1.00001
	ly=vy
	yn=0.
	yx=yn
	d=.5*y4
	yd=yc+d
	yp=0.
c
c  determine minimum and maximum of the surface
c
	tn=f(1,1)
	tx=tn
	do 50 j=1,ny
	do 50 i=1,nx
	tn=amin1(tn,f(i,j))
50	tx=amax1(tx,f(i,j))
c
c  compute scaling factor
c
	s=tx-tn
	if(s.eq.0.)s=tx
	s=amax1(xb,yb)/s
	if(scale.ne.0.)s=scale
c
c  compute projected x and y coordinates of surface
c
	do 120 j=1,ny
	xp=0.
	c=yp*sa+xc
	do 100 i=1,nx
	t=f(i,j)*s
	hl=xp*sa+yp*ca
	x(i,j)=(xp*ca-c)*yd/(hl+yc)+small
	y(i,j)=((hl-d)*h+t*yd)/(hl+yc)+small
	yn=amin1(yn,y(i,j))
	yx=amax1(yx,y(i,j))
100	xp=xp+dx
120	yp=yp+dy
c
c  x component of the vertices of the grid
c
	x1=x(1,1)
	x2=x(nx,1)
	x3=x(1,ny)
	x4=x(nx,ny)
	tn=tn*yd*s*1.1
c
c  y component of the vertices of the projected grid
c
	y1=(tn-d*h)/yc
	y2=((y2-d)*h+tn)/(y2+yc)
	y3=((y3-d)*h+tn)/(y3+yc)
	y4=((y4-d)*h+tn)/(y4+yc)
c
	xn=amin1(x1,x2,x3,x4)
	xx=amax1(x1,x2,x3,x4)
	yn=amin1(yn,y1,y2,y3,y4)
	xscale=scrx*xw/(xx-xn)
	yscale=scry*yw/(yx-yn)
	xn=xscale*xn
	yn=yscale*yn
c
c  convert to screen cordinates
c
	do 200 j=1,ny
	do 200 i=1,nx
	x(i,j)=aint(xscale*x(i,j)-xn+pt5)
200	y(i,j)=aint(yscale*y(i,j)-yn+pt5)
c
	x1=aint(xscale*x1-xn+pt5)
	y1=aint(yscale*y1-yn+pt5)
	x2=aint(xscale*x2-xn+pt5)
	y2=aint(yscale*y2-yn+pt5)
	x3=aint(xscale*x3-xn+pt5)
	y3=aint(yscale*y3-yn+pt5)
	x4=aint(xscale*x4-xn+pt5)
	y4=aint(yscale*y4-yn+pt5)
	xscale=xscale*yd
c
c  begin hidden line removal
c
	do 300 j=1,ny
	do 300 i=1,nx
	ix=isign(1,lx-i)
	iy=isign(1,ly-j)
	ih=i
	iplot=1
	if(0.eq.lx-i)goto 250
215	ih=ih+ix
	if(ih.lt.1.or.ih.gt.nx)goto 250
	xp=(ih-1)*dx
	xm=dx/dy*(dy*(j-1)-yv)/(dx*(i-1)-xv)
	yp=xm*(ih-i)+j
	jj=ifix(yp)
	jh=jj+1
	if(jj.lt.1.or.jh.gt.ny)goto 250
	yp=(yp-1)*dy
	hl=yc+xp*sa+yp*ca
	xi=aint(xscale*(xp*ca-yp*sa-xc)/hl-xn+pt5)
	yi=(y(ih,jj)-y(ih,jh))*(xi-aint(x(ih,jh)))/
	1(aint(x(ih,jj))-aint(x(ih,jh)))+y(ih,jh)
	if(aint(yi+pt5).lt.y(i,j))goto 215
	iplot=0
	goto 290
250	if(ly-j.eq.0)goto 290
	jh=j
260	jh=jh+iy
	if(jh.lt.1.or.jh.gt.ny)goto 290
	yp=(jh-1)*dy
	ym=dy/dx*(dx*(i-1)-xv)/(dy*(j-1)-yv)
	xp=ym*(jh-j)+i
	ii=ifix(xp)
	ih=ii+1
	if(ii.lt.1.or.ih.gt.nx)goto 290
	xp=(xp-1)*dx
	hl=yc+xp*sa+yp*ca
	xi=aint(xscale*(xp*ca-yp*sa-xc)/hl-xn+pt5)
	yi=(y(ii,jh)-y(ih,jh))*(xi-aint(x(ih,jh)))/
	1(aint(x(ii,jh))-aint(x(ih,jh)))+y(ih,jh)
	if(aint(yi+pt5).lt.y(i,j))goto 260
	iplot=0
290	x(i,j)=x(i,j)+.101*iplot
300	continue
c
c  plot x grid lines
c
305	if(line.eq.2)goto 350
	isubch=1
	do 310 j=1,ny
	iplot0=ifix(10*(x(1,j)-ifix(x(1,j))))
	if(iplot0.eq.1)call plot(ifix(x(1,j)),ifix(y(1,j)),0)
	do 310 i=2,nx
	iplot=ifix(10*(x(i,j)-aint(x(i,j))))
	call insect(i,j,iplot,iplot0,x,y,nx,ny,0)
310	iplot0=iplot
c
c  plot visible verticle bars
c
	if(lx.ge.1.and.lx.le.nx)goto 350
	if(lx.lt.1)goto 320
	ii=nx
	aa=(y2-y4)/(x2-x4)
	bb=y4-x4*aa
	ix=x2
	iy=y2
	goto 330
320	ii=1
	aa=(y3-y1)/(x3-x1)
	bb=y1-x1*aa
	ix=x1
	iy=y1
330	do 340 j=1,ny
	call plot(ifix(x(ii,j)),ifix(y(ii,j)),0)
340	call plot(ifix(x(ii,j)),ifix(aa*x(ii,j)+bb),1)
	call plot(ix,iy,1)
c
c  plot y grid lines
c
350	if(line.eq.1)return
	isubch=2
	do 400 i=1,nx
	iplot0=ifix(10*(x(i,1)-ifix(x(i,1))))
	if(iplot0.eq.1)call plot(ifix(x(i,1)),ifix(y(i,1)),0)
	do 400 j=2,ny
	iplot=ifix(10*(x(i,j)-aint(x(i,j))))
	call insect(i,j,iplot,iplot0,x,y,nx,ny,1)
400	iplot0=iplot
c
c  plot visible verticle bars
c
	if(ly.ge.1.and.ly.le.ny)return
	if(ly.lt.1)goto 420
	jj=ny
	aa=(y3-y4)/(x3-x4)
	bb=y4-x4*aa
	ix=x3
	iy=y3
	goto 430
420	jj=1
	aa=(y2-y1)/(x2-x1)
	bb=y1-x1*aa
	ix=x1
	iy=y1
430	do 440 i=1,nx
	call plot(ifix(x(i,jj)),ifix(y(i,jj)),0)
440	call plot(ifix(x(i,jj)),ifix(aa*x(i,jj)+bb),1)
	call plot(ix,iy,1)
	return
	end

ccc
c  SUBROUTINE INSECT
c
c  This routine computes the intersection point of a grid line which is 
c  partially hidden.
c
c
	subroutine insect(i1,j1,iplot,iplot0,x,y,nx,ny,l)
c
	dimension x(nx,ny),y(nx,ny)
	byte eq1,eq2
	common/pscale/xv,yv,dx,dy,ca,sa,ly,lx,xscale,xn,yscale,yn
	data pt5/.5/
	if(iplot.eq.iplot0)goto 290
	i=i1
	j=j1
	i0=i1
	j0=j1
	ix=isign(1,lx-i)
	iy=isign(1,ly-j)
	insct=0
	if(l.eq.0)goto 10
	j=j1-iplot0
	j0=j1-iplot
	xf=aint(x(i0,j0))
	yf=y(i0,j0)
	xl=aint(x(i,j))
	yl=y(i,j)
	imin=i1
	jmin=iy*min0(iy*j,iy*j0)
	if(jmin.eq.j0.and.(y(i+ix,j)-yl)/(x(i+ix,j)-xl).lt.
	1(yf-yl)/(xf-xl))goto 270
	goto 20
10	i=i1-iplot0
	i0=i1-iplot
	xf=aint(x(i0,j0))
	yf=y(i0,j0)
	xl=aint(x(i,j))
	yl=y(i,j)
	imin=ix*min0(ix*i,ix*i0)
	jmin=j1
	if(imin.eq.i0.and.(y(i,j+iy)-yl)/(x(i,j+iy)-xl).gt.
	1(yf-yl)/(xf-xl))goto 270
20	xint=xf-xl
	yint=yf-yl
	ih=imin
	if(lx-i.eq.0)goto 250
	xm1=dx/dy*(dy*(j-1)-yv)/(dx*(i-1)-xv)
	xm2=dx/dy*(dy*(j0-1)-yv)/(dx*(i0-1)-xv)
215	ih=ih+ix
	if(ih.lt.1.or.ih.gt.nx)goto 250
	yp1=xm1*(ih-i)+j
	if(i.eq.ih)yp1=yp1+iy
	yp2=xm2*(ih-i0)+j0
	if(i0.eq.ih)yp2=yp2+iy
	jn=ifix(amin1(yp1,yp2))
	jx=ifix(amax1(yp1,yp2))
	do 220 jj=jn,jx
	jh=jj+1
	if(jj.lt.1.or.jh.gt.ny)goto 220
	if((x(ih,jh)-aint(x(ih,jh)))+(x(ih,jj)-aint(x(ih,jj))).eq.0)
	1goto 220
	s1=(yf-yl)/(xf-xl)
	s2=(y(ih,jj)-y(ih,jh))/(aint(x(ih,jj))-aint(x(ih,jh)))
	if(s1.eq.s2)goto 220
	xi=(y(ih,jh)-yl+xl*s1-aint(x(ih,jh))*s2)/(s1-s2)
	if(abs(aint(xi+pt5)-.5*(xl+xf)).gt..5*abs(xf-xl))goto 220
	xmid=.5*aint(x(ih,jj)+x(ih,jh))
	xdif=.5*abs(aint(x(ih,jj))-aint(x(ih,jh)))
	if(abs(aint(xi+pt5)-xmid).gt.xdif)goto 220
	jh=(jj+jh+iy*(jj-jh))/2
	yi=s1*(xi-xl)+yl
	if(xint*(y(ih,jh+iy)-y(ih,jh))*ix*iy.lt.
	1yint*(aint(x(ih,jh+iy))-aint(x(ih,jh)))*ix*iy)goto 220
	xf=aint(xi+pt5)
	yf=aint(yi+pt5)
	insct=1
	if(xf-xl.eq.0)goto 270
220	continue
	goto 215
250	if(ly-j.eq.0)goto 270
	ym1=dy/dx*(dx*(i-1)-xv)/(dy*(j-1)-yv)
	ym2=dy/dx*(dx*(i0-1)-xv)/(dy*(j0-1)-yv)
	jh=jmin
260	jh=jh+iy
	if(jh.lt.1.or.jh.gt.ny)goto 270
	xp1=ym1*(jh-j)+i
	if(j.eq.jh)xp1=xp1+ix
	xp2=ym2*(jh-j0)+i0
	if(j0.eq.jh)xp2=xp2+ix
	jn=ifix(amin1(xp1,xp2))
	jx=ifix(amax1(xp1,xp2))
	do 265 ii=jn,jx
	ih=ii+1
	if(ii.lt.1.or.ih.gt.nx)goto 265
	if((x(ih,jh)-aint(x(ih,jh)))+(x(ii,jh)-aint(x(ii,jh))).eq.0)
	1goto 265
	eq1=i.eq.ih.or.i.eq.ii
	if(eq1.and.eq2)goto 265
	s1=(yf-yl)/(xf-xl)
	s2=(y(ii,jh)-y(ih,jh))/(aint(x(ii,jh))-aint(x(ih,jh)))
	if(s1.eq.s2)goto 265
	xi=(y(ih,jh)-yl+xl*s1-aint(x(ih,jh))*s2)/(s1-s2)
	if(abs(aint(xi+pt5)-.5*(xl+xf)).gt..5*abs(xf-xl))goto 265
	xmid=.5*aint(x(ii,jh)+x(ih,jh))
	xdif=.5*abs(aint(x(ii,jh))-aint(x(ih,jh)))
	if(abs(aint(xi+pt5)-xmid).gt.xdif)goto 265
	ih=(ih+ii+ix*(ii-ih))/2
	yi=s1*(xi-xl)+yl
	if(xint*(y(ih+ix,jh)-y(ih,jh))*ix*iy.gt.
	1yint*(aint(x(ih+ix,jh))-aint(x(ih,jh)))*ix*iy)goto 265
	xf=aint(xi+pt5)
	yf=aint(yi+pt5)
	insct=1
	if(xf.eq.xl)goto 270
265	continue
	goto 260
270	if(insct.eq.1)goto 280
	if(l.eq.0.and.i.eq.1.or.i.eq.nx)return
	if(l.eq.1.and.j.eq.1.or.j.eq.ny)return
	call plot(ifix(xl),ifix(yl),0)
	return
280	ix=ifix(xf)
	iy=ifix(yf)
	call plot(ix,iy,iplot0)
	if(xf.eq.xl)return
290	if(iplot.eq.0)return
	ix=ifix(x(i1,j1))
	iy=ifix(y(i1,j1))
	call plot(ix,iy,1)
	return
	end


	subroutine frame (xmin, xmax, ymin, ymax)
c
c	this routine frames a plot, i.e., no plotting of any kind
c	may take place outside of this frame without issuing a new
c	call to frame
c
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
c
	common /pagsiz/ xbond, ybond
c
	xn = amax1 ( 0., xmin )
	xx = amin1 ( float(xbond), xmax )
	yn = amax1 ( 0., ymin )
	yx = amin1 ( float(ybond), ymax )
	xmo = xmid
	ymo = ymid
	xmid = .5 * scrx * ( xx + xn )
	ymid = .5 * scry * ( yx + yn )
	dxb = .5 * scrx * ( xx - xn )
	dyb = .5 * scry * ( yx - yn )
c
	dxo = xmid - xmo
	dyo = ymid - ymo
	ixorig = scrx * xn
	iyorig = scry * yn
c
c	move window into frame
c
	iascr = iascr + dxo
	ibscr = ibscr + dxo
	icscr = icscr + dyo
	idscr = idscr + dyo
	xm = xm + dxo
	ym = ym + dyo
c
c	move view port into frame
c
	xconst = xconst + dxo
	yconst = yconst + dyo
c
	return
	end



	subroutine window(xwin0, xwin1, ywin0, ywin1)
c
c	this subroutine sets the dimensions of the plotting window.
c	xwin0, xwin1, ywin0, and ywin1 are the bounds in inches of
c	the plotting window. the actual plot less labels are contained
c	within this window. scaling is performed by subroutine setwin.
c
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
c
	common /world/ xmin, xmax, ymin, ymax
c
c	convert from inches to plot coordinates
c
	iascr = scrx * xwin0 + ixorig
	ibscr = scrx * xwin1 + ixorig
	icscr = scry * ywin0 + iyorig
	idscr = scry * ywin1 + iyorig
c
	dx = .5 * (ibscr - iascr)
	dy = .5 * (idscr - icscr)
	xm = .5 * (ibscr + iascr)
	ym = .5 * (idscr + icscr)
c
	call vuport ( xmin, xmax, ymin, ymax )
	return
	end


	subroutine vuport(awrld, bwrld, cwrld, dwrld)
c
c	this subroutine sets up the correspondence between the
c	world coordinates and the screen coordinates by initializing
c	the variables xslope, yslope, xconst and yconst which
c	implicitly define the affine map from world to screen
c	coordinates through the formulas
c
c		xscr = xwrld * xslope + xconst
c		yscr = ywrld * yslope + yconst
c
c
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
c
	common /world/ xmin, xmax, ymin, ymax
c
	xmin = awrld
	xmax = bwrld
	ymin = cwrld
	ymax = dwrld
c
	xslope = (ibscr - iascr) / (bwrld - awrld)
	yslope = (idscr - icscr) / (dwrld - cwrld)
c
	xconst = - xslope * awrld + float(iascr)
	yconst = - yslope * cwrld + float(icscr)
c
	return
	end


	subroutine axis(xtic, ytic, ixtit, ixlen, ixsiz, ixfmt, 
	1			    iytit, iylen, iysiz, iyfmt)
c
c  draws axes for two dimensional plots. xtic and ytic are the 
c  distances between tick marks on the x and y axes. ixtit is the
c  label for the x axis. ixlen is the length of the x label. ixsiz 
c  is the size (1 - 5) of the x label. iytit, iylen and iysiz have
c  obvious cooresponding definitions.
c
	byte ixtit(ixlen), iytit(iylen), label(9)
	integer fmt(2) , emt, num(5), form(3)
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
c
	data fmt / '(F', '9.' /
	data emt / '(E' /
	data num / '0)', '1)', '2)', '3)', '2)' /
c
c  draw window (plot) boundary
c
	call coltyp(0)
	call plot(iascr, icscr, 0)
	call plot(ibscr, icscr, 1)
	call plot(ibscr, idscr, 1)
	call plot(iascr, idscr, 1)
	call plot(iascr, icscr, 1)
c
	idx = ifix(xtic * xslope)
	idy = ifix(ytic * yslope)
	iy = icscr
	n = (ibscr - iascr) / idx - 1
ccc
c
c  The following code puts tick marks and labels on the horizontal axis 
c
	is = iascr
c
c  set up format of labels for horizontal axis
c
	iform = max0 ( 0 , ixfmt )
	form(1) = fmt(1)
	form(2) = fmt(2)
	form(3) = num( min0( 5 , iform + 1) )
	if ( iform .gt. 3 ) form(1) = emt
c
c  label tick mark on left end of horizontal axis
c
c    compute value, translate and trim leading blanks
c
	xmin = xworld( iascr )
	encode ( 9 , form, label ) xmin
	call triml(label, 9, len)
c
c    remove decimal point for integer label
c
	if ( iform .eq. 0 ) len = len - 1
c
c    determine size of label and center
c
	ixl = len * ichar(ixsiz)
	ix = is - ixl / 2
	iy = icscr - (ichar(ixsiz) * 7) / 3
c
c    write label
c
	call pltstr(ix, iy, label , len , 1 , ixsiz)
c
c  draw ticks and write labels on the interior of horizontal axis
c
	lentic = scrx / 4
c
c  if n < = 0 then there will be no tick marks drawn nor any labels written
c  on the horizontal axis
c
	if(n .le. 0)goto 20
	do 10 i = 1, n
		is = is + idx
c
c    draw tick mark
c
		call plot(is, icscr, 0)
		call plot(is, icscr + lentic, 1)
c
c    label tick mark
c
		xmin = xmin + xtic
		encode ( 9 , form, label ) xmin
		call triml(label, 9, len)
		if ( iform .eq. 0 ) len = len - 1
		ixl = len * ichar(ixsiz)
		ix = is - ixl / 2
c
		call pltstr(ix, iy, label , len , 1 , ixsiz)
c
10	continue
c
c  label tick mark on right end of horizontal axis
c
20	continue
	is = is + idx
	xmin = xmin + xtic
	encode ( 9 , form, label ) xmin
	call triml(label , 9, len)
	if ( iform .eq. 0 ) len = len - 1
	ixl = len * ichar(ixsiz)
	ix = is - ixl / 2
c
	call pltstr(ix, iy, label , len , 1 , ixsiz)
c
c	trim leading and trailing blanks
c
	do 23 i = 1, ixlen
	if ( ixtit(i) .ne. ' ' ) goto 25
23	continue
	i = ixlen
c
25	continue
c
	do 27 j = ixlen, i, -1
	if ( ixtit(j) .ne. ' ' ) goto 29
27	continue
	j = i
c
29	continue
c
c  center horizontal title and plot
c
	ixln = j + 1 - i
	ixl = ixln * ichar(ixsiz)
	ix = (iascr + ibscr - ixl) / 2
	iy = iy - (7 * ichar(ixsiz)) / 3
c
	call pltstr(ix, iy, ixtit(i), ixln, 1, ixsiz)
c
	ix = iascr
	n = ( idscr - icscr ) / idy - 1
ccc
c
c  The following code puts tick marks and labels on the vertical axis 
c
	call plot(iascr, icscr, 0)
	is = icscr
c
c  set up format of labels for vertical axis
c
	iform = max0 ( 0 , iyfmt )
	form(1) = fmt(1)
	form(2) = fmt(2)
	form(3) = num( min0( 5 , iform + 1) )
	if ( iform .gt. 3 ) form(1) = emt
c
c  label tick mark on bottom end of vertical axis
c
c    compute value, translate and trim leading blanks
c
	ymin = yworld( icscr )
	encode ( 9 , form, label ) ymin
	call triml(label, 9, len)
c
c    remove decimal point for integer label
c
	if ( iform .eq. 0 ) len = len - 1
c
c    determine size of label and center
c
	iyl = (len + 1) * ichar(iysiz )
	ix = iascr - iyl
	ixmin = ix
	ihgt = (7 * ichar(iysiz )) / 12
	iy = is - ihgt
c
c    write label
c
	call pltstr(ix, iy, label , len , 1 , iysiz)
c
c  draw ticks and write labels on the interior of vertical axis
c
	lentic = scry / 4
c
c  if n < = 0 then there will be no tick marks drawn nor any labels written
c  on the vertical axis
c
	if(n .le. 0)goto 40
	do 30 i = 1, n
		is = is + idy
c
c    draw tick mark
c
		call plot(iascr, is, 0)
		call plot(iascr + lentic, is, 1)
c
c    label tick mark
c
		ymin = ymin + ytic
		encode ( 9 , form, label ) ymin
		call triml ( label, 9, len )
		if ( iform .eq. 0 ) len = len - 1
		ix = iascr - ( len + 1 ) * ichar (iysiz)
		ixmin = min0 ( ixmin, ix )
		iy = is - ihgt
c
		call pltstr(ix, iy, label , len , 1 , iysiz)
c
30	continue
c
c  label tick mark on top end of vertical axis
c
40	continue
	is = is + idy
	ymin = ymin + ytic
	encode ( 9 , form, label ) ymin
	call triml ( label, 9, len )
	if ( iform .eq. 0 ) len = len - 1
	ix = iascr - ( len + 1 ) * ichar (iysiz)
	ixmin = min0 ( ixmin, ix )
	iy = is - ihgt
c
	call pltstr(ix, iy, label , len , 1 , iysiz)
c
c  center vertical title and plot
c
	ix = ixmin - (7 * ichar(ixsiz)) / 6
	iyl = iylen * ichar(iysiz)
	iy = (icscr + idscr - iyl) / 2
c
	call pltstr(ix, iy, iytit, iylen, 4, iysiz)
c
	return
	end


	subroutine line ( x, y, n, icol, isym, isize, inum )
	dimension x(n), y(n)
	call dashln(x, y, n, icol, isym, isize, inum, 0 )
	return
	end


	subroutine dashln(x, y, n, icol, isym, isize, inum, ityp)
c
c	this routine plots the n world coordinates, (x, y).
c	icol is the color of line which connects the points.
c	isym is the symbol to be plotted at the points.
c	isize is the size of the symbol.
c	ityp is the type of line which connects the points.
c	a symbol is plotted every abs(inum) points. if inum is
c	negative, only symbols are plotted.
c	if any of the following are out of range, only a line is 
c	drawn: isym, isize, inum
c
	byte mark
	dimension x(n), y(n)
	common/chkbon/ibond, istrt, linntp
c
	mark = isize .le. 0.or.isize .gt. 5
	mark = mark.or.isym .lt. 0.or.isym .gt. 5
	mark = mark.or.inum .eq. 0
	ibond = 0
c
	if(inum .lt. 0)goto 20
c
c	draw line
c
	call coltyp(icol)
	istrt = 1
	call moveto(x(1), y(1), 0, ityp)
c
	do 10 i = 2, n
		call moveto(x(i), y(i), 1, ityp)
10	continue
c
20	if(mark)return
c
c	draw symbols
c
	ibond = 0
	itemp = iabs(inum)
	do 30 i = 1, n, itemp
		call moveto(x(i), y(i), 0, ityp)
		if(ibond .eq. 0)call marker(isym, isize)
30	continue
c
	return
	end


	subroutine moveto (xwrld, ywrld, iud, lintyp)
c
c	moves the pen to the point with world coordinates (xwrld, ywrld)
c	draws a line only if iud = 1
c	the type of line drawn is set by the lintyp
c	move is made only to window boundary
c
	dimension ix1(2), iy1(2)
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
	common/chkbon/ibond, istrt, linntp
c
	if ( istrt .eq. 0 .and. linntp .eq. lintyp ) goto 500
	call mov1st(xwrld, ywrld, lintyp)
	return
c
500	ix = iscrx( xwrld )
	iy = iscry( ywrld )
c
c	(ix,iy) is screen coordinate inside window
c	(ixo,iyo) is current screen coordinate
c	(ixcur,iycur) is previous screen coordinate
c	(ix0,iy0) is screen coordinate outside window
c
	ixo = ix
	iyo = iy
c
c	check to see if in window
c	
	inx = (abs(ix - xm) .le. dx) 
	iny = (abs(iy - ym) .le. dy)
	in = inx * iny
c
c	in is one if inside window, zero if outside window
c	ibond is zero if previous point inside window, minus one if outside
c	in + ibond will be positive only if both previous point and current
c		point are inside window
c	in + ibond will be zero only if one point is inside and one point 
c		is outside
c	in + ibond will be negative only if both points are outside
c
	if ( in + ibond ) 300, 100, 200
ccc
c	this section finds the intersection with the window boundary for
c	consecutive points on opposite sides of the boundary
c
100	ix0 = ixcur
	iy0 = iycur
	if ( in ) goto 120
c
c	current point outside window
c
	ibond = - 1
c
c	find intersection (if any) with horizontal boundaries
c
	if ( inx ) goto 110
	xsl = (iy - iy0) / float(ix - ix0)
	xsc = xm + sign(dx, ix - xm)
	iy = xsl * (xsc - ix) + iy
	ix = xsc
c
c	find intersection (if any) with vertical boundaries
c
110	if ( abs( iy - ym ) .le. dy ) goto 200
	ysl = (ix - ixcur) / float(iy - iycur)
	ysc = ym + sign(dy, iy - ym)
	ix = ysl * (ysc - iy) + ix
	iy = ysc
	goto 200
c
c	current point inside window
c
120	ibond = 0
c
c	find intersection (if any) with horizontal boundaries
c
	if ( abs( ix0 - xm ) .le. dx ) goto 130
	xsl = (iy - iy0) / float(ix - ix0)
	xsc = xm + sign(dx, ix0 - xm)
	iy0 = xsl * (xsc - ix0) + iy0
	ix0 = xsc
c
c	find intersection (if any) with vertical boundaries
c
130	if ( abs( iy0 - ym ) .le. dy ) goto 140
	ysl = (ix - ix0) / float(iy - iy0)
	ysc = ym + sign(dy, iy0 - ym)
	ix0 = ysl * (ysc - iy0) + ix0
	iy0 = ysc
c
c	move pen to boundary
c
140	call dash(ix0, iy0, 0, lintyp)
ccc
c	this section plots points within the window
c
200	call dash(ix, iy, iud, lintyp)
c
c	save current coordinates
c
	ixcur = ixo
	iycur = iyo
c
	return
ccc
c	this section handels points outside window
c
300	insect = 0
	if ( ix .eq. ixcur ) goto 320
	if( abs(xm - dx - .5*(ix+ixcur)) .gt. .5*iabs(ix-ixcur) )goto 310
	xsl = ( iy - iycur ) / float( ix - ixcur )
	y0 = xsl * ( xm - dx - ix ) + iy
	if ( ( ( y0 - iy ) * ( y0 - iycur ) .gt. 0. ) 
	1    .or. ( abs( y0 - ym ) .gt. dy ) ) goto 310
	insect = 1
	ix1(insect) = xm - dx
	iy1(insect) = y0
c
310	if( abs(xm + dx - .5*(ix+ixcur)) .gt. .5*iabs(ix-ixcur) )goto 320
	y0 = xsl * ( xm + dx - ix ) + iy
	if ( ( ( y0 - iy ) * ( y0 - iycur ) .gt. 0. ) 
	1    .or. ( abs( y0 - ym ) .gt. dy ) ) goto 320
	insect = insect + 1
	ix1(insect) = xm + dx
	iy1(insect) = y0
	if ( insect .eq. 2 ) goto 350
c
320	if ( iy .eq. iycur ) goto 390
	if( abs(ym - dy - .5*(iy+iycur)) .gt. .5*iabs(iy-iycur) )goto 330
	ysl = ( ix - ixcur ) / float( iy - iycur )
	x0 = ysl * ( ym - dy - iy ) + ix
	if ( ( ( x0 - ix ) * ( x0 - ixcur ) .gt. 0. ) 
	1    .or. ( abs( x0 - xm ) .gt. dx ) ) goto 330
	insect = insect + 1
	ix1(insect) = x0
	iy1(insect) = ym - dy
	if ( insect .eq. 2 .and. ix1(2) .ne. ix1(1) ) goto 350
	insect = 1
c
330	if( abs(ym + dy - .5*(iy+iycur)) .gt. .5*iabs(iy-iycur) )goto 390
	x0 = ysl * ( ym + dy - iy ) + ix
	if ( ( ( x0 - ix ) * ( x0 - ixcur ) .gt. 0. ) 
	1    .or. ( abs( x0 - xm ) .gt. dx ) ) goto 390
	insect = insect + 1
	ix1(insect) = x0
	iy1(insect) = ym + dy
c
	if ( insect .lt. 2 ) goto 390
350	call dash ( ix1(1), iy1(1), 0, lintyp)
	call dash ( ix1(2), iy1(2), iud, lintyp)
390	ibond = -1
	ixcur = ixo
	iycur = iyo
	return
	end


	subroutine mov1st (xwrld, ywrld, lintyp)
c
c	moves the pen to the point with world coordinates (xwrld, ywrld)
c	draws a line only if iud = 1
c	the type of line drawn is set by the lintyp
c	move is made only to window boundary
c
	dimension ix1(2), iy1(2)
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
	common/chkbon/ibond, istrt, linntp
c
	ix = iscrx( xwrld )
	iy = iscry( ywrld )
	ibond = 0
	istrt = 1
c
c	(ix,iy) is screen coordinate 
c	check to see if in window
c	
	inx = (abs(ix - xm) .le. dx) 
	iny = (abs(iy - ym) .le. dy)
	in = inx * iny
c
c	in is one if inside window, zero if outside window
c
	if ( in ) goto 100
c
c	initial point outside of plotting window
	ibond = -1
	ixcur = ix
	iycur = iy
c
100	call dash(ix, iy, 0, lintyp)
c
c	save current coordinates
c
	ixcur = ix
	iycur = iy
	istrt = 0
c
	return
	end


	subroutine dash( ix, iy, iud, lintyp )
c
	integer lind(11)
	real idash(11), ispace(11)
	data lind / 1, 2, 3, 4, 5, 6, 9, 10, 7, 11, 8 /
	data idash/.02, .04, .08, .16, .24, .32, .16, .16, .04, .04, .04/
	data ispace/.08, .08, .08, .08, .12, .16, .08, .08, .08, .08, .08/
c
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
c
	common /chkbon/ ibond, istrt, lin
	data lin, lin1 / 0, 0 /
c
	if ( istrt .ne. 1 ) goto 50
	if ( lintyp .gt. 0 ) goto 60
	lin = 0
	goto 50
60	lin = max0 ( lintyp, 1 )
	lin = min0 ( lin, 8 )
	lin1 = lin
	dash = scrx * idash(lin1)
	space = 0
50	if ( lin .le. 0 ) goto 300
	if ( iud .ne. 1 ) goto 300
	dx1 = ix - ixcur
	dy1 = iy - iycur
	ds = sqrt ( dx1 * dx1 + dy1 * dy1 )
	if ( ds .eq. 0. ) return
	c = dx1 / ds
	s = dy1 / ds
	x0 = ixcur
	y0 = iycur
	if ( dash .eq. 0. ) goto 200
100	if ( ds .ge. dash ) goto 150
	dash = dash - ds
	ds = 0.
	call plot ( ix, iy, 1 )
	return
150	ds = ds - dash
	x0 = dash * c + x0
	y0 = dash * s + y0
	ix0 = x0 +.5
	iy0 = y0 +.5
	call plot ( ix0, iy0, 1 )
	dash = 0.
	space = scrx * ispace(lin1)
	if ( space .eq. 0. ) goto 270
200	if ( ds .ge. space ) goto 250
	space = space - ds
	ds = 0.
	call plot ( ix, iy, 0 )
	return
250	ds = ds - space
	x0 = space * c + x0
	y0 = space * s + y0
	ix0 = x0 +.5
	iy0 = y0 +.5
	call plot ( ix0, iy0, 0 )
	space = 0.
	lin1 = lind(lin1)
270	dash = scrx * idash(lin1)
	goto 100
300	call plot ( ix, iy, iud )
	return
	end


	subroutine pltstr( ix, iy, str, len, irot, isize )
c
c	this routine moves to the point (ix,iy) and writes a string
c	only if the entire string is in the boundary
c
	byte str(len)
c
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
c
c	attempt move to point and write string
c
	call plot ( ix, iy, 0 )
	if ( (ix .eq. ixcur) .and. (iy .eq. iycur) )
	1	call wrtstr ( str, len, irot, isize )
	return
	end


	subroutine triml( label, n, len )
	byte label(n)
c
	do 10 i = 1, n
	if ( label(i) .ne. ' ' ) goto 20
10	continue
c
20	len = 0
	do 30 j = i, n
	len = len + 1
	label(len) = label(j)
30	continue
c
	return
	end


	function iscrx( xwrld )
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
c
	x = amin1( xwrld * xslope + xconst, 32767. )
	x = amax1 ( x, -32767. )
	iscrx = ifix( x )
c
	return
	end


	function iscry( ywrld )
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
c
	y = amin1 ( ywrld * yslope + yconst, 32767. )
	y = amax1 ( y, -32767. )
	iscry = ifix ( y )
c
	return
	end


	function xworld( ix )
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
c
	xworld = ( ix - xconst ) / xslope
c
	return
	end


	function yworld( iy )
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
c
	yworld = ( iy - yconst ) / yslope
c
	return
	end


	subroutine curpos( icurx, icury )
	common /zgraph/ lunplt, iascr, ibscr, icscr, idscr, 
	1 xm, ym, dx, dy, dxb, dyb, ixcur, iycur, xmid, ymid, 
	1 xslope, yslope, xconst, yconst, scrx, scry, ichar(5),
	1 ixorig, iyorig
c
	icurx = ixcur
	icury = iycur
	return
	end