Usenet 1994 January/usenetsourcesnewsgroupsinfomagicjanuary1994.iso/sources/unix/volume10/crc_plot/part01

Text File | 1987-07-07 | 26.4 KB | 1,002 lines

Path: uunet!rs
From: [e-mail] (Rich Salz)
Newsgroups: comp.sources.unix
Subject: v10i045:  CRC Plotting Package, Part01/06
Message-ID: <[e-mail]>
Date: 9 Jul 87 01:14:13 GMT
Organization: UUNET Communications Services, Arlington, VA
Lines: 991
Approved: [e-mail]

Submitted-by: "Wombat" <[e-mail]>
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, [e-mail], 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    [e-mail]
Moderator, comp.sources.unix    [e-mail]