subroutine bv(x,n,ab,nbin,nc,nf,face,out,idir,tx,ty,eye,dist,iptr,square)
integer nf(3), nbin(2), nc(nbin(1),nbin(2)), idir(5*nbin(1)*nbin(2))
integer skip(5), iptr(nf(3))
dimension x(n,2), ab(2,2), tx(nbin(1)+1), ty(nbin(2)+1)
dimension face(nf(1),nf(2),nf(3)), out(nf(1),nf(2),nf(3))
dimension eye(3), v1(3), v2(3), v3(3), vdir(5,3), dist(nf(3))
dimension square(4,nf(3))

##### S function to plot bivariate histogram perspective view
##### Copyright 1988 David W. Scott	7/2/88
##### New S version 1/6/91

#		x	data matrix -  n by 2
#		ab	bin support rectangle (half-open)
#				x-axis  [ ab[1,1] , ab[1,2] )
#				y-axis	[ ab[2,1] , ab[2,2] )
#		nbin	number of x and y axis bins
#		opt	optional message control T=suppress F=default

#		nc	bin count matrix (see "persp" function)
#		nskip	number of points outside binning rectange

call realpr("x axis binned over interval",27,ab(1,1),1)
call realpr("x axis binned over interval",27,ab(1,2),1)
call realpr("y axis binned over interval",27,ab(2,1),1)
call realpr("y axis binned over interval",27,ab(2,2),1)

#### bin data

lx=nbin(1);  ly=nbin(2);  
ncmax=0;  nskip=0;  
ax = ab(1,1); bx = ab(1,2); ay = ab(2,1); by = ab(2,2)
do j=1,ly{ do i=1,lx { nc(i,j)=0 }}
dx = (bx-ax) / lx; dy = (by-ay) / ly
do i=1,n {
	kx = (x(i,1)-ax) / dx + 1.0
	ky = (x(i,2)-ay) / dy + 1.0
	if (kx>=1 & kx<=lx & ky>=1 & ky<=ly) nc(kx,ky) = nc(kx,ky)+1
		else nskip = nskip+1
	}
do ky=1,ly { do kx=1,lx { ncmax = max1(ncmax,nc(kx,ky)) }}
call intpr("max bin count=",14,ncmax,1)

if( nskip>0 ) call intpr("points not in open rectangle ab",31,nskip,1)

##### set up data structures for histogram faces ( cube (-1,1)^3 ) 
#		from (-1,1)
dx=2./lx; dy=2./ly
do kx=1,lx+1{ tx(kx)=-1.+dx*(kx-1)}; do ky=1,ly+1{ ty(ky)=-1.+dy*(ky-1)}

ic = 0	# pointer to face
do ky=1,ly { do kx=1,lx { ht=2*nc(kx,ky)/float(ncmax)-1; mx=kx+1; my=ky+1
	ic=ic+1; idir(ic)=1
		face(1,1,ic)=tx(kx); face(1,2,ic)=ty(ky); face(1,3,ic)=ht
		face(2,1,ic)=tx(mx); face(2,2,ic)=ty(ky); face(2,3,ic)=ht
		face(3,1,ic)=tx(mx); face(3,2,ic)=ty(my); face(3,3,ic)=ht
		face(4,1,ic)=tx(kx); face(4,2,ic)=ty(my); face(4,3,ic)=ht
		face(5,1,ic)=tx(kx); face(5,2,ic)=ty(ky); face(5,3,ic)=ht
   if(nc(kx,ky)>0) {	#	NEED indicator which case
	ic=ic+1; idir(ic)=2
		face(1,1,ic)=tx(kx); face(1,2,ic)=ty(ky); face(1,3,ic)=-1.
		face(2,1,ic)=tx(mx); face(2,2,ic)=ty(ky); face(2,3,ic)=-1.
		face(3,1,ic)=tx(mx); face(3,2,ic)=ty(ky); face(3,3,ic)=ht
		face(4,1,ic)=tx(kx); face(4,2,ic)=ty(ky); face(4,3,ic)=ht
		face(5,1,ic)=tx(kx); face(5,2,ic)=ty(ky); face(5,3,ic)=-1.
	ic=ic+1; idir(ic)=3
		face(1,1,ic)=tx(mx); face(1,2,ic)=ty(ky); face(1,3,ic)=-1.
		face(2,1,ic)=tx(mx); face(2,2,ic)=ty(my); face(2,3,ic)=-1.
		face(3,1,ic)=tx(mx); face(3,2,ic)=ty(my); face(3,3,ic)=ht
		face(4,1,ic)=tx(mx); face(4,2,ic)=ty(ky); face(4,3,ic)=ht
		face(5,1,ic)=tx(mx); face(5,2,ic)=ty(ky); face(5,3,ic)=-1.
	ic=ic+1; idir(ic)=4
		face(1,1,ic)=tx(mx); face(1,2,ic)=ty(my); face(1,3,ic)=-1.
		face(2,1,ic)=tx(kx); face(2,2,ic)=ty(my); face(2,3,ic)=-1.
		face(3,1,ic)=tx(kx); face(3,2,ic)=ty(my); face(3,3,ic)=ht
		face(4,1,ic)=tx(mx); face(4,2,ic)=ty(my); face(4,3,ic)=ht
		face(5,1,ic)=tx(mx); face(5,2,ic)=ty(my); face(5,3,ic)=-1.
	ic=ic+1; idir(ic)=5
		face(1,1,ic)=tx(kx); face(1,2,ic)=ty(my); face(1,3,ic)=-1.
		face(2,1,ic)=tx(kx); face(2,2,ic)=ty(ky); face(2,3,ic)=-1.
		face(3,1,ic)=tx(kx); face(3,2,ic)=ty(ky); face(3,3,ic)=ht
		face(4,1,ic)=tx(kx); face(4,2,ic)=ty(my); face(4,3,ic)=ht
		face(5,1,ic)=tx(kx); face(5,2,ic)=ty(my); face(5,3,ic)=-1.
		}
	}
}
icmax = ic
call intpr("number of faces",16,icmax,1)
call intpr("dimensioned",11,nf(3),1)

