src/util/contour.c

#ifndef lint
static const char       RCSid[] = "$Id$";
#endif
/*
 *  contour.c - program to make contour plots, mappings from 3 to
 *              2 dimenstions.
 *
 *     8/22/86
 */

#include  <stdio.h>

#include  <stdlib.h>

#include  <ctype.h>

#define  MAXPTS         2048            /* maximum number of input points */

typedef float  DATAPT[3];

DATAPT  xyz[MAXPTS];                    /* the input data */
int  xyzsiz;

double  zmin = 1e20;                    /* z minimum */
double  zmax = -1e20;                   /* z maximum */
int  minset = 0;                        /* user set minimum? */
int  maxset = 0;                        /* user set maximum? */

int  ncurves = 6;                       /* number of contours */


main(argc, argv)
int  argc;
char  *argv[];
{
        extern int  xycmp();
        FILE  *fp;
        int  i;

        for (i = 1; i < argc && argv[i][0] == '-'; i++)
                switch (argv[i][1]) {
                case 'n':
                        ncurves = atoi(argv[++i]);
                        break;
                case 'l':
                        zmin = atof(argv[++i]);
                        minset = 1;
                        break;
                case 'u':
                        zmax = atof(argv[++i]);
                        maxset = 1;
                        break;
                default:
                        fprintf(stderr, "%s: unknown option: %s\n",
                                        argv[0], argv[i]);
                        exit(1);
                }

        if (i == argc)
                getdata(stdin);
        else if (i < argc-1) {
                fprintf(stderr, "%s: too many file names\n", argv[0]);
                exit(1);
        } else if ((fp = fopen(argv[i], "r")) == NULL) {
                fprintf(stderr, "%s: file not found\n", argv[i]);
                exit(1);
        } else
                getdata(fp);

        qsort(xyz, xyzsiz, sizeof(DATAPT), xycmp);      /* sort data */

        for (i = 0; i < ncurves; i++)                   /* do contours */
                contour(i);

        exit(0);
}


getdata(fp)                     /* read input data */
FILE  *fp;
{
        register int  i;

        while (xyzsiz < MAXPTS) {
                for (i = 0; i < 3; i++)
                        if (fscanf(fp, "%f", &xyz[xyzsiz][i]) != 1)
                                break;
                if (i != 3)
                        break;
                if (xyz[xyzsiz][2] < zmin)
                        if (minset)
                                continue;
                        else
                                zmin = xyz[xyzsiz][2];
                if (xyz[xyzsiz][2] > zmax)
                        if (maxset)
                                continue;
                        else
                                zmax = xyz[xyzsiz][2];
                xyzsiz++;
        }
        if (!feof(fp)) {
                fprintf(stderr, "Error reading input data\n");
                exit(1);
        }
}


int
xycmp(p1, p2)                   /* -1 if p1 < p2, 0 if equal, 1 if p1 > p2 */
DATAPT  p1, p2;
{
        if (p1[0] > p2[0])
                return(1);
        if (p1[0] < p2[0])
                return(-1);
        if (p1[1] > p2[1])
                return(1);
        if (p1[1] < p2[1])
                return(-1);
        return(0);
}


contour(n)                      /* make contour n */
int  n;
{
        double  z;

        z = (n+.5)*(zmax-zmin)/ncurves + zmin;
        printf("%clabel=\"%f\";\n", n+'A', z);
        printf("%cdata=\n", n+'A');
        crossings(z);
        printf(";\n");
}


crossings(z)                    /* find crossings for z */
double  z;
{
        register DATAPT  *p0, *p1;

        p0 = p1 = xyz;
        while (p0 < xyz+xyzsiz-2) {
                while (p1 < xyz+xyzsiz-2) {
                        if (p0[0][0] < p1[0][0] && p0[0][1] <= p1[0][1])
                                break;
                        p1++;
                }
                if (p0[0][0] == p0[1][0])
                        cross(p0[0], p0[1], z);
                if (p0[0][1] == p1[0][1])
                        cross(p0[0], p1[0], z);
                p0++;
        }
}


