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
#	Content
1	#ifndef lint
2	static const char RCSid[] = "$Id$";
3	#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	#include <stdio.h>
12
13	#include <stdlib.h>
14
15	#include <ctype.h>
16
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	}