# outward normals of 5 histogram bin faces
	vdir(1,1)= 0; vdir(1,2)= 0; vdir(1,3)= 1	# 1 top
	vdir(2,1)= 0; vdir(2,2)=-1; vdir(2,3)= 0	# 2 front
	vdir(3,1)= 1; vdir(3,2)= 0; vdir(3,3)= 0	# 3 right
	vdir(4,1)= 0; vdir(4,2)= 1; vdir(4,3)= 0	# 4 back
	vdir(5,1)=-1; vdir(5,2)= 0; vdir(5,3)= 0	# 5 left

e1 = eye(1); e2 = eye(2); e3 = eye(3)
rho = e1**2 + e2**2 + e3**2; rho=sqrt(rho)
if(rho.le.0.) {call realpr("eye=0 quitting",14,eye,3); goto 102}
if(e3<=0.){call realpr("eye3 zaxis must be positive");goto 102}

phi = acos(e3/rho); pi = 3.14159265
if( e1!=0. | e2!=0. ) theta = atan2(eye(2),eye(1)) # atan2 handles special cases
if( e1==0. & e2==0. & e3>0. ) theta=-pi/2
if( e1==0. & e2==0. & e3<0. ) theta= pi/2  # don't have RATFOR manual
call realpr("rho",3,rho,1)
call realpr("phi",3,phi,1)
call realpr("theta",5,theta,1)
sp=sin(phi);  cp=cos(phi);  st=sin(theta);  ct=cos(theta)
v1(1)=-st;    v1(2)=ct;     v1(3)=0.0
v2(1)=-ct*cp; v2(2)=-st*cp; v2(3)=sp
v3(1)=ct*sp;  v3(2)=st*sp;  v3(3)=cp  #right hand system (not graphics LHS)
call realpr("rv1=",4,v1,3)
call realpr("rv2=",4,v2,3)
call realpr("rv3=",4,v3,3)

do i=1,3 {eye(i)=eye(i)/rho};  idskip=0
do i=1,5 {tot=0.; do j=1,3 {tot=tot+eye(j)*vdir(i,j)}	# cosine between vectors
	if (tot>0.) skip(i)=0		# between -pi/2 and pi/2
	else {skip(i)=1; idskip=idskip+1}
}
call intpr("number of faces skipped",23,nskip,1)

# project rectangles by eye

