| 1 | greg | 1.1 | #ifndef lint | 
| 2 | schorsch | 1.4 | static const char       RCSid[] = "$Id: contour.c,v 1.3 2003/02/22 02:07:30 greg Exp $"; | 
| 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 | schorsch | 1.4 | static void getdata(FILE  *fp); | 
| 32 |  |  | static int xycmp(const void*  p1, const void* p2); | 
| 33 |  |  | static void contour(int  n); | 
| 34 |  |  | static void crossings(double  z); | 
| 35 |  |  | static void cross(DATAPT  p0, DATAPT  p1, double  z); | 
| 36 | greg | 1.1 |  | 
| 37 | schorsch | 1.4 |  | 
| 38 |  |  | int | 
| 39 |  |  | main( | 
| 40 |  |  | int  argc, | 
| 41 |  |  | char  *argv[] | 
| 42 |  |  | ) | 
| 43 | greg | 1.1 | { | 
| 44 |  |  | FILE  *fp; | 
| 45 |  |  | int  i; | 
| 46 |  |  |  | 
| 47 |  |  | for (i = 1; i < argc && argv[i][0] == '-'; i++) | 
| 48 |  |  | switch (argv[i][1]) { | 
| 49 |  |  | case 'n': | 
| 50 |  |  | ncurves = atoi(argv[++i]); | 
| 51 |  |  | break; | 
| 52 |  |  | case 'l': | 
| 53 |  |  | zmin = atof(argv[++i]); | 
| 54 |  |  | minset = 1; | 
| 55 |  |  | break; | 
| 56 |  |  | case 'u': | 
| 57 |  |  | zmax = atof(argv[++i]); | 
| 58 |  |  | maxset = 1; | 
| 59 |  |  | break; | 
| 60 |  |  | default: | 
| 61 |  |  | fprintf(stderr, "%s: unknown option: %s\n", | 
| 62 |  |  | argv[0], argv[i]); | 
| 63 |  |  | exit(1); | 
| 64 |  |  | } | 
| 65 |  |  |  | 
| 66 |  |  | if (i == argc) | 
| 67 |  |  | getdata(stdin); | 
| 68 |  |  | else if (i < argc-1) { | 
| 69 |  |  | fprintf(stderr, "%s: too many file names\n", argv[0]); | 
| 70 |  |  | exit(1); | 
| 71 |  |  | } else if ((fp = fopen(argv[i], "r")) == NULL) { | 
| 72 |  |  | fprintf(stderr, "%s: file not found\n", argv[i]); | 
| 73 |  |  | exit(1); | 
| 74 |  |  | } else | 
| 75 |  |  | getdata(fp); | 
| 76 |  |  |  | 
| 77 |  |  | qsort(xyz, xyzsiz, sizeof(DATAPT), xycmp);      /* sort data */ | 
| 78 |  |  |  | 
| 79 |  |  | for (i = 0; i < ncurves; i++)                   /* do contours */ | 
| 80 |  |  | contour(i); | 
| 81 |  |  |  | 
| 82 |  |  | exit(0); | 
| 83 |  |  | } | 
| 84 |  |  |  | 
| 85 |  |  |  | 
| 86 | schorsch | 1.4 | static void | 
| 87 |  |  | getdata(                        /* read input data */ | 
| 88 |  |  | FILE  *fp | 
| 89 |  |  | ) | 
| 90 | greg | 1.1 | { | 
| 91 |  |  | register int  i; | 
| 92 |  |  |  | 
| 93 |  |  | while (xyzsiz < MAXPTS) { | 
| 94 |  |  | for (i = 0; i < 3; i++) | 
| 95 |  |  | if (fscanf(fp, "%f", &xyz[xyzsiz][i]) != 1) | 
| 96 |  |  | break; | 
| 97 |  |  | if (i != 3) | 
| 98 |  |  | break; | 
| 99 | schorsch | 1.4 | if (xyz[xyzsiz][2] < zmin) { | 
| 100 |  |  | if (minset) { | 
| 101 | greg | 1.1 | continue; | 
| 102 | schorsch | 1.4 | } else { | 
| 103 | greg | 1.1 | zmin = xyz[xyzsiz][2]; | 
| 104 | schorsch | 1.4 | } | 
| 105 |  |  | } | 
| 106 |  |  | if (xyz[xyzsiz][2] > zmax) { | 
| 107 |  |  | if (maxset) { | 
| 108 | greg | 1.1 | continue; | 
| 109 | schorsch | 1.4 | } else { | 
| 110 | greg | 1.1 | zmax = xyz[xyzsiz][2]; | 
| 111 | schorsch | 1.4 | } | 
| 112 |  |  | } | 
| 113 | greg | 1.1 | xyzsiz++; | 
| 114 |  |  | } | 
| 115 |  |  | if (!feof(fp)) { | 
| 116 |  |  | fprintf(stderr, "Error reading input data\n"); | 
| 117 |  |  | exit(1); | 
| 118 |  |  | } | 
| 119 |  |  | } | 
| 120 |  |  |  | 
| 121 |  |  |  | 
| 122 | schorsch | 1.4 | static int | 
| 123 |  |  | xycmp(                  /* -1 if p1 < p2, 0 if equal, 1 if p1 > p2 */ | 
| 124 |  |  | const void*  p1, | 
| 125 |  |  | const void*     p2 | 
| 126 |  |  | ) | 
| 127 |  |  | { | 
| 128 |  |  | const float* pp1 = p1; | 
| 129 |  |  | const float* pp2 = p2; | 
| 130 |  |  | if (pp1[0] > pp2[0]) | 
| 131 | greg | 1.1 | return(1); | 
| 132 | schorsch | 1.4 | if (pp1[0] < pp2[0]) | 
| 133 | greg | 1.1 | return(-1); | 
| 134 | schorsch | 1.4 | if (pp1[1] > pp2[1]) | 
| 135 | greg | 1.1 | return(1); | 
| 136 | schorsch | 1.4 | if (pp1[1] < pp2[1]) | 
| 137 | greg | 1.1 | return(-1); | 
| 138 |  |  | return(0); | 
| 139 |  |  | } | 
| 140 |  |  |  | 
| 141 |  |  |  | 
| 142 | schorsch | 1.4 | static void | 
| 143 |  |  | contour(                        /* make contour n */ | 
| 144 |  |  | int  n | 
| 145 |  |  | ) | 
| 146 | greg | 1.1 | { | 
| 147 |  |  | double  z; | 
| 148 |  |  |  | 
| 149 |  |  | z = (n+.5)*(zmax-zmin)/ncurves + zmin; | 
| 150 |  |  | printf("%clabel=\"%f\";\n", n+'A', z); | 
| 151 |  |  | printf("%cdata=\n", n+'A'); | 
| 152 |  |  | crossings(z); | 
| 153 |  |  | printf(";\n"); | 
| 154 |  |  | } | 
| 155 |  |  |  | 
| 156 |  |  |  | 
| 157 | schorsch | 1.4 | static void | 
| 158 |  |  | crossings(                      /* find crossings for z */ | 
| 159 |  |  | double  z | 
| 160 |  |  | ) | 
| 161 | greg | 1.1 | { | 
| 162 |  |  | register DATAPT  *p0, *p1; | 
| 163 |  |  |  | 
| 164 |  |  | p0 = p1 = xyz; | 
| 165 |  |  | while (p0 < xyz+xyzsiz-2) { | 
| 166 |  |  | while (p1 < xyz+xyzsiz-2) { | 
| 167 |  |  | if (p0[0][0] < p1[0][0] && p0[0][1] <= p1[0][1]) | 
| 168 |  |  | break; | 
| 169 |  |  | p1++; | 
| 170 |  |  | } | 
| 171 |  |  | if (p0[0][0] == p0[1][0]) | 
| 172 |  |  | cross(p0[0], p0[1], z); | 
| 173 |  |  | if (p0[0][1] == p1[0][1]) | 
| 174 |  |  | cross(p0[0], p1[0], z); | 
| 175 |  |  | p0++; | 
| 176 |  |  | } | 
| 177 |  |  | } | 
| 178 |  |  |  | 
| 179 |  |  |  | 
| 180 | schorsch | 1.4 | static void | 
| 181 |  |  | cross(          /* mark crossing between p0 and p1 */ | 
| 182 |  |  | register DATAPT  p0, | 
| 183 |  |  | register DATAPT  p1, | 
| 184 |  |  | double  z | 
| 185 |  |  | ) | 
| 186 | greg | 1.1 | { | 
| 187 |  |  | if (p1[2] - p0[2] == 0.0) | 
| 188 |  |  | return; | 
| 189 |  |  | z = (z - p0[2])/(p1[2] - p0[2]); | 
| 190 |  |  | if (z < 0.0 || z >= 1.0) | 
| 191 |  |  | return; | 
| 192 |  |  | printf("%e\t", p0[0]*(1.0 - z) + p1[0]*z); | 
| 193 |  |  | printf("%e\n", p0[1]*(1.0 - z) + p1[1]*z); | 
| 194 |  |  | } |