[comp.sources.amiga] v89i091: plplot - scientific plotting library, Part06/07

page%rishathra@Sun.COM (Bob Page) (04/25/89)

Submitted-by: amr@dukee.egr.duke.edu (Tony Richardson)
Posting-number: Volume 89, Issue 91
Archive-name: applications/plplot

# This is a shell archive.
# Remove anything above and including the cut line.
# Then run the rest of the file through 'sh'.
# Unpacked files will be owned by you and have default permissions.
#----cut here-----cut here-----cut here-----cut here----#
#!/bin/sh
# shar: SHell ARchive
# Run the following text through 'sh' to create:
#	plbeg.c
#	plbin.c
#	plbox.c
#	plbox3.c
#	plccal.c
#	plclr.c
#	plcntr.c
#	plcol.c
#	plcont.c
#	plcvec.c
#	pldeco.c
#	pldtik.c
#	plend.c
#	plenv.c
#	plerrx.c
#	plerry.c
#	plerx1.c
#	plery1.c
#	plfont.c
#	plform.c
#	plgra.c
#	plgrid3.c
#	plgspa.c
#	plhist.c
#	plhrsh.c
#	pljoin.c
#	pllab.c
#	pllclp.c
#	plline.c
#	plmtex.c
#	plnxtv.c
# This is archive 6 of a 7-part kit.
# This archive created: Thu Apr 20 13:47:05 1989
echo "extracting plbeg.c"
sed 's/^X//' << \SHAR_EOF > plbeg.c
X/* Sets up the device "dev" for plotting, dividing the page into "nx" */
X/* by "ny" subpages. */
X
X#include "plplot.h"
X#include <math.h>
X
Xvoid plbeg(dev,nx,ny)      
Xint dev, nx, ny;
X{
X      float scale, def, ht;
X
X      if ((nx <= 0) || (ny <= 0 )) 
X         fatal("Cannot have negative number of subpages in PLSTAR");
X
X      scale = 1.0/sqrt((double)ny);
X      grbeg(dev); 
X      ssub(nx,ny,0);
X
X      /* Set up character, symbol and tick sizes for requested number of */
X      /* subpages */
X
X      gchr(&def,&ht);
X      schr(def*scale,def*scale);
X      gsym(&def,&ht);
X      ssym(def*scale,def*scale);
X      gmaj(&def,&ht);
X      smaj(def*scale,def*scale);
X      gmin(&def,&ht);
X      smin(def*scale,def*scale);
X      slev(1);
X}
SHAR_EOF
echo "extracting plbin.c"
sed 's/^X//' << \SHAR_EOF > plbin.c
X/* Plot a histogram using the arrays x and y to represent data       */
X/* values and frequencies respectively. If centre is false, x values     */
X/* denote the lower edge of the bin, and if centre is true, they        */
X/* they denote the centre of the bin                                */
X
X#include "plplot.h"
X#include <math.h>
X
Xvoid plbin(nbin,x,y,centre)
Xint nbin;
Xint centre;
Xfloat x[], y[];
X{
X
X      int i;
X      float xmin, xmax, vpwxmi, vpwxma, vpwymi, vpwyma;
X
X      int level;
X      glev(&level);
X      if (level<3)fatal("Please set up window before calling PLBIN.");
X
X/*    Check x[i] are in ascending order */
X
X      for(i=0;i<nbin-1;i++)
X        if (x[i] >= x[i+1])
X           fatal("Elements of X(*) must be increasing in PLBIN.");
X
X      gvpw(&vpwxmi,&vpwxma,&vpwymi,&vpwyma);
X      if (!centre) {
X        for(i=0;i<nbin;i++){
X          pljoin(x[i],vpwymi,x[i],y[i]);
X          pljoin(x[i],y[i],x[i+1],y[i]);
X          pljoin(x[i+1],y[i],x[i+1],vpwymi);
X        }
X        if (x[nbin-1] < vpwxma) {
X          pljoin(x[nbin-1],vpwymi,x[nbin-1],y[nbin-1]);
X          pljoin(x[nbin-1],y[nbin-1],vpwxma,y[nbin-1]);
X          pljoin(vpwxma,y[nbin-1],vpwxma,vpwymi);
X        }
X      }
X      else {
X        if (nbin < 2) return;
X        xmin = vpwxmi;
X        xmax = max(0.5*(x[1]+x[2]),vpwxmi);
X        if (xmin < xmax) {
X          pljoin(xmin,vpwymi,xmin,y[1]);
X          pljoin(xmin,y[1],xmax,y[1]);
X          pljoin(xmax,y[1],xmax,vpwymi);
X        }
X        for (i=1;i<nbin-1;i++) {
X          xmin = xmax;
X          xmax = min(0.5*(x[i]+x[i+1]),vpwxma);
X          pljoin(xmin,vpwymi,xmin,y[i]);
X          pljoin(xmin,y[i],xmax,y[i]);
X          pljoin(xmax,y[i],xmax,vpwymi);
X        }
X        xmin = xmax;
X        xmax = vpwxma;
X        if (xmin < xmax) {
X          pljoin(xmin,vpwymi,xmin,y[nbin]);
X          pljoin(xmin,y[nbin],xmax,y[nbin]);
X          pljoin(xmax,y[nbin],xmax,vpwymi);
X        }
X      }
X}
SHAR_EOF
echo "extracting plbox.c"
sed 's/^X//' << \SHAR_EOF > plbox.c
X/* This draws a box around the current viewport. XOPT and YOPT are */
X/* character strings which define the box as follows: */
X
X/* A: Draw axis (X is horizontal line Y=0, Y is vertical line X=0) */
X/* B: Draw bottom (X) or left (Y) edge of frame */
X/* C: Draw top (X) or right (Y) edge of frame */
X/* G: Draws a grid at the major tick interval */
X/* I: Inverts tick marks */
X/* L: Logarithmic axes, major ticks at decades, minor ticks at units */
X/* N: Write numeric label at conventional location */
X/* M: Write numeric label at unconventional location */
X/* T: Draw major tick marks */
X/* S: Draw minor tick marks */
X/* V: (for Y only) Label vertically */
X
X/* xtick, ytick are the major tick intervals required, zero for  */
X/*              automatic selection */
X/* nxsub, nysub are the number of subtick intervals in a major tick */
X/*              interval */
X
X#include "plplot.h"
X#include <stdio.h>
X#include <math.h>
X#include <string.h>
X
X#define betw(c,a,b) ((a<=c && c<=b) || (b<=c && c<=a))
X
Xstatic float xlog[8] = {0.301030, 0.477121, 0.602060, 0.698970,
X                        0.778151, 0.845098, 0.903090, 0.954243};
X
Xvoid plbox(xopt,xtick,nxsub,yopt,ytick,nysub)
Xchar xopt[],yopt[];
Xfloat xtick, ytick;
Xint nxsub, nysub;
X{
X      char string[40];
X      char strtmp[4];
X      int lax,lbx,lcx,lgx,lix,llx,lmx,lnx,lsx,ltx;
X      int lay,lby,lcy,lgy,liy,lly,lmy,lny,lsy,lty,lvy;
X      int xmajor, xminor, ymajor, yminor, xmode, xprec;
X      int ymode, yprec;
X      int i, i1x, i2x, i3x, i4x, i1y, i2y, i3y, i4y, it0;
X      int nxsub1, nysub1;
X      int lxmin, lxmax, lymin, lymax;
X      int pxmin, pxmax, pymin, pymax;
X      int vppxmi, vppxma, vppymi, vppyma;
X      int level;
X      float xpmm, ypmm, defmaj, defmin, htmaj, htmin;
X      float xtick1, ytick1, vpwxmi, vpwxma, vpwymi, vpwyma;
X      float pos, tn, tp, temp;
X
X      glev(&level);
X      if (level<3) fatal("Please set up window before calling PLBOX.");
X
X/* Open  the clip limits to the subpage limits */
X
X      gclp(&lxmin,&lxmax,&lymin,&lymax);
X      gphy(&pxmin,&pxmax,&pymin,&pymax);
X      sclp(pxmin,pxmax,pymin,pymax);
X
X      gvpp(&vppxmi,&vppxma,&vppymi,&vppyma);
X
X/* Tick and subtick sizes in device coords */
X
X      gpixmm(&xpmm,&ypmm);
X      gmaj(&defmaj,&htmaj);
X      gmin(&defmin,&htmin);
X      
X      xmajor=max(round(htmaj*ypmm),1);
X      ymajor=max(round(htmaj*xpmm),1);
X      xminor=max(round(htmin*ypmm),1);
X      yminor=max(round(htmin*xpmm),1);
X
X      xtick1=xtick;
X      nxsub1=nxsub;
X      ytick1=ytick;
X      nysub1=nysub;
X
X      lax=strpos(xopt,'A')!=-1 || strpos(xopt,'a')!=-1;
X      lbx=strpos(xopt,'B')!=-1 || strpos(xopt,'b')!=-1;
X      lcx=strpos(xopt,'C')!=-1 || strpos(xopt,'c')!=-1;
X      lgx=strpos(xopt,'G')!=-1 || strpos(xopt,'g')!=-1;
X      lix=strpos(xopt,'I')!=-1 || strpos(xopt,'i')!=-1;
X      llx=strpos(xopt,'L')!=-1 || strpos(xopt,'l')!=-1;
X      lmx=strpos(xopt,'M')!=-1 || strpos(xopt,'m')!=-1;
X      lnx=strpos(xopt,'N')!=-1 || strpos(xopt,'n')!=-1;
X      lsx=strpos(xopt,'S')!=-1 || strpos(xopt,'s')!=-1;
X      ltx=strpos(xopt,'T')!=-1 || strpos(xopt,'t')!=-1;
X
X      lay=strpos(yopt,'A')!=-1 || strpos(yopt,'a')!=-1;
X      lby=strpos(yopt,'B')!=-1 || strpos(yopt,'b')!=-1;
X      lcy=strpos(yopt,'C')!=-1 || strpos(yopt,'c')!=-1;
X      lgy=strpos(yopt,'G')!=-1 || strpos(yopt,'g')!=-1;
X      liy=strpos(yopt,'I')!=-1 || strpos(yopt,'i')!=-1;
X      lly=strpos(yopt,'L')!=-1 || strpos(yopt,'l')!=-1;
X      lmy=strpos(yopt,'M')!=-1 || strpos(yopt,'m')!=-1;
X      lny=strpos(yopt,'N')!=-1 || strpos(yopt,'n')!=-1;
X      lsy=strpos(yopt,'S')!=-1 || strpos(yopt,'s')!=-1;
X      lty=strpos(yopt,'T')!=-1 || strpos(yopt,'t')!=-1;
X      lvy=strpos(yopt,'V')!=-1 || strpos(yopt,'v')!=-1;
X
X      gvpw(&vpwxmi,&vpwxma,&vpwymi,&vpwyma);
X      lax=lax && (vpwymi*vpwyma<0.0) && !llx;
X      lay=lay && (vpwxmi*vpwxma<0.0) && !lly;
X      if (llx) xtick1=1.0;
X      if (lly) ytick1=1.0;
X      if (ltx || lgx) 
X          pldtik(vpwxmi,vpwxma,&xtick1,&nxsub1,&xmode,&xprec);
X      if (lty || lgy) 
X          pldtik(vpwymi,vpwyma,&ytick1,&nysub1,&ymode,&yprec);
X
X/* Set up tick variables */
X
X      if (lix) {
X        i1x=xminor;
X        i2x=0;
X        i3x=xmajor;
X        i4x=0;
X      }
X      else {
X        i1x=0;
X        i2x=xminor;
X        i3x=0;
X        i4x=xmajor;
X      }
X
X      if (liy) {
X        i1y=yminor;
X        i2y=0;
X        i3y=ymajor;
X        i4y=0;
X      }
X      else {
X        i1y=0;
X        i2y=yminor;
X        i3y=0;
X        i4y=ymajor;
X      }
X
X/* Draw the bottom edge of the box */
X
X      if (lbx) {
X        movphy(vppxmi,vppymi);
X        if (ltx) {
X          tp=xtick1*floor(vpwxmi/xtick1);
Xbedge:
X          tn=tp+xtick1;
X          if (lsx) {
X            if (llx) {
X              for(i=0; i<=7;i++){
X                temp=tp+xlog[i];
X                if (betw(temp,vpwxmi,vpwxma))
X                     plxtik(wcpcx(temp),vppymi,i1x,i2x);
X              }
X            }
X            else {
X              for ( i=1;i<=nxsub1-1;i++) {
X                temp=tp+i*(tn-tp)/nxsub1;
X                if (betw(temp,vpwxmi,vpwxma))
X                     plxtik(wcpcx(temp),vppymi,i1x,i2x);
X              }
X            }
X          }
X          temp=tn;
X          if (betw(temp,vpwxmi,vpwxma)) {
X            plxtik(wcpcx(temp),vppymi,i3x,i4x);
X            tp=tn;
X            goto bedge;
X          }
X        }
X        draphy(vppxma,vppymi);
X      }
X                                           
X/* Draw right-hand edge of box */
X
X      if (lcy) {
X        movphy(vppxma,vppymi);
X        if (lty) {
X          tp=ytick1*floor(vpwymi/ytick1);
Xredge:
X          tn=tp+ytick1;
X          if (lsy) {
X            if (lly) {
X              for(i=0;i<=7;i++) {
X                temp=tp+xlog[i];
X                if (betw(temp,vpwymi,vpwyma))
X                     plytik(vppxma,wcpcy(temp),i2y,i1y);
X              }
X            }
X            else {
X              for(i=1;i<=nysub1-1;i++) {
X                temp=tp+i*(tn-tp)/nysub1;
X                if (betw(temp,vpwymi,vpwyma))
X                     plytik(vppxma,wcpcy(temp),i2y,i1y);
X              }
X            }
X          }
X          temp=tn;
X          if (betw(temp,vpwymi,vpwyma)) {
X            plytik(vppxma,wcpcy(temp),i4y,i3y);
X            tp=tn;
X            goto redge;
X          }
X        }
X        draphy(vppxma,vppyma);
X      }
X
X/* Draw the top edge of the box */
X
X      if (lcx) {
X        movphy(vppxma,vppyma);
X        if (ltx) {
X          tp=xtick1*(floor(vpwxma/xtick1)+1);
Xtedge:
X          tn=tp-xtick1;
X          if (lsx) {
X            if (llx) {
X              for(i=7;i>=0;i--) {
X                temp=tn+xlog[i];
X                if (betw(temp,vpwxmi,vpwxma))
X                     plxtik(wcpcx(temp),vppyma,i2x,i1x);
X              }
X            }
X            else {
X              for(i=nxsub1-1;i>=1;i--) {
X                temp=tn+i*(tp-tn)/nxsub1;
X                if (betw(temp,vpwxmi,vpwxma))
X                     plxtik(wcpcx(temp),vppyma,i2x,i1x);
X              }
X            }
X          }
X          temp=tn;
X          if (betw(temp,vpwxmi,vpwxma)) {
X            plxtik(wcpcx(temp),vppyma,i4x,i3x);
X            tp=tn;
X            goto tedge;
X          }
X        }
X        draphy(vppxmi,vppyma);
X      }
X
X/* Draw left-hand edge of box */
X
X      if (lby) {
X        movphy(vppxmi,vppyma);
X        if (lty) {
X          tp=ytick1*(floor(vpwyma/ytick1)+1);
Xledge:
X          tn=tp-ytick1;
X          if (lsy) {
X            if (lly) {
X              for(i=7;i>=0;i--) {
X                temp=tn+xlog[i];
X                if (betw(temp,vpwymi,vpwyma))
X                     plytik(vppxmi,wcpcy(temp),i1y,i2y);
X              }
X            }
X            else {
X              for(i=nysub1-1;i>=1;i--) {
X                temp=tn+i*(tp-tn)/nysub1;
X                if (betw(temp,vpwymi,vpwyma))
X                     plytik(vppxmi,wcpcy(temp),i1y,i2y);
X              }
X            }
X          }
X          temp=tn;
X          if (betw(temp,vpwymi,vpwyma)) {
X            plytik(vppxmi,wcpcy(temp),i3y,i4y);
X            tp=tn;
X            goto ledge;
X          }
X        }
X        draphy(vppxmi,vppymi);
X      }
X
X/* Draw the horizontal axis */
X
X      if (lax) {
X        it0=wcpcy(0.0);
X        movphy(vppxmi,it0);
X        if (ltx) {
X          tp=xtick1*floor(vpwxmi/xtick1);
Xhaxis:
X          tn=tp+xtick1;
X          if (lsx) {
X            if (llx) {
X              for(i=0;i<=7;i++) {
X                temp=tp+xlog[i];
X                if (betw(temp,vpwxmi,vpwxma))
X                     plxtik(wcpcx(temp),it0,xminor,xminor);
X              }
X            }
X            else {
X              for(i=1;i<=nxsub1-1;i++) {
X                temp=tp+i*(tn-tp)/nxsub1;
X                if (betw(temp,vpwxmi,vpwxma))
X                     plxtik(wcpcx(temp),it0,xminor,xminor);
X              }
X            }
X          }
X          temp=tn;
X          if (betw(temp,vpwxmi,vpwxma)) {
X            plxtik(wcpcx(temp),it0,xmajor,xmajor);
X            tp=tn;
X            goto haxis;
X          }
X        }
X        draphy(vppxma,it0);
X      }
X
X/* Draw the vertical axis */
X
X      if (lay) {
X        it0=wcpcx(0.0);
X        movphy(it0,vppymi);
X        if (lty) {
X          tp=ytick1*floor(vpwymi/ytick1);
Xvaxis:
X          tn=tp+ytick1;
X          if (lsy) {
X            if (lly) {
X              for(i=0;i<=7;i++) {
X                temp=tp+xlog[i];
X                if (betw(temp,vpwymi,vpwyma))
X                   plytik(it0,wcpcy(temp),yminor,yminor);
X              }
X            }
X            else {
X              for(i=1;i<=nysub1-1;i++) {
X                temp=tp+i*(tn-tp)/nysub1;
X                if (betw(temp,vpwymi,vpwyma))
X                     plytik(it0,wcpcy(temp),yminor,yminor);
X              }
X            }
X          }
X          temp=tn;
X          if (betw(temp,vpwymi,vpwyma)) {
X            plytik(it0,wcpcy(temp),ymajor,ymajor);
X            tp=tn;
X            goto vaxis;
X          }
X        }
X        draphy(it0,vppyma);
X      }
X
X/* Draw grid in x direction */
X
X      if (lgx) {
X        tp=xtick1*floor(vpwxmi/xtick1);
Xxgrid:
X          tn=tp+xtick1;
X          if (betw(tn,vpwxmi,vpwxma)) {
X            pljoin(tn,vpwymi,tn,vpwyma);
X            tp=tn;
X            goto xgrid;
X          }
X      }
X
X/* Draw grid in y direction */
X
X      if (lgy) {
X        tp=ytick1*floor(vpwymi/ytick1);
Xygrid:
X          tn=tp+ytick1;
X          if (betw(tn,vpwymi,vpwyma)) {
X            pljoin(vpwxmi,tn,vpwxma,tn);
X            tp=tn;
X            goto ygrid;
X          }
X      }
X
X/* Write horizontal label(s) */
X 
X      if ((lmx || lnx) && ltx) {
X        tp=xtick1*floor(vpwxmi/xtick1);
Xhlabel:
X        tn=tp+xtick1;
X        if (betw(tn,vpwxmi,vpwxma)) {
X          if (!llx)
X            plform(tn,xmode,xprec,string);
X          else {
X            sprintf(strtmp,"%-d",round(tn));
X            strcpy(string,"10\\u");
X            strcat(string,strtmp);
X          }
X          pos=(tn-vpwxmi)/(vpwxma-vpwxmi);
X          if (lnx) plmtex("b",1.5,pos,0.5,string);
X          if (lmx) plmtex("t",1.5,pos,0.5,string);
X          tp=tn;
X          goto hlabel;
X        }
X      }
X
X/* Write vertical label(s) */
X 
X      if ((lmy || lny) && lty) {
X        tp=ytick1*floor(vpwymi/ytick1);
Xvlabel:
X        tn=tp+ytick1             ;
X        if (betw(tn,vpwymi,vpwyma)) {
X          if (!lly) 
X            plform(tn,ymode,yprec,string);
X          else {
X            sprintf(strtmp,"%-d",round(tn));
X            strcpy(string,"10\\u");
X            strcat(string,strtmp);
X          }
X          pos=(tn-vpwymi)/(vpwyma-vpwymi);
X          if (lny) {
X            if (lvy)
X              plmtex("lv",0.5,pos,1.0,string);
X            else
X              plmtex("l",1.5,pos,0.5,string);
X          }
X          if (lmy) {
X            if (lvy) 
X              plmtex("rv",0.5,pos,0.0,string);
X            else
X              plmtex("r",1.5,pos,0.5,string);
X          }
X          tp=tn;
X          goto vlabel;
X        }
X      }
X
X/* Restore the clip limits to viewport edge */
X
X      sclp(lxmin,lxmax,lymin,lymax);
X}
SHAR_EOF
echo "extracting plbox3.c"
sed 's/^X//' << \SHAR_EOF > plbox3.c
X/* Draws axes and axis labels for 3-d plots */
X
X#include "plplot.h"
X
Xvoid plbox3(xopt,xlabel,xtick,nsubx,yopt,ylabel,ytick,nsuby,
X            zopt,zlabel,ztick,nsubz)
Xchar *xopt, *xlabel, *yopt, *ylabel, *zopt, *zlabel;
Xint nsubx, nsuby, nsubz;
Xfloat xtick, ytick, ztick;
X{      
X      float dx,dy,tx,ty,ux,uy;
X      float xmin,xmax,ymin,ymax,zmin,zmax,zscale;
X      float cxx,cxy,cyx,cyy,cyz;
X      int ln;
X
X      int level;
X      glev(&level);
X      if (level < 3) fatal("Please set up window before calling PLBOX3");
X
X      gw3wc(&cxx,&cxy,&cyx,&cyy,&cyz);
X      gdom(&xmin,&xmax,&ymin,&ymax);
X      grange(&zscale,&zmin,&zmax);
X
X      if (cxx >= 0.0 && cxy <= 0.0) {
X        ln=strpos(xopt,'N') != -1 || strpos(xopt,'n') != -1;
X        tx=w3wcx(xmin,ymin,zmin);
X        ty=w3wcy(xmin,ymin,zmin);
X        ux=w3wcx(xmax,ymin,zmin);
X        uy=w3wcy(xmax,ymin,zmin);
X        plxybx(xopt,xlabel,tx,ty,ux,uy,xmin,xmax,xtick,nsubx,0);
X        dx = ux - tx;
X        dy = uy - ty;
X        plzbx(zopt,zlabel,1,dx,dy,ux,uy,
X             w3wcy(xmax,ymin,zmax),zmin,zmax,ztick,nsubz);
X        tx=w3wcx(xmin,ymax,zmin);
X        ty=w3wcy(xmin,ymax,zmin);
X        ux=w3wcx(xmin,ymin,zmin);
X        uy=w3wcy(xmin,ymin,zmin);
X        plxybx(yopt,ylabel,tx,ty,ux,uy,ymax,ymin,ytick,nsuby,ln);
X        dx = ux - tx;
X        dy = uy - ty;
X        plzbx(zopt,zlabel,0,dx,dy,tx,ty,
X             w3wcy(xmin,ymax,zmax),zmin,zmax,ztick,nsubz);
X
X      }
X      else if (cxx <= 0.0 && cxy <= 0.0) {
X        ln=strpos(yopt,'N') != -1 || strpos(yopt,'n') != -1;
X        tx=w3wcx(xmin,ymax,zmin);
X        ty=w3wcy(xmin,ymax,zmin);
X        ux=w3wcx(xmin,ymin,zmin);
X        uy=w3wcy(xmin,ymin,zmin);
X        plxybx(yopt,ylabel,tx,ty,ux,uy,ymax,ymin,ytick,nsuby,0);
X        dx = ux - tx;
X        dy = uy - ty;
X        plzbx(zopt,zlabel,1,dx,dy,ux,uy,
X             w3wcy(xmin,ymin,zmax),zmin,zmax,ztick,nsubz);
X        tx=w3wcx(xmax,ymax,zmin);
X        ty=w3wcy(xmax,ymax,zmin);
X        ux=w3wcx(xmin,ymax,zmin);
X        uy=w3wcy(xmin,ymax,zmin);
X        plxybx(xopt,xlabel,tx,ty,ux,uy,xmax,xmin,xtick,nsubx,ln);
X        dx = ux - tx;
X        dy = uy - ty;
X        plzbx(zopt,zlabel,0,dx,dy,tx,ty,
X             w3wcy(xmax,ymax,zmax),zmin,zmax,ztick,nsubz);
X
X      }
X      else if (cxx <= 0.0 && cxy >= 0.0) {
X        ln=strpos(xopt,'N') != -1 || strpos(xopt,'n') != -1;
X        tx=w3wcx(xmax,ymax,zmin);
X        ty=w3wcy(xmax,ymax,zmin);
X        ux=w3wcx(xmin,ymax,zmin);
X        uy=w3wcy(xmin,ymax,zmin);
X        plxybx(xopt,xlabel,tx,ty,ux,uy,xmax,xmin,xtick,nsubx,0);
X        dx = ux - tx;
X        dy = uy - ty;
X        plzbx(zopt,zlabel,1,dx,dy,ux,uy,
X             w3wcy(xmin,ymax,zmax),zmin,zmax,ztick,nsubz);
X        tx=w3wcx(xmax,ymin,zmin);
X        ty=w3wcy(xmax,ymin,zmin);
X        ux=w3wcx(xmax,ymax,zmin);
X        uy=w3wcy(xmax,ymax,zmin);
X        plxybx(yopt,ylabel,tx,ty,ux,uy,ymin,ymax,ytick,nsuby,ln);
X        dx = ux - tx;
X        dy = uy - ty;
X        plzbx(zopt,zlabel,0,dx,dy,tx,ty,
X             w3wcy(xmax,ymin,zmax),zmin,zmax,ztick,nsubz);
X      }
X      else if (cxx >= 0.0 && cxy >= 0.0) {
X        ln=strpos(yopt,'N') != -1 || strpos(yopt,'n') != -1;
X        tx=w3wcx(xmax,ymin,zmin);
X        ty=w3wcy(xmax,ymin,zmin);
X        ux=w3wcx(xmax,ymax,zmin);
X        uy=w3wcy(xmax,ymax,zmin);
X        plxybx(yopt,ylabel,tx,ty,ux,uy,ymin,ymax,ytick,nsuby,0);
X        dx = ux - tx;
X        dy = uy - ty;
X        plzbx(zopt,zlabel,1,dx,dy,ux,uy,
X             w3wcy(xmax,ymax,zmax),zmin,zmax,ztick,nsubz);
X        tx=w3wcx(xmin,ymin,zmin);
X        ty=w3wcy(xmin,ymin,zmin);
X        ux=w3wcx(xmax,ymin,zmin);
X        uy=w3wcy(xmax,ymin,zmin);
X        plxybx(xopt,xlabel,tx,ty,ux,uy,xmin,xmax,xtick,nsubx,ln);
X        dx = ux - tx;
X        dy = uy - ty;
X        plzbx(zopt,zlabel,0,dx,dy,tx,ty,
X             w3wcy(xmin,ymin,zmax),zmin,zmax,ztick,nsubz);
X
X      }
X}
X
SHAR_EOF
echo "extracting plccal.c"
sed 's/^X//' << \SHAR_EOF > plccal.c
X/* Subroutine to interpolate the position of a contour which is known */
X/* to be next to ix,iy in the direction ixg,iyg. The unscaled distance */
X/* along ixg,iyg is returned as dist */
X
X#include "plplot.h"
X
Xvoid plccal(pts,nx,ny,zlev,ix,iy,ixg,iyg,dist)
Xint nx,ny,ix,iy,ixg,iyg;
Xfloat *pts,zlev,*dist;
X{
X      int ia,ib;
X      float dbot, dtop, pmid, qmid, zmid;
X
X      ia = ix+ixg;
X      ib = iy+iyg;
X      if (ixg == 0 || iyg == 0) {
X        dtop = zlev - *(pts+(ix-1)*ny+iy-1);
X        dbot = *(pts+(ia-1)*ny+ib-1) - *(pts+(ix-1)*ny+iy-1);
X        *dist = 0.0;
X        if (dbot != 0.0) *dist = dtop/dbot;
X      }
X      else {
X        pmid = *(pts+(ix-1)*ny+iy-1) + *(pts+(ia-1)*ny+ib-1);
X        qmid = *(pts+(ix-1)*ny+ib-1) + *(pts+(ia-1)*ny+iy-1);
X        zmid = (pmid+qmid)/4.0;
X        if (zmid >= zlev) {
X          dtop = zlev - *(pts+(ix-1)*ny+iy-1);
X          dbot = zmid - *(pts+(ix-1)*ny+iy-1);
X          *dist = 0.0;
X          if (dbot != 0.0) *dist = 0.5*dtop/dbot;
X        }
X        else {
X          dtop = zlev-zmid;
X          dbot = *(pts+(ia-1)*ny+ib-1) - zmid;
X          *dist = 0.5;
X          if (dbot != 0.0) *dist = 0.5 + 0.5*dtop/dbot;
X        }
X      }
X}
SHAR_EOF
echo "extracting plclr.c"
sed 's/^X//' << \SHAR_EOF > plclr.c
X#include "plplot.h"
X
Xvoid plclr()
X{
X      int level;
X      glev(&level);
X      if (level < 1) fatal("Please call PLSTAR before calling PLCLR.");
X
X      grclr();
X}
SHAR_EOF
echo "extracting plcntr.c"
sed 's/^X//' << \SHAR_EOF > plcntr.c
X/* Contour "ny" by "nx" real array "points" at the level "zlev" */
X/* points is a pointer to a 2d array of nx by ny points. */
X/* iscan has nx elements. ixstor and iystor each have nstor elements. */
X
X#include "plplot.h"
X
Xvoid plcntr(points,nx,ny,kx,lx,ky,ly,zlev,iscan,ixstor,iystor,nstor,tr)
Xint nx, ny, ky, ly, kx, lx, nstor;
Xfloat zlev, *points;
Xint iscan[], ixstor[], iystor[];
Xvoid (*tr)();
X{
X      int kcol, krow, kstor, kscan, iwbeg, ixbeg, iybeg, izbeg;
X      int iboun, iw, ix, iy, iz, ifirst, istep, ixgo, iygo;
X      int l, ixg, iyg, ia, ib, ixt, iyt, jstor, next;
X      float dist, dx, dy, xnew, ynew, x, y;
X      float xlas=0., ylas=0., tpx, tpy, xt, yt;
X
X      /* Initialize memory pointers */
X
X      kstor = 0;
X      kscan = 0;
X
X      for (krow=ky; krow<=ly; krow++) {
X        for (kcol=kx+1; kcol <= lx; kcol++)  {
X
X          /* Check if a contour has been crossed */
X
X          x = *(points + (kcol-2)*ny + krow-1);
X          y = *(points + (kcol-1)*ny + krow-1);
X          if (x < zlev && y >= zlev) {
X            ixbeg = kcol-1;
X            iwbeg = kcol;
X          }
X          else if (y < zlev && x >= zlev) {
X            ixbeg = kcol;
X            iwbeg = kcol-1;
X          }
X          else
X            goto lab70;
X
X          iybeg = krow;
X          izbeg = krow;
X
X          /* Yes, a contour has been crossed. Check to see if it */
X          /* is a new one. */
X
X          for(l=0;l<kscan;l++)
X            if (ixbeg == iscan[l]) goto lab70;
X
X          /* Here is the section which follows and draws a contour */
X
X          for (iboun=1; iboun>= -1; iboun -= 2) {
X
X            /* Set up starting point and initial search directions */
X
X            ix = ixbeg;
X            iy = iybeg;
X            iw = iwbeg;
X            iz = izbeg;
X            ifirst = 1;
X            istep = 0;
X            ixgo = iw - ix;
X            iygo = iz - iy;
X
Xlab20:
X              plccal(points,nx,ny,zlev,ix,iy,ixgo,iygo,&dist);
X              dx = dist * ixgo;
X              dy = dist * iygo;
X              xnew = ix+dx;
X              ynew = iy+dy;
X
X             /* Has a step occured in search? */
X
X              if (istep != 0) {
X                if (ixgo*iygo == 0) {
X
X                  /* This was a diagonal step, so interpolate missed */
X                  /* point, rotating 45 degrees to get it */
X
X                  ixg = ixgo;
X                  iyg = iygo;
X                  plr45(&ixg,&iyg,iboun);
X                  ia = iw-ixg;
X                  ib = iz-iyg;
X                  plccal(points,nx,ny,zlev,ia,ib,ixg,iyg,&dist);
X                  (*tr)(xlas,ylas,&tpx,&tpy);
X                  drawor(tpx,tpy);
X                  dx = dist*ixg;
X                  dy = dist*iyg;
X                  xlas = ia+dx;
X                  ylas = ib+dy;
X                }
X                else {
X                  if (dist > 0.5) {
X                    xt = xlas;
X                    xlas = xnew;
X                    xnew = xt;
X                    yt = ylas;
X                    ylas = ynew;
X                    ynew = yt;
X                  }
X                }
X              }
X              if (ifirst != 1) {
X                (*tr)(xlas,ylas,&tpx,&tpy);
X                drawor(tpx,tpy);
X              }
X              else {
X                (*tr)(xnew,ynew,&tpx,&tpy);
X                movwor(tpx,tpy);
X              }
X              xlas = xnew;
X              ylas = ynew;
X
X              /* Check if the contour is closed */
X
X              if (ifirst != 1 && ix == ixbeg && iy == iybeg
X                  && iw == iwbeg && iz == izbeg) {
X                (*tr)(xlas,ylas,&tpx,&tpy);
X                drawor(tpx,tpy);
X                goto lab70;
X              }
X              ifirst = 0;
X
X              /* Now the rotation */
X
X              istep = 0;
X              plr45(&ixgo,&iygo,iboun);
X              iw = ix+ixgo;
X              iz = iy+iygo;
X
X              /* Check if out of bounds */
X
X              if (iw<kx || iw>lx || iz<ky || iz>ly) goto  lab50;
X
X              /* Has contact been lost with the contour? */
X
X              if (*(points+(iw-1)*ny+iz-1) < zlev) {
X
X                /* Yes, lost contact => step to new centre */
X
X                istep = 1;
X                ix = iw;
X                iy = iz;
X                plr135(&ixgo,&iygo,iboun);
X                iw = ix+ixgo;
X                iz = iy+iygo;
X
X                /* And do the contour memory */
X
X                if (iy == krow) {
X                  kscan = kscan+1;
X                  iscan[kscan-1] = ix;
X                }
X                else if (iy>krow) {
X                  kstor = kstor+1;
X                  if (kstor>nstor) {
X                    fatal("Heap exhausted in PLCONT.");
X                    goto lab70;
X                  }
X                  ixstor[kstor-1] = ix;
X                  iystor[kstor-1] = iy;
X                }
X              }
X            goto lab20;
Xlab50:
X            /* Reach here only if boundary encountered - Draw last bit */
X
X            (*tr)(xnew,ynew,&tpx,&tpy);
X            drawor(tpx,tpy);
X          }
Xlab70:
X          ;    /* Null statement to carry label */
X        }
X
X        /* Search of row complete - set up memory of next row in iscan and */
X        /* edit ixstor and iystor */
X
X        if (krow<ny) {
X          jstor = 0;
X          kscan = 0;
X          next = krow+1;
X          for (l=1; l<=kstor; l++) {
X            ixt = ixstor[l-1];
X            iyt = iystor[l-1];
X
X            /* Memory of next row into iscan */
X
X            if (iyt == next) {
X              kscan = kscan+1;
X              iscan[kscan-1] = ixt;
X      
X            /* Retain memory of rows to come, and forget rest */
X
X            }
X            else if (iyt>next) {
X              jstor = jstor+1;
X              ixstor[jstor-1] = ixt;
X              iystor[jstor-1] = iyt;
X            }
X          }
X          kstor = jstor;
X        }
X      }
X}
SHAR_EOF
echo "extracting plcol.c"
sed 's/^X//' << \SHAR_EOF > plcol.c
X/* Sets line colour */
X
X#include "plplot.h"
X
Xvoid plcol(colour)
Xint colour;
X{
X      int font,col;
X
X      int level;
X      glev(&level);
X      if (level < 1)  fatal("Please call PLSTAR before calling PLCOL.");
X      if (colour < 0)  fatal("Invalid colour in PLCOL.");
X
X      gatt(&font,&col);
X      satt(font,colour);
X      grcol();
X}
SHAR_EOF
echo "extracting plcont.c"
sed 's/^X//' << \SHAR_EOF > plcont.c
X/* Draws a contour plot from data in z(nx,ny), using the subarray */
X/* from kx to lx in the x direction and from ky to ly in the y */
X/* direction. The array of contour levels is clevel(nlevel), and */
X/* "tr" is the name of a subroutine which transforms array indicies */
X/* into world coordinates */
X
X#include <stdio.h>
X#include <stdlib.h>
X#include "plplot.h"
X#include "declare.h"
X
Xvoid plcont(z,nx,ny,kx,lx,ky,ly,clevel,nlevel,tr)
Xint nx, ny, kx, lx, ky, ly, nlevel;
Xfloat *z, clevel[];
Xvoid (*tr)();
X{
X      int i, mx, my, nstor;
X
X      mx = lx - kx + 1;
X      my = ly - ky + 1;
X
X      if (kx < 1 || lx > nx || kx >= lx || ky < 1 || ky > ny || ky >= ly)
X               fatal("Argument error in PLCONT");
X      
X      nstor = mx*my/5;
X
X      if( heapc != NULL)
X         free(heapc);
X      if(( heapc = (int *)malloc((mx+2*nstor)*sizeof(int))) == NULL)
X         fatal("Not enough heap space in PLCONT.");
X
X      for (i=0; i<nlevel; i++)
X        plcntr(z,nx,ny,kx,lx,ky,ly,clevel[i],&heapc[0],
X               &heapc[nx],&heapc[nx+nstor],nstor,tr);
X
X}
SHAR_EOF
echo "extracting plcvec.c"
sed 's/^X//' << \SHAR_EOF > plcvec.c
X/* Gets the character digitisation of Hershey table entry "char" */
X/* Returns 1 if there is a valid entry */
X
Xextern short int *buffer[];
Xextern int *findex[];
X
X#include "plplot.h"
X
Xint plcvec(ch,xygrid)
Xint ch;
Xshort int xygrid[];
X{
X
X      int nc1, nc2, k, ib, ix, iy;
X
X      nc1=1;
X      nc2=3000;
X      if (ch < nc1 || ch > nc2) return(0);
X      ib = *(findex[(ch-nc1)/100] + (ch-nc1)%100);
X      if (ib == 0) return(0);
X
X      k=1;
X      xygrid[k-1]= *(buffer[(ib-1)/100] + (ib-1)%100);
X
X      k=k+1;
X      do {
X        ib=ib+1;
X        ix= *(buffer[(ib-1)/100]+(ib-1)%100)/128 - 64;
X        xygrid[k-1]=ix;
X        k=k+1;
X        iy = *(buffer[(ib-1)/100]+(ib-1)%100)%128 - 64;
X        xygrid[k-1]=iy;
X        k=k+1;
X      } while (ix != -64 || iy != -64);
X
X      return(1);
X
X} 
SHAR_EOF
echo "extracting pldeco.c"
sed 's/^X//' << \SHAR_EOF > pldeco.c
X/* Decode a character string, and return an array of integer symbol */
X/* numbers. This routine is responsible for interpreting all escape */
X/* sequences. At present the following escape sequences are defined */
X/* (the letter following the \ may be either upper or lower case): */
X
X/* \u       :      up one level (returns -1) */
X/* \d       :      down one level (returns -2) */
X/* \b       :      backspace (returns -3) */
X/* \+       :      toggles overline mode (returns -4) */
X/* \-       :      toggles underline mode (returns -5) */
X/* \\       :      \, returns the code for backslash */
X/* \gx      :      greek letter corresponding to roman letter x */
X/* \fn      :      switch to Normal font */
X/* \fr      :      switch to Roman font */
X/* \fi      :      switch to Italic font */
X/* \fs      :      switch to Script font */
X/* \(nnn)   :      Hershey symbol number nnn (any number of digits) */
X
X#include "plplot.h"
X#include <stdio.h>    /* include for definition of NULL */
X#include <string.h>
X
Xstatic char font[] = "nrisnris";
Xstatic char greek[] = "ABGDEZYHIKLMNCOPRSTUFXQWabgdezyhiklmncoprstufxqw";
X
Xextern short int *hersh[];
X
Xvoid pldeco(symbol,length,text)
Xint symbol[],*length;
Xchar text[];
X{
X      int ch,icol,ifont,ig,j,lentxt;
X      char test;
X
X      /* Initialize parameters. */
X
X      lentxt=strlen(text);
X      *length=0;
X      j=0;
X      gatt(&ifont,&icol);
X
X      /* Get next character; treat non-printing characters as spaces. */
X
Xlab100:
X      j=j+1;
X      if (j>lentxt) return;
X      test=text[j-1];
X      ch=test;
X      if (ch<0)   ch = 32;
X
X      if (ch>175) ch = 32;
X
X      /* Test for escape sequence (\) */
X
X      if (ch=='\\') {
X        if ((lentxt-j)>=1) {
X          test=text[j];
X          if (test=='\\')
X            j = j+1;
X          else if (test=='U' || test=='u') {
X            *length = *length + 1;
X            symbol[*length-1] = -1;
X            j = j+1;
X            goto lab100;
X          }
X          else if (test=='D' || test=='d') {
X            *length = *length + 1;
X            symbol[*length-1] = -2;
X            j = j+1;
X            goto lab100;
X          }
X          else if (test=='B' || test=='b') {
X            *length = *length + 1;
X            symbol[*length-1] = -3;
X            j = j+1;
X            goto lab100;
X          }
X          else if (test=='+') {
X            *length = *length + 1;
X            symbol[*length-1] = -4;
X            j = j+1;
X            goto lab100;
X          }
X          else if (test=='-') {
X            *length = *length + 1;
X            symbol[*length-1] = -5;
X            j = j+1;
X            goto lab100;
X          }
X          else if (test=='(') {
X            *length = *length + 1;
X            symbol[*length-1] = 0;
X            j = j+2;
Xlab10:
X            if ('0'<=text[j-1] && text[j-1]<='9') {
X              symbol[*length-1] = symbol[*length-1]*10 + text[j-1] - '0';
X              j = j+1;
X              goto lab10;
X            }
X            if (text[j-1]!=')') j = j-1;
X            goto lab100;
X          }
X          else if (test=='F' || test=='f') {
X            test=text[j+1];
X            ifont = strpos(font,test) + 1;
X            if (ifont>4) ifont = ifont-4;
X            if (ifont==0) ifont = 1;
X            j = j+2;
X            goto lab100;
X          }
X          else if (test=='G' || test=='g') {
X            test=text[j+1];
X            ig = strpos(greek,test) + 1;
X            *length = *length + 1;
X            symbol[*length-1] = *(hersh[ifont-1] + 127 + ig);
X            j = j+2;
X            goto lab100;
X          }
X        }
X      }
X
X      /* Decode character. */
X
X      *length = *length + 1;
X      symbol[*length-1] = *(hersh[ifont-1]+ch);
X      goto lab100;
X}
SHAR_EOF
echo "extracting pldtik.c"
sed 's/^X//' << \SHAR_EOF > pldtik.c
X/* If tick == 0, this works out a "nice" interval, so that there */
X/* are between 3 and 7.5 major tick intervals in the input range */
X/* "vmin" to "vmax". Using this value for the tick interval or   */
X/* supplied value, it also computes "prec" which specifies       */
X/* the number of places that should be written after the decimal */
X/* point. The recommended number of subticks is returned in      */
X/* "nsubt" unless the routine is entered with a non-zero value   */
X/* of "nsubt". The output variable "mode" is set to 0 if         */
X/* labels are to be written in floating-point format, or to 1 if */
X/* they are to be written in fixed-point format.                 */
X                                                                
X#include "plplot.h"
X#include <math.h>
X
Xvoid pldtik(vmin, vmax, tick, nsubt, mode, prec)
Xfloat vmin, vmax, *tick;
Xint *nsubt, *mode, *prec;
X{
X      
X      float t1, t2, vmod;
X      int msd, np, ns;
X
X      vmod = max(fabs(vmin),fabs(vmax));
X      *mode = 0;
X      if (vmod < 1.0e-2 || vmod > 1.0e6) *mode = 1;
X      t1 = (float)log10(vmod);
X      msd = (int)t1;
X
X      t1 = (float)log10(fabs(vmax-vmin));
X      np = (int)t1;
X      t1 = t1 - np;
X      
X      if (t1 > 0.7781512503) {
X        t2 = 2.0 ;
X        ns = 4;
X      }
X      else if (t1 > 0.4771212549) {
X        t2 = 1.0 ;
X        ns = 5;
X      }
X      else if (t1 > 0.1760912591) {
X        t2 = 5.0;
X        ns = 5;
X        np = np-1;
X      }
X      else  {
X        t2 = 2.0;
X        ns = 4;
X        np = np-1;
X      }
X
X      if (*tick == 0.0) {
X        *tick = t2 * pow(10.0,(double)np);
X        if (vmin > vmax) *tick = -*tick;
X        if (*nsubt == 0) *nsubt = ns;
X        if (*mode != 0) 
X          *prec = msd - np;
X        else
X          *prec = max(-np,0);
X      }
X      else  {
X        t1 = (float)log10(fabs(*tick));
X        np = (int)t1;
X        if (*mode != 0) 
X          *prec = msd - np + 1;
X        else
X          *prec = max(-np+1,0);
X      }
X}
X
SHAR_EOF
echo "extracting plend.c"
sed 's/^X//' << \SHAR_EOF > plend.c
X/* Ends a plotting session */
X
X#include <stdio.h>
X#include <stdlib.h>
X#include "plplot.h"
X#include "declare.h"
X
Xvoid plend()
X{
X      int dev, term, gra, level;
X      
X      glev(&level);
X      if (level == 0) return;
X      gdev(&dev,&term,&gra);
X      if( heap3 != NULL) free((void *)heap3);
X      if( heapc != NULL) free((void *)heapc);
X      if (gra != 0) pltext();
X      grtidy();
X      slev(0);
X}
SHAR_EOF
echo "extracting plenv.c"
sed 's/^X//' << \SHAR_EOF > plenv.c
X/* Simple interface for defining viewport and window. If "just"=1,  */
X/* X and Y scales will be the same, otherwise they are scaled       */
X/* independently. The "axis" parameter is interpreted as follows:   */
X/* axis=-2 : draw no box, axis or labels                            */
X/* axis=-1 : draw box only                                          */
X/* axis= 0 : Draw box and label with coordinates                    */
X/* axis= 1 : Also draw the coordinate axes                          */
X/* axis= 2 : Draw a grid at major tick positions                    */
X/* axis=10 : Logarithmic X axis, L!=r Y axis, No X=0 axis         */
X/* axis=11 : Logarithmic X axis, L!=r Y axis, X=0 axis            */
X/* axis=20 : L!=r X axis, Logarithmic Y axis, No Y=0 axis         */
X/* axis=21 : L!=r X axis, Logarithmic Y axis, Y=0 axis            */
X/* axis=30 : Logarithmic X and Y axes                               */
X
X#include "plplot.h"
X#include <math.h>
X
Xvoid plenv(xmin,xmax,ymin,ymax,just,axis)
Xint just,axis;
Xfloat xmin, xmax, ymin, ymax;
X{
X    int level;
X    float chrdef, chrht;
X    float lb, rb, tb, bb, dx, dy;
X    float xsize, ysize, xscale, yscale;
X    float spxmin, spxmax, spymin, spymax;
X    float vpxmin, vpxmax, vpymin, vpymax;
X    float scale;
X
X    glev(&level);
X    if (level < 1) fatal("Please call PLSTAR before PLENV.");
X    
X    if (xmin == xmax) fatal("Invalid xmin and xmax arguments in PLENV");
X    if (ymin == ymax) fatal("Invalid ymin and ymax arguments in PLENV");
X    if ((just != 0) && (just != 1)) fatal("Invalid just option in PLENV");
X
X
X    pladv(0);
X    if (just == 0)
X      plvsta();
X    else {
X      gchr(&chrdef,&chrht);
X      lb = 7.0 * chrht;
X      rb = 4.0 * chrht;
X      tb = 4.0 * chrht;
X      bb = 4.0 * chrht;
X      dx = fabs(xmax-xmin);
X      dy = fabs(ymax-ymin);
X      plgspa(&spxmin,&spxmax,&spymin,&spymax);
X      xsize = spxmax - spxmin;
X      ysize = spymax - spymin;
X      xscale = dx/(xsize - lb - rb);
X      yscale = dy/(ysize - tb - bb);
X      scale = max(xscale,yscale);
X      vpxmin = max(lb,0.5*(xsize - dx/scale));
X      vpxmax = vpxmin + (dx/scale);
X      vpymin = max(bb,0.5*(ysize - dy/scale));
X      vpymax = vpymin + (dy/scale);
X      plsvpa(vpxmin,vpxmax,vpymin,vpymax);
X    }
X    plwind(xmin,xmax,ymin,ymax);
X    if (axis == -2)
X       ;
X    else if (axis == -1)
X       plbox("bc",0.0,0,"bc",0.0,0);
X    else if (axis == 0) 
X       plbox("bcnst",0.0,0,"bcnstv",0.0,0);
X    else if (axis == 1) 
X       plbox("abcnst",0.0,0,"abcnstv",0.0,0);
X    else if (axis == 2) 
X       plbox("abcgnst",0.0,0,"abcgnstv",0.0,0);
X    else if (axis == 10) 
X       plbox("bclnst",0.0,0,"bcnstv",0.0,0);
X    else if (axis == 11) 
X       plbox("bclnst",0.0,0,"abcnstv",0.0,0);
X    else if (axis == 20) 
X       plbox("bcnst",0.0,0,"bclnstv",0.0,0);
X    else if (axis == 21) 
X       plbox("bcnst",0.0,0,"abclnstv",0.0,0);
X    else if (axis == 30) 
X       plbox("bclnst",0.0,0,"bclnstv",0.0,0);
X    else
X       fatal("Invalid axis argument in plenv");
X}
SHAR_EOF
echo "extracting plerrx.c"
sed 's/^X//' << \SHAR_EOF > plerrx.c
X/* Plots horizontal error bars (xmin(i),y(i)) to (xmax(i),y(i)) */
X
X#include "plplot.h"
X
Xvoid plerrx(n,xmin,xmax,y)
Xint n;
Xfloat xmin[], xmax[], y[];
X{
X      int i, level;
X
X      glev(&level);
X      if (level < 3) fatal("Please set up window before calling PLERRX.");
X      for (i=0;i<n;i++)
X        plerx1(xmin[i],xmax[i],y[i]);
X}
SHAR_EOF
echo "extracting plerry.c"
sed 's/^X//' << \SHAR_EOF > plerry.c
X/* Plots vertical error bars (x,ymin(i)) to (x(i),ymax(i)) */
X
X#include "plplot.h"
X
Xvoid plerry(n,x,ymin,ymax)
Xint n;
Xfloat x[], ymin[], ymax[];
X{
X      int i, level;
X
X      glev(&level);
X      if (level < 3) fatal("Please set up window before calling PLERRY.");
X      for(i=0; i<n; i++)
X        plery1(x[i],ymin[i],ymax[i]);
X}
SHAR_EOF
echo "extracting plerx1.c"
sed 's/^X//' << \SHAR_EOF > plerx1.c
X/* Plots single horizontal error bar */
X
X#include "plplot.h"
X#include <math.h>
X
Xvoid plerx1(xmin,xmax,y)
Xfloat xmin, xmax, y;
X{
X      float mindef, minht, xpmm, ypmm;
X      int yminor;
X      
X      gmin(&mindef,&minht);
X      gpixmm(&xpmm,&ypmm);
X      yminor = max(1.0,minht*ypmm);
X      movwor(xmin,y);
X      plxtik(wcpcx(xmin),wcpcy(y),yminor,yminor);
X      drawor(xmax,y);
X      plxtik(wcpcx(xmax),wcpcy(y),yminor,yminor);
X}
SHAR_EOF
echo "extracting plery1.c"
sed 's/^X//' << \SHAR_EOF > plery1.c
X/* Plots single vertical error bar */
X
X#include "plplot.h"
X#include <math.h>
X
Xvoid plery1(x,ymin,ymax)
Xfloat x, ymin, ymax;
X{
X      float mindef, minht, xpmm, ypmm;
X      int xminor;
X      
X      gmin(&mindef,&minht);
X      gpixmm(&xpmm,&ypmm);
X      xminor = max(1.0,minht*xpmm);
X      movwor(x,ymin);
X      plytik(wcpcx(x),wcpcy(ymin),xminor,xminor);
X      drawor(x,ymax);
X      plytik(wcpcx(x),wcpcy(ymax),xminor,xminor);
X}
SHAR_EOF
echo "extracting plfont.c"
sed 's/^X//' << \SHAR_EOF > plfont.c
X#include "plplot.h"
X
Xvoid plfont(ifont)
Xint ifont;
X{
X      int ifnt,icol;
X
X      int level;
X      glev(&level);
X      if (level < 1) fatal("Please call PLSTAR before calling PLFONT.");
X
X      if (ifont < 1 || ifont > 4) fatal("Invalid font in PLFONT");
X      gatt(&ifnt,&icol);
X      satt(ifont,icol);
X}
SHAR_EOF
echo "extracting plform.c"
sed 's/^X//' << \SHAR_EOF > plform.c
X/* Formats a floating point value in one of the following formats   */
X/* (i)  If mode == 0, use floating point format with "precision"    */
X/*      places after the decimal point.                             */
X/* (ii) If mode == 1, use scientific notation with one place before */
X/*      the decimal point and "precision" places after.             */
X
X#include "plplot.h"
X#include <stdio.h>
X#include <string.h>
X
Xvoid plform(value,mode,prec,result)
Xfloat value;
Xint mode, prec;
Xchar *result;
X{
X      int j, expon;
X      char form[10];
X      char temp[30];
X
X      if (mode == 0) {
X        sprintf(form,"%%-.%df",prec);
X        sprintf(temp,form,value);
X        strcpy(result,temp);
X      }
X      else {
X        sprintf(form,"%%-.%dE",prec);
X        sprintf(temp,form,value);
X        j = strpos(temp,'E') + 1;
X        if (j == 0) fatal("Unrecognized scientific notation");
X        sscanf(&temp[j],"%d",&expon);
X        sprintf(form,"%-d",expon);
X        strcpy(&temp[j-1],"x10\\u");
X        strcat(temp,form);
X        strcpy(result,temp);
X      }
X}
X
SHAR_EOF
echo "extracting plgra.c"
sed 's/^X//' << \SHAR_EOF > plgra.c
X/* Switches to graphics mode */
X
X#include "plplot.h"
X
Xvoid plgra()
X{
X      int level;
X      glev(&level);
X      if (level < 1) fatal("Please call PLSTAR before PLGRA.");
X      grgra();
X}
SHAR_EOF
echo "extracting plgrid3.c"
sed 's/^X//' << \SHAR_EOF > plgrid3.c
X/* Routine to draw a grid around the back side of the 3d plot */
X/* wih hidden line removal */
X
X#include "plplot.h"
X#include <math.h>
X
Xvoid plgrid3(tick)
Xfloat tick;
X{
X    float xmin, ymin, zmin, xmax, ymax, zmax, zscale;
X    float cxx, cxy, cyx, cyy, cyz;
X    int u[3], v[3];
X    int nsub, mode, prec;
X    float tp;
X
X    gw3wc(&cxx,&cxy,&cyx,&cyy,&cyz);
X    gdom(&xmin,&xmax,&ymin,&ymax);
X    grange(&zscale,&zmin,&zmax);
X
X    nsub = 0;
X    pldtik(zmin,zmax,&tick,&nsub,&mode,&prec);
X    tp = tick * floor(zmin/tick) + tick;
X
X    if (cxx >= 0.0 && cxy <= 0.0)  {
X        while ( tp <= zmax ) {
X            u[0] = wcpcx(w3wcx(xmin,ymax,tp));
X            v[0] = wcpcy(w3wcy(xmin,ymax,tp));
X            u[1] = wcpcx(w3wcx(xmax,ymax,tp));
X            v[1] = wcpcy(w3wcy(xmax,ymax,tp));
X            u[2] = wcpcx(w3wcx(xmax,ymin,tp));
X            v[2] = wcpcy(w3wcy(xmax,ymin,tp));
X            plnxtv(u,v,3,2);
X
X            tp += tick;
X        }
X        u[0] = wcpcx(w3wcx(xmax,ymax,zmin));
X        v[0] = wcpcy(w3wcy(xmax,ymax,zmin));
X        u[1] = wcpcx(w3wcx(xmax,ymax,zmax));
X        v[1] = wcpcy(w3wcy(xmax,ymax,zmax));
X        plnxtv(u,v,2,2);
X    }
X    else if(cxx <= 0.0  && cxy <= 0.0)  {
X        while ( tp <= zmax ) {
X            u[0] = wcpcx(w3wcx(xmax,ymax,tp));
X            v[0] = wcpcy(w3wcy(xmax,ymax,tp));
X            u[1] = wcpcx(w3wcx(xmax,ymin,tp));
X            v[1] = wcpcy(w3wcy(xmax,ymin,tp));
X            u[2] = wcpcx(w3wcx(xmin,ymin,tp));
X            v[2] = wcpcy(w3wcy(xmin,ymin,tp));
X            plnxtv(u,v,3,2);
X
X            tp += tick;
X        }
X        u[0] = wcpcx(w3wcx(xmax,ymin,zmin));
X        v[0] = wcpcy(w3wcy(xmax,ymin,zmin));
X        u[1] = wcpcx(w3wcx(xmax,ymin,zmax));
X        v[1] = wcpcy(w3wcy(xmax,ymin,zmax));
X        plnxtv(u,v,2,2);
X    }
X    else if(cxx <= 0.0  && cxy >= 0.0)  {
X        while ( tp <= zmax ) {
X            u[0] = wcpcx(w3wcx(xmax,ymin,tp));
X            v[0] = wcpcy(w3wcy(xmax,ymin,tp));
X            u[1] = wcpcx(w3wcx(xmin,ymin,tp));
X            v[1] = wcpcy(w3wcy(xmin,ymin,tp));
X            u[2] = wcpcx(w3wcx(xmin,ymax,tp));
X            v[2] = wcpcy(w3wcy(xmin,ymax,tp));
X            plnxtv(u,v,3,2);
X
X            tp += tick;
X        }
X        u[0] = wcpcx(w3wcx(xmin,ymin,zmin));
X        v[0] = wcpcy(w3wcy(xmin,ymin,zmin));
X        u[1] = wcpcx(w3wcx(xmin,ymin,zmax));
X        v[1] = wcpcy(w3wcy(xmin,ymin,zmax));
X        plnxtv(u,v,2,2);
X    }
X    else if(cxx >= 0.0  && cxy >= 0.0)  {
X        while ( tp <= zmax ) {
X            u[0] = wcpcx(w3wcx(xmin,ymin,tp));
X            v[0] = wcpcy(w3wcy(xmin,ymin,tp));
X            u[1] = wcpcx(w3wcx(xmin,ymax,tp));
X            v[1] = wcpcy(w3wcy(xmin,ymax,tp));
X            u[2] = wcpcx(w3wcx(xmax,ymax,tp));
X            v[2] = wcpcy(w3wcy(xmax,ymax,tp));
X            plnxtv(u,v,3,2);
X
X            tp += tick;
X        }
X        u[0] = wcpcx(w3wcx(xmin,ymax,zmin));
X        v[0] = wcpcy(w3wcy(xmin,ymax,zmin));
X        u[1] = wcpcx(w3wcx(xmin,ymax,zmax));
X        v[1] = wcpcy(w3wcy(xmin,ymax,zmax));
X        plnxtv(u,v,2,2);
X    }
X}
SHAR_EOF
echo "extracting plgspa.c"
sed 's/^X//' << \SHAR_EOF > plgspa.c
X/* Gets subpage boundaries in absolute coordinates (mm from bottom */
X/* left-hand corner of page) */
X
X#include "plplot.h"
X
Xvoid plgspa(xmin,xmax,ymin,ymax)
Xfloat *xmin, *xmax, *ymin, *ymax;
X{
X      float spdxmi, spdxma, spdymi, spdyma;
X
X      int level;
X
X      glev(&level);
X      if (level < 1) fatal("Please call PLSTAR before PLGSPA.");
X      gspd(&spdxmi,&spdxma,&spdymi,&spdyma);
X      *xmin = dcmmx(spdxmi);
X      *xmax = dcmmx(spdxma);
X      *ymin = dcmmy(spdymi);
X      *ymax = dcmmy(spdyma);
X}
SHAR_EOF
echo "extracting plhist.c"
sed 's/^X//' << \SHAR_EOF > plhist.c
X/* Draws a histogram of n values of a variable in array data(1..n) */
X/* in the range datmin to datmax using nbin bins. If "oldwin" is   */
X/* true, the histogram is plotted in the current window. If not,   */
X/* the routine calls "plenv" to set up the graphics environment.   */
X
X#include "plplot.h"
X#include <math.h>
X
Xvoid plhist(n,data,datmin,datmax,nbin,oldwin)
Xint n, nbin;
Xint oldwin;
Xfloat data[], datmin, datmax;
X{
X      int i, bin, level;
X      float x[201], y[201], dx, ymax;
X      
X      glev(&level);
X      if (level<1) fatal("Please call PLSTAR before PLHIST.");
X      if (level<3 && oldwin)
X              fatal("Please set up window before calling PLHIST.");
X      
X      if (nbin < 1 || nbin > 200)
X              fatal("Cannot have <1 or >200 bins in PLHIST.");
X      if (datmin >= datmax)
X              fatal("Data range invalid in PLHIST.");
X      
X      dx = (datmax-datmin) / nbin;
X      for (i=0; i<=nbin; i++) {
X        x[i] = datmin + i*dx;
X        y[i] = 0.0;
X      }
X
X      for (i=0; i<n; i++)  {
X        bin = (data[i] - datmin)/dx;
X        if (bin >= 0 && bin < nbin) y[bin]++;
X      }
X
X      if (!oldwin) {
X        ymax = 0.0;
X        for (i=0; i<nbin; i++) ymax = max(ymax,y[i]);
X        plenv(datmin,datmax,0.0,1.1*ymax,0,0);
X      }
X         
X      plbin(nbin+1,x,y,0);
X}
SHAR_EOF
echo "extracting plhrsh.c"
sed 's/^X//' << \SHAR_EOF > plhrsh.c
X/* Writes the Hershey symbol "ch" centred at the physical */
X/* coordinate (x,y) */
X
X#include "plplot.h"
X
Xvoid plhrsh(ch,x,y)
Xint ch, x, y;
X{
X      int cx, cy, k, penup;
X      short int xygrid[300];
X      float symdef, symht, scale, xscale, yscale, xpmm, ypmm;
X
X      gsym(&symdef,&symht);
X      gpixmm(&xpmm,&ypmm);
X      k = 5;
X      penup = 1;
X      scale = 0.05 * symht;
X
X      if (!plcvec(ch,xygrid)) {
X        movphy(x,y);
X        return;
X      }
X      
X      /* Compute how many physical pixels correspond to a character pixel */
X
X      xscale = scale * xpmm;
X      yscale = scale * ypmm;
X
X      for(;;) {
X        cx = xygrid[k];
X        k = k+1;
X        cy = xygrid[k];
X        k = k+1;
X        if (cx == -64 && cy == -64) {
X          movphy(x,y);
X          return;
X        }
X        else if (cx == -64 && cy == 0) 
X          penup = 1;
X        else  {
X          if (penup != 0) {
X            movphy(round(x+xscale*cx),round(y+yscale*cy));
X            penup = 0;
X          }
X          else
X            draphy(round(x+xscale*cx),round(y+yscale*cy));
X        }
X      }
X}
SHAR_EOF
echo "extracting pljoin.c"
sed 's/^X//' << \SHAR_EOF > pljoin.c
X#include "plplot.h"
X
Xvoid pljoin(x1,y1,x2,y2)
Xfloat x1, x2, y1, y2;
X{
X   movwor(x1,y1);
X   drawor(x2,y2);
X}
SHAR_EOF
echo "extracting pllab.c"
sed 's/^X//' << \SHAR_EOF > pllab.c
X/* Simple routine for labelling graphs */
X
X#include "plplot.h"
X
Xvoid pllab(xlabel, ylabel, tlabel)
Xchar *xlabel, *ylabel, *tlabel;
X{
X      int level;
X      glev(&level);
X      if (level < 2) fatal("Please set up viewport before calling PLLAB.");
X
X      plmtex("t",2.0,0.5,0.5,tlabel);
X      plmtex("b",3.2,0.5,0.5,xlabel);
X      plmtex("l",5.0,0.5,0.5,ylabel);
X}
SHAR_EOF
echo "extracting pllclp.c"
sed 's/^X//' << \SHAR_EOF > pllclp.c
X/* Draws a (x1,y1) to (x2,y2) within the clip limits */
X
X#include "plplot.h"
X
X#define betw(ix,ia,ib) ((ix<=ia && ix>=ib) || (ix>=ia && ix<=ib))
X#define inside(ix,iy) (betw(ix,clpxmi,clpxma) && betw(iy,clpymi,clpyma))
X#define ixcut(iy) round((float)((x2-x1)*(iy-y1)+x1*(y2-y1))/(float)(y2-y1))
X#define iycut(ix) round((float)((y2-y1)*(ix-x1)+y1*(x2-x1))/(float)(x2-x1))
X
Xvoid pllclp(x1,y1,x2,y2)
Xint x1,y1,x2,y2;
X{
X      int xt[2],yt[2];
X      int it, k;
X      int clpxmi,clpxma,clpymi,clpyma;
X
X      gclp(&clpxmi,&clpxma,&clpymi,&clpyma);
X      k=0;
X      if (inside(x1,y1)) {
X        xt[k]=x1;
X        yt[k]=y1;
X        k=k+1;
X      }
X      if (inside(x2,y2)) {
X        xt[k]=x2;
X        yt[k]=y2;
X        k=k+1;
X      }
X      if (k == 2) goto draw;
X      if (y1 != y2) {
X        it=ixcut(clpymi);
X        if (inside(it,clpymi) && betw(clpymi,y1,y2)) {
X          xt[k]=it;
X          yt[k]=clpymi;
X          k=k+1;
X        }
X        if (k == 2) goto draw;
X        it=ixcut(clpyma);
X        if (inside(it,clpyma) && betw(clpyma,y1,y2)) {
X          xt[k]=it;
X          yt[k]=clpyma;
X          k=k+1;
X        }
X        if (k == 2) goto draw;
X      }
X      if (x1 != x2) {
X        it=iycut(clpxmi);
X        if (inside(clpxmi,it) && betw(clpxmi,x1,x2)) {
X          xt[k]=clpxmi;
X          yt[k]=it;
X          k=k+1;
X        }
X        if (k == 2) goto draw;
X        it=iycut(clpxma);
X        if (inside(clpxma,it) && betw(clpxma,x1,x2)) {
X          xt[k]=clpxma;
X          yt[k]=it;
X          k=k+1;
X        }
X        if (k == 2) goto draw;
X      }
Xdraw:
X      if (k == 2) genlin(xt[0],yt[0],xt[1],yt[1]);
X      scurr(x2,y2);
X} 
SHAR_EOF
echo "extracting plline.c"
sed 's/^X//' << \SHAR_EOF > plline.c
X#include "plplot.h"
X
Xvoid plline(n,x,y)
Xint n;
Xfloat x[],y[];
X{
X      int i,level;
X      glev(&level);
X      if (level<3) fatal("Please set up window before calling PLLINE.");
X      movwor(x[0],y[0]);
X      for(i=0; i<n; i++)
X        drawor(x[i],y[i]);
X}
SHAR_EOF
echo "extracting plmtex.c"
sed 's/^X//' << \SHAR_EOF > plmtex.c
X/* Prints out "text" at specified position relative to viewport       */
X/* (may be inside or outside)                                         */
X
X/* side   String which is one of the following:                       */
X/*        B or b  :  Bottom of viewport                               */
X/*        T or t  :  Top of viewport                                  */
X/*        L or l  :  Left of viewport                                 */
X/*        R or r  :  Right of viewport                                */
X/*        LV or lv : Left of viewport, vertical text                  */
X/*        RV or rv : Right of viewport, vertical text                 */
X/* disp   Displacement from specified edge of viewport, measured      */
X/*         outwards from the viewport in units of the current         */
X/*         character height. The CENTRELINES of the characters are    */
X/*         aligned with the specified position.                       */
X/* pos    Position of the reference point of the string relative      */
X/*         to the viewport edge, ranging from 0.0 (left-hand edge)    */
X/*         to 1.0 (right-hand edge)                                   */
X/* just   Justification of string relative to reference point         */
X/*         just = 0.0 => left hand edge of string is at reference     */
X/*         just = 1.0 => right hand edge of string is at reference    */
X/*         just = 0.5 => centre of string is at reference             */
X
X#include "plplot.h"
X
Xvoid plmtex(side,disp,pos,just,text)
Xchar side[], text[];
Xfloat disp, pos, just;
X{
X      int clpxmi, clpxma, clpymi, clpyma;
X      int sppxmi, sppxma, sppymi, sppyma;
X      int vert, refx, refy;
X      float shift, xform[4];
X      float vpdxmi, vpdxma, vpdymi, vpdyma;
X      float chrdef, chrht;
X      float mpxscl, mpxoff, mpyscl, mpyoff;
X
X      int level;
X      glev(&level);
X      if (level < 2) fatal("Please set up viewport before calling PLMTEX.");
X      
X      /* Open clip limits to subpage limits */
X
X      gclp(&clpxmi,&clpxma,&clpymi,&clpyma);
X      gspp(&sppxmi,&sppxma,&sppymi,&sppyma);
X      sclp(sppxmi,sppxma,sppymi,sppyma);
X
X      gvpd(&vpdxmi,&vpdxma,&vpdymi,&vpdyma);
X      gmp(&mpxscl,&mpxoff,&mpyscl,&mpyoff);
X      gchr(&chrdef,&chrht);
X
X      shift = 0.0;
X      if (just!=0.0) shift = just * plstrl(text);
X
X      if (strpos(side,'B')!=-1 || strpos(side,'b')!=-1) {
X        vert = 0;
X        refx = dcpcx(vpdxmi + (vpdxma-vpdxmi) * pos) - shift*mpxscl;
X        refy = mmpcy(dcmmy(vpdymi) - disp * chrht);
X      }
X      else if (strpos(side,'T')!=-1 || strpos(side,'t')!=-1) {
X        vert = 0;
X        refx = dcpcx(vpdxmi + (vpdxma-vpdxmi) * pos) - shift*mpxscl;
X        refy = mmpcy(dcmmy(vpdyma) + disp * chrht);
X      }
X      else if (stindex(side,"LV")!=-1 || stindex(side,"lv")!=-1) {   
X        vert = 0;
X        refy = dcpcy(vpdymi + (vpdyma-vpdymi) * pos);
X        refx = mmpcx(dcmmx(vpdxmi) - disp * chrht - shift);
X      }
X      else if (stindex(side,"RV")!=-1 || stindex(side,"rv")!=-1) {
X        vert = 0;
X        refy = dcpcy(vpdymi + (vpdyma-vpdymi) * pos);
X        refx = mmpcx(dcmmx(vpdxma) + disp * chrht - shift);
X      }
X      else if (strpos(side,'L')!=-1 || strpos(side,'l')!=-1) {   
X        vert = 1;
X        refy = dcpcy(vpdymi + (vpdyma-vpdymi) * pos) - shift*mpyscl;
X        refx = mmpcx(dcmmx(vpdxmi) - disp * chrht);
X      }
X      else if (strpos(side,'R')!=-1 || strpos(side,'r')!=-1) {
X        vert = 1;
X        refy = dcpcy(vpdymi + (vpdyma-vpdymi) * pos) - shift*mpyscl;
X        refx = mmpcx(dcmmx(vpdxma) + disp * chrht);
X      }
X      else {
X        sclp(clpxmi,clpxma,clpymi,clpyma);
X        return;
X      }
X
X      if (vert != 0) {
X        xform[0] = 0.0;
X        xform[1] = -1.0;
X        xform[2] = 1.0;
X        xform[3] = 0.0;
X      }
X      else {
X        xform[0] = 1.0;
X        xform[1] = 0.0;
X        xform[2] = 0.0;
X        xform[3] = 1.0;
X      }
X      plstr(0,xform,refx,refy,text);
X      sclp(clpxmi,clpxma,clpymi,clpyma);
X}
SHAR_EOF
echo "extracting plnxtv.c"
sed 's/^X//' << \SHAR_EOF > plnxtv.c
X/* Draws the next view of a 3-d plot. The physical coordinates of */
X/* the points for the next view are placed in the n points of arrays  */
X/* u and v. The silhouette found so far is stored in the heap as a */
X/* set of m peak points. */
X
X/* Define the heap size.  This can be increased if necessary.  Someday */
X/* I hope to just allocate a heap space of minimum required size.      */
X#define  heaps    4000
X
X#include <stdio.h>
X#include <stdlib.h>
X#include "plplot.h"
X#include "declare.h"
X
Xstatic int m, x, y;
X
Xvoid plnxtv(u,v,n,init)
Xint u[], v[], n, init;
X{
X      int xx;
X      int i, j, first;
X      int sx1, sx2, sy1, sy2;
X      int su1, su2, sv1, sv2;
X      int cx, cy, px, py;
X      int seg, ptold, lstold, pthi, pnewhi, newhi, change, ochange;
X      
X      first = 1;
X      pnewhi = 0;
X      
X      /* For the initial set of points, just display them and store them as */
X      /* the peak points. */
X
X      if (init == 1) {
X         if( heap3 == NULL)  /* heap not yet allocated so ... */
X            if(( heap3 = (int *)malloc(heaps*sizeof(int))) == NULL)
X               fatal("Not enough heap space in PLNXTV.");
X         y = heaps - n;
X         x = y - n;
X         movphy(u[0],v[0]);
X         heap3[x] = u[0];
X         heap3[y] = v[0];
X         for (i=1; i<n; i++){
X           draphy(u[i],v[i]);
X           heap3[x+i] = u[i];
X           heap3[y+i] = v[i];
X         }
X         m = n;
X         return;
X      }
X
X      /* Otherwise, we need to consider hidden-line removal problem. We scan */
X      /* over the points in both the old (i.e. heap3[x...] and heap3[y...]) and */
X      /* new (i.e. u[] and v[]) arrays in order of increasing x coordinate. */
X      /* At each stage, we find the line segment in the other array (if one */
X      /* exists) that straddles the x coordinate of the point. We */
X      /* have to determine if the point lies above or below the line segment, */
X      /* and to check if the below/above status has changed since the last */
X      /* point. */
X
X      /* If init = 2 we do not update the view, this is useful for drawing */
X      /* lines on the graph after we are done plotting points.  Hidden line */
X      /* removal is still done, but the view is not updated. */
X      xx = 0;
X      i = 0;
X      j = 0;
X
X      /* (heap3[x+i], heap3[y+i]) is the i'th point in the old array */
X      /* (u[j], v[j])           is the j'th point in the new array */
X
X      while (i < m || j < n) {
X
X        /* The coordinates of the point under consideration are (px,py). */
X        /* The line segment joins (sx1,sy1) to (sx2,sy2). */
X
X        /* "ptold" is true if the point lies in the old array. We set it */
X        /*  by comparing the x coordinates of the i'th old point */
X        /*  and the j'th new point, being careful if we have fallen past */
X        /*  the edges. Having found the point, load up the point and */
X        /*  segment coordinates appropriately. */
X 
X        ptold = ((heap3[x+i] < u[j] && i<m) || j>= n);
X        if (ptold) {
X          px = heap3[x+i];
X          py = heap3[y+i];
X          seg = j>0 && j<n;
X          if (seg) {
X            sx1 = u[j-1];
X            sy1 = v[j-1];
X            sx2 = u[j];
X            sy2 = v[j];
X          }
X        }
X        else {
X          px = u[j];
X          py = v[j];
X          seg = i>0 && i<m;
X          if (seg) {
X            sx1 = heap3[x+i-1];
X            sy1 = heap3[y+i-1];
X            sx2 = heap3[x+i];
X            sy2 = heap3[y+i];
X          }
X        }
X
X        /* Now determine if the point is higher than the segment, using the */
X        /* logical function "above". We also need to know if it is */
X        /* the old view or the new view that is higher. "newhi" is set true */
X        /* if the new view is higher than the old */
X
X        if(seg)
X            pthi = plabv(px,py,sx1,sy1,sx2,sy2);
X        else
X            pthi = 1;
X
X        newhi = (ptold && !pthi) || (!ptold && pthi);
X        change = (newhi && !pnewhi) || (!newhi && pnewhi);
X
X        /* There is a new intersection point to put in the peak array if the */
X        /* state of "newhi" changes */
X
X        if (first) {
X          movphy(px,py);
X          first = 0;
X          lstold = ptold;
X          if (init != 2) {
X            heap3[xx] = px;
X            heap3[xx+1] = py;
X            xx = xx + 2;
X          }
X          pthi = 0;
X          ochange = 0;
X        }
X        else if (change) {
X          /* Take care of special cases at end of arrays.  If init is 2 the */
X          /* endpoints are not connected to the old view. */
X          if (init==2 && ((!ptold && j==0)||(ptold && i==0))) {
X            movphy(px,py);
X            lstold = ptold;
X            pthi = 0;
X            ochange = 0;
X          }
X          else if (init==2 && ((!ptold && i>=m)||(ptold && j>=n))) {
X            movphy(px,py);
X            lstold = ptold;
X            pthi = 0;
X            ochange = 0;
X          }
X          /* If init is not 2 then we do want to connect the current line */
X          /* with the previous view at the endpoints.  */
X          /* Also find intersection point with old view. */
X          else  {
X            if (i == 0) {
X              sx1 = heap3[x];
X              sy1 = -1;
X              sx2 = heap3[x];
X              sy2 = heap3[y];
X            }
X            else if (i >= m) {
X              sx1 = heap3[x+m-1];
X              sy1 = heap3[y+m-1];
X              sx2 = heap3[x+m-1];
X              sy2 = -1;
X            }
X            else {
X              sx1 = heap3[x+i-1];
X              sy1 = heap3[y+i-1];
X              sx2 = heap3[x+i];
X              sy2 = heap3[y+i];
X            }
X
X            if (j == 0) {
X              su1 = u[0];
X              sv1 = -1;
X              su2 = u[0];
X              sv2 = v[0];
X            }
X            else if (j >= n) {
X              su1 = u[n-1];
X              sv1 = v[n-1];
X              su2 = u[n];
X              sv2 = -1;
X            }
X            else  {
X              su1 = u[j-1];
X              sv1 = v[j-1];
X              su2 = u[j];
X              sv2 = v[j];
X            }
X
X            /* Determine the intersection */
X            pl3cut(sx1,sy1,sx2,sy2,su1,sv1,su2,sv2,&cx,&cy);
X            if (cx == px && cy == py)  {
X              if (lstold && !ochange)
X                movphy(px,py);
X              else 
X                draphy(px,py); 
X
X              if (init != 2) {
X                heap3[xx] = px;
X                heap3[xx+1] = py;
X                xx = xx+2;
X              }
X              lstold = 1;
X              pthi = 0;
X            }
X            else {
X              if (lstold && !ochange)
X                 movphy(cx,cy);
X              else
X                 draphy(cx,cy);
X
X              lstold = 1;
X              if (init != 2) {
X                heap3[xx] = cx;
X                heap3[xx+1] = cy;
X                xx = xx+2;
X              }
X            }
X            ochange =1;
X          }
X        }
X
X        /* If point is high then draw plot to point and update view. */
X        if (pthi) {
X          if (lstold && ptold)
X            movphy(px,py);
X          else 
X            draphy(px,py); 
X
X          if (init != 2) {
X            heap3[xx] = px;
X            heap3[xx+1] = py;
X            xx = xx+2;
X          }
X          lstold = ptold;
X          ochange = 0;
X        }
X
X        pnewhi = newhi;
X
X        if (ptold)
X          i = i+1;
X        else
X          j = j+1;
X
X        if (xx>=x) fatal("Heap overflow in plnxtv.");
X      }
X
X      /* Transfer the peak points to the RHS of the heap */
X
X      if(init != 2) {
X        m = xx/2;
X        xx = 0;
X        y = heaps - m;
X        x = y - m;
X        for(i=0; i<m; i++) {
X          heap3[x+i] = heap3[xx];
X          heap3[y+i] = heap3[xx+1];
X          xx = xx+2;
X        }
X      }
X}
X
SHAR_EOF
echo "End of archive 6 (of 7)"
# if you want to concatenate archives, remove anything after this line
exit