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 /hwindo/ xscr,yscr,ixorig,iyorig,iascr,ibscr,icscr,idscr 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 scrx = xscr scry = yscr 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 call window( 0., xw, 0., yw ) call vuport( 0., scrx * xw, 0., scry * yw ) 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 moveto(aint(x(1,j)),aint(y(1,j)),0,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 moveto(aint(x(ii,j)),aint(y(ii,j)),0,0) 340 call moveto(aint(x(ii,j)),aint(aa*x(ii,j)+bb),1,0) call moveto(float(ix),float(iy),1,0) 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 moveto(aint(x(i,1)),aint(y(i,1)),0,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 moveto(aint(x(i,jj)),aint(y(i,jj)),0,0) 440 call moveto(aint(x(i,jj)),aint(aa*x(i,jj)+bb),1,0) call moveto(float(ix),float(iy),1,0) return end subroutine insect(i1,j1,iplot,iplot0,x,y,nx,ny,l) c dimension x(nx,ny),y(nx,ny) byte eq1,eq2 common /hwindo/ scrx,scry,ixorig,iyorig,iascr,ibscr,icscr,idscr 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 moveto(aint(xl),aint(yl),0,0) return 280 ix=ifix(xf) iy=ifix(yf) call moveto(float(ix),float(iy),iplot0,0) if(xf.eq.xl)return 290 if(iplot.eq.0)return ix=ifix(x(i1,j1)) iy=ifix(y(i1,j1)) call moveto(float(ix),float(iy),1,0) return end