jcmax = 0	# number of faces "visible"
do jz=1,icmax{					# icmax <= nf(3)
   if(skip(idir(jz))==0) { jcmax=jcmax+1;  iptr(jcmax)=jcmax
	do jx=1,5 {				# 5 == nf(1)
		tot=0.; do j=1,3 {tot=tot+v1(j)*face(jx,j,jz)}; out(jx,1,jcmax)=tot
		tot=0.; do j=1,3 {tot=tot+v2(j)*face(jx,j,jz)}; out(jx,2,jcmax)=tot
		tot=0.; do j=1,3 {tot=tot+v3(j)*face(jx,j,jz)}; out(jx,3,jcmax)=tot
	}
	if (idir(jz)==1) {	# compute closest point as if face has no height (ie z=-1)
		di=-1.e9
		do jx=1,4{ tot=-v3(3); do j=1,2 {tot=tot+v3(j)*face(jx,j,jz)}; di=amax1(di,tot)}
		dist(jcmax) = -di }
	else {			# case of a side --- just first coordinates
		di=-1.e9
		do jx=1,2{ tot=-v3(3); do j=1,2 {tot=tot+v3(j)*face(jx,j,jz)}; di=amax1(di,tot)}
		dist(jcmax) = -di }
   }
}

# sort by closest point on each face's bottom

call sorts ( dist, iptr, 1, jcmax )	# iptr = pointer to original array front-to-back
					#   notice -di so that closest to eye is smallest

# smallest square covering each face (xl,xr,yl,yr)

do jz=1,jcmax{
	di=0.; do jx=1,4{di=amin1(di,out(jx,1,jz))}; square(1,jz)= di
	di=0.; do jx=1,4{di=amax1(di,out(jx,1,jz))}; square(2,jz)= di
	di=0.; do jx=1,4{di=amin1(di,out(jx,2,jz))}; square(3,jz)= di
	di=0.; do jx=1,4{di=amax1(di,out(jx,2,jz))}; square(4,jz)= di
}

# plot histogram pieces

do jz=1,jcmax{
	do jx=1,4{
		call inter( jz, jx, iptr, out, square)
	}
}

102	return
end

#------------------- internew.r ---------------------

subroutine inter( jz,jx,iptr,out,sq )
dimension out(5,3,1),sq(4,1),xy(4,20),xout(2),yout(2)
integer iptr(1),status(20)

# working on side "jx" on sorted face "jz"
#   check part of side jx visible behind faces 1,...,jz-1
#  xy = (x1,x2,y1,y2)  (same as sq==square)

lenxy=1;  status(1)=1		# visible part of segment(s)

ip=iptr(jz)
xy(1,1)=out(jx,1,ip); xy(2,1)=out(jx+1,1,ip); xy(3,1)=out(jx,2,ip); xy(4,1)=out(jx+1,2,ip)

do ii=1,jz-1 {		## go through list of closer faces (excluding current one)
   ip=iptr(ii)
			##  compute intersections
   jmax = lenxy

   do j=1,jmax{		## loop on all current subsegments of this side
	if(status(j)==0) goto 101	# skip (instead of popping stack)

#--check if edge not covered by face

	if( amax1(xy(1,j),xy(2,j)) < sq(1,ip) | amin1(xy(1,j),xy(2,j)) > sq(2,ip) |
	    amax1(xy(3,j),xy(4,j)) < sq(3,ip) | amin1(xy(3,j),xy(4,j)) > sq(4,ip) ) {
		goto 101 }	# disjoint so continue search

#--check if edge completely covered by face (later do all segments)
	is1 = insider( xy(1,j),xy(3,j),ip,out)
	is2 = insider( xy(2,j),xy(4,j),ip,out)

	if(is1+is2==2) {
		status(j) = 0		# endpoints covered --- invisible segment --- delete
		}

	if(is1==1 & is2==0) {	# shorten segment by replacing xy(1 and 3) since covered
		ind = 1		# want only one (preferably nonzero) cut
		call cuts( xy(1,j), ip, out, xout, yout, ind)
		if(ind==1) {xy(1,j) = xout(1);  xy(3,j) = yout(1)}
		}

	if(is1==0 & is2==1) {	# shorten segment by replacing xy(2 and 4) since covered
				### in either case   status==1 still
		ind = 1		# want only one (preferably nonzero) cut
		call cuts( xy(1,j), ip, out, xout, yout, ind)
		if(ind==1) {xy(2,j) = xout(1);  xy(4,j) = yout(1)}
		}

	if(is1+is2==0) {	# either two segments or entirely not covered (ie don't change)
		### if add segment   lenxy=lenxy+1;  status(lenxy)=1
		ind = 2		# want two (preferably nonzero) cuts  --- or none
		call cuts( xy(1,j), ip, out, xout, yout, ind)
			# ind=0 usually  but ind=1 if one end on face edge

		if(ind==2) {
			lenxy = lenxy+1; if(lenxy>20) call intpr("???opppss???????lenxy",21,lenxy,1)
			xy(1,lenxy) = xout(2); xy(2,lenxy) = xy(2,j)
			xy(3,lenxy) = yout(2); xy(4,lenxy) = xy(4,j); status(lenxy) = 1
			xy(2,j)     = xout(1); xy(4,j)     = yout(1)
			}
		}

	if(is1+is2<0 | is1+is2>2) {
			call intpr("*******iside error",18,is1,1)
			call intpr("*******iside error",18,is2,1) }

101	continue
   }
}

