rs@uunet.UU.NET (Rich Salz) (07/09/87)
Submitted-by: "Wombat" <rsk@j.cc.purdue.edu>
Posting-Number: Volume 10, Issue 45
Archive-name: crc_plot/Part01
[ No more inconsistant headers; I've got my posting tool back. --r$ ]
# This is a shell archive.
# Remove everything above and including the cut line.
# Then run the rest of the file through sh.
#----cut here-----cut here-----cut here-----cut here----#
#!/bin/sh
# shar: Shell Archiver
# Run the following text with /bin/sh to create:
# README
# Makefile
# crc.h
# examples
cat << \SHAR_EOF > README
The CRC plotting package is a device independent graphics system.
It includes subroutines for generating graphics which may be called
from Fortran or C, a two-dimensional plotting utility, and a
three-dimensional plotting utility.
The CRC package was originally developed at the Purdue University School
of Electrical Engineering by Carl Crawford; additional work has been
contributed by Mani Azimi and Malcolm Slaney, notably "plot3d". This
software has been in use locally for several years, and so most of the
obvious bugs have hopefully been caught and fixed. Although nobody's
willing to promise to fix future bugs immediately, it is not unreasonable
to assume that this package will continue to be supported, so please
do report bugs. (If you like, send them to me, and I'll forward them
to the folks at EE.) HOWEVER, no guarantees, folks.
This software has been developed on Vaxen running 4.XBSD; it's likely that
it will work on most machines running some variant of 4.XBSD. The two
user programs contained herein (qplot and plot3d) are probably of some use
to folks who need quick plots with reasonable resolution and labels
and annotation and so on without a lot of bother. Nice features of qplot
and plot3d include the ability to overlay multiple plots, tolerance of
a lot of different data formats, automatic or explicit scaling, logarithmic
plotting, ability to plot "bar graphs", and adjustable surface tilt (plot3d).
--
Rich Kulawiec, rsk@j.cc.purdue.edu, j.cc.purdue.edu!rsk
SHAR_EOF
cat << \SHAR_EOF > Makefile
#
# crc Makefile
#
# $Header: /usr/src/unsup/bin/crc/RCS/Makefile,v 2.0 87/02/28 13:47:25 ksb Exp $
#
# Richard S. Kulawiec, Purdue University Computing Center
# 9/26/86
#
SUBDIR= src lib
LOOP= for i in ${SUBDIR}; do\
echo $$i:;\
cd $$i;\
[ -f Makefile ] || co Makefile;\
make ${MFLAGS} DESTDIR=${DESTDIR} $@;\
cd ..;\
done
all clean depend install lint print source spotless tags: FRC
${LOOP}
FRC:
SHAR_EOF
cat << \SHAR_EOF > crc.h
#
/*
crc.h - include file for the CRC graphics package
Carl Crawford
Purdue University
W. Lafayette, IN 47907
Jan. 1981
*/
#include <stdio.h>
#include <math.h>
#include <signal.h>
unsigned short *_pic; /* pointer to bit plane */
int _xp,_yp; /* integer position */
float _axp,_ayp; /* real position */
float _xo,_yo; /* current origin */
int _ud; /* indicates up/down for pen */
int _error; /* indicates error in plotting */
float _fac; /* scale factor */
float _ipsz; /* size of the internal file - 1 */
float _ipsz10; /* ipsize / 10.0 */
int DEV; /* major device number */
char DEVN; /* minor device number */
int BLANK; /* 1 = don't blank device before plotting */
char *STORE; /* default storage file */
char *PLOTFILT; /* Plot Filter Name */
float TICDIS; /* distance between tic marks on the axis */
float HEIGHT; /* char height in axis routines */
int DIGITS; /* number of dec. digits + 1 in axis annotation */
unsigned _bufsize; /* size of point buffer */
char _abuf[100]; /* char buffer for anyone */
char *SITE; /* site for gplp */
FILE *_pipe_fd; /* file descriptor for pipes and pseudo pipes */
int (*_isig)(); /* save SIGINT signal */
int (*_qsig)(); /* save SIGQUIT signal */
int (*_hsig)(); /* save SIGHUP signal */
int _intty[3]; /* save current tty modes in here */
/* control characters */
#define NUL 0 /* <nul> */
#define SOH 1 /* <soh> */
#define STX 2 /* <stx> */
#define ETC 3 /* <etc> */
#define ETX 3 /* <etx> */
#define EOT 4 /* <eot> */
#define ENQ 5 /* <enq> */
#define ACK 6 /* <ack> */
#define BEL 7 /* <bel> */
#define BS 8 /* <bs> */
#define HT 9 /* <ht> */
#define LF 10 /* <lf> */
#define VT 11 /* <vt> */
#define FF 12 /* <ff> */
#define CR 13 /* <cr> */
#define SO 14 /* <so> */
#define SI 15 /* <si> */
#define DLE 16 /* <dle> */
#define DC1 17 /* <dc1> */
#define DC2 18 /* <dc2> */
#define DC3 19 /* <dc3> */
#define DC4 20 /* <dc4> */
#define NAK 21 /* <nak> */
#define SYN 22 /* <syn> */
#define ETB 23 /* <etb> */
#define CAN 24 /* <can> */
#define EM 25 /* <em> */
#define SUB 26 /* <sub> */
#define ESC 27 /* <esc> */
#define FS 28 /* <fs> */
#define GS 29 /* <gs> */
#define RS 30 /* <rs> */
#define US 31 /* <us> */
/* variables for HP and TEK */
int _CM; /* current mode */
int _X; /* x position */
int _Y; /* y position */
int _FILL; /* number of fill characters */
#define BINARY_FONT_FILE "/usr/unsup/lib/crc/font.5x7"
#define PLOTBIN "/usr/bin/plot"
#define BIT 0 /* major device table */
#define GOV 1
#define IMAGE 2
#define GGOV 3
#define GIMAGE 4
#define PLOT 5
#define TEK 6
#define HP 7
#define MBIT 4
/* maximum device in bit plane mode */
/*
Major and minor device tables
DEV DEVN dev OUTPUT
0 0 0 file or standard output
1 8 Versatec through gp (I)
2 16 Printronix through gplp (I) and opr (I)
1 0 1 Comtal graphics overlay 0(*)
1 9 Comtal graphics overlay 1(*)
2 17 Comtal graphics overlay 2(*)
2 0 2 Comtal image image displayed(*)
1 10 Comtal image 0(*)
2 18 Comtal image 1(*)
3 26 Comtal image 2(*)
3 0 3 Grinnell graphics overlay 0(*)
1 11 Grinnell graphics overlay 1(*)
2 19 Grinnell graphics overlay 2(*)
3 27 Grinnell graphics overlay 3(*)
4 0 4 Grinnell Image being Displayed (*)
1 12 Grinnell Image Plane 0(*)
2 20 Grinnell Image Plane 1(*)
3 28 Grinnell Image Plane 2(*)
4 36 Grinnell Image Plane 3(*)
5 44 Grinnell Image Plane 4(*)
5 0 5 Plot Subroutines
6 0 6 Tektronix through standard output
1 14 Retro-Graphics through standard output
2 22 Tektronix 4113
7 0 7 HP through /u/lib/graphics/hpd
(*) - through /u/lib/graphics/gd
*/
SHAR_EOF
mkdir examples
chdir examples
cat << \SHAR_EOF > ex1.c
/*
* C Example 1 - Compute and plot a sinc(r) function.
*
* Compile with
* cc ex1.c -lm -o ex1
*
* Run with
* ex1
*
* Look at the data by typing the command
* od -f data
*/
#include <stdio.h>
#include <math.h>
#define N 64
float z[N][N];
main(){
int i, j;
double x, y, r;
FILE *output;
output = fopen("data","w"); /* Open the output file */
if (!output){ /* And make sure the open succeeded */
fprintf(stderr,"Can't open data file for output.\n");
exit(1);
}
for (i=0;i<N;i++){ /* Increment the x-direction */
x = i - N/2;
for (j=0;j<N;j++){ /* Increment the y-direction */
y = j - N/2;
r = sqrt(x*x+y*y);
if ( r < .0001)
z[i][j] = 1.0;
else
z[i][j] = sin(r)/r;
}
}
/* Write the data out in binary
* format.
*/
fwrite(z,sizeof(z[0][0]),N*N,output);
}
SHAR_EOF
cat << \SHAR_EOF > ex2.c
/*
* C Example 2 - Compute and plot a sinc(x)*sinc(y) function.
*
* Compile with
* cc ex2.c -lm -o ex2
*
* Run with
* ex2 > data
*/
#include <stdio.h>
#include <math.h>
#define N 64
main(){
int i, j;
double x, y, z, sinc();
for (i=0;i<N;i++){ /* For each x */
x = i - N/2;
for (j=0;j<N;j++){ /* and for each y */
y = j - N/2;
z = sinc(x)*sinc(y);
printf("%f\n",z); /* Write out the value */
}
}
}
double sinc(x) /* Compute the sinc(x) */
double x;
{
if (fabs(x) < .0001)
return(1.0);
else
return(sin(x)/x);
}
SHAR_EOF
cat << \SHAR_EOF > ex3.c
/*
* Qplot C Example - Compute and plot a Gaussian Random Variable
*
* Compile with
* cc ex3.c -lm -o ex3
*
* Run with
* ex3
*/
#include <stdio.h>
#include <math.h>
main(){
int i;
double x, y, Gauss();
FILE *xfile, *yfile;
xfile = fopen("x","w"); /* Open the x and y files */
yfile = fopen("y","w");
if (!xfile || !yfile){ /* Check for Errors */
printf("Can't open files for output.\n");
exit(1);
}
for (i=0;i<100;i++){ /* Now compute 100 RV's */
x = Gauss()*3.0;
y = Gauss();
/* Check for out of bounds */
if (x < -1 || x > 1 || y < -1 || y > 1)
continue;
fprintf(xfile,"%f\n",x); /* Print out the values */
fprintf(yfile,"%f\n",y);
}
fclose(xfile); /* Close the files */
fclose(yfile);
}
#define NUM 25
/*
* Compute a Gaussian random
* variable by summing a number
* of uniformly distributed
* variables.
*
* The returned value will
* have a mean of 0.
*/
double
Gauss(){
int i;
float x;
x = 0;
for (i=0;i<NUM;i++)
x += (float)random();
/*
* Scale the sum by the
* maximum value from the
* random() subroutine and
* the number of RV's summed.
*/
return(x/((float)0x7fffffff*NUM/2) - 1.0);
}
SHAR_EOF
cat << \SHAR_EOF > a.f
c
c Plot3d Fortran Example - Compute and plot a semi-Gaussian function
c
c Compile with
c f77 ex.3 -lU77 -o ex3
c
c Run with
c ex3
c
c
c Set the different damping constants
c of the Gaussian function.
sigx1 = 7.0
sigx2 = 15.0
sigy = 5.0
c Open the z file.
open(unit=2,file="z",status="unknown",form="formatted")
c The damping constant is different for
c negative and positive values of x while
c it is constant along the y axis.
c Compute the function.
do 30 j=1 , 32
yfact = exp( - ( abs(j-7.0) / sigy ) ** 2 )
c Compute the first half of the Gaussian function.
do 10 i=1 , 32
tmp = exp( - ( abs(i-33.0) / sigx1 ) ** 2 ) * yfact
write(2,*)tmp
10 continue
c Compute the second half of the Gaussian
c function with a different damping constant.
do 20 i=33 , 64
tmp = exp( - ( abs(i-33.0) / sigx2 ) ** 2 ) * yfact
write(2,*)tmp
20 continue
30 continue
c Close the z file.
close(2)
stop
end
SHAR_EOF
cat << \SHAR_EOF > b.f
c
c Plot3d Fortran Example - Compute and plot a semi-Gaussian function
c for non-uniform values of x and y.
c
c Compile with
c f77 ex.4 -lU77 -o ex4
c
c Run with
c ex4
c
real tmp(32)
integer ucreat
c The damping constant is different for
c negative and positive values of x while
c it is constant along the y axis.
c Set the different damping constants
c of the Gaussian function.
sigx1 = 3.5
sigx2 = 7.5
sigy = 12.5
c Create the z file.
ifd = ucreat("z",420)
c Compute the function.
do 30 j=1 , 16
yfact = exp( - ( abs(j-3.5) / sigy ) ** 2 )
c Compute the first half of the Gaussian function.
do 10 i=1 , 16
tmp(i) = exp( - ( abs(i-16.5) / sigx1 ) ** 2 ) * yfact
10 continue
c Compute the second half of the Gaussian
c function with a different damping constant.
do 20 i=17 , 32
tmp(i) = exp( - ( abs(i-16.5) / sigx2 ) ** 2 ) * yfact
20 continue
call uwrite(ifd,tmp,4*32)
30 continue
c Open the file x.
open(unit=2,file="x",status="unknown",form="formatted")
c Compute the values of x.
x = - 10.0 * 16.0 * 0.5 * 0.5
c To make the sampling non-uniform,
c use the random function generater.
do 40 i=1 , 32
x = x + 10.0 * rand(13*i) * rand(17*i)
write(2,*)x
40 continue
c Close the file x.
close(2)
c Open the file x.
open(unit=3,file="y",status="unknown",form="formatted")
c Compute the values of y.
y = - 30.0 * 8.0 * 0.5 * 0.5
c To make the sampling non-uniform,
c use the random function generater.
do 50 i=1 , 16
y = y + 30.0 * rand(13*i) * rand(17*i)
write(3,*)y
50 continue
c Close the file y.
close(3)
stop
end
SHAR_EOF
cat << \SHAR_EOF > hide.f
subroutine hide(x,y,xg,g,xh,h,ng,maxdim,n1,nfns,
&xlnth,ylnth,xmin,deltax,xlabel,ymin,deltay,ylabel)
c
c n1 = the number of points plotted in a given call. If
c n1 < 0 Y vs X will be plotted in reverse order.
c x = a real array containing the horizontal coodinates.
c The contents of x must be in increasing order.
c x has dimension n1.
c y = a real array containing the vertical coordinates. y has
c dimension n1.
c g & xg = two real arrays that hold the current visual
c maxima.
c h & xh = two work arrays.
c ng = a non-positive integer.
c ng = 0 - draw the 8 1/2 x 11 border and plot visual
c maxima.
c ng = -1 - don't draw the 8 1/2 x 11 border but plot
c visual maxima.
c ng = -2 - draw the border and plot visual minima. This
c results in the "bottom view" of the graph.
c ng = -3 - don't draw the border but plot the visual
c minima.
c maxdim = the dimension of g, xg, h, xh. If the program is
c about to go out of bounds in these arrays maxdim will
c be returned as its negative. When the subroutine is
c called with maxdim < 0 it will immediately return.
c nfns = the total number of curves to be plotted. If a plot
c is desired with no shift then nfns is the negative of
c this number. nfns = 0 will plot the curve with the
c same ammount of shift as in the last call.
c xlnth = length (in inches) of the horizontal axis.
c ylnth = length (in inches) of the vertical axis.
c xmin = the minimum value of x.
c ymin = the minimun value of y.
c deltax = the x increment per inch.
c deltay = the y increment per inch.
c
dimension x(1), y(1), xg(1), g(1), h(1), xh(1)
character xlabel(1), ylabel(1)
c
c the only purpose of the following equivalence statement
c is to save storage.
c
equivalence (k1,iwhich), (k2,slope), (fnsm1,z1),
+ (iggp1,k1), (k1,n2)
c
c eps1 is the relative abcissa increment used to simulate
c discontinuities in the maximum function.
c
data eps1 /1.0e-3/
c
c the following statement function computes the ordinate on
c the line joining (xi,yi) and (xip1,yip1) corresponding to
c the abcissa xx.
c
f(xx,xi,yi,xip1,yip1) = yi + (xx - xi) * (yip1 - yi) /
+ (xip1 - xi)
if (maxdim.le.0) return
ifplot = 1
if (n1.gt.0) go to 76
n1 = -n1
ifplot = 0
76 do 71 i=2,n1
if (x(i-1).lt.x(i)) go to 71
maxdim = 0
write(6,1020)
1020 format('abcissa array not in increasing order')
go to 75
71 continue
if (ng.gt.0) go to 5000
if (n1 + 4.0.le.maxdim) go to 74
maxdim = -maxdim
75 return
c
c we want sign = 1 if we are plotting maximum, = -1 if
c minimum
c
74 sign = 1.0
if (ng.lt.-1) sign = -1.0
c
c the kth curve to be plotted will (optionally) be
c translated by the vector (-dxin,dyin) * (k - 1) to
c simulate stepping in the depth dimension.
c
fnsm1 = 0.0
if (nfns.le.0) go to 46
fnsm1 = nfns - 1
dxin = (9.0 - abs(xlnth)) * deltax / fnsm1
dyin = (6.0 - abs(ylnth)) * deltay / fnsm1
c
c systems routine plot moves the pen to a point whose
c coordinates are specified in inches by the first two
c parameters. the pen is picked up if the absolute value of
c the third parameter is 3, is put down if 2, and is left as
c after last call if 1. if the third parameter is negative,
c a new reference point will be established.
c
46 if (ng.eq.-1.or.ng.eq.-3) go to 41
c
c draw 8 1/2 by 11 inch border.
c
call plot(0.0,0.0,3)
call plot (11.0,0.0,2)
call plot (11.0,8.5,1)
call plot (00.0,8.5,1)
call plot (00.0,0.0,1)
call plot (01.625,2.0,-3)
c
c call systems routine to plot the 80-character title.
c the first two arguments are the coordinates in inches
c relative to the reference point of the lower left-hand
c corner of the first character. the third argument
c determines the height in inches of the characters. the
c fifth argument gives the angle relative to horizontal of
c the plotted characters.
c
41 if (xlnth.lt.0.0) go to 42
c
c call systems routine to draw the horizontal axis. the
c left end is specified in inches relative to the reference
c point by the first two arguments.
c
call axis (9.0 - xlnth, 0.0, xlabel, 0, xlnth,
*xmin, xlnth*deltax+xmin, 1)
if (ylnth.lt.0.0) go to 43
c
c draw the depth axis.
c
call plot (9.0 - xlnth, 0.0, 3)
call plot (0.0, 6.0 - ylnth, 2)
42 if (ylnth.lt.0.0) go to 43
c
c draw the vertical axis. the bottom point is specified in
c inches relative to the reference point by the first two
c arguments.
c
call axis (0.0, 6.0 - ylnth, ylabel, 1, ylnth, ymin,
*ylnth*deltay+ymin, 0)
c
c curves successively farther in the background will be
c plotted where they are not hidden by g vs xg. g vs xg
c will be updated each time a new curve is drawn and will be
c the visual maximum (or minimum) function of the curves
c already plotted.
c
43 indext = 3
do 3 j=1,n1
xg(indext) = x(j)
g(indext) = sign * y(j)
3 indext = indext + 1
c
c the following precautionary step is used in place of a
c test in subroutine lookup to see if the value for which we
c want an index is outside the table.
c the last xg value will be set equal to the last abcissa
c of the curve to be plotted in the next call to hide.
c
eps = eps1 * (abs(xmin) + abs(deltax))
ng = n1 + 4
xg(1) = -fnsm1 * dxin + xmin - abs(xmin) - abs(xg(3)) - 1.0
xg(2) = xg(3) - eps
xg(n1 + 3) = xg(n1 + 2) + eps
zz = ymin
if (sign.lt.0.0) zz = -ymin - 50.0 * deltay
g(1) = zz
g(2) = zz
g(n1+3) = zz
g(ng) = zz
c call systems routine to produce a line plot of
c (x(i), y(i), i=1,n1) - this is the curve farthest in the
c foreground.
c xstart is the x value at the reference point.
c
xstart = xmin - (9.0 - abs(xlnth)) * deltax
c
if(ifplot.eq.1) call pdatax(x,y,n1,xstart,deltax,ymin,deltay)
dxkk = 0.0
dykk = 0.0
relinc = deltax / deltay
xg(ng) = sign
return
c
c statement 5000 is reached if any except the curve farthest
c in the foreground is to be plotted.
c
5000 sign = xg(ng)
xg(ng) = x(n1)
c
c translate the arrays before plotting to simulate stepping
c in the depth dimension.
c
if (nfns) 52, 48, 49
49 dxkk = dxkk + dxin
dykk = dykk + dyin
48 do 4 j=1,n1
y(j) = sign * (y(j) + dykk)
4 x(j) = x(j) - dxkk
52 call lookup (x(1), xg(1), jj)
if (jj.ge.maxdim) go to 700
do 31 j=1,jj
xh(j) = xg(j)
31 h(j) = g(j)
ig = jj + 1
xh(ig) = x(1)
h(ig) = f(x(1), xg(jj), g(jj), xg(ig), g(ig))
c
c we will be making table lookups for an increasing sequence
c of numbers - therefore, we do not have to search from the
c first of the (xg and x) tables each time. hence indexg
c and indext.
c
indexg = jj
indext = 1
z1 = x(1)
f1 = h(ig) - y(1)
it = 2
jj = ig
if (h(ig).ge.y(1)) go to 32
if (jj.ge.maxdim) go to 700
jj = ig + 1
h(jj) = y(1)
xh(jj) = z1 + eps
32 last = 0
x1 = z1
c
c find the first zero, z2, of the function g-y to the right
c of z1.
c
1100 if (xg(ig).lt.x(it)) go to 1001
c
c do not jump if we are to look for a zero between x1 and
c x(i).
c
iwhich = 0
x2 = x(it)
f2 = f(x2, xg(ig - 1), g(ig - 1), xg(ig), g(ig)) - y(it)
it = it + 1
go to 1002
c
c come to 1001 if we are to look for a zero between x1 and
c xg(ig).
c
1001 x2 = xg(ig)
iwhich = 1
f2 = g(ig) - f(x2, x(it - 1), y(it - 1), x(it), y(it))
ig = ig + 1
c
c the function (g - y) has a zero z2 such that x1 le z2 le x2
c if and only if (g - y at x1) * (g - y at x2) le 0.
c (g - y is assumed, for plotting purposes, to be linear on
c each interval (x1, x2).)
c
1002 if (f1 * f2.gt.0.) go to 1005
if (f1.eq.f2) go to 1005
slope = (f2 - f1) / (x2 - x1)
igg = ig - 1 - iwhich
itt = it - 2 + iwhich
if (abs(slope * relinc) .gt. eps1) goto 1007
c
c if g and y differ imperceptibly (for plotting purposes)
c on the interval (x1, x2), set z2 = x2. this step prevents
c division by zero.
c
z2 = x2
go to 1006
c
c otherwise, compute the zero z2.
c
1007 z2 = x1 - f1 / slope
go to 1006
c
c if no zero was found between x1 and x2, continue the
c search for zeroes.
c
1005 x1 = x2
f1 = f2
if (it.le.n1) go to 1100
c
c if the end of the x table has been reached, consider the
c interval from the last zero found to the end of the x
c table (plot, update maximum function as indicated).
c
1008 last = 1
z2 = x(n1)
call lookup (z2, xg(indexg), igg)
igg = indexg + igg - 1
itt = n1 - 1
c
c it is necessary to plot y vs x on the interval (z1, z2)
c only if y is unhidden at each zz such that z1 lt zz lt z2.
c we choose zz near the left end of the interval for
c efficiency in the table lookup.
c note that it is more efficient to choose this value for zz
c than, say, 0.99 * x(indext) + 0.01 * x(indext + 1), which
c would eliminate one of the two table lookups, but would
c necessitate a test to determine if zz was between z1 and z2.
c
1006 zz = 0.99 * z1 + 0.01 * z2
call lookup (zz, x(indext), k1)
call lookup (zz, xg(indexg), k2)
k1 = k1 + indext - 1
k2 = k2 + indexg - 1
if (f(zz, x(k1), y(k1), x(k1 + 1), y(k1 + 1)).gt.
+ f(zz, xg(k2), g(k2), xg(k2 + 1), g(k2 + 1))) go to 7
c
c if y is hidden between z1 and z2, update the maximum
c function.
c for generality, the maximum function is updated even if
c this is the (nfns)th curve.
c
if (jj + igg - indexg.ge.maxdim) go to 700
if (indexg.eq.igg) go to 712
j1 = indexg + 1
do 12 i=j1,igg
jj = jj + 1
xh(jj) = xg(i)
12 h(jj) = g(i)
712 jj = jj +1
xh(jj) = z2
h(jj) = f(z2, xg(igg), g(igg), xg(igg + 1), g(igg + 1))
indexg = igg
indext = itt
go to 60
c
c if y is not hidden between z1 and z2, update the maximum
c function and plot.
c
7 ngraph = itt - indext + 2
if (jj + ngraph - 1.gt.maxdim) go to 700
n2 = jj
if (ngraph.eq.2) go to 9
j1 = indext + 1
do 11 i=j1,itt
jj = jj + 1
xh(jj) = x(i)
11 h(jj) = y(i)
9 jj = jj + 1
xh(jj) = z2
h(jj) = f(z2, x(itt), y(itt), x(itt + 1), y(itt + 1))
c
c call systems routine to produce line plot of
c (xh(i), h(i), i=n2, n2 + ngraph - 1).
c
if(ifplot.eq.1) call pdatax(xh(n2),h(n2),ngraph,xstart,deltax,
1 sign*ymin,sign*deltay)
c
indext = itt
indexg = igg
60 if (last.eq.1) go to 61
x1 = x2
f1 = f2
z1 = z2
c
c after plotting and/or updating the maximum function on the
c interval (z1, z2), search for the next zero if the end of
c the abcissa table xt has not been reached.
c
if (it.le.n1) go to 1100
go to 1008
c
c after y vs x has been plotted, finish updating and store
c the new maximum function.
c allow for the possibility that the previous maximum
c function extends to the right of the function just
c plotted.
c
61 if (xg(ng).le.xg(ng - 1)) ng = ng - 1
if (xg(ng).le.x(n1)) go to 33
if (jj + 3 + ng - igg.gt.maxdim) go to 700
xh(jj + 1) = xh(jj) + eps
jj = jj + 1
h(jj) = f(x(n1), xg(igg), g(igg), xg(igg + 1), g(igg + 1))
iggp1 = igg + 1
do 34 j = iggp1, ng
jj = jj + 1
xh(jj) = xg(j)
34 h(jj) = g(j)
33 ng = jj + 2
if (ng.gt.maxdim) go to 700
do 13 i=1,jj
g(i) = h(i)
13 xg(i) = xh(i)
c
c the following precautionary step is used in place of a
c test in subroutine lookup to see if the value for which we
c want an index is outside the table.
c the last xg value will be set equal to the last abcissa
c of the next curve to be plotted.
c
xg(jj + 1) = xg(jj) + eps
g(jj + 1) = ymin + dykk
if (sign.lt.0.) g(jj + 1) = -ymin - 50.0 * deltay + dykk
g(ng) = g(jj + 1)
c
c restore arrays x and y before returning.
c
66 if (nfns.lt.0) go to 53
do 82 i=1,n1
x(i) = x(i) + dxkk
82 y(i) = sign * y(i) - dykk
53 xg(ng) = sign
return
c
c if statement 700 is reached, dimensions would have been
c exceeded. see comments on calling sequence for hide.
c
700 maxdim = -maxdim
write(6,1030)
1030 format('visual maximum exceeds maxdim')
go to 66
end
subroutine lookup (x, xtbl, j)
c
c this subroutine is called by hide to perform a table
c lookup. because of precautions taken in hide, a test to
c see if x is outside the table is unnecessary.
c
dimension xtbl(1)
j = 2
4 if (xtbl(j) - x) 1, 2, 3
1 j = j + 1
go to 4
2 return
3 j = j - 1
return
end
subroutine pdatax(x,y,n,xm,dx,ym,dy)
c
c purdue calcomp/dipl compatable data ploting routine
c
dimension x(n),y(n)
data cx,cy/2*0.0/
px(i)=(x(i)-xm)/dx
py(i)=(y(i)-ym)/dy
c
i1 = 1
i2 = 1
if(amax1(abs(cx-px(1)),abs(cy-py(1))).lt.
1 amax1(abs(cx-px(n)),abs(cy-py(n)))) goto 1
i1 = n
i2 = -i2
1 call plot(px(i1),py(i1),3)
do 2 i3=2,n
i1=i1+i2
2 call plot(px(i1),py(i1),2)
cx = px(i1)
cy = py(i1)
return
end
SHAR_EOF
cat << \SHAR_EOF > pl3d.f
dimension xg(1400),g(1400),xh(1400),h(1400),x(1400)
integer uread,uopen
dimension pic(32,32),p(1024)
equivalence (pic(1,1),p(1))
maxdim = 1000
write(6,1000)
1000 format('hi carl')
npic = 32
npic2 = npic * npic
ng = 0
xl = 6.5
yl = 4.0
ifd = uopen('pic.r',0)
write(6,1000)
if(ifd .ne. -1) goto 20
write(6,10)
10 format(' can not open: pic.r')
stop
20 ic = uread(ifd,pic,4 * npic2)
write(6,1000)
if(ic .eq. 4 * npic2) goto 40
write(6,30)
30 format(' bad data structure')
stop
40 amax = p(1)
amin = amax
call plots(3,0)
do 50 i=1,npic2
a = p(i)
if(a .gt. amax) amax = a
if(a .lt. amin) amin = a
50 continue
dx = float(npic) / xl
rmin = amin
dy = (amax - amin) / yl
l = 1
do 60 i=1,npic
60 x(i) = i
do 70 i=1,npic
call hide(x,p(l),xg,g,xh,h,ng,maxdim,npic,npic,xl,yl
* ,1.0,dx,rmin,dy)
l = l + npic
70 continue
call plot(0.0,0.0,999)
stop
end
SHAR_EOF
chdir ..
# End of shell archive
exit 0
--
Rich $alz "Anger is an energy"
Cronus Project, BBN Labs rsalz@bbn.com
Moderator, comp.sources.unix sources@uunet.uu.net