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