do j=1,lenxy { if(status(j)==1) call segmtz( xy(1,j),xy(3,j),xy(2,j),xy(4,j),1 ) }

return
end

#------------------cutsnew.r-------------------------

subroutine cuts	( xy, ip, out, xout, yout, ind)
dimension xy(4),out(5,3,1),xout(2),yout(2)	# just a column of xy
dimension iptr(10),xl1(10)	# a large margin of error (should only need 2)

# find intersection of edge with 4 face edges
# if more than one, order by xl1

#	ind=1 ==> try to find one nonzero cut
#	ind=2 ==> try to find two nonzero cuts
#		put 0 or 1 in as last resort

ixl1=0;  do i=1,10{iptr(i)=i}
e1=xy(1); e2=xy(2); f1=xy(3); f2=xy(4)

do i=1,4{		# check four sides
	a1=out(i,1,ip);  a2=out(i+1,1,ip)
	b1=out(i,2,ip);  b2=out(i+1,2,ip)

	det = (e2-e1)*(b1-b2) - (f2-f1)*(a1-a2)

	if (det!=0.) {		# != should work here
		x1 = ( (a1-e1)*(b1-b2) - (b1-f1)*(a1-a2) ) / det
		x2 = ( (e2-e1)*(b1-f1) - (f2-f1)*(a1-e1) ) / det
		if (x1>=0. & x1<=1. & x2>=0. & x2<=1.) {
			ixl1=ixl1+1; xl1(ixl1)=x1
			}
		}
	else {		# parallel --- check if coincident
		if(e1==e2) {		# vertical parallel lines
			if (abs(a1-e1) > 1.e-6) goto 101	# not coincident a1!=e1
			if ((b1-f1)*(b1-f2)<=0.){
				ixl1=ixl1+1; xl1(ixl1)=(b1-f1)/(f2-f1)}
			if ((b2-f1)*(b2-f2)<=0.){
				ixl1=ixl1+1; xl1(ixl1)=(b2-f1)/(f2-f1)}
			}
		if(f1==f2) {		# horizontal parallel lines
			if (abs(b1-f1) > 1.e-6) goto 101	# not coincident b1!=f1
			if((a1-e1)*(a1-e2)<=0.){
				ixl1=ixl1+1; xl1(ixl1)=(a1-e1)/(e2-e1)}
			if((a2-e1)*(a2-e2)<=0.){
				ixl1=ixl1+1; xl1(ixl1)=(a2-e1)/(e2-e1)}
			}
		if(e1!=e2 & f1!=f2) {	# oblique parallel lines    != works here
			c1 = (a1-e1)/(e2-e1)
			if( abs( c1-(b1-f1)/(f2-f1) )<1.e-6 & c1>=0. & c1<=1.)
				{ixl1=ixl1+1; xl1(ixl1)=c1}
			c2 = (a2-e1)/(e2-e1)
			if( abs( c2-(b2-f1)/(f2-f1) )<1.e-6 & c2>=0. & c2<=1.)
				{ixl1=ixl1+1; xl1(ixl1)=c2}
			}
		}
101	continue
}

# wrapup

if(ixl1==0) { ind=0; return }

call sorts( xl1, iptr, 1, ixl1 )		# sort cut values

iup=ixl1-1
do i=iup,1,-1 {		# get rid of duplicates (0's and 1's)
	if( abs( xl1(i+1)-xl1(i)) < 1.e-7 ) {	# == doesn't work in middle
		ixl1=ixl1-1
		do j=i,iup { xl1(j)=xl1(j+1) }
	}
}

k=iup+1-ixl1

if(ind==ixl1) goto 201		# got 1 or 2 as desired

if(ind==1){	# find cuts closest to 1/2
	di=1.e9; do i=1,ixl1{ dt=abs(xl1(i)-.5); if(dt<di) {k=i;di=dt} }
	xl1(1)=xl1(k)		# put in first spot
	}

if(ind==2 & ixl1==1) {ind=1; goto 201}	# on edge I suppose

