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