adrian@milton.u.washington.edu (Adrian Mariano) (05/02/91)
Submitted-by: Adrian Mariano <adrian@milton.u.washington.edu>
Posting-number: Volume 19, Issue 12
Archive-name: plot3d/part01
This program plots 3d surfaces as wire meshes. It supports cartesian,
cylindrical and spherical coordinates. This version compiles under
MSDOS with Turbo C. It shouldn't be very difficult to port to other
environments, however.
Adrian Mariano
----------
#! /bin/sh
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
# calc.c
# graph.c
# makefile
# plot.c
# plot3d.c
# plot3d.h
# readme
# This archive created: Tue Apr 30 14:21:27 1991
export PATH; PATH=/bin:/usr/bin:$PATH
if test -f 'calc.c'
then
echo shar: "will not over-write existing file 'calc.c'"
else
cat << \SHAR_EOF > 'calc.c'
/*
Plot3d -- calc.c, the function evaluation code
By Karl Crary -- kc2m+@andrew.cmu.edu
Copyright (c) 1991 by Karl Crary
You may use and distribute this program as much as you like so long as you do
not charge for this service. I am not liable for failure of this program to
perform in any way.
*/
#include "plot3d.h"
int convert(char string[],char **errormsg,char stackno); /* convert infix to RPN in fstack */
void process(char operator,unsigned char *findex,char cstack[],char *cstacktop,term fstack[]);
char precedence(char operator, char push);
void simplify(char stackno);
double gamma(double arg);
extern char stop;
char ferr;
double ff(double arg1, double arg2, char stackno)
{
static unsigned char begin = 0;
long double stack[25], intermediate;
unsigned char findex;
char height, ii, recurse;
if ((arg1 < -MAXDOUBLE) || (arg1 > MAXDOUBLE)) {
ferr = 1;
return(0);
}
if (range(-MINDOUBLE,arg1,MINDOUBLE)) arg1 = 0;
ferr = 0;
height = recurse = 0;
if (begin) recurse = 1;
for ((findex = begin),(begin = 0); fstack[stackno][findex].operator != '='; findex++) {
if (!stop && kbhit()) {
ii = getch();
if (ii == ' ') stop = 2;
if (ii == 24) terminate();
if (ii == 27) {
stop = 1;
return(0);
}
}
switch (fstack[stackno][findex].operator) {
case '#' :
if (height >= 25) {
ferr = 1;
break;
}
for (ii = height; ii > 0; ii--) stack[ii] = stack[ii-1];
height++;
stack[0] = fstack[stackno][findex].operand;
break;
case 'y' :
if (height >= 25) {
ferr = 1;
break;
}
for (ii = height; ii > 0; ii--) stack[ii] = stack[ii-1];
height++;
stack[0] = arg2;
break;
case 'x' :
if (height >= 25) {
ferr = 1;
break;
}
for (ii = height; ii > 0; ii--) stack[ii] = stack[ii-1];
height++;
stack[0] = arg1;
break;
case '^' :
if (((stack[1] == 0) && (stack[0] <= 0)) || ((stack[1] < 0) && (stack[0] != floor(stack[0])))) ferr = 1;
else intermediate = pow(stack[1],stack[0]);
if (intermediate == HUGE_VAL) ferr = 1;
else stack[0] = intermediate;
goto popstack;
case '\\' : /* root */
if ((stack[0] == 0) && (stack[1] <= 0)) ferr = 1; /* division by zero */
if (stack[1] == 0) ferr = 1; /* zeroth root */
if ((stack[0] < 0) && (floor((stack[1]+1)/2) != (stack[1]+1)/2)) ferr = 1;
if (!ferr) stack[0] = (stack[0] < 0 ? -1 : 1)*pow(fabs(stack[0]),1/stack[1]);
goto popstack;
case '*' :
stack[0] = stack[0] * stack[1];
goto popstack;
case '/' :
if (stack[0] == 0) ferr = 1;
else stack[0] = stack[1] / stack[0];
goto popstack;
case '+' :
stack[0] = stack[0] + stack[1];
goto popstack;
case '-' :
stack[0] = stack[1] - stack[0];
goto popstack;
case 'm' :
stack[0] = -stack[0];
break;
case 's' :
stack[0] = sin(stack[0]);
break;
case 'c' :
stack[0] = cos(stack[0]);
break;
case 't' :
stack[0] = tan(stack[0]);
break;
case 'S' :
if (!range(-1,stack[0],1)) ferr = 1;
else stack[0] = asin(stack[0]);
break;
case 'C' :
if (!range(-1,stack[0],1)) ferr = 1;
else stack[0] = acos(stack[0]);
break;
case 'T' :
stack[0] = atan(stack[0]);
break;
case 1 : /* sinh */
stack[0] = sinh(stack[0]);
break;
case 2 : /* cosh */
stack[0] = cosh(stack[0]);
break;
case 3 : /* tanh */
stack[0] = tanh(stack[0]);
break;
case 4 : /* invsinh */
stack[0] = log(stack[0] + sqrt(pow(stack[0],2) + 1));
break;
case 5 : /* invcosh */
intermediate = pow(stack[0],2);
if (intermediate < 1) ferr = 1;
else {
intermediate = stack[0] + sqrt(intermediate - 1);
if (intermediate <= 0) ferr = 1;
else stack[0] = log(stack[0] + sqrt(pow(stack[0],2) - 1));
}
break;
case 6 : /* invtanh */
if (stack[0] == 1) ferr = 1;
else {
intermediate = (1 + stack[0])/(1 - stack[0]);
if (intermediate <= 0) ferr = 1;
else stack[0] = log(intermediate)/2;
}
break;
case 'l' :
if (stack[0] <= 0) ferr = 1;
else stack[0] = log10(stack[0]);
break;
case 'n' :
if (stack[0] <= 0) ferr = 1;
else stack[0] = log(stack[0]);
break;
case 'a' : /* absolute value */
stack[0] = fabs(stack[0]);
break;
case 'g' : /* gamma */
stack[0] = gamma(stack[0]);
break;
case 'f' : /* user function */
stack[0] = ff(stack[0],0,1);
break;
popstack:
height--;
for (ii = 1; ii < height; ii++) stack[ii] = stack[ii+1];
break;
}
if ((stack[0] < -2 * MAXDOUBLE) || (stack[0] > 2 * MAXDOUBLE)) ferr = 1; /* roundoff makes MAXDOUBLE > MAXDOUBLE */
if (ferr) return(0);
}
if (recurse) begin = findex; /* return where ended */
return(stack[0]);
}
int matherr(struct exception *e)
{
ferr = 1;
return(1);
}
double gamma(double arg)
{
double xx, total;
char ii;
/* magic gamma generating coefficients */
static double coef[7] = { 0.988205891,-0.897056937, 0.918206857,
-0.756704078, 0.482199394,-0.193527818,
0.035868343};
if (!stop && kbhit()) {
ii = getch();
if (ii == ' ') stop = 2;
if (ii == 24) terminate();
if (ii == 27) {
stop = 1;
return(0);
}
}
if (arg <= 0) {
if (arg == floor(arg)) {
ferr = 1;
return(0);
}
else {
arg = -arg + 1; /* recurse if x < 0 */
return(M_PI/(gamma(arg)*sin(M_PI*arg)));
}
}
else if (arg < 1) return(gamma(arg+1)/arg); /* recurse if x < 1 */
else if (arg <= 2) {
arg--;
total = 1 - 0.577191652 * arg;
xx = arg;
for (ii = 0; ii <= 6; ii++) {
xx *= arg;
total += coef[ii] * xx;
}
return(total);
}
else { /* arg > 2 */
xx = arg * arg; /* another approximation for larger x */
total = (arg - 0.5) * log(arg);
total -= arg;
total += log(2 * M_PI) / 2;
total += 1 / (12 * arg);
total -= 1 / (360 * arg * xx);
total += 1 / (1260 * arg * xx * xx);
total -= arg / (1680 * xx * xx * xx * xx);
return(exp(total));
}
}
/* errormessages */
char errorlist[6][45] = {"Values must precede each binary operator.",\
"Unmatched right parenthesis.",\
"Expression must end with value.",\
"Stack overflow.",\
"Function ends in operator.",\
"Unbalanced parentheses."};
int convert(char string[], char **errormsg, char stackno)
{
unsigned char findex;
char *cursor, value, inverse;
char cstack[100], cstacktop;
#define PROCESS(operator) process((operator),&findex,cstack,&cstacktop,fstack[stackno])
cursor = string;
findex = value = inverse = 0;
cstack[0] = '('; /* initialize cstack */
cstacktop = 0;
while (*cursor != 0) {
switch (*cursor) {
case '=' : /* explicit definition */
cursor++;
findex = value = inverse = 0; /* discard everything and start over */
cstack[0] = '(';
cstacktop = 0;
break;
case 'h' : /* hyperbolics already taken care of */
cursor++;
break;
case '0' : /* numbers */
case '1' :
case '2' :
case '3' :
case '4' :
case '5' :
case '6' :
case '7' :
case '8' :
case '9' :
case '.' :
if (value) PROCESS('$'); /* process implicit multiplication */
fstack[stackno][findex].operator = '#'; /* term is a number */
fstack[stackno][findex].operand = strtod(cursor,&cursor); /* store number */
findex++;
value = 1;
break;
case 'x' : /* variables */
case 'y' :
case 'o' : /* theta */
case 'r' :
case 'u' :
case 'p' : /* pi */
case 'e' : /* e */
case '!' : /* infinity */
if (value) PROCESS('$');
switch (*cursor) {
case 'p' :
fstack[stackno][findex].operator = '#';
fstack[stackno][findex].operand = M_PI;
break;
case 'e' :
fstack[stackno][findex].operator = '#';
fstack[stackno][findex].operand = M_E;
break;
case '!' :
fstack[stackno][findex].operator = '#';
fstack[stackno][findex].operand = MAXDOUBLE;
break;
case 'y' :
fstack[stackno][findex].operator = 'y';
break;
default : /* variables */
fstack[stackno][findex].operator = 'x';
break;
}
findex++;
cursor++;
value = 1;
break;
case '^' : /* operators */
case '\\' :
case '*' :
case '/' :
case '+' :
case '-' :
if (!value) { /* process unary operators */
if (*cursor == '+') { /* unary plus, do nothing */
cursor++;
break;
}
else if (*cursor == '-') { /* unary minus */
PROCESS('m');
cursor++;
break;
}
if (errormsg != NULL) (*errormsg) = errorlist[0];
return(cursor-string); /* return location of error */
}
PROCESS(*cursor); /* process operator */
cursor++;
value = 0;
break;
case '(' :
if (value) PROCESS('$'); /* process implicit multiplication */
PROCESS('('); /* process paren */
cursor++;
value = 0;
break;
case ')' :
if (!value) { /* expression ends in operator */
if (errormsg != NULL) *errormsg = errorlist[2];
return(cursor-string);
}
PROCESS(')');
cstacktop--; /* zap open paren */
if (cstacktop < 0) {
if (errormsg != NULL) (*errormsg) = errorlist[1];
return(cursor-string); /* unbalanced parens */
}
cursor++;
value = 1;
break;
case 'd' :
if (value) PROCESS('$'); /* process implicit multiplication */
fstack[stackno][findex].operator = 'd'; /* mark start of derivative */
findex++;
PROCESS('='); /* process end of derivative */
cursor++; /* will be followed by parenthesis */
value = 0;
break;
case 'A' :
cursor++;
/* inverse of non-trig prevented in IO */
inverse = 32;
case 's' :
case 'c' :
case 't' :
if (value) PROCESS('$'); /* process implicit multiplication */
if (cursor[1] != 'h') PROCESS(*cursor - inverse); /* inverse is lowercase */
else switch (*cursor) {
case 's' :
PROCESS(inverse ? 4 : 1);
break;
case 'c' :
PROCESS(inverse ? 5 : 2);
break;
case 't' :
PROCESS(inverse ? 6 : 3);
break;
} /* advance cursor on next pass */
cursor++;
value = inverse = 0;
break;
case ' ' :
cursor++;
break;
default : /* functions */
if (value) PROCESS('$'); /* process implicit multiplication */
PROCESS(*cursor);
cursor++;
value = 0;
break;
}
if ((cstacktop >= 100) || (findex >= 99)) {
if (errormsg != NULL) *errormsg = errorlist[3];
return(cursor-string); /* stack overflow */
}
}
if (!value) {
if (errormsg != NULL) *errormsg = errorlist[4];
return(cursor-1-string); /* function ends with operator, return location error */
}
PROCESS(')'); /* flush conversion stack */
if (cstacktop != 0) { /* unbalanced parens */
if (errormsg != NULL) *errormsg = errorlist[5];
return(cursor-string);
}
if (findex >= 99) {
if (errormsg != NULL) *errormsg = errorlist[3];
return(cursor-string); /* stack overflow or unbalanced parens */
}
fstack[stackno][findex].operator = '=';
return(-1);
}
void process(char operator,unsigned char *findex,char cstack[],char *cstacktop,term fstack[])
{
char pushprec;
pushprec = precedence(operator,1);
while (precedence(cstack[*cstacktop],0) >= pushprec) {
if (cstack[*cstacktop] == '$') cstack[*cstacktop] = '*';
fstack[*findex].operator = cstack[*cstacktop];
if (++(*findex) >= 99) break; /* stack overflow */
(*cstacktop)--;
}
if ((operator != ')') && (*cstacktop < 99)) {
(*cstacktop)++;
cstack[*cstacktop] = operator;
}
}
char precedence(char operator, char push)
{
switch (operator) {
case '(' :
return(push ? 10 : 0);
case ')' :
return(1);
case '=' : /* end of derivative */
return(9);
case '$' : /* implicit multiplication */
return(6);
case '\\' :
case '^' :
return(push ? 7 : 5);
case '*' :
case '/' :
return(3);
case '+' :
case '-' :
return(2);
default : /* functions */
return(push ? 8 : 4);
}
}
void simplify(char stackno)
{
int ii, jj, restart;
char last, prevlast;
restart = -1;
last = prevlast = 0;
for (ii = 0; fstack[stackno][ii].operator != '='; ii++) switch (fstack[stackno][ii].operator) {
case '+' : /* binary operators */
case '-' :
case '/' :
case '*' :
case '\\' :
case '^' :
if (restart == ii - 1) {
restart = ii;
break;
}
if (last && prevlast) {
/* load auxiliary function to calculate simplified value */
fstack[maxffs][0].operator = fstack[stackno][ii-2].operator; /* might not be # */
fstack[maxffs][0].operand = fstack[stackno][ii-2].operand;
fstack[maxffs][1].operator = fstack[stackno][ii-1].operator; /* might not be # */
fstack[maxffs][1].operand = fstack[stackno][ii-1].operand;
fstack[maxffs][2].operator = fstack[stackno][ii].operator;
fstack[maxffs][3].operator = '=';
/* simplify stack */
ii -= 2; /* update stack position */
fstack[stackno][ii].operand = ff(0,0,maxffs); /* calculate new value */
fstack[stackno][ii].operator = '#';
for (jj = ii+1; fstack[stackno][jj+1].operator != '='; jj++) { /* shift stack left by 2 */
fstack[stackno][jj].operator = fstack[stackno][jj+2].operator;
fstack[stackno][jj].operand = fstack[stackno][jj+2].operand;
}
last = prevlast = 0;
ii = restart; /* restart */
continue;
}
else last = 0; /* prevlast irrelevant */
break;
case '#' :
case 'p' :
case 'e' :
case '!' :
prevlast = last;
last = 1;
break;
case 'x' :
case 'y' :
restart = ii; /* no need to restart before unsimplifyable character */
last = 0; /* prevlast irrelevant */
break;
case 'm' :
if (fstack[stackno][ii-1].operator == 'm') { /* nuke double negate */
for (jj = ii-1; fstack[stackno][jj+1].operator != '='; jj++) { /* shift stack left by 2 */
fstack[stackno][jj].operator = fstack[stackno][jj+2].operator;
fstack[stackno][jj].operand = fstack[stackno][jj+2].operand;
}
last = prevlast = 0;
ii = restart;
continue;
}
/* no break */
default : /* unary operators */
if (restart == ii - 1) {
restart = ii;
break;
}
if (last) {
/* load auxiliary function to calculate simplified value */
fstack[maxffs][0].operator = fstack[stackno][ii-1].operator; /* might not be # */
fstack[maxffs][0].operand = fstack[stackno][ii-1].operand;
fstack[maxffs][1].operator = fstack[stackno][ii].operator;
fstack[maxffs][2].operator = '=';
/* simplify stack */
ii -= 1; /* update stack position */
fstack[stackno][ii].operand = ff(0,0,maxffs); /* calculate new value */
fstack[stackno][ii].operator = '#';
for (jj = ii+1; fstack[stackno][jj].operator != '='; jj++) { /* shift stack left */
fstack[stackno][jj].operator = fstack[stackno][jj+1].operator;
fstack[stackno][jj].operand = fstack[stackno][jj+1].operand;
}
last = prevlast = 0;
ii = restart; /* restart */
continue;
}
else last = 0; /* prevlast irrelevant */
break;
}
}
SHAR_EOF
fi
if test -f 'graph.c'
then
echo shar: "will not over-write existing file 'graph.c'"
else
cat << \SHAR_EOF > 'graph.c'
/*
Plot3d -- graph.c, the graphics code
By Adrian Mariano -- adrian@milton.u.washington.edu
Copyright (c) 1991 by Adrian Mariano
You may use and distribute this program as much as you like so long as you do
not charge for this service. I am not liable for failure of this program to
perform in any way.
*/
#include <graphics.h>
#include <conio.h>
#include <dos.h>
int returnmode, maxx, maxy, mode;
/* This function should return to text mode from graphics mode */
void backtotext()
{
restorecrtmode();
textmode(returnmode);
}
/* Enter graphics mode from text mode */
void entergraphics()
{
setgraphmode(mode);
cleardevice();
setfillstyle(1, 0); /* Sections will be black filled */
}
/* Do graphics initialization, including setting maxx and maxy to the
maximum x and y values */
void initgraphics()
{
int mm, driver = DETECT;
char *path;
struct text_info ti;
gettextinfo(&ti);
if (ti.screenheight == 50 || ti.screenheight == 43)
returnmode = C4350;
else
returnmode = C80;
path = "";
detectgraph((int far *) (&driver), (int far *) (&mode));
initgraph((int far *) (&driver), (int far *) (&mode), (char far *) (path));
if ((mm = graphresult()) != grOk) {
printf("Error: %d\n", mm);
exit(0);
}
maxx = getmaxx();
maxy = getmaxy();
backtotext();
}
/* Draw the quadrilateral specified by the four coordinate pairs, filled
in with a different color than the boundaries are drawn in */
#pragma argsused
void fill(int x1, int y1, int x2, int y2, int x3, int y3, int x4, int y4)
{
fillpoly(4, MK_FP(_SS, &x1));
}
/* Draw a line */
void drawline(int x1, int x2, int x3, int x4)
{
line(x1, x2, x3, x4);
}
SHAR_EOF
fi
if test -f 'makefile'
then
echo shar: "will not over-write existing file 'makefile'"
else
cat << \SHAR_EOF > 'makefile'
# Makefile for Turbo C++ 1.0
# Place your library directory below
LIB=c:\prog\c\lib
plot3d.exe: plot.obj calc.obj plot3d.obj graph.obj
tlink /c /x $(LIB)\c0h plot3d plot calc graph,plot3d,,\
$(LIB)\graphics $(LIB)\emu $(LIB)\mathh $(LIB)\ch
.c.obj:
bcc -mh -c -G $*.c
SHAR_EOF
fi
if test -f 'plot.c'
then
echo shar: "will not over-write existing file 'plot.c'"
else
cat << \SHAR_EOF > 'plot.c'
/*
Plot3d -- plot.c, the plotting code
By Adrian Mariano -- adrian@milton.u.washington.edu
Copyright (c) 1991 by Adrian Mariano
You may use and distribute this program as much as you like so long as you do
not charge for this service. I am not liable for failure of this program to
perform in any way.
*/
#include "plot3d.h"
/******************************************************************/
/* */
/* Parameters to plot3d() */
float xmin=-3; /* Dimensions of display area */
float xmax=3;
float ymin=-3;
float ymax=3;
float zmin=-3;
float zmax=3;
float aspect=1.0; /* Aspect ratio */
float distance=3; /* Perspective distance */
float tip=-30.0; /* Rotation around y axis */
float spin=0.0; /* Rotation around z axis */
float var1start;
float var2start;
float var1end;
float var2end;
long var1stepcount=30; /* Number of steps over var 1 */
long var2stepcount=30; /* Number of steps over var 2 */
/* */
/******************************************************************/
#define TRUE 1
#define FRONT 1
#define FALSE 0
#define BACK 0
#define PI 3.14159265
/*************************************
Internal variables
Global for convenience
*************************************/
struct sectiontype {
int x1,y1,x2,y2,x3,y3,x4,y4;
float x;
/* char ferr;*/
};
float xrange,yrange,zrange,sinspin,cosspin,sintip,costip;
float multx,multy,addx,addy;
/* The main function */
int plot3d(void (*transform)(float *,float *,float *,float,float,float));
/* Takes (x,y,z) and returns (x,y,z) in screen coordinates */
void realscale(float *xx,float *yy,float *zz);
/* Compare to sections -- used by qsort */
int sectioncmp(struct sectiontype *a,struct sectiontype *b);
/* plots the data contained the section */
void plotsections(struct sectiontype *section,int count);
/* Get screen coordinates from float (x,y,z) */
void scale(int *xint,int *yint,float x,float y,float z);
/* Draw 3d line */
void line3d(float a,float b,float c,float d,float e,float f);
/* Initialize multx, multy, addx, and addy for use in scaling */
void calcscale(float xmin,float ymin,float zmin,float xmax,float ymax,float zmax);
/* Draw the axes */
void drawaxes(char flag);
/* Is this face (of the cube) visible ? */
int visible(int i);
/* Function to be plotted */
#define FIX(flag,min,max) ( flag==1 ? max : min )
#define round(x) floor((x)+0.5)
typedef int(*function)(const void *,const void *);
int plot3d(void (*transform)(float *,float *,float *,float,float,float)) /* transform points to transform to cartesian */
{
int i,edgecount,steps,steps2; /* Counters */
float x,y,z; /* Cartesian coordinates */
float var1step,var2step,var1,var2; /* Values for independent vars */
struct sectiontype *section; /* Table of sections */
/* Get some memory */
section=farcalloc(var1stepcount+var2stepcount+3+(var1stepcount+1)*(var2stepcount+1),sizeof(struct sectiontype));
if (!section){
printf("Insufficient memory\n\n");
return 0;
}
xrange=xmax-xmin; /* Calculate widths for x, y, and z */
yrange=ymax-ymin;
zrange=zmax-zmin;
costip=cos(tip*PI/180); /* Precalculate trig functions for angles */
sintip=sin(tip*PI/180);
sinspin=sin(spin*PI/180);
cosspin=cos(spin*PI/180);
calcscale(xmin,ymin,zmin,xmax,ymax,zmax); /* Calculate scaling factors */
var1step=(var1end-var1start)/(var1stepcount);/* Set stepsizes for the */
var2step=(var2end-var2start)/(var2stepcount);/* independent variables */
i=0;
edgecount=0;
var1=var1start; /* Initialize var1 */
/* First make a partial table of needed boundary values. These are given */
/* very large x distances so that they will get sorted to the end of the */
/* section list. A second counter, 'edgecount' is kept. When the */
/* sections are plotted, the this value will be subtracted from the */
/* total counter so these edge sections won't be plotted */
/* Create a table of the edge values along the "top" of the function */
for(var2=var2start,steps=var2stepcount+1; steps ;steps--, i++,edgecount++,var2+=var2step) {
(*transform)(&x,&y,&z,ff(var1,var2-var2step,0),var1,var2-var2step);
realscale(&x,&y,&z);
(section+i)->x4=round(addx+multx*x); /* Bottom left point */
(section+i)->y4=round(addy+multy*y);
(section+i)->x=1e37; /* This section will get sorted to the end */
(*transform)(&x,&y,&z,ff(var1,var2,0),var1,var2);
realscale(&x,&y,&z);
(section+i)->x3=round(addx+multx*x); /* Bottom right point */
(section+i)->y3=round(addy+multy*y);
}
/* Main loop, cycles through var 1 */
for (var1=var1start+var1step,steps2=var1stepcount; steps2 ;steps2--, var1+=var1step) {
/* Calculate value for left edge of function */
var2=var2start;
(*transform)(&x,&y,&z,ff(var1,var2 ,0),var1,var2);
realscale(&x,&y,&z);
(section+i)->x3=round(addx+multx*x); /* Bottom right point */
(section+i)->y3=round(addy+multy*y);
(section+i)->x=1e37; /* Will get sorted to the end */
(*transform)(&x,&y,&z,ff(var1-var1step,var2,0),var1-var1step,var2);
realscale(&x,&y,&z);
(section+i)->x2=round(addx+multx*x); /* Top right point */
(section+i)->y2=round(addy+multy*y);
edgecount++;
i++;
/* Main loop. Calculate sections for one cycle of var 2 */
for(var2=var2start+var2step,steps=var2stepcount; steps ; i++ , steps--, var2+=var2step) {
(*transform)(&x,&y,&z,ff(var1,var2,0),var1,var2);
realscale(&x,&y,&z);
(section+i)->x3=round(addx+multx*x); /* Bottom right point */
(section+i)->y3=round(addy+multy*y);
(section+i)->x=z;
/* These points have already been calculated */
(section+i)->x4=(section+i-1)->x3; /* Bottom left point */
(section+i)->y4=(section+i-1)->y3;
(section+i)->x2=(section+i-var2stepcount-1)->x3; /* Top right */
(section+i)->y2=(section+i-var2stepcount-1)->y3;
(section+i)->x1=(section+i-var2stepcount-1)->x4; /* Top left */
(section+i)->y1=(section+i-var2stepcount-1)->y4;
}
printf("."); /* The it hasn't crashed indicator */
}
printf("\nSorting...");
qsort(section,i,sizeof(struct sectiontype),(function)sectioncmp);
printf("\done");
entergraphics();
i-=edgecount; /* Subtract the edges out */
drawaxes(BACK); /* Draw back part of axes */
plotsections(section,i);
drawaxes(FRONT); /* And plot those sections */
farfree(section);
return 1;
}
void realscale(float *xx,float *yy,float *zz) {
float x,y,z,dist;
x=(*xx/xrange)-xmin/xrange/2-xmax/xrange/2;
y=(*yy/yrange)-ymin/yrange/2-ymax/yrange/2;
z=(*zz/zrange)-zmin/zrange/2-zmax/zrange/2;
*zz=(x*cosspin-y*sinspin)*costip-z*sintip;
if (distance==0.0) dist=1;
else dist=distance+0.5-*zz;
*xx=(y*cosspin+x*sinspin)*aspect/dist;
*yy=(z*costip+(x*cosspin-y*sinspin)*sintip)/dist;
}
/* Compare to sections for qsort() */
int sectioncmp(struct sectiontype *a,struct sectiontype *b)
{
if (a->x < b->x) return(-1);
if (a->x == b->x) return(0);
return(1);
}
void plotsections(struct sectiontype *section,int count) /* Plot count sections */
{
int i;
for(i=0;i<count;i++)
fill((section+i)->x1,(section+i)->y1,(section+i)->x2,(section+i)->y2,
(section+i)->x3,(section+i)->y3,(section+i)->x4,(section+i)->y4);
}
void scale(int *xint,int *yint,float x,float y,float z) {
realscale(&x,&y,&z);
*xint=round(addx+multx*x);
*yint=round(addy+multy*y);
}
void line3d(float a,float b,float c,float d,float e,float f) {
int x1,y1,x2,y2;
scale(&x1,&y1,a,b,c);
scale(&x2,&y2,d,e,f);
drawline(x1,y1,x2,y2);
}
/* Determine scaling constants */
void calcscale(float xmin,float ymin,float zmin,float xmax,float ymax,float zmax)
{
float xs[2],ys[2],zs[2],yb=-1e37,ym=1e37,xm=1e37,xb=-1e37;
char i,j,k;
multx=1;
addx=0;
multy=1;
addy=0;
xs[0]=xmin;
xs[1]=xmax;
ys[0]=ymin;
ys[1]=ymax;
zs[0]=zmin;
zs[1]=zmax;
for(i=0;i<2;i++)
for(j=0;j<2;j++)
for(k=0;k<2;k++){
xmin=xs[i];
ymin=ys[j];
zmin=zs[k];
realscale(&xmin,&ymin,&zmin);
if (xmin>xb) xb=xmin;
else if (xmin<xm) xm=xmin;
if (ymin>yb) yb=ymin;
else if (ymin<ym) ym=ymin;
}
multx=((float)maxx)/(xb-xm);
multy=((float)maxy)/(yb-ym);
if (multx>multy) multx=multy; else multy=multx;
multy=-multy;
addx=-xm*multx+((float)maxx-multx*(xb-xm))/2;
addy=(float)maxy-ym*multy-((float)maxy+multy*(yb-ym))/2;
}
static int faces[6][3]={{0,0,-1},{-1,0,0},{0,1,0},{1,0,0},{0,-1,0},{0,0,1}};
static int edges[12][8]={
{0,1,-1,-1,-1,-1,1,-1},
{0,2,-1,1,-1,1,1,-1},
{0,3,1,-1,-1,1,1,-1},
{0,4,-1,-1,-1,1,-1,-1},
{4,1,-1,-1,-1,-1,-1,1},
{1,2,-1,1,-1,-1,1,1},
{2,3,1,1,-1,1,1,1},
{3,4,1,-1,-1,1,-1,1},
{5,1,-1,-1,1,-1, 1,1},
{5,2,-1, 1,1, 1, 1,1},
{5,3, 1,-1,1, 1, 1,1},
{5,4,-1,-1,1, 1,-1,1}};
void drawaxes(char flag) /* Draw the axes */
{
int i;
for(i=0;i<12;i++)
if ((visible(edges[i][0]) && visible(edges[i][1]))==flag)
line3d(FIX(edges[i][2],xmin,xmax), FIX(edges[i][3],ymin,ymax),
FIX(edges[i][4],zmin,zmax), FIX(edges[i][5],xmin,xmax),
FIX(edges[i][6],ymin,ymax), FIX(edges[i][7],zmin,zmax));
}
int visible(int i) /* determine if a face of the coordinate cube is visible */
{
return(faces[i][2] * sintip-(faces[i][0]*cosspin-faces[i][1]*sinspin)*costip<0);
}
SHAR_EOF
fi
if test -f 'plot3d.c'
then
echo shar: "will not over-write existing file 'plot3d.c'"
else
cat << \SHAR_EOF > 'plot3d.c'
/*
Plot3d -- plot3d.c, the main program
By Adrian Mariano -- adrian@milton.u.washington.edu
Copyright (c) 1991 by Adrian Mariano
You may use and distribute this program as much as you like so long as you do
not charge for this service. I am not liable for failure of this program to
perform in any way.
*/
#include "plot3d.h"
#include <string.h>
extern unsigned _stklen = 32000;/* stack space to make qsort happy */
/* status variables */
char stop;
void terminate(void)
{
backtotext();
exit(0);
}
term fstack[maxffs + 1][100]; /* function stack for ff evaluation */
typedef struct {
float v1start, v1end, v2start, v2end;
char *name;
} bound;
bound *currentbound;
/* Functions to be passed to plot3d() for different plot types */
void sphererho(float *x, float *y, float *z, float rho, float phi, float theta);
void spherephi(float *x, float *y, float *z, float phi, float rho, float theta);
void spheretheta(float *x, float *y, float *z, float theta, float phi, float rho);
void cartz(float *x, float *y, float *z, float zval, float xval, float yval);
void cartx(float *x, float *y, float *z, float xval, float yval, float zval);
void carty(float *x, float *y, float *z, float yval, float xval, float zval);
void cylinderz(float *x, float *y, float *z, float zval, float rval, float theta);
void cylinderr(float *x, float *y, float *z, float rval, float theta, float zval);
void cylindertheta(float *x, float *y, float *z, float theta, float rval, float zval);
void parametric(float *x, float *y, float *z, float dummy, float u, float v);
bound sphererho_bound = {0, M_PI, 0, PI2, "sphererho"};
bound spherephi_bound = {0, 3, 0, PI2, "spherephi"};
bound spheretheta_bound = {0, M_PI, 0, 3, "spheretheta"};
bound cartz_bound = {-3, 3, -3, 3, "cartz"};
bound cartx_bound = {-3, 3, -3, 3, "cartx"};
bound carty_bound = {-3, 3, -3, 3, "carty"};
bound cylinderz_bound = {0, 3, 0, PI2, "cylinderz"};
bound cylinderr_bound = {0, PI2, -3, 3, "cylinderr"};
bound cylindertheta_bound = {0, 3, -3, 3, "cylindertheta"};
char *xstr = "x";
char *ystr = "y";
char *zstr = "z";
char *thetastr = "theta";
char *rhostr = "rho";
char *phistr = "phi";
char *rstr = "r";
char *var1str;
char *var2str;
void (*plottype) (float *, float *, float *, float, float, float);
struct {
char *st, *en;
} replaces[] = {
{"sinh", "sh"}, {"cosh", "ch"}, {"tanh", "th"},{"sin", "s"},
{"cos", "c"}, {"arc", "A"}, {"ln", "n"}, {"log", "l"},{"abs","a"}
};
void setbounds(bound * newbound)
{
var1start = newbound->v1start;
var2start = newbound->v2start;
var1end = newbound->v1end;
var2end = newbound->v2end;
currentbound = newbound;
}
void substitute(char *str, char *start, char *end)
{
while (str = strstr(str, start)) {
strnset(str, ' ', strlen(start));
strncpy(str, end, strlen(end));
}
}
void fixfun(char *fun)
{
int i;
substitute(fun, var1str, "X");
substitute(fun, var2str, "Y");
substitute(fun, "X", "x");
substitute(fun, "Y", "y");
for (i = 0; i < sizeof(replaces) / (2 * sizeof(char *)); i++)
substitute(fun, replaces[i].st, replaces[i].en);
}
void doreplot()
{
if (plot3d(plottype)) {
getch();
backtotext();
}
}
void doplot(char *params)
{
char *errmesg;
int where;
fixfun(params);
if ((where = convert(params, &errmesg, 0)) == -1) {
simplify(0);
if (plot3d(plottype)) {;
getch();
backtotext();
}
} else {
while (where--)
putchar(' ');
printf(" ^\n%s\n", errmesg);
}
}
float getval(char *expr, float def)
{
int where;
char *errmesg;
if ((where = convert(expr, &errmesg, maxffs)) == -1)
return ff(def, def, maxffs);
else {
while (where--)
putchar(' ');
printf(" ^\n%s\n", errmesg);
return def;
}
}
void dotype(char *newtype)
{
if (!strcmp(newtype, "cartx")) {
plottype = cartx;
setbounds(&cartx_bound);
var1str = ystr;
var2str = zstr;
} else if (!strcmp(newtype, "carty")) {
plottype = carty;
setbounds(&carty_bound);
var1str = xstr;
var2str = zstr;
} else if (!strcmp(newtype, "cartz")) {
plottype = cartz;
setbounds(&cartz_bound);
var1str = xstr;
var2str = ystr;
} else if (!strcmp(newtype, "sphererho")) {
plottype = sphererho;
setbounds(&sphererho_bound);
var1str = phistr;
var2str = thetastr;
} else if (!strcmp(newtype, "spherephi")) {
plottype = spherephi;
setbounds(&spherephi_bound);
var1str = rhostr;
var2str = thetastr;
} else if (!strcmp(newtype, "spheretheta")) {
plottype = spheretheta;
setbounds(&spheretheta_bound);
var1str = phistr;
var2str = rhostr;
} else if (!strcmp(newtype, "cylinderz")) {
plottype = cylinderz;
setbounds(&cylinderz_bound);
var1str = rstr;
var2str = thetastr;
} else if (!strcmp(newtype, "cylinderr")) {
plottype = cylinderr;
setbounds(&cylinderr_bound);
var1str = thetastr;
var2str = zstr;
} else if (!strcmp(newtype, "cylindertheta")) {
plottype = cylindertheta;
setbounds(&cylindertheta_bound);
var1str = rstr;
var2str = zstr;
} else
printf("Unknown coordinate system: %s\n\n", newtype);
}
void doshow()
{
printf("Coordinate system: %s\n",currentbound->name);
printf("%s range: [%.17g,%.17g]\n%s range: [%.17g,%.17g]\n",
var1str, var1start, var1end, var2str, var2start, var2end);
printf("X: [%.17g,%.17g]\nY: [%.17g,%.17g]\nZ: [%.17g,%.17g]\n",
xmin, xmax, ymin, ymax, zmin, zmax);
printf("%s steps: %d\t%s steps: %d\n", var1str,
(int) var1stepcount, var2str, (int) var2stepcount);
printf("Aspect: %.17g\tDistance: %.17g\n", aspect, distance);
printf("Spin: %.17g \tTip: %.17g\n", spin, tip);
}
void dohelp()
{
printf(
"Commands:\n\n"
"plot <function> -- function is in terms of the current coordinate system\n"
"replot -- plot the same function again.\n"
"type <cartz|cartx|carty|sphererho|spheretheta|spherephi|cylinderz|\n"
" cylindertheta|cylinderr> -- Selects coordinate system for plotting\n"
" where the variable name specified (sphereRHO) is the dependent variable\n"
"show -- displays current plotting parameters\n"
"xmin, xmax, ymin, ymax, zmin, zmax <value> -- set region to display on screen\n"
"%sstart, %send, %sstart, %send, %ssteps, %ssteps <value> --\n"
" Set the range for the independent variables, and the number of steps at\n"
" which to sample the function\n"
"aspect <value> -- set aspect ratio of screen\n"
"distance <value> -- set perspect distance, 0 for no perspective\n"
"spin, tip <value> -- set view angle in degrees\n\n"
"Functions supported:\n"
" [arc]sin, [arc]cos, [arc]tan, [arc]sinh, [arc]cosh, [arc]tanh\n"
" ln, log, abs\n",var1str,var1str,var2str,var2str,var1str,var2str);
}
#define IS(XX) (!strcmp(command,XX))
void main()
{
char input[100], command[100], params[100];
char xcount[15], ycount[15], xstart[15], xend[15], ystart[15], yend[15];
initgraphics();
plottype = cartz;
setbounds(&cartz_bound);
var1str = xstr;
var2str = ystr;
printf(" -- Plot3d Version 1.0 by Adrian Mariano --\n"
" -- Type 'help' for help --\n");
do {
printf("> ");
gets(input);
*command = *params = 0;
sscanf(input, "%s %[^\n\r]", command, params);
/* printf("'%s' '%s'\n",command,params); */
strcpy(xcount, var1str);
strcat(xcount, "steps");
strcpy(ycount, var2str);
strcat(ycount, "steps");
strcpy(xstart, var1str);
strcat(xstart, "start");
strcpy(ystart, var2str);
strcat(ystart, "start");
strcpy(xend, var1str);
strcat(xend, "end");
strcpy(yend, var2str);
strcat(yend, "end");
if IS ("quit") break;
else if IS ("plot") doplot(params);
else if IS ("type") dotype(params);
else if IS ("show") doshow();
else if IS ("help") dohelp();
else if IS ("zmin") zmin = getval(params, zmin);
else if IS ("zmax") zmax = getval(params, zmax);
else if IS ("ymin") ymin = getval(params, ymin);
else if IS ("xmin") xmin = getval(params, xmin);
else if IS ("ymax") ymax = getval(params, ymax);
else if IS ("xmax") xmax = getval(params, xmax);
else if IS (xcount) var1stepcount = getval(params, var1stepcount);
else if IS (ycount) var2stepcount = getval(params, var2stepcount);
else if IS (xstart) {
currentbound->v1start = getval(params, var1start);
setbounds(currentbound);
}
else if IS (ystart) {
currentbound->v2start = getval(params, var2start);
setbounds(currentbound);
}
else if IS (xend) {
currentbound->v1end = getval(params, var1end);
setbounds(currentbound);
}
else if IS (yend) {
currentbound->v2end = getval(params, var2end);
setbounds(currentbound);
}
else if IS ("replot") doreplot();
else if IS ("aspect") aspect = getval(params, aspect);
else if IS ("tip") tip = getval(params, tip);
else if IS ("spin") spin = getval(params, spin);
else if IS ("distance") distance = getval(params, distance);
else if (!*command);
else
printf("Unknown command: %s\n", command);
} while (1);
clrscr();
}
/* A transform function for sphere notation */
void sphererho(float *x, float *y, float *z, float rho, float phi, float theta)
{
*x = rho * sin(phi);
*y = *x;
*x *= cos(theta);
*y *= sin(theta);
*z = rho * cos(phi);
}
void spherephi(float *x, float *y, float *z, float phi, float rho, float theta)
{
*x = rho * sin(phi);
*y = *x;
*x *= cos(theta);
*y *= sin(theta);
*z = rho * cos(phi);
}
void spheretheta(float *x, float *y, float *z, float theta, float phi, float rho)
{
*x = rho * sin(phi);
*y = *x;
*x *= cos(theta);
*y *= sin(theta);
*z = rho * cos(phi);
}
void cartz(float *x, float *y, float *z, float zval, float xval, float yval)
{
*x = xval;
*y = yval;
*z = zval;
}
void cartx(float *x, float *y, float *z, float xval, float yval, float zval)
{
*x = xval;
*y = yval;
*z = zval;
}
void carty(float *x, float *y, float *z, float yval, float xval, float zval)
{
*x = xval;
*y = yval;
*z = zval;
}
void cylinderz(float *x, float *y, float *z, float zval, float rval, float theta)
{
*x = rval * cos(theta);
*y = rval * sin(theta);
*z = zval;
}
void cylinderr(float *x, float *y, float *z, float rval, float theta, float zval)
{
*x = rval * cos(theta);
*y = rval * sin(theta);
*z = zval;
}
void cylindertheta(float *x, float *y, float *z, float theta, float rval, float zval)
{
*x = rval * cos(theta);
*y = rval * sin(theta);
*z = zval;
}
#pragma argsused
void parametric(float *x, float *y, float *z, float dummy, float u, float v)
{
/* *x =1.7*cos(2*PI*u); y =1.7*sin(2*PI*u); z =v;
*
* x = 0.5*v*cos(2*PI*u); y = 0.5*v*sin(2*PI*u); z = 0.866*v;
*
* x = cos(2*PI*u)*sin(PI*v); y = sin(2*PI*u)*sin(PI*v); z = cos(PI*v);
*
* x = cos(PI*(2*u-1))*sin(PI*v); y = sin(PI*(2*u-1))*sin(PI*v); z = cos(PI*v);
*
* x = exp(u+v); y = exp(u*v); z = u*v*v;
*
* x = -v + 1; y = u - v; z = v; */
float a = 2.0;
float b = 0.8;
*x = (a + b * cos(u)) * cos(v);
*y = (a + b * cos(u)) * sin(v);
*z = b * sin(u);
}
SHAR_EOF
fi
if test -f 'plot3d.h'
then
echo shar: "will not over-write existing file 'plot3d.h'"
else
cat << \SHAR_EOF > 'plot3d.h'
/* header files */
#include <conio.h>
#include <dos.h>
#include <stdio.h>
#include <math.h>
#include <alloc.h>
#include <stdlib.h>
#include <values.h>
/* constants */
#define maxffs 3
#define PI2 6.28318530
/* macros */
#define range(min,value,max) (((min) <= (value)) && ((value) <= (max)))
/* typedefs */
typedef struct {
char operator;
double operand;
} term;
extern char ferr;
extern term fstack[maxffs+1][100]; /* function stack for ff evaluation */
extern int maxx,maxy,mode;
extern float xmin; /* Dimensions of display area */
extern float xmax;
extern float ymin;
extern float ymax;
extern float zmin;
extern float zmax;
extern float aspect; /* Aspect ratio */
extern float distance; /* Perspective distance */
extern float tip; /* Rotation around y axis */
extern float spin; /* Rotation around z axis */
extern float var1start;
extern float var2start;
extern float var1end;
extern float var2end;
extern long var1stepcount; /* Number of steps over var 1 */
extern long var2stepcount; /* Number of steps over var 2 */
extern double ff(double arg1,double arg2, char stackno); /* function evaluator */
extern int plot3d(void (*transform)(float *,float *,float *,float,float,float));
extern void backtotext();
extern void entergraphics();
extern void initgraphics();
extern void fill(int x1,int y1,int x2,int y2,int x3,int y3,int x4,int y4);
extern void drawline(int x1,int x2,int x3,int x4);
SHAR_EOF
fi
if test -f 'readme'
then
echo shar: "will not over-write existing file 'readme'"
else
cat << \SHAR_EOF > 'readme'
Plot3d
by Adrian Mariano
and Karl Crary
Copyright 1991 by Adrian Mariano
You may use and distribute this program as much as you like so long as you do
not charge for this service. I am not liable for failure of this program to
perform in any way.
Please send any comments, patches, or improvements to me at:
adrian@milton.u.washington.edu
Plot3d is a program for plotting 3 dimensional surfaces specified via one of
the three standard coordinate systems.
It uses a command driven interface with the following commands:
type <type><dependent var>
Selects coordinate system for plotting. The three available coordinate
systems are: cartesian, a right handed orthoginal coordinate system with
x, y and z as the variables; cylindrical, polar coordinate system with r
the orthogonal distance from the z axis, theta the angle around the z
axis, and z the distance along the z axis; and spherical, where rho is
the distance from the origin, theta is the angle around the z axis, and
phi is the angle between the z axis and the point.
The <type> above should be one of cart, cylinder, or sphere.
The dependent variable is the one which depends on the other two. It
should be specified immediately after the type (without spaces).
The default is cartz, which allows you to plot surfaces of the form
z=f(x,y).
To plot phi=f(rho, theta) in spherical coordinates, use the command
'type spherephi'
plot <function>
Plots the function with the current settings. The function should be
specified WITHOUT an equals sign.
So x^2+y^2 is ok. z=x^2+y^2 is not ok.
Functions use the standard operations, +, -, /, *, and ^ (for exponents).
e give 2.71828..., and pi produces the value pi.
The '*' for multiplication may be omitted.
The following functions are supported:
[arc]sin, [arc]cos, [arc]tan, [arc]sinh, [arc]cosh, [arc]tanh
ln, log, abs
replot
plots the same function again with the current settings (which may be
different from when the function was last plotted).
show
displays current plotting parameters
xmin, xmax, ymin, ymax, zmin, zmax <number>
set the region to display on screen.
<var>start, <var>end, <var>steps
The program maintains separate values for the ranges over when the
independent variables should vary, and the number of steps to take over
this range for each coordinate system. You can set these values for the
current coordinate system only by using these commands, where <var> is
the variable you want to affect.
aspect <number>
Set aspect ratio of your monitor
distance <number>
Set perspective distance from image. Use 0 for no perspective. Distance
is measured in multiples of the image depth (direction in the dimension
perpendicular to your monitor).
spin, tip <number>
Set view angle in degrees. Spin rotates the image around the z axis, tip
rotates the image around the y axis.
LIMITATIONS
Selecting more than about 55 x 55 resolution results in weird behavior for
reasons I don't understand.
EXAMPLES
Here is an example session that produces a few interesting graphs:
plot sin(x)
plot sin(y)
plot sin(x)*sin(y)
type cylinderr
plot 2.5
plot sin(2z)
plot sin(3z)+.5
type sphererho
plot 2.5sin(3theta)
thetasteps 50
phisteps 50
replot
SHAR_EOF
fi
exit 0
# End of shell archive
exit 0 # Just in case...
--
Kent Landfield INTERNET: kent@sparky.IMD.Sterling.COM
Sterling Software, IMD UUCP: uunet!sparky!kent
Phone: (402) 291-8300 FAX: (402) 291-4362
Please send comp.sources.misc-related mail to kent@uunet.uu.net.