if(ind==2){	# find two cuts closest to 1/2 (ixl1 is > 2)
	k=0
	do i=1,ixl1{		# find nonzero nonone cuts
		dt=xl1(i);  if(dt>0.&dt<1.) {k=k+1; iptr(k)=i}
	}
	if(k==1){	# return either (0,dt) or (dt,1)
		if(xl1(1)==0.) { xl1(2)=xl1(iptr(1)) }	# == should work here
		else { xl1(1)=xl1(iptr(1)); xl1(2)=1. }
		}
	if(k==2){	# put in positions 1 and 2
		xl1(1)=xl1(iptr(1)); xl1(2)=xl1(iptr(2))
		}
	}

201	continue
do i=1,ind{ xout(i) = e1 + xl1(i)*(e2-e1); yout(i) = f1 + xl1(i)*(f2-f1) }

return
end

#------------------------insidernew.r---------------------

integer function insider(x,y,ip,out)
dimension out(5,3,1)

#			draw vertical line at x

fact=1
insider=0	# 0 ==> no     1 ==> yes    2 ==> on edge
icount=0	# used to count number of face intersections with vertical line thru point

do i=1,4{		# check four sides
	x1=out(i,1,ip);  x2=out(i+1,1,ip)
	y1=out(i,2,ip);  y2=out(i+1,2,ip)
	if( x1==x2 ) {			# edge is vertical
		if ( abs(x-x1) < 1.e-6 ) {	# x==x1 point in same vertical plane
			if( (y-y1)*(y-y2) <= 0. ) {
				insider=1	###### or 2 at corner ######
				return
				}	# else no intersection
			}
		}

	else {				# x1 != x2  edge not vertical
		if( (x-x1)*(x-x2) <= 0. ) {
			xl2 = (x1-x)/(x1-x2)		# lambda intersection factors
			xl1 = (y1-y)-xl2*(y1-y2)	# intersection with vertical line

			if (abs(xl1) < 1.e-6) {		# xl1 == 0. approximation
				insider=1	#   point on face edge
				return
				}

			if( xl2 != 0. ) {
				icount=icount+1
				fact = fact*xl1		# trying to avoid having
				}			# same fact twice at corner
							# once with xl2=1 then xl2=0

			}		# should happen twice
		}			#   if insider, xl1(1)*xl1(2) will be negative (or zero)
}
if(icount==0) return
if(icount==1) {			# point in line through extreme corner of face
	insider=0		# special corner condition
	return}
if(fact	< 0.) {			# == 0. handled above:  strictly inside
	insider=1
	return}

return
end

      subroutine sorts (v,a,ii,jj)
      dimension a(jj),v(1),iu(20),il(20)
      integer t,tt
      integer a
      real v
      m=1
      i=ii
      j=jj
 10   if (i.ge.j) go to 80
 20   k=i
      ij=(j+i)/2
      t=a(ij)
      vt=v(ij)
      if (v(i).le.vt) go to 30
      a(ij)=a(i)
      a(i)=t
      t=a(ij)
      v(ij)=v(i)
      v(i)=vt
      vt=v(ij)
 30   l=j
      if (v(j).ge.vt) go to 50
      a(ij)=a(j)
      a(j)=t
      t=a(ij)
      v(ij)=v(j)
      v(j)=vt
      vt=v(ij)
      if (v(i).le.vt) go to 50
      a(ij)=a(i)
      a(i)=t
      t=a(ij)
      v(ij)=v(i)
      v(i)=vt
      vt=v(ij)
      go to 50
 40   a(l)=a(k)
      a(k)=tt
      v(l)=v(k)
      v(k)=vtt
 50   l=l-1
      if (v(l).gt.vt) go to 50
      tt=a(l)
      vtt=v(l)
 60   k=k+1
      if (v(k).lt.vt) go to 60
      if (k.le.l) go to 40
      if (l-i.le.j-k) go to 70
      il(m)=i
      iu(m)=l
      i=k
      m=m+1
      go to 90
 70   il(m)=k
      iu(m)=j
      j=l
      m=m+1
      go to 90
 80   m=m-1
      if (m.eq.0) return
      i=il(m)
      j=iu(m)
 90   if (j-i.gt.10) go to 20
      if (i.eq.ii) go to 10
      i=i-1
 100  i=i+1
      if (i.eq.j) go to 80
      t=a(i+1)
      vt=v(i+1)
      if (v(i).le.vt) go to 100
      k=i
 110  a(k+1)=a(k)
      v(k+1)=v(k)
      k=k-1
      if (vt.lt.v(k)) go to 110
      a(k+1)=t
      v(k+1)=vt
      go to 100
      end