cross(p0, p1, z)                /* mark crossing between p0 and p1 */
register DATAPT  p0, p1;
double  z;
{
        if (p1[2] - p0[2] == 0.0)
                return;
        z = (z - p0[2])/(p1[2] - p0[2]);
        if (z < 0.0 || z >= 1.0)
                return;
        printf("%e\t", p0[0]*(1.0 - z) + p1[0]*z);
        printf("%e\n", p0[1]*(1.0 - z) + p1[1]*z);
}
Revision:	1.3
Committed:	Sat Feb 22 02:07:30 2003 UTC (21 years, 2 months ago) by greg
Content type:	text/plain
Branch:	MAIN
CVS Tags:	rad3R5
Changes since 1.2:	+5 -7 lines
Log Message:	Changes and check-in for 3.5 release Includes new source files and modifications not recorded for many years See ray/doc/notes/ReleaseNotes for notes between 3.1 and 3.5 release
#	User	Rev	Content
1	greg	1.1	#ifndef lint
2	greg	1.3	static const char RCSid[] = "$Id$";
3	greg	1.1	#endif
4			/*
5			* contour.c - program to make contour plots, mappings from 3 to
6			* 2 dimenstions.
7			*
8			* 8/22/86
9			*/
10
11	greg	1.3	#include <stdio.h>
12
13			#include <stdlib.h>
14	greg	1.1
15	greg	1.3	#include <ctype.h>
16	greg	1.1
17			#define MAXPTS 2048 /* maximum number of input points */
18
19			typedef float DATAPT[3];
20
21			DATAPT xyz[MAXPTS]; /* the input data */
22			int xyzsiz;
23
24			double zmin = 1e20; /* z minimum */
25			double zmax = -1e20; /* z maximum */
26			int minset = 0; /* user set minimum? */
27			int maxset = 0; /* user set maximum? */
28
29			int ncurves = 6; /* number of contours */
30
31
32			main(argc, argv)
33			int argc;
34			char *argv[];
35			{
36			extern int xycmp();
37			FILE *fp;
38			int i;
39
40			for (i = 1; i < argc && argv[i][0] == '-'; i++)
41			switch (argv[i][1]) {
42			case 'n':
43			ncurves = atoi(argv[++i]);
44			break;
45			case 'l':
46			zmin = atof(argv[++i]);
47			minset = 1;
48			break;
49			case 'u':
50			zmax = atof(argv[++i]);
51			maxset = 1;
52			break;
53			default:
54			fprintf(stderr, "%s: unknown option: %s\n",
55			argv[0], argv[i]);
56			exit(1);
57			}
58
59			if (i == argc)
60			getdata(stdin);
61			else if (i < argc-1) {
62			fprintf(stderr, "%s: too many file names\n", argv[0]);
63			exit(1);
64			} else if ((fp = fopen(argv[i], "r")) == NULL) {
65			fprintf(stderr, "%s: file not found\n", argv[i]);
66			exit(1);
67			} else
68			getdata(fp);
69
70			qsort(xyz, xyzsiz, sizeof(DATAPT), xycmp); /* sort data */
71
72			for (i = 0; i < ncurves; i++) /* do contours */
73			contour(i);
74
75			exit(0);
76			}
77
78
79			getdata(fp) /* read input data */
80			FILE *fp;
81			{
82			register int i;
83
84			while (xyzsiz < MAXPTS) {
85			for (i = 0; i < 3; i++)
86			if (fscanf(fp, "%f", &xyz[xyzsiz][i]) != 1)
87			break;
88			if (i != 3)
89			break;
90			if (xyz[xyzsiz][2] < zmin)
91			if (minset)
92			continue;
93			else
94			zmin = xyz[xyzsiz][2];
95			if (xyz[xyzsiz][2] > zmax)
96			if (maxset)
97			continue;
98			else
99			zmax = xyz[xyzsiz][2];
100			xyzsiz++;
101			}
102			if (!feof(fp)) {
103			fprintf(stderr, "Error reading input data\n");
104			exit(1);
105			}
106			}
107
108
109			int
110			xycmp(p1, p2) /* -1 if p1 < p2, 0 if equal, 1 if p1 > p2 */
111			DATAPT p1, p2;
112			{
113			if (p1[0] > p2[0])
114			return(1);
115			if (p1[0] < p2[0])
116			return(-1);
117			if (p1[1] > p2[1])
118			return(1);
119			if (p1[1] < p2[1])
120			return(-1);
121			return(0);
122			}
123
124
125			contour(n) /* make contour n */
126			int n;
127			{
128			double z;
129
130			z = (n+.5)*(zmax-zmin)/ncurves + zmin;
131			printf("%clabel=\"%f\";\n", n+'A', z);
132			printf("%cdata=\n", n+'A');
133			crossings(z);
134			printf(";\n");
135			}
136
137
138			crossings(z) /* find crossings for z */
139			double z;
140			{
141			register DATAPT p0, p1;
142
143			p0 = p1 = xyz;
144			while (p0 < xyz+xyzsiz-2) {
145			while (p1 < xyz+xyzsiz-2) {
146			if (p0[0][0] < p1[0][0] && p0[0][1] <= p1[0][1])
147			break;
148			p1++;
149			}
150			if (p0[0][0] == p0[1][0])
151			cross(p0[0], p0[1], z);
152			if (p0[0][1] == p1[0][1])
153			cross(p0[0], p1[0], z);
154			p0++;
155			}
156			}
157
158
159			cross(p0, p1, z) /* mark crossing between p0 and p1 */
160			register DATAPT p0, p1;
161			double z;
162			{
163			if (p1[2] - p0[2] == 0.0)
164			return;
165			z = (z - p0[2])/(p1[2] - p0[2]);
166			if (z < 0.0 \|\| z >= 1.0)
167			return;
168			printf("%e\t", p0[0](1.0 - z) + p1[0]z);
169			printf("%e\n", p0[1](1.0 - z) + p1[1]z);
170			}