1 /* pdnmesh - a 2D finite element solver
2 Copyright (C) 2001-2005 Sarod Yatawatta <sarod@users.sf.net>
3 This program is free software; you can redistribute it and/or modify
4 it under the terms of the GNU General Public License as published by
5 the Free Software Foundation; either version 2 of the License, or
6 (at your option) any later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16 $Id: contour.c,v 1.36 2005/02/18 18:27:10 sarod Exp $
17 */
18
19 /* Werner Hoch provided code to produce color in EPS files
20 */
21
22
23 #ifdef HAVE_CONFIG_H
24 #include <config.h>
25 #endif
26
27 #include <stdio.h>
28 #include <stdlib.h>
29 #include <math.h>
30 #include "types.h"
31
32
33 #define LEGEND_WIDTH 40
34 #define LEGEND_HEIGHT DEFAULT_HEIGHT
35 #define LEGEND_BORDER 5
36
37 int contour_levels; /* no of contour levels */
38 int current_plotting_contour; /* number in z array */
39 /* current max. value of gradient */
40 MY_DOUBLE max_gradient=0.01; /* initially a low one */
41
42
43 /* contour plotting */
44 static int
plot_contour(triangle * t,double potential,Mesh MM)45 plot_contour(triangle *t, double potential, Mesh MM)
46 {
47 /* variables for contour plotting */
48 int H,L,E;
49
50 MY_INT p1,p2,p3;
51 MY_DOUBLE x1,x2,y1,y2;
52
53 #ifdef DEBUG
54 printf("contour for triangle %d\n",t->n);
55 #endif
56 /* initialize the following to ignore compiler warnings */
57 p1=p2=p3=0;
58 x1=x2=y1=y2=0;
59
60 /* now determine the contours if any */
61 /* first determine point potentials */
62 E=H=L=0;
63 if ( Mzz(t->p0,current_plotting_contour,MM) == potential )
64 E++;
65 else if ( Mzz(t->p0,current_plotting_contour,MM) < potential )
66 L++;
67 else
68 H++;
69 if ( Mzz(t->p1,current_plotting_contour,MM) == potential )
70 E++;
71 else if ( Mzz(t->p1,current_plotting_contour,MM) < potential )
72 L++;
73 else
74 H++;
75 if ( Mzz(t->p2,current_plotting_contour,MM) == potential )
76 E++;
77 else if ( Mzz(t->p2,current_plotting_contour,MM) < potential )
78 L++;
79 else
80 H++;
81
82 /* now according to the values of E,H and L */
83 /* fit the general case */
84 if ( E == 3 ) { /* 3 contours */
85 /* we do not plot this */
86 #ifdef DEBUG
87 glBegin(GL_LINE_LOOP);
88 /* get first point */
89 glVertex2f( (GLfloat)Mx(t->p0,MM), (GLfloat)My(t->p0,MM));
90 /* get second point */
91 glVertex2f( (GLfloat)Mx(t->p1,MM), (GLfloat)My(t->p1,MM));
92 /* get third point */
93 glVertex2f( (GLfloat)Mx(t->p2,MM), (GLfloat)My(t->p2,MM));
94 glEnd();
95 #endif
96 #ifdef DEBUG
97 printf("3 lines\n");
98 #endif
99 return(3);
100 }
101 if ( E == 2 ) { /* 1 line */
102 /* generalize to p1,p2 */
103 if ( Mzz(t->p0,current_plotting_contour,MM) != potential ) {
104 p1 = t->p1; p2 = t->p2;
105 }
106 else if ( Mzz(t->p1,current_plotting_contour,MM) != potential ) {
107 p1 = t->p2; p2 = t->p0;
108 }
109 else if ( Mzz(t->p2,current_plotting_contour,MM) != potential ) {
110 p1 = t->p0; p2 = t->p1;
111 }
112 else
113 return(-1); /*error*/
114 /* now plot the general line */
115 glBegin(GL_LINE_LOOP);
116 /* get first point */
117 glVertex2f( (GLfloat)Mx(p1,MM), (GLfloat)My(p1,MM));
118 /* get second point */
119 glVertex2f( (GLfloat)Mx(p2,MM), (GLfloat)My(p2,MM));
120 glEnd();
121 #ifdef DEBUG
122 printf("1 line boundary\n");
123 #endif
124 return(1);
125 }
126 if ( (E == 1) && ( H == 1 ) && ( L == 1 ) ) {
127 /* againe one line generalize */
128 if ( Mzz(t->p0,current_plotting_contour,MM) == potential ) {
129 p1 = t->p1; p2 = t->p2; p3 = t->p0;
130 }
131 else if ( Mzz(t->p1,current_plotting_contour,MM) == potential ) {
132 p1 = t->p2; p2 = t->p0; p3 = t->p1;
133 }
134 else if ( Mzz(t->p2,current_plotting_contour,MM) == potential ) {
135 p1 = t->p0; p2 = t->p1; p3 = t->p2;
136 }
137 else
138 return(-1);
139 /* now interpolate p1,p2 to find the correct point */
140 if ( Mzz(p1,current_plotting_contour,MM) == Mzz(p2,current_plotting_contour,MM) )
141 return(-1); /* error */
142 x1 = Mx(p2,MM)
143 +( Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
144 y1 = My(p2,MM)
145 +( My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
146
147 /* now draw line between p3,(x1,y1) */
148 glBegin(GL_LINE_LOOP);
149 /* get first point */
150 glVertex2f( (GLfloat)Mx(p3,MM), (GLfloat)My(p3,MM));
151 glVertex2f( (GLfloat)x1,(GLfloat) y1 );
152 glEnd();
153 #ifdef DEBUG
154 printf("1 line equal\n");
155 #endif
156 return(1);
157 }
158
159 if (( ( H == 2 ) && ( L ==1 ) ) || ( ( H == 1 ) && ( L == 2 ))) {
160 /* two lines, need to generalize */
161
162 if ( ( H == 2 ) && ( L ==1 ) ) {
163 /* find the lowest point, assign it to p1 */
164 if ( Mzz(t->p0,current_plotting_contour,MM) < potential ) {
165 p1 = t->p0; p2 = t->p1; p3 = t->p2;
166 } else
167 if ( Mzz(t->p1,current_plotting_contour,MM) < potential ) {
168 p1 = t->p1; p2 = t->p2; p3 = t->p0;
169 } else
170 if ( Mzz(t->p2,current_plotting_contour,MM) < potential ) {
171 p1 = t->p2; p2 = t->p0; p3 = t->p1;
172 } else
173 return(-1); /* error */
174 #ifdef DEBUG
175 printf("1 line case1\n");
176 #endif
177 }
178 if ( ( H == 1 ) && ( L == 2 ) ) {
179 /* find the highest point, assign it to p1 */
180 if ( Mzz(t->p0,current_plotting_contour,MM) > potential ) {
181 p1 = t->p0; p2 = t->p1; p3 = t->p2;
182 } else
183 if ( Mzz(t->p1,current_plotting_contour,MM) > potential ) {
184 p1 = t->p1; p2 = t->p2; p3 = t->p0;
185 } else
186 if ( Mzz(t->p2,current_plotting_contour,MM) > potential ) {
187 p1 = t->p2; p2 = t->p0; p3 = t->p1;
188 } else
189 return(-1);
190 #ifdef DEBUG
191 printf("1 line case2\n");
192 #endif
193 }
194 /* now interpolate between (p1,p2) and (p1,p3) */
195
196 /* first p1,p2 interpolation */
197 if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p2,current_plotting_contour,MM) ) {
198 x1 = Mx(p2,MM)
199 + (Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
200 y1 = My(p2,MM)
201 + (My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
202 }
203 /* now p1,p3 interpolation */
204 if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p3,current_plotting_contour,MM) ) {
205 x2 = Mx(p3,MM)
206 + (Mx(p1,MM)-Mx(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
207 y2 = My(p3,MM)
208 + (My(p1,MM)-My(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
209
210 }
211 /* now draw line (x1,y1) (x2,y2) */
212 glBegin(GL_LINE_LOOP);
213 /* get first point */
214 glVertex2f( (GLfloat)x1, (GLfloat) y1 );
215 glVertex2f( (GLfloat)x2, (GLfloat) y2 );
216 glEnd();
217 return(1);
218
219 }
220
221 /* if we reach here, no contour */
222 return(0);
223 }
224
225 /* 3d plot with potential */
226 static int
plot_contour_in_3d(triangle * t,double potential,Mesh MM)227 plot_contour_in_3d(triangle *t, double potential,Mesh MM)
228 {
229 /* variables for contour plotting */
230 int H,L,E;
231
232 unsigned long int p1,p2,p3;
233 long double x1,x2,y1,y2;
234 double elevation;/*=potential/g_maxpot;*/
235 if ( ABS(g_maxpot)<ABS(g_minpot) ) {
236 elevation=potential/ABS(1.1*g_minpot);
237 } else {
238 elevation=potential/ABS(1.1*g_maxpot);
239 }
240 #ifdef DEBUG
241 printf("contour for triangle %d\n",t->n);
242 #endif
243 p1=p2=p3=0;
244 x1=x2=y1=y2=0;
245
246
247 /* now determine the contours if any */
248 /* first determine point potentials */
249 E=H=L=0;
250 if ( Mzz(t->p0,current_plotting_contour,MM) == potential )
251 E++;
252 else if ( Mzz(t->p0,current_plotting_contour,MM) < potential )
253 L++;
254 else
255 H++;
256 if ( Mzz(t->p1,current_plotting_contour,MM) == potential )
257 E++;
258 else if ( Mzz(t->p1,current_plotting_contour,MM) < potential )
259 L++;
260 else
261 H++;
262 if ( Mzz(t->p2,current_plotting_contour,MM) == potential )
263 E++;
264 else if ( Mzz(t->p2,current_plotting_contour,MM) < potential )
265 L++;
266 else
267 H++;
268
269 /* now according to the values of E,H and L */
270 /* fit the general case */
271 if ( E == 3 ) { /* 3 contours */
272 /* we do not plot this */
273 #ifdef DEBUG
274 glBegin(GL_LINE_LOOP);
275 /* get first point */
276 glVertex3f( (GLfloat)Mx(t->p0,MM), (GLfloat)My(t->p0,MM),(GLfloat)elevation);
277 /* get second point */
278 glVertex3f( (GLfloat)Mx(t->p1,MM), (GLfloat)My(t->p1,MM),(GLfloat)elevation);
279 /* get third point */
280 glVertex3f( (GLfloat)Mx(t->p2,MM), (GLfloat)My(t->p2,MM),(GLfloat)elevation);
281 glEnd();
282 #endif
283 #ifdef DEBUG
284 printf("3 lines\n");
285 #endif
286 return(3);
287 }
288 if ( E == 2 ) { /* 1 line */
289 /* generalize to p1,p2 */
290 if ( Mzz(t->p0,current_plotting_contour,MM) != potential ) {
291 p1 = t->p1; p2 = t->p2;
292 }
293 else if ( Mzz(t->p1,current_plotting_contour,MM) != potential ) {
294 p1 = t->p2; p2 = t->p0;
295 }
296 else if ( Mzz(t->p2,current_plotting_contour,MM) != potential ) {
297 p1 = t->p0; p2 = t->p1;
298 }
299 else
300 return(-1); /*error*/
301 /* now plot the general line */
302 glBegin(GL_LINE_LOOP);
303 /* get first point */
304 glVertex3f( (GLfloat)Mx(p1,MM), (GLfloat)My(p1,MM),(GLfloat)elevation);
305 /* get second point */
306 glVertex3f( (GLfloat)Mx(p2,MM), (GLfloat)My(p2,MM),(GLfloat)elevation);
307 glEnd();
308 #ifdef DEBUG
309 printf("1 line boundary\n");
310 #endif
311 return(1);
312 }
313 if ( (E == 1) && ( H == 1 ) && ( L == 1 ) ) {
314 /* againe one line generalize */
315 if ( Mzz(t->p0,current_plotting_contour,MM) == potential ) {
316 p1 = t->p1; p2 = t->p2; p3 = t->p0;
317 }
318 else if ( Mzz(t->p1,current_plotting_contour,MM) == potential ) {
319 p1 = t->p2; p2 = t->p0; p3 = t->p1;
320 }
321 else if ( Mzz(t->p2,current_plotting_contour,MM) == potential ) {
322 p1 = t->p0; p2 = t->p1; p3 = t->p2;
323 }
324 else
325 return(-1);
326 /* now interpolate p1,p2 to find the correct point */
327 if ( Mzz(p1,current_plotting_contour,MM) == Mzz(p2,current_plotting_contour,MM) )
328 return(-1); /* error */
329 x1 = Mx(p2,MM)
330 +( Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
331 y1 = My(p2,MM)
332 +( My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
333
334 /* now draw line between p3,(x1,y1) */
335 glBegin(GL_LINE_LOOP);
336 /* get first point */
337 glVertex3f( (GLfloat)Mx(p3,MM), (GLfloat)My(p3,MM),(GLfloat)elevation);
338 glVertex3f( (GLfloat)x1,(GLfloat) y1,(GLfloat)elevation );
339 glEnd();
340 #ifdef DEBUG
341 printf("1 line equal\n");
342 #endif
343 return(1);
344 }
345
346 if (( ( H == 2 ) && ( L ==1 ) ) || ( ( H == 1 ) && ( L == 2 ))) {
347 /* two lines, need to generalize */
348
349 if ( ( H == 2 ) && ( L ==1 ) ) {
350 /* find the lowest point, assign it to p1 */
351 if ( Mzz(t->p0,current_plotting_contour,MM) < potential ) {
352 p1 = t->p0; p2 = t->p1; p3 = t->p2;
353 } else
354 if ( Mzz(t->p1,current_plotting_contour,MM) < potential ) {
355 p1 = t->p1; p2 = t->p2; p3 = t->p0;
356 } else
357 if ( Mzz(t->p2,current_plotting_contour,MM) < potential ) {
358 p1 = t->p2; p2 = t->p0; p3 = t->p1;
359 } else
360 return(-1); /* error */
361 #ifdef DEBUG
362 printf("1 line case1\n");
363 #endif
364 }
365 if ( ( H == 1 ) && ( L == 2 ) ) {
366 /* find the highest point, assign it to p1 */
367 if ( Mzz(t->p0,current_plotting_contour,MM) > potential ) {
368 p1 = t->p0; p2 = t->p1; p3 = t->p2;
369 } else
370 if ( Mzz(t->p1,current_plotting_contour,MM) > potential ) {
371 p1 = t->p1; p2 = t->p2; p3 = t->p0;
372 } else
373 if ( Mzz(t->p2,current_plotting_contour,MM) > potential ) {
374 p1 = t->p2; p2 = t->p0; p3 = t->p1;
375 } else
376 return(-1);
377 #ifdef DEBUG
378 printf("1 line case2\n");
379 #endif
380 }
381 /* now interpolate between (p1,p2) and (p1,p3) */
382
383 /* first p1,p2 interpolation */
384 if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p2,current_plotting_contour,MM) ) {
385 x1 = Mx(p2,MM)
386 + (Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
387 y1 = My(p2,MM)
388 + (My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
389 }
390 /* now p1,p3 interpolation */
391 if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p3,current_plotting_contour,MM) ) {
392 x2 = Mx(p3,MM)
393 + (Mx(p1,MM)-Mx(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
394 y2 = My(p3,MM)
395 + (My(p1,MM)-My(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
396
397 }
398 /* now draw line (x1,y1) (x2,y2) */
399 glBegin(GL_LINE_LOOP);
400 /* get first point */
401 glVertex3f( (GLfloat)x1, (GLfloat) y1,(GLfloat)elevation );
402 glVertex3f( (GLfloat)x2, (GLfloat) y2,(GLfloat)elevation );
403 glEnd();
404 return(1);
405
406 }
407
408 /* if we reach here, no contour */
409 return(0);
410 }
411
412 /* to get a color to plot */
413 int
get_colour(float * redp,float * greenp,float * bluep,int level,int max_levels)414 get_colour(float* redp, float* greenp, float *bluep, int level, int max_levels)
415 {
416 float cl = (float)level/(float)max_levels;
417 /*
418 if ( cl < 0.166666667 ) {
419 *redp = 0.0f;
420 *greenp = 1.0f;
421 *bluep = cl*6.0;
422 } else if ( cl < 0.333333333 ) {
423 *redp = 0.0f;
424 *greenp = 2.0f - 6.0*cl;
425 *bluep = 1.0f;
426 } else if ( cl < 0.5 ) {
427 *redp = cl*6.0-3.0;
428 *greenp = 0.0f;
429 *bluep = 1.0f;
430 } else if ( cl < 0.6666666667 ) {
431 *redp = 1.0f;
432 *greenp = 0.0f;
433 *bluep = 4.0f-cl*6.0;
434 } else if ( cl < 0.833333333 ) {
435 *redp = 1.0f;
436 *greenp = cl*6.0 -4.0f;
437 *bluep = 0.0f;
438 } else {
439 *redp = 1.0f;
440 *greenp = 1.0f;
441 *bluep = cl*6.0-5.0f;
442 }
443 */
444 if ( cl < 0.25 ) {
445 *redp = 0.0f;
446 *greenp = cl*4.0f;
447 *bluep = 1.0;
448 } else if ( cl < 0.5 ) {
449 *redp = 0.0f;
450 *greenp = 1.0f;
451 *bluep = 2.0f -cl*4.0f;
452 } else if ( cl < 0.75 ) {
453 *redp = cl*4.0f-2.0f;
454 *greenp = 1.0f;
455 *bluep = 0.0f;
456 } else {
457 *redp = 1.0f;
458 *greenp = 4.0f - cl*4.0f;
459 *bluep = 0.0f;
460 }
461
462 /* blue - green */
463 /*
464 if ( cl < 1.5 ) {
465 *redp = 0.5;
466 *greenp = cl*0.5;
467 *bluep = 1.0 - *greenp;
468 */
469 /* green - red */
470 /*
471 } else {
472 *redp = 0.25 + (cl-1.5)*0.5;
473 *bluep = 0.0;
474 *greenp = 1.0 - *redp;
475 }
476 */
477
478 /*
479 *redp = (float)level/(float)max_levels;
480 *greenp = 0.7f;
481 *bluep = 1.0f - *redp;
482 */
483 return(0);
484
485 }
486
487 /* function to plot the contour */
488 void
plot_contour_all(Mesh MM)489 plot_contour_all(Mesh MM)
490 {
491 triangle *t;
492 int j; /* for contour calculation */
493 float r,g,b; /* colours */
494 double step = (g_maxpot-g_minpot)/(float)contour_levels; /* contour step */
495
496
497
498 #ifdef DEBUG
499 printf("plotting from %lf to %lf at %lf intervals\n",g_maxpot,g_minpot,step);
500 #endif
501
502 DAG_traverse_list_reset(&dt);
503 t=DAG_traverse_prune_list(&dt);
504 while(t) {
505 if ( t && t->status==LIVE ) {
506 for ( j = 0; j < contour_levels ; j++ ) {
507 /* get the colour */
508 get_colour(&r,&g,&b,j, contour_levels);
509 /* set the colour */
510 glColor3f((GLfloat)r, (GLfloat)g, (GLfloat)b);
511 plot_contour(t,g_minpot+j*step,MM);
512 }
513
514 }
515
516 t=DAG_traverse_prune_list(&dt);
517 }
518
519 }
520
521 void
plot_contour_all_in_3d(Mesh MM)522 plot_contour_all_in_3d(Mesh MM)
523 {
524 triangle *t;
525 int j; /* for contour calculation */
526 float r,g,b; /* colours */
527 double step = (g_maxpot-g_minpot)/(float)contour_levels; /* contour step */
528
529
530
531 #ifdef DEBUG
532 printf("plotting from %lf to %lf at %lf intervals\n",g_maxpot,g_minpot,step);
533 #endif
534
535 DAG_traverse_list_reset(&dt);
536 t=DAG_traverse_prune_list(&dt);
537 while(t) {
538 if ( t && t->status==LIVE ) {
539 for ( j = 0; j < contour_levels ; j++ ) {
540 /* get the colour */
541 /* get_colour(&r,&g,&b,j, contour_levels); */
542 r=g=b=0.0;
543 glColor3f((GLfloat)r, (GLfloat)g, (GLfloat)b);
544 plot_contour_in_3d(t,g_minpot+j*step,MM);
545 }
546
547 }
548
549 t=DAG_traverse_prune_list(&dt);
550 }
551
552 }
553
554
555 /* a function to print the contour of a triangle to a postscript file */
556 static int
print_contour(FILE * plotf,triangle * t,MY_DOUBLE potential,MY_DOUBLE myscale,float r,float g,float b,int color_flag,Mesh MM)557 print_contour(FILE *plotf, triangle *t, MY_DOUBLE potential, MY_DOUBLE myscale,
558 float r, float g, float b, int color_flag, Mesh MM)
559 {
560 /* variables for contour plotting */
561 int H,L,E;
562
563 MY_INT p1,p2,p3;
564 MY_DOUBLE x1,x2,y1,y2;
565
566 /* note: we use shorthand here as follows:
567 /l {lineto} bind def
568 /m {moveto} bind def
569 /n {newpath} bind def
570 /s {stroke} bind def
571 /slw {setlinewidth} bind def
572 */
573
574
575 #ifdef DEBUG
576 printf("contour for triangle %d\n",t->n);
577 #endif
578 /* initialize the following to ignore compiler warnings */
579 p1=p2=p3=0;
580 x1=x2=y1=y2=0;
581
582 /* now determine the contours if any */
583 /* first determine point potentials */
584 E=H=L=0;
585 if ( Mzz(t->p0,current_plotting_contour,MM) == potential )
586 E++;
587 else if ( Mzz(t->p0,current_plotting_contour,MM) < potential )
588 L++;
589 else
590 H++;
591 if ( Mzz(t->p1,current_plotting_contour,MM) == potential )
592 E++;
593 else if ( Mzz(t->p1,current_plotting_contour,MM) < potential )
594 L++;
595 else
596 H++;
597 if ( Mzz(t->p2,current_plotting_contour,MM) == potential )
598 E++;
599 else if ( Mzz(t->p2,current_plotting_contour,MM) < potential )
600 L++;
601 else
602 H++;
603
604 /* now according to the values of E,H and L */
605 /* fit the general case */
606 if ( E == 3 ) { /* 3 contours */
607 /* we do not plot this */
608 #ifdef DEBUG
609 printf("3 lines\n");
610 #endif
611 return(3);
612 }
613 if ( E == 2 ) { /* 1 line */
614 /* generalize to p1,p2 */
615 if ( Mzz(t->p0,current_plotting_contour,MM) != potential ) {
616 p1 = t->p1; p2 = t->p2;
617 }
618 else if ( Mzz(t->p1,current_plotting_contour,MM) != potential ) {
619 p1 = t->p2; p2 = t->p0;
620 }
621 else if ( Mzz(t->p2,current_plotting_contour,MM) != potential ) {
622 p1 = t->p0; p2 = t->p1;
623 }
624 else
625 return(-1); /*error*/
626 /* set color */
627 if ( color_flag )
628 fprintf(plotf,"%f %f %f srgb\n",r,g,b);
629 /* now plot the general line */
630 fprintf(plotf,"n ");
631 /* get first point */
632 fprintf(plotf,"%lf %lf m ",(double) (g_xrange+Mx(p1,MM)*g_xscale)*myscale, (double)(g_yrange + My(p1,MM)*g_yscale)*myscale);
633 /* get second point */
634 fprintf(plotf,"%lf %lf l ",(double) (g_xrange+Mx(p2,MM)*g_xscale)*myscale, (double)(g_yrange + My(p2,MM)*g_yscale)*myscale);
635 fprintf(plotf,"cp ");
636 fprintf(plotf,"s\n");
637 #ifdef DEBUG
638 printf("1 line boundary\n");
639 #endif
640 return(1);
641 }
642 if ( (E == 1) && ( H == 1 ) && ( L == 1 ) ) {
643 /* againe one line generalize */
644 if ( Mzz(t->p0,current_plotting_contour,MM) == potential ) {
645 p1 = t->p1; p2 = t->p2; p3 = t->p0;
646 }
647 else if ( Mzz(t->p1,current_plotting_contour,MM) == potential ) {
648 p1 = t->p2; p2 = t->p0; p3 = t->p1;
649 }
650 else if ( Mzz(t->p2,current_plotting_contour,MM) == potential ) {
651 p1 = t->p0; p2 = t->p1; p3 = t->p2;
652 }
653 else
654 return(-1);
655 /* now interpolate p1,p2 to find the correct point */
656 if ( Mzz(p1,current_plotting_contour,MM) == Mzz(p2,current_plotting_contour,MM) )
657 return(-1); /* error */
658 x1 = Mx(p2,MM)
659 +( Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
660 y1 = My(p2,MM)
661 +( My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
662
663 /* set color */
664 if ( color_flag )
665 fprintf(plotf,"%f %f %f srgb\n",r,g,b);
666
667 /* now draw line between p3,(x1,y1) */
668 fprintf(plotf,"n ");
669 /* get first point */
670 fprintf(plotf,"%lf %lf m ",(double) (g_xrange+Mx(p3,MM)*g_xscale)*myscale, (double)(g_yrange + My(p3,MM)*g_yscale)*myscale);
671 /* get second point */
672 fprintf(plotf,"%lf %lf l ",(double) (g_xrange+x1*g_xscale)*myscale, (double)(g_yrange + y1*g_yscale)*myscale);
673 fprintf(plotf,"cp ");
674 fprintf(plotf,"s\n");
675
676 #ifdef DEBUG
677 printf("1 line equal\n");
678 #endif
679 return(1);
680 }
681
682 if (( ( H == 2 ) && ( L ==1 ) ) || ( ( H == 1 ) && ( L == 2 ))) {
683 /* two lines, need to generalize */
684
685 if ( ( H == 2 ) && ( L ==1 ) ) {
686 /* find the lowest point, assign it to p1 */
687 if ( Mzz(t->p0,current_plotting_contour,MM) < potential ) {
688 p1 = t->p0; p2 = t->p1; p3 = t->p2;
689 } else
690 if ( Mzz(t->p1,current_plotting_contour,MM) < potential ) {
691 p1 = t->p1; p2 = t->p2; p3 = t->p0;
692 } else
693 if ( Mzz(t->p2,current_plotting_contour,MM) < potential ) {
694 p1 = t->p2; p2 = t->p0; p3 = t->p1;
695 } else
696 return(-1); /* error */
697 #ifdef DEBUG
698 printf("1 line case1\n");
699 #endif
700 }
701 if ( ( H == 1 ) && ( L == 2 ) ) {
702 /* find the highest point, assign it to p1 */
703 if ( Mzz(t->p0,current_plotting_contour,MM) > potential ) {
704 p1 = t->p0; p2 = t->p1; p3 = t->p2;
705 } else
706 if ( Mzz(t->p1,current_plotting_contour,MM) > potential ) {
707 p1 = t->p1; p2 = t->p2; p3 = t->p0;
708 } else
709 if ( Mzz(t->p2,current_plotting_contour,MM) > potential ) {
710 p1 = t->p2; p2 = t->p0; p3 = t->p1;
711 } else
712 return(-1);
713 #ifdef DEBUG
714 printf("1 line case2\n");
715 #endif
716 }
717 /* now interpolate between (p1,p2) and (p1,p3) */
718
719 /* first p1,p2 interpolation */
720 if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p2,current_plotting_contour,MM) ) {
721 x1 = Mx(p2,MM)
722 + (Mx(p1,MM)-Mx(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
723 y1 = My(p2,MM)
724 + (My(p1,MM)-My(p2,MM))*(potential-Mzz(p2,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p2,current_plotting_contour,MM));
725 }
726 /* now p1,p3 interpolation */
727 if ( Mzz(p1,current_plotting_contour,MM) != Mzz(p3,current_plotting_contour,MM) ) {
728 x2 = Mx(p3,MM)
729 + (Mx(p1,MM)-Mx(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
730 y2 = My(p3,MM)
731 + (My(p1,MM)-My(p3,MM))*(potential-Mzz(p3,current_plotting_contour,MM))/(Mzz(p1,current_plotting_contour,MM)-Mzz(p3,current_plotting_contour,MM));
732
733 }
734 /* set color */
735 if ( color_flag )
736 fprintf(plotf,"%f %f %f srgb\n",r,g,b);
737
738 /* now draw line (x1,y1) (x2,y2) */
739 fprintf(plotf,"n ");
740 /* get first point */
741 fprintf(plotf,"%lf %lf m ",(double) (g_xrange+x1*g_xscale)*myscale, (double)(g_yrange + y1*g_yscale)*myscale);
742 /* get second point */
743 fprintf(plotf,"%lf %lf l ",(double) (g_xrange+x2*g_xscale)*myscale, (double)(g_yrange + y2*g_yscale)*myscale);
744 fprintf(plotf,"cp ");
745 fprintf(plotf,"s\n");
746
747 }
748
749 /* if we reach here, no contour */
750 return(0);
751
752 }
753
754 static void
print_polygon(polygon * poly,double myscale,FILE * plotf,Mesh MM)755 print_polygon(polygon *poly, double myscale, FILE *plotf, Mesh MM)
756 {
757
758 if (poly&&poly->head) {
759 int i=poly->nedges;
760 poly_elist *tempe=poly->head;
761 for (;i;i--) {
762 fprintf(plotf,"n ");
763 fprintf(plotf,"%lf %lf m ",(double)(Mx(tempe->data.p1,MM)*g_xscale+g_xrange)*myscale, (double)(My(tempe->data.p1,MM)*g_yscale+g_yrange)*myscale);
764 fprintf(plotf,"%lf %lf l ",(double)(Mx(tempe->data.p2,MM)*g_xscale+g_xrange)*myscale, (double)(My(tempe->data.p2,MM)*g_yscale+g_yrange)*myscale);
765 fprintf(plotf,"s\n");
766 tempe=tempe->next;
767 }
768 }
769
770 }
771
772
773 /* printing contour plot */
774 void
print_contour_all(const char * filename,int color_flag,Mesh MM)775 print_contour_all(const char* filename,int color_flag, Mesh MM)
776 {
777 /* color_flag: 0, BW contours, 1: color contours, 2: color contours+legend */
778 FILE *plotf;
779 MY_DOUBLE myscale,delta;
780 int i;
781 float r,g,b;
782 triangle *t;
783 double step = (g_maxpot-g_minpot)/(float)contour_levels; /* contour step */
784
785 if ((plotf=fopen(filename,"w"))==NULL){
786 fprintf(stderr,"%s: %d: could not open file %s\n",__FILE__,__LINE__,filename);
787 return;
788 }
789
790
791 /* determine the right bound of bounding box */
792 if ( g_xrange > g_yrange )
793 myscale = 400.0/g_xrange;
794 else
795 myscale = 400.0/g_yrange;
796
797 /* print the header */
798 fprintf(plotf,"%%!PS-Adobe-2.0 EPSF-3.0\n");
799 fprintf(plotf,"%%%%Title: %s\n",filename);
800 fprintf(plotf,"%%%%Generated By: "PACKAGE" "VERSION" contour levels: %d, min= "MDF", max= "MDF"\n", contour_levels, g_minpot, g_maxpot);
801 /* now bounding box -add border of 20 */
802 /* add x margin of 60 for legend
803 the myscale is selected such that max dimension is 820+80(x) of 820 (y)
804 */
805 if ( color_flag==2 ) {
806 fprintf(plotf,"%%%%BoundingBox: 0 0 %4.2lf %4.2lf\n",(double)(2*g_xrange*myscale+20+100), (double)(2*g_yrange*myscale+20) );
807 } else {
808 fprintf(plotf,"%%%%BoundingBox: 0 0 %4.2lf %4.2lf\n",(double)(2*g_xrange*myscale+20), (double)(2*g_yrange*myscale+20) );
809 }
810 fprintf(plotf,"%%%%Magnification: 1.0000\n");
811
812
813
814 /* start plotting from 10 10 */
815 fprintf(plotf,"10 10 translate\n");
816
817
818 /* now print some shorter definitions */
819 fprintf(plotf,"/l {lineto} bind def\n"); /* lineto */
820 fprintf(plotf,"/f {fill} bind def\n"); /* lineto */
821 fprintf(plotf,"/m {moveto} bind def\n"); /* moveto */
822 fprintf(plotf,"/n {newpath} bind def\n"); /* newpath */
823 fprintf(plotf,"/s {stroke} bind def\n"); /* stroke */
824 fprintf(plotf,"/slw {setlinewidth} bind def\n"); /* setlinewidth */
825 fprintf(plotf,"/cp {closepath} bind def\n"); /* closepath */
826 if ( color_flag )
827 fprintf(plotf,"/srgb {setrgbcolor} bind def\n"); /* rgb color */
828
829 /* now print contours */
830 fprintf(plotf,"1 slw\n");
831 DAG_traverse_list_reset(&dt);
832 t=DAG_traverse_prune_list(&dt);
833 while(t) {
834 if ( t && t->status==LIVE ) {
835 for ( i = 0; i < contour_levels ; i++ ) {
836 /* get the colour */
837 get_colour(&r,&g,&b,i, contour_levels);
838 /* set the colour */
839 print_contour(plotf,t,g_minpot+i*step,myscale,r,g,b,color_flag,MM);
840 }
841
842 }
843
844 t=DAG_traverse_prune_list(&dt);
845 }
846
847 /* then plot the boundaries */
848 /* set the line thikness */
849 fprintf(plotf,"2 slw\n");
850
851 if ( color_flag )
852 fprintf(plotf,"0 0 0 srgb\n");
853 for ( i = 0; i < bound.npoly; i++ ) {
854 print_polygon(&bound.poly_array[i],myscale,plotf,MM);
855 }
856
857 if ( color_flag==2 ) {
858 /* draw legend */
859 fprintf(plotf,"1 slw\n");
860 /* box */
861 fprintf(plotf,"n %lf %lf m %lf %lf l %lf %lf l %lf %lf l %lf %lf l s\n",(double)(2*g_xrange*myscale+20),0.0, (double)(2*g_xrange*myscale+60),0.0, (double)(2*g_xrange*myscale+60),(double)(2*g_yrange*myscale), (double)(2*g_xrange*myscale+20), (double)(2*g_yrange*myscale), (double)(2*g_xrange*myscale+20),0.0);
862 /* tick marks */
863 if ( contour_levels > 1 ) {
864 delta=(double)(2*g_yrange*myscale)/(contour_levels-1);
865 for (i=1; i< contour_levels-1; i++) {
866 fprintf(plotf,"n %lf %lf m %lf %lf l %lf %lf l %lf %lf l cp\n",(double)(2*g_xrange*myscale+20),(i-1)*delta, (double)(2*g_xrange*myscale+60),(i-1)*delta,(double)(2*g_xrange*myscale+60),(i)*delta,(double)(2*g_xrange*myscale+20),(i)*delta);
867 get_colour(&r,&g,&b,i, contour_levels);
868 fprintf(plotf,"%f %f %f srgb f\n",r,g,b);
869 fprintf(plotf,"n %lf %lf m %lf %lf l s\n",(double)(2*g_xrange*myscale+20),i*delta,(double)(2*g_xrange*myscale+20+5),i*delta);
870 }
871 fprintf(plotf,"n %lf %lf m %lf %lf l %lf %lf l %lf %lf l cp\n",(double)(2*g_xrange*myscale+20),(i-1)*delta, (double)(2*g_xrange*myscale+60),(i-1)*delta,(double)(2*g_xrange*myscale+60),(i)*delta,(double)(2*g_xrange*myscale+20),(i)*delta);
872 get_colour(&r,&g,&b,i, contour_levels);
873 fprintf(plotf,"%f %f %f srgb f\n",r,g,b);
874
875 /* now text for potentials */
876 fprintf(plotf,"/Times-Roman findfont\n");
877 fprintf(plotf,"0.5 slw\n");
878 fprintf(plotf,"9 scalefont setfont\n");
879 fprintf(plotf,"0 0 0 srgb\n");
880 for (i=0; i< contour_levels; i++) {
881 fprintf(plotf,"n %lf %lf m (%5.2lf) show\n",(double)(2*g_xrange*myscale+62),i*delta,g_minpot+i*step);
882 }
883 }
884 }
885
886 fprintf(plotf,"%%%%EOF\n");
887 fclose(plotf);
888 }
889
890 /* printing mesh as eps*/
891 void
print_mesh_eps(const char * filename,Mesh MM)892 print_mesh_eps(const char* filename,Mesh MM)
893 {
894
895 FILE *plotf;
896 MY_DOUBLE myscale;
897 int i;
898 triangle *t;
899
900 if ((plotf=fopen(filename,"w"))==NULL){
901 fprintf(stderr,"%s: %d: could not open file %s\n",__FILE__,__LINE__,filename);
902 return;
903 }
904
905
906 /* determine the right bound of bounding box */
907 if ( g_xrange > g_yrange )
908 myscale = 400.0/g_xrange;
909 else
910 myscale = 400.0/g_yrange;
911
912 /* print the header */
913 fprintf(plotf,"%%!PS-Adobe-2.0 EPSF-3.0\n");
914 fprintf(plotf,"%%%%Title: %s\n",filename);
915 fprintf(plotf,"%%%%Generated By: "PACKAGE" "VERSION" contour levels: %d, min= "MDF", max= "MDF"\n", contour_levels, g_minpot, g_maxpot);
916 /* now bounding box -add border of 20 */
917 fprintf(plotf,"%%%%BoundingBox: 0 0 %4.2lf %4.2lf\n",(double)(2*g_xrange*myscale+20), (double)(2*g_yrange*myscale+20) );
918 fprintf(plotf,"%%%%Magnification: 1.0000\n");
919
920
921
922 /* start plotting from 10 10 */
923 fprintf(plotf,"10 10 translate\n");
924
925
926 /* now print some shorter definitions */
927 fprintf(plotf,"/l {lineto} bind def\n"); /* lineto */
928 fprintf(plotf,"/m {moveto} bind def\n"); /* moveto */
929 fprintf(plotf,"/n {newpath} bind def\n"); /* newpath */
930 fprintf(plotf,"/s {stroke} bind def\n"); /* stroke */
931 fprintf(plotf,"/slw {setlinewidth} bind def\n"); /* setlinewidth */
932 fprintf(plotf,"/cp {closepath} bind def\n"); /* closepath */
933
934 /* now print contours */
935 fprintf(plotf,"1 slw\n");
936 DAG_traverse_list_reset(&dt);
937 t=DAG_traverse_prune_list(&dt);
938 while(t) {
939 if ( t && t->status==LIVE ) {
940 /* get the colour */
941 /*get_colour(&r,&g,&b,j, contour_levels); */
942 /* set the colour */
943 fprintf(plotf,"n ");
944 fprintf(plotf,"%lf %lf m ",(double)(Mx(t->p0,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p0,MM)*g_yscale+g_yrange)*myscale);
945 fprintf(plotf,"%lf %lf l ",(double)(Mx(t->p1,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p1,MM)*g_yscale+g_yrange)*myscale);
946 fprintf(plotf,"s\n");
947 fprintf(plotf,"n ");
948 fprintf(plotf,"%lf %lf m ",(double)(Mx(t->p1,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p1,MM)*g_yscale+g_yrange)*myscale);
949 fprintf(plotf,"%lf %lf l ",(double)(Mx(t->p2,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p2,MM)*g_yscale+g_yrange)*myscale);
950 fprintf(plotf,"s\n");
951 fprintf(plotf,"n ");
952 fprintf(plotf,"%lf %lf m ",(double)(Mx(t->p2,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p2,MM)*g_yscale+g_yrange)*myscale);
953 fprintf(plotf,"%lf %lf l ",(double)(Mx(t->p0,MM)*g_xscale+g_xrange)*myscale, (double)(My(t->p0,MM)*g_yscale+g_yrange)*myscale);
954 fprintf(plotf,"s\n");
955
956 }
957
958 t=DAG_traverse_prune_list(&dt);
959 }
960
961 /* then plot the boundaries */
962 /* set the line thikness */
963 fprintf(plotf,"2 slw\n");
964 for ( i = 0; i < bound.npoly; i++ ) {
965 print_polygon(&bound.poly_array[i],myscale,plotf,MM);
966 }
967 fprintf(plotf,"%%%%EOF\n");
968
969 fclose(plotf);
970 }
971
972
973 /* printing mesh as ACII file*/
974 void
print_mesh_ascii(const char * filename,Mesh MM)975 print_mesh_ascii(const char* filename,Mesh MM)
976 {
977 FILE *outf;
978 unsigned int i,j;
979 triangle *t;
980 unsigned int *renumb,*re_renumb;
981 unsigned MY_INT total,unknowns,triangles;
982
983 if ((outf=fopen(filename,"w"))==NULL){
984 fprintf(stderr,"%s: %d: could not open file %s\n",__FILE__,__LINE__,filename);
985 return;
986 }
987
988 /* allocate memory for arrays */
989 if ((renumb = (MY_INT *)calloc(M.count,sizeof(MY_INT)))==0){
990 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
991 exit(1);
992 }
993
994 if((re_renumb = (MY_INT *)calloc(M.count,sizeof(MY_INT)))==0){
995 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
996 exit(1);
997 }
998
999 total=unknowns=triangles=0;
1000 /* need to solve for only points referred by triangles */
1001 DAG_traverse_list_reset(&dt);
1002 t=DAG_traverse_prune_list(&dt);
1003 while ( t ) {
1004 if ( t && t->status==LIVE ) {
1005 /* mark the points in this triangle */
1006 /* use re_renumber[] array for this */
1007 re_renumb[t->p0] = 1;
1008 re_renumb[t->p1] = 1;
1009 re_renumb[t->p2] = 1;
1010
1011 triangles++;
1012 }
1013 t=DAG_traverse_prune_list(&dt);
1014 }
1015
1016 /* now find unknowns, or INSIDE points */
1017 j = 0;
1018 for ( i = 0; i < M.count; i++ )
1019 if ( Mval(i,MM) == VAR && re_renumb[i] == 1) {
1020 renumb[i] = j++;
1021 unknowns++;
1022 }
1023 /* now renumber the known points */
1024 total = unknowns;
1025 for ( i = 0; i < M.count; i++ )
1026 if ( Mval(i,MM) == FX && re_renumb[i]==1 ) {
1027 renumb[i] = j++;
1028 total++;
1029 }
1030
1031 /* finally give numbers to points 0,1,2 */
1032 /* we will not use these points */
1033 renumb[0]=j++;
1034 renumb[1]=j++;
1035 renumb[2]=j++;
1036
1037 /* now construct re-renmbering array */
1038 i = 0;
1039 /* do an insertion sort */
1040 while ( i < M.count ) {
1041 for ( j = 0; j < M.count; j++ ) {
1042 if ( renumb[j] == i )
1043 break;
1044 }
1045 re_renumb[i] = j;
1046 i++;
1047 }
1048
1049 fprintf(outf,"## Generated by "PACKAGE" "VERSION"\n");
1050 fprintf(outf,"## Total points %d, (known points %d, unknown points %d) Total triangles %d of %d\n",total,unknowns,total-unknowns,triangles,ntriangles);
1051
1052 /* now print the unknown point coords, with the number */
1053 fprintf(outf,"## unknown points: point number | x | y\n");
1054
1055 for ( i = 0; i < unknowns; i++ )
1056 fprintf(outf,"%d "MDF" "MDF"\n",i, (Mx(re_renumb[i],MM))*g_xscale-g_xoff, (My(re_renumb[i],MM))*g_yscale-g_yoff );
1057 /* now the known coordinates with potentials */
1058 fprintf(outf,"## known points: point number | x | y | potential\n");
1059 for ( i = unknowns; i < total; i++ )
1060 fprintf(outf,"%d "MDF" "MDF" "MDF"\n",i, (Mx(re_renumb[i],MM))*g_xscale-g_xoff, (My(re_renumb[i],MM))*g_yscale-g_yoff, Mz(re_renumb[i],MM));
1061 /* now the triangles data with rho and epsilon */
1062 fprintf(outf,"## triangles : triangle number | point 1 | point 2 | point 3 | rho| mu or epsilon\n");
1063
1064
1065 DAG_traverse_list_reset(&dt);
1066 i=0;
1067 t=DAG_traverse_prune_list(&dt);
1068 while(t) {
1069 if ( t && t->status==LIVE ) {
1070 fprintf(outf,MIF" "MIF" "MIF" "MIF" ",i,renumb[t->p0],renumb[t->p1],renumb[t->p2]);
1071 fprintf(outf,MDF" ",(bound.poly_array[t->boundary].rhoexp?
1072 (ex(bound.poly_array[t->boundary].rhoexp,t->p0)+ex(bound.poly_array[t->boundary].rhoexp,t->p0)+ex(bound.poly_array[t->boundary].rhoexp,t->p0))*ONE_OVER_THREE:
1073 bound.poly_array[t->boundary].rho));
1074 fprintf(outf,MDF"\n",bound.poly_array[t->boundary].mu);
1075 i++;
1076 }
1077 t=DAG_traverse_prune_list(&dt);
1078 }
1079
1080 fclose(outf);
1081 free(renumb);
1082 free(re_renumb);
1083
1084 }
1085
1086 /* printing potentials as ACII file*/
1087 void
print_potential(const char * filename,Mesh MM)1088 print_potential(const char* filename, Mesh MM)
1089 {
1090 /* file format
1091 first line: no. of points
1092 other lines: point number, x, y, potential
1093 */
1094
1095 FILE *outf;
1096 unsigned int i;
1097 if ((outf=fopen(filename,"w"))==NULL){
1098 fprintf(stderr,"%s: %d: could not open file %s\n",__FILE__,__LINE__,filename);
1099 return;
1100 }
1101
1102 fprintf(outf,"%d\n",M.count);
1103 for ( i = 0; i < M.count; i++ ) {
1104 fprintf(outf,"%d "MDF" "MDF" "MDF"\n",i, (Mx(i,MM))*g_xscale-g_xoff, (My(i,MM))*g_yscale-g_yoff, Mzz(i,current_plotting_contour,MM));
1105 }
1106
1107 fclose(outf);
1108 }
1109
1110
1111 /* a function to plot the gradient of potential */
1112 void
plot_gradient(Mesh MM)1113 plot_gradient(Mesh MM)
1114 {
1115 triangle *t;
1116 MY_INT n1,n2,n3;
1117 float r,g,bl;
1118 MY_DOUBLE del, a,b, yy12, yy23, yy31, x0,x1,x2,y0,y1,y2;
1119
1120 max_gradient=0.01;
1121 DAG_traverse_list_reset(&dt);
1122 t=DAG_traverse_prune_list(&dt);
1123 while(t) {
1124 if ( t && t->status==LIVE ) {
1125 #ifdef DEBUG
1126 printf("plotting gradient for "MIF" contour %d\n",t->n,current_plotting_contour);
1127 #endif
1128 /* calculate gradient */
1129 /* phi =a*x+b*y+c for each triangle */
1130 /* need to find a and b only */
1131 /* using cramer's rule for this */
1132 yy23 = My(t->p1,MM)-My(t->p2,MM);
1133 yy31 = My(t->p2,MM)-My(t->p0,MM);
1134 yy12 = My(t->p0,MM)-My(t->p1,MM);
1135 del = ABS(Mx(t->p0,MM)*yy23 + Mx(t->p1,MM)*yy31 + Mx(t->p2,MM)*yy12);
1136 if ( del < TOL ) {
1137 t=DAG_traverse_prune_list(&dt);
1138 continue; /* do not plot */
1139 }
1140 /* else */
1141 del = 0.001/del;
1142 a = del*(Mzz(t->p0,current_plotting_contour,MM)*yy23 + Mzz(t->p1,current_plotting_contour,MM)*yy31 + Mzz(t->p2,current_plotting_contour,MM)*yy12);
1143 b = del*(Mx(t->p0,MM)*(Mzz(t->p1,current_plotting_contour,MM)-Mzz(t->p2,current_plotting_contour,MM))
1144 +Mx(t->p1,MM)*(Mzz(t->p2,current_plotting_contour,MM)-Mzz(t->p0,current_plotting_contour,MM))
1145 +Mx(t->p2,MM)*(Mzz(t->p0,current_plotting_contour,MM)-Mzz(t->p1,current_plotting_contour,MM)));
1146
1147 /* now find the centroid of triangle */
1148 /* x coord */
1149 x0 = (Mx(t->p0,MM)+Mx(t->p1,MM)+Mx(t->p2,MM))*ONE_OVER_THREE;
1150 /* y coord */
1151 y0 = (My(t->p0,MM)+My(t->p1,MM)+My(t->p2,MM))*ONE_OVER_THREE;
1152
1153 /* now find the two intersection points of the gradient line */
1154 /* p2 */
1155 /* /\ */
1156 /* p1/ \ */
1157 /* /` \ */
1158 /* / ` \ */
1159 /* / ` \ */
1160 /* p0 -------`----p1 */
1161 /* p2 */
1162
1163 /* at each corner, evaluate line value */
1164 /* check to see if n1 on line */
1165 yy12= a*(My(t->p0,MM)-y0) - b*(Mx(t->p0,MM)-x0);
1166 /* check to see if n2 on line */
1167 yy23= a*(My(t->p1,MM)-y0) - b*(Mx(t->p1,MM)-x0);
1168 /* check to see if n2 on line */
1169 yy31= a*(My(t->p2,MM)-y0) - b*(Mx(t->p2,MM)-x0);
1170
1171 /* now 2 of the above values will have opposite signes from the other one */
1172 /* so generalize */
1173 if ( ((yy12 >= TOL) && (yy23 < 0) && (yy31 < 0)) || ((yy12 <= TOL) && (yy23 > 0) && (yy31 > 0))) {
1174 n1 = t->p0;
1175 n2 = t->p1;
1176 n3 = t->p2;
1177 } else if ( ((yy23 >= TOL) && (yy12 < 0) && (yy31 < 0)) || ((yy23 <= TOL) && (yy12 > 0) && (yy31 > 0))) {
1178 n1 = t->p1;
1179 n2 = t->p2;
1180 n3 = t->p0;
1181 } else if ( ((yy31 >= TOL) && (yy12 < 0) && (yy23 < 0)) || ((yy31 <= TOL) && (yy12 > 0) && (yy23 > 0)) ) {
1182 n1 = t->p2;
1183 n2 = t->p0;
1184 n3 = t->p1;
1185 } else {
1186 t=DAG_traverse_prune_list(&dt);
1187 continue; /* error */
1188 }
1189 /* now check intersection between (n1,n2) and (n1,n3) */
1190 /* first (n1,n2) */
1191 yy12= a*(My(n1,MM)-y0) - b*(Mx(n1,MM)-x0);
1192 yy23= a*(My(n2,MM)-y0) - b*(Mx(n2,MM)-x0);
1193 if ( (yy23 < TOL) && (yy23 > -TOL) ) {
1194 x1 = Mx(n2,MM);
1195 y1 = My(n2,MM);
1196 } else {
1197 del = -yy12/yy23;
1198 yy23 = 1/(1+del);
1199 x1 = (Mx(n1,MM)+del*Mx(n2,MM))*yy23;
1200 y1 = (My(n1,MM)+del*My(n2,MM))*yy23;
1201 }
1202 /* now (n1,n3) */
1203 yy23= a*(My(n3,MM)-y0) - b*(Mx(n3,MM)-x0);
1204 if ( (yy23 < TOL) && (yy23 > -TOL) ) {
1205 x2 = Mx(n3,MM);
1206 y2 = My(n3,MM);
1207 } else {
1208 del = -yy12/yy23;
1209 yy12 = 1/(1+del);
1210 x2 = (Mx(n1,MM)+del*Mx(n3,MM))*yy12;
1211 y2 = (My(n1,MM)+del*My(n3,MM))*yy12;
1212 }
1213
1214 /* now we know the two endpoints */
1215 /* we draw an arrowhead in the direction point (x1,y1) */
1216 /* need to determine that now */
1217 if ( a > 0 ) {
1218 if ( x1 > x2 )
1219 n1 = 1; /* correct point */
1220 else
1221 n1 =2; /* need to switch */
1222 } else if ( a < 0 ) {
1223 if ( x1 < x2 )
1224 n1 = 1; /* correct point */
1225 else
1226 n1 =2; /* need to switch */
1227 } else { /* a == 0 */
1228 if ( b > 0 ) {
1229 if ( y1 > y2 )
1230 n1 = 1; /* correct point */
1231 else
1232 n1 =2; /* need to switch */
1233 } else if ( b < 0 ) {
1234 if ( y1 < y2 )
1235 n1 = 1; /* correct point */
1236 else
1237 n1 =2; /* need to switch */
1238 }
1239 }
1240
1241 /* need to switch points? */
1242 if ( n1 == 2 ) {
1243 x0 = x1;
1244 y0 = y1;
1245 x1 = x2;
1246 y1 = y2;
1247 x2 = x0;
1248 y2 = y0;
1249 }
1250
1251
1252 /* arrow will look like this */
1253 /* (? , ?) */
1254 /* |\ */
1255 /* | \ */
1256 /* (x2,y2)----------------- \ (x1,y1) */
1257 /* (x0,y0)| / */
1258 /* |/ */
1259 /* (? ?) */
1260
1261 /* keep 5:1 ratio */
1262 del=5;
1263 x0 = (x1*del+x2)/(del+1);
1264 y0 = (y1*del+y2)/(del+1);
1265
1266 /* calculate value of gradient */
1267 a = sqrt(a*a + b*b);
1268 /* if this is higher than the maximum, update it */
1269 if ( max_gradient < a )
1270 max_gradient = a;
1271 /* set colour */
1272 get_colour(&r,&g,&bl,(int)(100*a), (int)(100*max_gradient));
1273 glColor3f(r, g, bl);
1274 /* arrow size */
1275 a = 0.2*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
1276 /* gradient of normal */
1277 /* temporary variables */
1278 if ( y2 != y1 ) {
1279 del = -(x2-x1)/(y2-y1);
1280 yy12 = sqrt(1+del*del);
1281 yy23 = a/yy12;
1282 yy12 = del*a/yy12;
1283 } else {
1284 yy23 = 0;
1285 yy12 = a;
1286 }
1287
1288
1289 /* now plot gradient line (x1,y1) (x2,y2) */
1290 glBegin(GL_LINE_LOOP);
1291 /* draw line */
1292 glVertex2d( (GLfloat)(x2), (GLfloat)(y2) );
1293 glVertex2d( (GLfloat)(x1), (GLfloat)(y1) );
1294 glEnd();
1295 /* now arrow head*/
1296 glBegin(GL_TRIANGLES);
1297 glVertex2d( (GLfloat)(x0+yy23), (GLfloat)(y0+yy12) );
1298 glVertex2d( (GLfloat)(x0-yy23), (GLfloat)(y0-yy12) );
1299 glVertex2d( (GLfloat)(x1), (GLfloat)(y1) );
1300 glEnd();
1301 }
1302
1303 t=DAG_traverse_prune_list(&dt);
1304 }
1305 }
1306 /**********************************************/
1307 /* legend plotting routines */
1308 #ifndef WIN32
1309 static void
realize(GtkWidget * widget,gpointer data)1310 realize(GtkWidget *widget, gpointer data)
1311 {
1312 GdkGLContext *glcontext =gtk_widget_get_gl_context(widget);
1313 GdkGLDrawable *gldrawable =gtk_widget_get_gl_drawable(widget);
1314 if (!gdk_gl_drawable_gl_begin(gldrawable,glcontext))
1315 return;
1316 gdk_gl_drawable_gl_end(gldrawable);
1317 }
1318
1319 static gboolean
configure_event(GtkWidget * widget,GdkEventConfigure * event,gpointer data)1320 configure_event(GtkWidget *widget, GdkEventConfigure *event, gpointer data)
1321 {
1322 GdkGLContext *glcontext=gtk_widget_get_gl_context(widget);
1323 GdkGLDrawable *gldrawable=gtk_widget_get_gl_drawable(widget);
1324 GLfloat w=widget->allocation.width;
1325 GLfloat h=widget->allocation.height;
1326 if (!gdk_gl_drawable_gl_begin(gldrawable,glcontext))
1327 return(FALSE);
1328 glViewport(0,0,w,h);
1329 glMatrixMode(GL_PROJECTION);
1330 glLoadIdentity();
1331 gluOrtho2D(-1,1,-1,1);
1332 glMatrixMode(GL_MODELVIEW);
1333 gdk_gl_drawable_gl_end(gldrawable);
1334 return(TRUE);
1335 }
1336
1337 static gboolean
expose_event(GtkWidget * widget,GdkEventExpose * event,gpointer data)1338 expose_event(GtkWidget *widget,
1339 GdkEventExpose *event,
1340 gpointer data)
1341 {
1342 int i;
1343 float r,g,b;
1344 GLfloat dy; /* y step */
1345 GLfloat y1,y2,x1,x2,scale;
1346 GdkGLDrawable *gldrawable =gtk_widget_get_gl_drawable(widget);
1347 GdkGLContext *glcontext=gtk_widget_get_gl_context(widget);
1348
1349 GLfloat w=(widget->allocation.width);
1350 GLfloat h=(widget->allocation.height);
1351
1352 if (!gdk_gl_drawable_gl_begin(gldrawable,glcontext))
1353 return(FALSE);
1354
1355 /* rescale */
1356 if ( h >w )
1357 scale=1/h;
1358 else
1359 scale=1/w;
1360 w=2;
1361 h=2;
1362
1363 /* top, bottom border =h/20 */
1364 /* left border = w/10, color block=9w/10, right border=w/10 */
1365 /* drawing */
1366 dy=(h*0.9)/(GLfloat)contour_levels;
1367
1368 y2=-0.45*h;
1369 y1=dy+y2;
1370 x1=-0.45*w;
1371 x2=-x1;
1372 glClear(GL_COLOR_BUFFER_BIT);
1373 #ifdef DEBUG
1374 printf("expose_ev: legend: rectangle (%f,%f)--(%f,%f)\n",x1,y1,x2,y2);
1375 #endif
1376
1377
1378 /* horizontal line */
1379 glColor3f(0.0f,0.0f,0.0f);
1380 glBegin(GL_LINES);
1381 glVertex2d(x1,y2);
1382 glVertex2d(x1+0.9*w,y2);
1383 glEnd();
1384
1385 for (i=0; i<contour_levels; i++) {
1386 get_colour(&r,&g,&b,i,contour_levels);
1387 glColor3f((GLfloat)r,(GLfloat)g,(GLfloat)b);
1388 glBegin(GL_QUADS);
1389 glVertex2d(x1,y1);
1390 glVertex2d(x2,y1);
1391 glVertex2d(x2,y2);
1392 glVertex2d(x1,y2);
1393 glEnd();
1394
1395 /* horizontal line */
1396 glColor3f(0.0f,0.7f,0.7f);
1397 glBegin(GL_LINES);
1398 glVertex2d(x1,y1);
1399 glVertex2d(x1+0.0*w,y1);
1400 glEnd();
1401
1402 y1+=dy;
1403 y2+=dy;
1404
1405 }
1406 if ( gdk_gl_drawable_is_double_buffered(gldrawable))
1407 gdk_gl_drawable_swap_buffers(gldrawable);
1408 else
1409 glFlush();
1410 gdk_gl_drawable_gl_end(gldrawable);
1411 return(TRUE);
1412 }
1413
1414 /* mouse motion */
1415 static gboolean
motion_notify_event(GtkWidget * widget,GdkEventMotion * event,gpointer data)1416 motion_notify_event(GtkWidget *widget,
1417 GdkEventMotion *event,
1418 gpointer data)
1419 {
1420 GtkWidget *label;
1421 GLfloat h,y;
1422 int i;
1423 gchar buf[128];
1424 h=widget->allocation.height;
1425
1426 y=-2*event->y/h+1.0;
1427 label=(GtkWidget*)data;
1428 /* determine exact contour we are in */
1429 i=(int)floor((double)(y+0.9)*(double)contour_levels/1.8);
1430 /* sanity check: i has to be in [0,contour_levels-1]*/
1431 if ( (i>= 0) && (i<contour_levels)) {
1432 g_snprintf(buf,128,"%4.2f",(float)i*(g_maxpot-g_minpot)/(float)contour_levels+g_minpot);
1433 /* printf("%4.3f\n",(float)i*(g_maxpot-g_minpot)/(float)contour_levels+g_minpot); */
1434 gtk_label_set_markup(GTK_LABEL(label),buf);
1435 return(TRUE);
1436 }
1437 return(FALSE);
1438 }
1439
1440 static gint
close_legend_window(GtkWidget * widget,GdkEvent * event,gpointer * data)1441 close_legend_window(GtkWidget *widget, GdkEvent *event, gpointer *data)
1442 {
1443 return(FALSE);
1444 }
1445
1446 /* plotting legend */
1447 void
display_legend_window(void)1448 display_legend_window(void)
1449 {
1450 GtkWidget *window;
1451 GtkWidget *vbox;
1452 GtkWidget *drawing_area;
1453 GtkWidget *legend_label;
1454 GdkGLConfig *glconfig;
1455
1456
1457 /* Try double-buffered visual */
1458 glconfig = gdk_gl_config_new_by_mode (GDK_GL_MODE_RGB |
1459 GDK_GL_MODE_DOUBLE);
1460 if (glconfig == NULL)
1461 {
1462 g_print ("\n*** Cannot find the double-buffered visual.\n");
1463 g_print ("\n*** Trying single-buffered visual.\n");
1464
1465 /* Try single-buffered visual */
1466 glconfig = gdk_gl_config_new_by_mode (GDK_GL_MODE_RGB);
1467 if (glconfig == NULL)
1468 {
1469 g_print ("*** No appropriate OpenGL-capable visual found.\n");
1470 exit(1);
1471 }
1472 }
1473
1474
1475 /*
1476 * Top-level window.
1477 */
1478
1479 window = gtk_window_new (GTK_WINDOW_TOPLEVEL);
1480 gtk_window_set_title (GTK_WINDOW (window), "Legend");
1481
1482 /* Get automatically redrawn if any of their children changed allocation. */
1483 gtk_container_set_reallocate_redraws (GTK_CONTAINER (window), TRUE);
1484
1485 /* Connect signal handlers to the window */
1486 g_signal_connect (G_OBJECT (window), "destroy",
1487 G_CALLBACK (gtk_widget_destroy), (gpointer)window);
1488 gtk_signal_connect (GTK_OBJECT (window), "destroy",
1489 GTK_SIGNAL_FUNC (close_legend_window), NULL);
1490
1491 /*
1492 * VBox.
1493 */
1494
1495 vbox = gtk_vbox_new (FALSE, 0);
1496 gtk_container_add (GTK_CONTAINER (window), vbox);
1497 gtk_widget_show (vbox);
1498
1499 /*
1500 * Drawing area to draw OpenGL scene.
1501 */
1502
1503 drawing_area = gtk_drawing_area_new ();
1504 gtk_widget_set_size_request (drawing_area, LEGEND_WIDTH, LEGEND_HEIGHT);
1505
1506 /* Set OpenGL-capability to the widget */
1507 gtk_widget_set_gl_capability (drawing_area,
1508 glconfig,
1509 NULL,
1510 TRUE,
1511 GDK_GL_RGBA_TYPE);
1512
1513 gtk_widget_add_events (drawing_area,
1514 GDK_POINTER_MOTION_MASK | /* any pointer motion */
1515 GDK_VISIBILITY_NOTIFY_MASK);
1516
1517 /* Connect signal handlers to the drawing area */
1518 g_signal_connect_after (G_OBJECT (drawing_area), "realize",
1519 G_CALLBACK (realize), NULL);
1520 g_signal_connect (G_OBJECT (drawing_area), "configure_event",
1521 G_CALLBACK (configure_event), NULL);
1522 g_signal_connect (G_OBJECT (drawing_area), "expose_event",
1523 G_CALLBACK (expose_event), NULL);
1524
1525 gtk_box_pack_start (GTK_BOX (vbox), drawing_area, TRUE, TRUE, 0);
1526 gtk_widget_show (drawing_area);
1527
1528 /* label to show potential */
1529 legend_label=gtk_label_new("000");
1530 gtk_label_set_justify(GTK_LABEL(legend_label),GTK_JUSTIFY_CENTER);
1531 gtk_box_pack_start(GTK_BOX(vbox),legend_label,FALSE,TRUE,0);
1532 gtk_widget_show(legend_label);
1533
1534 /* connect mouse motion */
1535 g_signal_connect (G_OBJECT(drawing_area),"motion_notify_event",
1536 G_CALLBACK(motion_notify_event), (gpointer)legend_label);
1537
1538
1539 gtk_widget_show(window);
1540
1541 }
1542 #endif/*WIN32*/
1543
1544 /* a function to print the gradient of potential */
1545 void
print_gradient(const char * filename,Mesh MM)1546 print_gradient(const char *filename, Mesh MM)
1547 {
1548 triangle *t;
1549 MY_INT n1,n2,n3;
1550 MY_DOUBLE del, a,b, yy12, yy23, yy31, x0,x1,x2,y0,y1,y2;
1551 FILE *plotf;
1552 MY_DOUBLE myscale;
1553 int i;
1554
1555 if ((plotf=fopen(filename,"w"))==NULL){
1556 fprintf(stderr,"%s: %d: could not open file %s\n",__FILE__,__LINE__,filename);
1557 return;
1558 }
1559
1560
1561 /* determine the right bound of bounding box */
1562 if ( g_xrange > g_yrange )
1563 myscale = 400.0/g_xrange;
1564 else
1565 myscale = 400.0/g_yrange;
1566
1567 /* print the header */
1568 fprintf(plotf,"%%!PS-Adobe-2.0 EPSF-3.0\n");
1569 fprintf(plotf,"%%%%Title: %s\n",filename);
1570 fprintf(plotf,"%%%%Generated By: "PACKAGE" "VERSION" contour levels: %d, min= "MDF", max= "MDF"\n", contour_levels, g_minpot, g_maxpot);
1571 fprintf(plotf,"%%%%BoundingBox: 0 0 %4.2lf %4.2lf\n",(double)(2*g_xrange*myscale+20), (double)(2*g_yrange*myscale+20) );
1572 fprintf(plotf,"%%%%Magnification: 1.0000\n");
1573
1574
1575
1576 /* start plotting from 10 10 */
1577 fprintf(plotf,"10 10 translate\n");
1578
1579 /* now print some shorter definitions */
1580 fprintf(plotf,"/l {lineto} bind def\n"); /* lineto */
1581 fprintf(plotf,"/f {fill} bind def\n"); /* lineto */
1582 fprintf(plotf,"/m {moveto} bind def\n"); /* moveto */
1583 fprintf(plotf,"/n {newpath} bind def\n"); /* newpath */
1584 fprintf(plotf,"/s {stroke} bind def\n"); /* stroke */
1585 fprintf(plotf,"/slw {setlinewidth} bind def\n"); /* setlinewidth */
1586 fprintf(plotf,"/cp {closepath} bind def\n"); /* closepath */
1587
1588 fprintf(plotf,"1 slw\n");
1589 DAG_traverse_list_reset(&dt);
1590 t=DAG_traverse_prune_list(&dt);
1591 while(t) {
1592 if ( t && t->status==LIVE ) {
1593 #ifdef DEBUG
1594 printf("plotting gradient for "MIF" contour %d\n",t->n,current_plotting_contour);
1595 #endif
1596 /* calculate gradient */
1597 /* phi =a*x+b*y+c for each triangle */
1598 /* need to find a and b only */
1599 /* using cramer's rule for this */
1600 yy23 = My(t->p1,MM)-My(t->p2,MM);
1601 yy31 = My(t->p2,MM)-My(t->p0,MM);
1602 yy12 = My(t->p0,MM)-My(t->p1,MM);
1603 del = ABS(Mx(t->p0,MM)*yy23 + Mx(t->p1,MM)*yy31 + Mx(t->p2,MM)*yy12);
1604 if ( del < TOL ) {
1605 t=DAG_traverse_prune_list(&dt);
1606 continue; /* do not plot */
1607 }
1608 /* else */
1609 del = 0.001/del;
1610 a = del*(Mzz(t->p0,current_plotting_contour,MM)*yy23 + Mzz(t->p1,current_plotting_contour,MM)*yy31 + Mzz(t->p2,current_plotting_contour,MM)*yy12);
1611 b = del*(Mx(t->p0,MM)*(Mzz(t->p1,current_plotting_contour,MM)-Mzz(t->p2,current_plotting_contour,MM))
1612 +Mx(t->p1,MM)*(Mzz(t->p2,current_plotting_contour,MM)-Mzz(t->p0,current_plotting_contour,MM))
1613 +Mx(t->p2,MM)*(Mzz(t->p0,current_plotting_contour,MM)-Mzz(t->p1,current_plotting_contour,MM)));
1614
1615 /* now find the centroid of triangle */
1616 /* x coord */
1617 x0 = (Mx(t->p0,MM)+Mx(t->p1,MM)+Mx(t->p2,MM))*ONE_OVER_THREE;
1618 /* y coord */
1619 y0 = (My(t->p0,MM)+My(t->p1,MM)+My(t->p2,MM))*ONE_OVER_THREE;
1620
1621 /* now find the two intersection points of the gradient line */
1622 /* p2 */
1623 /* /\ */
1624 /* p1/ \ */
1625 /* /` \ */
1626 /* / ` \ */
1627 /* / ` \ */
1628 /* p0 -------`----p1 */
1629 /* p2 */
1630
1631 /* at each corner, evaluate line value */
1632 /* check to see if n1 on line */
1633 yy12= a*(My(t->p0,MM)-y0) - b*(Mx(t->p0,MM)-x0);
1634 /* check to see if n2 on line */
1635 yy23= a*(My(t->p1,MM)-y0) - b*(Mx(t->p1,MM)-x0);
1636 /* check to see if n2 on line */
1637 yy31= a*(My(t->p2,MM)-y0) - b*(Mx(t->p2,MM)-x0);
1638
1639 /* now 2 of the above values will have opposite signes from the other one */
1640 /* so generalize */
1641 if ( ((yy12 >= TOL) && (yy23 < 0) && (yy31 < 0)) || ((yy12 <= TOL) && (yy23 > 0) && (yy31 > 0))) {
1642 n1 = t->p0;
1643 n2 = t->p1;
1644 n3 = t->p2;
1645 } else if ( ((yy23 >= TOL) && (yy12 < 0) && (yy31 < 0)) || ((yy23 <= TOL) && (yy12 > 0) && (yy31 > 0))) {
1646 n1 = t->p1;
1647 n2 = t->p2;
1648 n3 = t->p0;
1649 } else if ( ((yy31 >= TOL) && (yy12 < 0) && (yy23 < 0)) || ((yy31 <= TOL) && (yy12 > 0) && (yy23 > 0)) ) {
1650 n1 = t->p2;
1651 n2 = t->p0;
1652 n3 = t->p1;
1653 } else {
1654 t=DAG_traverse_prune_list(&dt);
1655 continue; /* error */
1656 }
1657 /* now check intersection between (n1,n2) and (n1,n3) */
1658 /* first (n1,n2) */
1659 yy12= a*(My(n1,MM)-y0) - b*(Mx(n1,MM)-x0);
1660 yy23= a*(My(n2,MM)-y0) - b*(Mx(n2,MM)-x0);
1661 if ( (yy23 < TOL) && (yy23 > -TOL) ) {
1662 x1 = Mx(n2,MM);
1663 y1 = My(n2,MM);
1664 } else {
1665 del = -yy12/yy23;
1666 yy23 = 1/(1+del);
1667 x1 = (Mx(n1,MM)+del*Mx(n2,MM))*yy23;
1668 y1 = (My(n1,MM)+del*My(n2,MM))*yy23;
1669 }
1670 /* now (n1,n3) */
1671 yy23= a*(My(n3,MM)-y0) - b*(Mx(n3,MM)-x0);
1672 if ( (yy23 < TOL) && (yy23 > -TOL) ) {
1673 x2 = Mx(n3,MM);
1674 y2 = My(n3,MM);
1675 } else {
1676 del = -yy12/yy23;
1677 yy12 = 1/(1+del);
1678 x2 = (Mx(n1,MM)+del*Mx(n3,MM))*yy12;
1679 y2 = (My(n1,MM)+del*My(n3,MM))*yy12;
1680 }
1681
1682 /* now we know the two endpoints */
1683 /* we draw an arrowhead in the direction point (x1,y1) */
1684 /* need to determine that now */
1685 if ( a > 0 ) {
1686 if ( x1 > x2 )
1687 n1 = 1; /* correct point */
1688 else
1689 n1 =2; /* need to switch */
1690 } else if ( a < 0 ) {
1691 if ( x1 < x2 )
1692 n1 = 1; /* correct point */
1693 else
1694 n1 =2; /* need to switch */
1695 } else { /* a == 0 */
1696 if ( b > 0 ) {
1697 if ( y1 > y2 )
1698 n1 = 1; /* correct point */
1699 else
1700 n1 =2; /* need to switch */
1701 } else if ( b < 0 ) {
1702 if ( y1 < y2 )
1703 n1 = 1; /* correct point */
1704 else
1705 n1 =2; /* need to switch */
1706 }
1707 }
1708
1709 /* need to switch points? */
1710 if ( n1 == 2 ) {
1711 x0 = x1;
1712 y0 = y1;
1713 x1 = x2;
1714 y1 = y2;
1715 x2 = x0;
1716 y2 = y0;
1717 }
1718
1719
1720 /* arrow will look like this */
1721 /* (? , ?) */
1722 /* |\ */
1723 /* | \ */
1724 /* (x2,y2)----------------- \ (x1,y1) */
1725 /* (x0,y0)| / */
1726 /* |/ */
1727 /* (? ?) */
1728
1729 /* keep 5:1 ratio */
1730 del=5;
1731 x0 = (x1*del+x2)/(del+1);
1732 y0 = (y1*del+y2)/(del+1);
1733
1734 /* calculate value of gradient */
1735 a = sqrtf((float)(a*a + b*b));
1736 /* arrow size */
1737 a = 0.2*sqrt((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
1738 /* gradient of normal */
1739 /* temporary variables */
1740 if ( y2 != y1 ) {
1741 del = -(x2-x1)/(y2-y1);
1742 yy12 = sqrt(1+del*del);
1743 yy23 = a/yy12;
1744 yy12 = del*a/yy12;
1745 } else {
1746 yy23 = 0;
1747 yy12 = a;
1748 }
1749
1750
1751 /* now plot gradient line (x1,y1) (x2,y2) */
1752 fprintf(plotf,"n ");
1753 fprintf(plotf,"%lf %lf m ",(double) (g_xrange+x2*g_xscale)*myscale, (double)(g_yrange + y2*g_yscale)*myscale);
1754 fprintf(plotf,"%lf %lf l ",(double) (g_xrange+x0*g_xscale)*myscale, (double)(g_yrange + y0*g_yscale)*myscale);
1755 fprintf(plotf,"cp s\n");
1756 /* now arrow head*/
1757 fprintf(plotf,"n ");
1758 fprintf(plotf,"%lf %lf m ",(double) (g_xrange+(x0+yy23)*g_xscale)*myscale, (double)(g_yrange +(y0+yy12)*g_yscale)*myscale);
1759 fprintf(plotf,"%lf %lf l ",(double) (g_xrange+(x0-yy23)*g_xscale)*myscale, (double)(g_yrange +(y0-yy12)*g_yscale)*myscale);
1760 fprintf(plotf,"%lf %lf l ",(double) (g_xrange+x1*g_xscale)*myscale, (double)(g_yrange + y1*g_yscale)*myscale);
1761 fprintf(plotf,"cp s\n");
1762 }
1763
1764 t=DAG_traverse_prune_list(&dt);
1765 }
1766
1767 fprintf(plotf,"2 slw\n");
1768 /* print boundaries */
1769 for ( i = 0; i < bound.npoly; i++ ) {
1770 print_polygon(&bound.poly_array[i],myscale,plotf,MM);
1771 }
1772
1773 fprintf(plotf,"%%%%EOF\n");
1774 fclose(plotf);
1775
1776 }
1777