1 /***************************************************************************
2 JSPICE3 adaptation of Spice3e2 - Copyright (c) Stephen R. Whiteley 1992
3 Copyright 1990 Regents of the University of California. All rights reserved.
4 Authors: 1987 UCB
5 1992 Stephen R. Whiteley
6 ****************************************************************************/
7
8 /*
9 * Routines to draw the various sorts of grids -- linear, log, polar.
10 */
11
12 #include "spice.h"
13 #include "ftedefs.h"
14 #include "plotdev.h"
15 #include "plotext.h"
16
17 typedef enum { x_axis, y_axis } Axis;
18
19 #ifdef __STDC__
20 static double *lingrid(GRAPH*,double,double,double,int,Axis);
21 static void drawlingrid(GRAPH*,char[16],int,int,int,int,int,bool,int,Axis);
22 static double *loggrid(GRAPH*,double,double,int,Axis);
23 static void drawloggrid(GRAPH*,int,int,int,int,int,Axis);
24 static void polargrid(GRAPH*);
25 static void drawpolargrid(GRAPH*);
26 static void adddeglabel(GRAPH*,int,int,int,int,int);
27 static void addradlabel(GRAPH*,int,double,int,int);
28 static void smithgrid(GRAPH*);
29 static void drawsmithgrid(GRAPH*);
30 static void arcset(int,int,GRAPH*,int,int,char*,char*);
31 static void cliparc(int,int,int,double,double,int,int,int);
32 #else
33 static double *lingrid();
34 static void drawlingrid();
35 static double *loggrid();
36 static void drawloggrid();
37 static void polargrid();
38 static void drawpolargrid();
39 static void adddeglabel();
40 static void addradlabel();
41 static void smithgrid();
42 static void drawsmithgrid();
43 static void arcset();
44 static void cliparc();
45 #endif
46
47
48 /* note: scaleunits is never changed in this file */
49 static bool scaleunits = true;
50
51
52 void
gr_fixgrid(graph,xdelta,ydelta,xtype,ytype)53 gr_fixgrid(graph, xdelta, ydelta, xtype, ytype)
54
55 GRAPH *graph;
56 double xdelta, ydelta;
57 int xtype, ytype;
58 {
59 double *dd;
60
61 if (graph->grid.gridtype == GRID_NONE) {
62 graph->grid.gridtype = GRID_LIN;
63 }
64
65 DevSetColor(1);
66 DevSetLinestyle(1);
67
68 if ((graph->datawindow.xmin > graph->datawindow.xmax)
69 || (graph->datawindow.ymin > graph->datawindow.ymax)) {
70 fprintf(cp_err,
71 "gr_fixgrid: Internal Error - bad limits: %lg, %lg, %lg, %lg\r\n",
72 graph->datawindow.xmin, graph->datawindow.xmax,
73 graph->datawindow.ymin, graph->datawindow.xmax);
74 return;
75 }
76
77 if (graph->grid.gridtype == GRID_POLAR) {
78 graph->grid.circular = true;
79 polargrid(graph);
80 return;
81 }
82 else if (graph->grid.gridtype == GRID_SMITH) {
83 graph->grid.circular = true;
84 smithgrid(graph);
85 return;
86 }
87 graph->grid.circular = false;
88
89 if ((graph->grid.gridtype == GRID_XLOG)
90 || (graph->grid.gridtype == GRID_LOGLOG))
91 dd = loggrid(graph, graph->datawindow.xmin,
92 graph->datawindow.xmax,
93 xtype, x_axis);
94 else
95 dd = lingrid(graph, graph->datawindow.xmin,
96 graph->datawindow.xmax,
97 xdelta, xtype, x_axis);
98 graph->datawindow.xmin = dd[0];
99 graph->datawindow.xmax = dd[1];
100 if ((graph->grid.gridtype == GRID_YLOG)
101 || (graph->grid.gridtype == GRID_LOGLOG))
102 dd = loggrid(graph, graph->datawindow.ymin,
103 graph->datawindow.ymax,
104 ytype, y_axis);
105 else
106 dd = lingrid(graph, graph->datawindow.ymin,
107 graph->datawindow.ymax,
108 ydelta, ytype, y_axis);
109 graph->datawindow.ymin = dd[0];
110 graph->datawindow.ymax = dd[1];
111
112 return;
113 }
114
115
116 void
gr_redrawgrid(graph)117 gr_redrawgrid(graph)
118
119 GRAPH *graph;
120 {
121 DevSetColor(1);
122 DevSetLinestyle(1);
123 /* draw labels */
124 if (graph->grid.xlabel) {
125 DevText(graph->grid.xlabel,
126 graph->absolute.width -
127 (strlen(graph->grid.xlabel) + 3) * graph->fontwidth,
128 graph->fontheight);
129 }
130 if (graph->grid.ylabel) {
131 if (graph->grid.gridtype == GRID_POLAR
132 || graph->grid.gridtype == GRID_SMITH ) {
133 DevText(graph->grid.ylabel,
134 graph->fontwidth,
135 graph->absolute.height * 3 / 4 );
136 }
137 else {
138 DevText(graph->grid.ylabel,
139 graph->fontwidth,
140 graph->absolute.height / 2 );
141 }
142 }
143
144 switch ( graph->grid.gridtype ) {
145
146 case GRID_POLAR:
147 drawpolargrid(graph);
148 break;
149
150 case GRID_SMITH:
151 drawsmithgrid(graph);
152 break;
153
154 case GRID_XLOG:
155 case GRID_LOGLOG:
156 drawloggrid(graph, graph->grid.xaxis.log.hmt,
157 graph->grid.xaxis.log.lmt,
158 graph->grid.xaxis.log.decsp,
159 graph->grid.xaxis.log.subs,
160 graph->grid.xaxis.log.pp, x_axis);
161 break;
162
163 default:
164 drawlingrid(graph,
165 graph->grid.xaxis.lin.units,
166 graph->grid.xaxis.lin.spacing,
167 graph->grid.xaxis.lin.numspace,
168 graph->grid.xaxis.lin.distance,
169 graph->grid.xaxis.lin.lowlimit,
170 graph->grid.xaxis.lin.highlimit,
171 graph->grid.xaxis.lin.onedec,
172 graph->grid.xaxis.lin.mult,
173 x_axis);
174 break;
175 }
176
177 switch ( graph->grid.gridtype ) {
178 case GRID_POLAR:
179 case GRID_SMITH:
180 break;
181
182 case GRID_YLOG:
183 case GRID_LOGLOG:
184 drawloggrid(graph, graph->grid.yaxis.log.hmt,
185 graph->grid.yaxis.log.lmt,
186 graph->grid.yaxis.log.decsp,
187 graph->grid.yaxis.log.subs,
188 graph->grid.yaxis.log.pp, y_axis);
189 break;
190
191 default:
192 drawlingrid(graph,
193 graph->grid.yaxis.lin.units,
194 graph->grid.yaxis.lin.spacing,
195 graph->grid.yaxis.lin.numspace,
196 graph->grid.yaxis.lin.distance,
197 graph->grid.yaxis.lin.lowlimit,
198 graph->grid.yaxis.lin.highlimit,
199 graph->grid.yaxis.lin.onedec,
200 graph->grid.yaxis.lin.mult,
201 y_axis);
202 break;
203 }
204
205 }
206
207
208 static double *
lingrid(graph,lo,hi,delta,type,axis)209 lingrid(graph, lo, hi, delta, type, axis)
210
211 /* Plot a linear grid. Returns the new hi and lo limits. */
212 GRAPH *graph;
213 double lo, hi;
214 double delta;
215 int type;
216 Axis axis;
217 {
218 int mag, hmt, lmt, nsp, spacing, dst;
219 double tenpowmag;
220 bool onedec = false;
221 int margin;
222 int max;
223 static double dd[2];
224 int mult = 1;
225 char buf[16], *s;
226
227 if (axis == x_axis) {
228 margin = graph->viewportxoff;
229 max = graph->viewport.width + graph->viewportxoff;
230 }
231 else {
232 margin = graph->viewportyoff;
233 max = graph->viewport.height + graph->viewportyoff;
234 }
235
236 if (delta < 0.0) {
237 fprintf(cp_err, "Warning: %cdelta is negative -- ignored\n",
238 (axis == x_axis) ? 'x' : 'y');
239 delta = 0;
240 }
241
242 /* Express the difference between the high and low values as
243 * diff = d * 10^mag. We know diff >= 0.0. If scaleunits is
244 * set then make sure that mag is a multiple of 3.
245 */
246 mag = floor(mylog10(hi - lo));
247 tenpowmag = pow(10.0, (double) mag);
248
249 if (scaleunits) {
250 if ((mag % 3) && !((mag - 1) % 3)) {
251 mult = 10;
252 mag--;
253 }
254 else if ((mag % 3) && !((mag + 1) % 3)) {
255 onedec = true;
256 mag++;
257 }
258 }
259
260 /* Now we have to round lo down, and hi up */
261 lmt = (int) floor(lo / tenpowmag);
262 hmt = (int) ceil(hi / tenpowmag);
263 if (fabs((hi - tenpowmag)/tenpowmag) < 1e-20) {
264 /* This is a hack to fix a strange case */
265 hmt = 1;
266 }
267 dst = hmt - lmt;
268
269 /* This is a strange case. */
270 if (dst == 1) {
271 dst = 10;
272 mag--;
273 tenpowmag /= 10;
274 hmt *= 10.0;
275 lmt *= 10.0;
276 }
277 dd[0] = lo = lmt * tenpowmag;
278 dd[1] = hi = hmt * tenpowmag;
279
280 /* We have to go from lmt to hmt, so think of some useful places
281 * to put grid lines. We will have a total of nsp lines, one
282 * every spacing pixels, which is every dst / nsp units.
283 */
284 if (delta == 0.0) {
285 for (nsp = 5; nsp < 10; nsp++)
286 if (!(dst % nsp))
287 break;
288 if (nsp == 10)
289 for (nsp = 2; nsp < 5; nsp++)
290 if (!(dst % nsp))
291 break;
292 /* Aaaaaghhh... */
293 if ((dst == 11) || (dst == 13) || (dst == 17) || (dst == 19))
294 nsp = dst;
295 spacing = (max - margin) / nsp;
296 }
297 else {
298 /* The user told us where to put the grid lines. They will
299 * not be equally spaced in this case (i.e, the right edge
300 * won't be a line).
301 */
302 nsp = (hi - lo) / delta;
303 spacing = (max - margin) * delta / (hi - lo);
304 }
305
306 /* In case we have a relatively small delta */
307 while (nsp > dst) {
308 dst *= 10;
309 mag--;
310 hmt *= 10.0;
311 lmt *= 10.0;
312 }
313
314 /* Reset the max X coordinate to deal with round-off error. */
315 if (nsp && (delta == 0.0)) {
316 if (axis == x_axis)
317 graph->viewport.width = spacing * nsp;
318 else
319 graph->viewport.height = spacing * nsp;
320 }
321 else if (!nsp) {
322 nsp = 1;
323 }
324
325 if (scaleunits) {
326 /* Figure out what to put here */
327 switch (mag) {
328 case -18: /* atto */
329 (void) strcpy(buf, "a");
330 break;
331 case -15: /* femto */
332 (void) strcpy(buf, "f");
333 break;
334 case -12: /* pico */
335 (void) strcpy(buf, "p");
336 break;
337 case -9: /* nano */
338 (void) strcpy(buf, "n");
339 break;
340 case -6: /* micro */
341 (void) strcpy(buf, "u");
342 break;
343 case -3: /* milli */
344 (void) strcpy(buf, "m");
345 break;
346 case 0: /* -- */
347 (void) strcpy(buf, "");
348 break;
349 case 3: /* kilo */
350 (void) strcpy(buf, "K");
351 break;
352 case 6: /* mega */
353 (void) strcpy(buf, "M");
354 break;
355 case 9: /* giga */
356 (void) strcpy(buf, "G");
357 break;
358 case 12: /* terra */
359 (void) strcpy(buf, "T");
360 break;
361 case 15: /* peta */
362 (void) strcpy(buf, "P");
363 break;
364 case 18: /* exa */
365 (void) strcpy(buf, "E");
366 break;
367 default:
368 (void) sprintf(buf, "e%d", mag);
369 break;
370 }
371 if (s = ft_typabbrev(type))
372 (void) strcat(buf, s);
373 else
374 /* Get rid of the prefix */
375 (void) sprintf(buf, "e%d", mag);
376 }
377 else
378 (void) sprintf(buf, "e%d", mag);
379
380 /* have to save non-intuitive variables left over
381 from old algorithms for redraws */
382 if (axis == x_axis) {
383 graph->grid.xaxis.lin.spacing = spacing;
384 graph->grid.xaxis.lin.numspace = nsp;
385 graph->grid.xaxis.lin.distance = dst;
386 graph->grid.xaxis.lin.lowlimit = lmt;
387 graph->grid.xaxis.lin.highlimit = hmt;
388 graph->grid.xaxis.lin.onedec = onedec;
389 graph->grid.xaxis.lin.mult = mult;
390 (void) strcpy(graph->grid.xaxis.lin.units, buf);
391 }
392 else {
393 graph->grid.yaxis.lin.spacing = spacing;
394 graph->grid.yaxis.lin.numspace = nsp;
395 graph->grid.yaxis.lin.distance = dst;
396 graph->grid.yaxis.lin.lowlimit = lmt;
397 graph->grid.yaxis.lin.highlimit = hmt;
398 graph->grid.yaxis.lin.onedec = onedec;
399 graph->grid.yaxis.lin.mult = mult;
400 (void) strcpy(graph->grid.yaxis.lin.units, buf);
401 }
402
403 return (dd);
404 }
405
406
407 static void
drawlingrid(graph,units,spacing,nsp,dst,lmt,hmt,onedec,mult,axis)408 drawlingrid(graph, units, spacing, nsp, dst, lmt, hmt, onedec, mult, axis)
409
410 GRAPH *graph;
411 char units[16];
412 bool onedec;
413 int hmt, lmt, nsp, spacing, dst, mult;
414 Axis axis;
415 {
416
417 int i, j;
418 char buf[16];
419
420 /* i counts how many pixels we have drawn, and j counts which unit
421 * we are at.
422 */
423 for (i = 0, j = lmt; j <= hmt; i += spacing, j += dst / nsp) {
424 if (j == 0)
425 DevSetLinestyle(0);
426 if (graph->grid.gridtype != GRID_NONE) {
427 /* note: get rid of this parameter and draw both axes at once */
428 if (axis == x_axis)
429 DevLine(graph->viewportxoff + i,
430 graph->viewportyoff, graph->viewportxoff + i,
431 graph->viewport.height + graph->viewportyoff);
432 else
433 DevLine(graph->viewportxoff,
434 graph->viewportyoff + i,
435 graph->viewport.width + graph->viewportxoff,
436 graph->viewportyoff + i);
437 }
438 if (j == 0)
439 DevSetLinestyle(1);
440 if (onedec)
441 (void) sprintf(buf, "%.1lf", (double) j / 10);
442 else
443 (void) sprintf(buf, "%d", j * mult);
444 if (axis == x_axis)
445 DevText(buf, graph->viewportxoff + i -
446 strlen(buf) / 2 * graph->fontwidth,
447 (int) (graph->fontheight * 2.5));
448 else
449 DevText(buf, graph->viewportxoff -
450 graph->fontwidth * (strlen(buf) + 1),
451 graph->viewportyoff + i -
452 graph->fontheight / 2);
453 /* This is to make sure things work when delta > hi - lo. */
454 if (nsp == 1)
455 j++;
456 }
457 if (axis == x_axis)
458 DevText(units,
459 graph->absolute.width -
460 (strlen(units) + 15) * graph->fontwidth,
461 graph->fontheight);
462 else
463 DevText(units, graph->fontwidth,
464 (int) (graph->absolute.height * 0.75));
465 DevUpdate();
466 }
467
468
469 /* ARGSUSED */
470 static double *
loggrid(graph,lo,hi,type,axis)471 loggrid(graph, lo, hi, type, axis)
472
473 /* Plot a log grid. Note that we pay no attention to x- and y-delta here. */
474 GRAPH *graph;
475 double lo, hi;
476 int type;
477 Axis axis;
478 {
479 static double dd[2];
480 int margin;
481 int max;
482 int decs, subs, pp, decsp, lmt, hmt;
483
484 if (axis == x_axis) {
485 margin = graph->viewportxoff;
486 max = graph->viewport.width + graph->viewportxoff;
487 }
488 else {
489 margin = graph->viewportyoff;
490 max = graph->viewport.height + graph->viewportyoff;
491 }
492
493 /* How many orders of magnitude. We are already guaranteed that hi
494 * and lo are positive. We want to have something like 8 grid lines
495 * total, so if there are few orders of magnitude put in intermediate
496 * lines.
497 */
498 lmt = floor(mylog10(lo));
499 hmt = ceil(mylog10(hi));
500
501 dd[0] = pow(10.0, (double) lmt);
502 dd[1] = pow(10.0, (double) hmt);
503
504 decs = hmt - lmt;
505 subs = 8 / decs;
506 if (decs > 10)
507 pp = decs / 8 + 1;
508 else
509 pp = 1;
510
511 decsp = (max - margin) * pp / decs;
512
513 if (axis == x_axis) {
514 graph->grid.xaxis.log.hmt = hmt;
515 graph->grid.xaxis.log.lmt = lmt;
516 graph->grid.xaxis.log.decsp = decsp;
517 graph->grid.xaxis.log.subs = subs;
518 graph->grid.xaxis.log.pp = pp;
519 }
520 else {
521 graph->grid.yaxis.log.hmt = hmt;
522 graph->grid.yaxis.log.lmt = lmt;
523 graph->grid.yaxis.log.decsp = decsp;
524 graph->grid.yaxis.log.subs = subs;
525 graph->grid.yaxis.log.pp = pp;
526 }
527
528 return (dd);
529 }
530
531
532 static void
drawloggrid(graph,hmt,lmt,decsp,subs,pp,axis)533 drawloggrid(graph, hmt, lmt, decsp, subs, pp, axis)
534 GRAPH *graph;
535 int hmt, lmt, decsp, subs, pp;
536 Axis axis;
537 {
538
539 int i, j, k, l, m;
540 char buf[16];
541
542 /* Now plot every pp'th decade line, with subs lines between them. */
543 for (i = 0, j = lmt; j <= hmt; i += decsp, j += pp) {
544 /* Draw the decade line */
545 if (graph->grid.gridtype != GRID_NONE) {
546 if (axis == x_axis)
547 DevLine(graph->viewportxoff + i,
548 graph->viewportyoff,
549 graph->viewportxoff + i,
550 graph->viewport.height
551 +graph->viewportyoff);
552 else
553 DevLine(graph->viewportxoff,
554 graph->viewportyoff + i,
555 graph->viewport.width
556 + graph->viewportxoff,
557 graph->viewportyoff + i);
558 }
559 (void) sprintf(buf, "e%d", j);
560 if (axis == x_axis)
561 DevText(buf, graph->viewportxoff + i - strlen(buf) / 2,
562 (int) (graph->fontheight * 2.5));
563 else
564 DevText(buf, graph->viewportxoff - graph->fontwidth *
565 (strlen(buf) + 1),
566 graph->viewportyoff + i -
567 graph->fontheight / 2);
568
569 if (j == hmt)
570 break;
571
572 /* Now draw the other lines */
573 for (k = 1; k <= subs; k++) {
574 l = ceil((double) k * 10 / subs);
575 if (l == 10)
576 break;
577 m = decsp * log10((double ) l) + i;
578 if (graph->grid.gridtype != GRID_NONE) {
579 if (axis == x_axis)
580 DevLine(graph->viewportxoff + m,
581 graph->viewportyoff,
582 graph->viewportxoff + m,
583 graph->viewport.height
584 + graph->viewportyoff);
585 else
586 DevLine(graph->viewportxoff,
587 graph->viewportyoff + m,
588 graph->viewport.width
589 + graph->viewportxoff,
590 graph->viewportyoff + m);
591 }
592 }
593 }
594 DevUpdate();
595 }
596
597
598 /* Polar grids */
599
600 #define minrad (graph->grid.xaxis.circular.mrad)
601 #define maxrad (graph->grid.yaxis.circular.mrad)
602
603
604 static void
polargrid(graph)605 polargrid(graph)
606
607 GRAPH *graph;
608 {
609 double rad, tenpowmag, theta;
610 int hmt, lmt, mag, step;
611 int x1, y1, x2, y2;
612 double xx, yy;
613 bool centered = true;
614
615 /* Make sure that our area is square. */
616 if (graph->viewport.width > graph->viewport.height) {
617 graph->viewport.width = graph->viewport.height;
618 }
619 else {
620 graph->viewport.height = graph->viewport.width;
621 }
622 /* Make sure that the borders are even */
623 if (graph->viewport.width & 1) {
624 graph->viewport.width += 1;
625 graph->viewport.height += 1;
626 }
627
628 graph->grid.xaxis.circular.center = (graph->viewport.width/2.0
629 + graph->viewportxoff);
630 graph->grid.yaxis.circular.center = (graph->viewport.height/2.0
631 + graph->viewportyoff);
632 graph->grid.xaxis.circular.radius = graph->viewport.width / 2;
633
634 /* Figure out the minimum and maximum radii we're dealing with. */
635 rad = sqrt(graph->datawindow.xmin * graph->datawindow.xmin
636 + graph->datawindow.ymin * graph->datawindow.ymin);
637 minrad = rad;
638 maxrad = rad;
639
640 rad = sqrt(graph->datawindow.xmin * graph->datawindow.xmin
641 + graph->datawindow.ymax * graph->datawindow.ymax);
642 if (rad > maxrad)
643 maxrad = rad;
644 if (rad < minrad)
645 minrad = rad;
646 rad = sqrt(graph->datawindow.xmax * graph->datawindow.xmax
647 + graph->datawindow.ymin * graph->datawindow.ymin);
648 if (rad > maxrad)
649 maxrad = rad;
650 if (rad < minrad)
651 minrad = rad;
652 rad = sqrt(graph->datawindow.xmax * graph->datawindow.xmax
653 + graph->datawindow.ymax * graph->datawindow.ymax);
654 if (rad > maxrad)
655 maxrad = rad;
656 if (rad < minrad)
657 minrad = rad;
658
659 if (maxrad == 0.0) {
660 fprintf(cp_err, "Error: 0 radius in polargrid\n");
661 return;
662 }
663 if ((graph->datawindow.xmin < 0) && (graph->datawindow.ymin < 0) &&
664 (graph->datawindow.xmax > 0) && (graph->datawindow.ymax > 0))
665 minrad = 0;
666 if ((graph->datawindow.xmin == - graph->datawindow.xmax)
667 && (graph->datawindow.ymin == -graph->datawindow.ymax)
668 && (graph->datawindow.xmin == graph->datawindow.ymin))
669 centered = true;
670
671 mag = floor(mylog10(maxrad));
672 tenpowmag = pow(10.0, (double) mag);
673 hmt = maxrad / tenpowmag;
674 lmt = minrad / tenpowmag;
675 if (hmt * tenpowmag < maxrad)
676 hmt++;
677 if (lmt * tenpowmag > minrad)
678 lmt--;
679 maxrad = hmt * tenpowmag;
680 minrad = lmt * tenpowmag;
681
682 xx = graph->datawindow.xmin + graph->datawindow.xmax;
683 yy = graph->datawindow.ymin + graph->datawindow.ymax;
684 graph->datawindow.xmin = xx - maxrad;
685 graph->datawindow.xmax = xx + maxrad;
686 graph->datawindow.ymin = yy - maxrad;
687 graph->datawindow.ymax = yy + maxrad;
688
689
690 if (ft_grdb)
691 printf("polar: maxrad = %g, center = (%g, %g)\n",
692 maxrad, xx, yy);
693
694 if ((minrad == 0) && ((hmt - lmt) > 5)) {
695 if (!((hmt - lmt) % 2))
696 step = 2;
697 else if (!((hmt - lmt) % 3))
698 step = 3;
699 else
700 step = 1;
701 }
702 else
703 step = 1;
704
705 graph->grid.xaxis.circular.lmt = lmt;
706 graph->grid.yaxis.circular.lmt = step;
707 }
708
709
710 static void
drawpolargrid(graph)711 drawpolargrid(graph)
712
713 GRAPH *graph;
714 {
715 double mx, my, rad, tenpowmag, theta;
716 int i, step, mag;
717 int relcx, relcy, relrad, dist, degs;
718 int x1, y1, x2, y2;
719 double pixperunit, xx, yy;
720 char buf[64];
721
722 /* The distance from the center of the plotting area to the center of
723 * the logical area.
724 */
725 step = graph->grid.yaxis.circular.lmt;
726 mag = floor(mylog10(maxrad));
727 tenpowmag = pow(10.0, (double) mag);
728 pixperunit = graph->grid.xaxis.circular.radius / (maxrad - minrad);
729
730 relcx = - (graph->datawindow.xmin + graph->datawindow.xmax) / 2
731 * pixperunit;
732 relcy = - (graph->datawindow.ymin + graph->datawindow.ymax) / 2
733 * pixperunit;
734 dist = sqrt((double) (relcx * relcx + relcy * relcy));
735
736 DevSetLinestyle(0);
737 DevArc(graph->grid.xaxis.circular.center,
738 graph->grid.yaxis.circular.center,
739 graph->grid.xaxis.circular.radius,
740 0.0, 0.0);
741 DevSetLinestyle(1);
742
743 /* Now draw the circles. */
744 for (i = graph->grid.xaxis.circular.lmt; (relrad = i * tenpowmag * pixperunit) <=
745 dist + graph->grid.xaxis.circular.radius; i += step) {
746 cliparc(graph->grid.xaxis.circular.center + relcx,
747 graph->grid.yaxis.circular.center + relcy,
748 relrad, 0.0, 0.0,
749 graph->grid.xaxis.circular.center,
750 graph->grid.yaxis.circular.center,
751 graph->grid.xaxis.circular.radius);
752 /* Toss on the label */
753 if (relcx || relcy)
754 theta = atan2((double) relcy, (double) relcx);
755 else
756 theta = M_PI;
757 if (i && (relrad > dist - graph->grid.xaxis.circular.radius))
758 addradlabel(graph, i, theta,
759 (int) (graph->grid.xaxis.circular.center -
760 (relrad - dist) * cos(theta)),
761 (int) (graph->grid.yaxis.circular.center
762 - (relrad - dist) * sin(theta)));
763 }
764
765 /* Now draw the spokes. We have two possible cases -- first, the
766 * origin may be inside the area -- in this case draw 12 spokes.
767 * Otherwise, draw several spokes at convenient places.
768 */
769 if ((graph->datawindow.xmin <= 0.0)
770 && (graph->datawindow.xmax >= 0.0)
771 && (graph->datawindow.ymin <= 0.0)
772 && (graph->datawindow.ymax >= 0.0)) {
773 for (i = 0; i < 12; i++) {
774 x1 = graph->grid.xaxis.circular.center + relcx;
775 y1 = graph->grid.yaxis.circular.center + relcy;
776 x2 = x1 + graph->grid.xaxis.circular.radius * 2
777 * cos(i * M_PI / 6);
778 y2 = y1 + graph->grid.xaxis.circular.radius * 2
779 * sin(i * M_PI / 6);
780 if (!clip_to_circle(&x1, &y1, &x2, &y2,
781 graph->grid.xaxis.circular.center,
782 graph->grid.yaxis.circular.center,
783 graph->grid.xaxis.circular.radius)) {
784 DevLine(x1, y1, x2, y2);
785 /* Add a label here */
786 adddeglabel(graph, i * 30, x2, y2, x1, y1);
787 }
788 }
789 }
790 else {
791 /* Figure out the angle that we have to fill up */
792 theta = 2 * asin((double) graph->grid.xaxis.circular.radius
793 / dist);
794 theta = theta * 180 / M_PI; /* Convert to degrees. */
795
796 /* See if we should put lines at 30, 15, 5, or 1 degree
797 * increments.
798 */
799 if (theta / 30 > 3)
800 degs = 30;
801 else if (theta / 15 > 3)
802 degs = 15;
803 else if (theta / 5 > 3)
804 degs = 5;
805 else
806 degs = 1;
807
808 /* We'll be cheap */
809 for (i = 0; i < 360; i+= degs) {
810 x1 = graph->grid.xaxis.circular.center + relcx;
811 y1 = graph->grid.yaxis.circular.center + relcy;
812 x2 = x1 + dist * 2 * cos(i * M_PI / 180);
813 y2 = y1 + dist * 2 * sin(i * M_PI / 180);
814 if (!clip_to_circle(&x1, &y1, &x2, &y2,
815 graph->grid.xaxis.circular.center,
816 graph->grid.yaxis.circular.center,
817 graph->grid.xaxis.circular.radius)) {
818 DevLine(x1, y1, x2, y2);
819 /* Put on the label */
820 adddeglabel(graph, i, x2, y2, x1, y1);
821 }
822 }
823 }
824
825 (void) sprintf(buf, "e%d", mag);
826 DevText(buf, graph->grid.xaxis.circular.center
827 + graph->grid.xaxis.circular.radius,
828 graph->grid.yaxis.circular.center
829 - graph->grid.xaxis.circular.radius);
830 DevUpdate();
831 return;
832 }
833
834
835 /* Put a degree label on the screen, with 'deg' as the label, near point (x, y)
836 * such that the perpendicular to (lx, ly) and (x, y) doesn't overwrite the
837 * label. Sort of... If the distance between the center and the point is
838 * too small, don't put the label on.
839 */
840
841 #define LOFF 5
842 #define MINDIST 10
843
844
845 static void
adddeglabel(graph,deg,x,y,lx,ly)846 adddeglabel(graph, deg, x, y, lx, ly)
847
848 GRAPH *graph;
849 int deg, x, y, lx, ly;
850 {
851 char buf[8];
852 int d, w, h;
853 double angle;
854
855 if (sqrt((double) (x - lx) * (x - lx) + (y - ly) * (y - ly)) < MINDIST)
856 return;
857 (void) sprintf(buf, "%d", deg);
858 w = graph->fontwidth * (strlen(buf) + 1);
859 h = graph->fontheight * 1.5;
860 angle = atan2((double) (y - ly), (double) (x - lx));
861 d = fabs(cos(angle)) * w / 2 + fabs(sin(angle)) * h / 2 + LOFF;
862
863 x = x + d * cos(angle) - w / 2;
864 y = y + d * sin(angle) - h / 2;
865
866 DevText(buf, x, y);
867 DevText("o", x + strlen(buf) * graph->fontwidth,
868 y + graph->fontheight / 2);
869 return;
870 }
871
872
873 /* This is kind of wierd. If dist = 0, then this is the normal case, where
874 * the labels should go along the positive X-axis. Otherwise, to make
875 * sure that all circles drawn have labels, put the label near the circle
876 * along the line from the logical center to the physical center.
877 */
878
879
880 static void
addradlabel(graph,lab,theta,x,y)881 addradlabel(graph, lab, theta, x, y)
882
883 GRAPH *graph;
884 double theta;
885 int lab, x, y;
886 {
887 char buf[32];
888
889 (void) sprintf(buf, "%d", lab);
890 if (theta == M_PI) {
891 y = y - graph->fontheight - 2;
892 x = x - graph->fontwidth * strlen(buf) - 3;
893 }
894 DevText(buf, x, y);
895 return;
896 }
897
898
899 /* Smith charts. */
900
901 static void
smithgrid(graph)902 smithgrid(graph)
903
904 GRAPH *graph;
905 {
906 double mx, my, d;
907 int mt, i, j, k;
908
909 /* Make sure that our area is square. */
910 if (graph->viewport.width > graph->viewport.height) {
911 graph->viewport.width = graph->viewport.height;
912 }
913 else {
914 graph->viewport.height = graph->viewport.width;
915 }
916
917 graph->grid.xaxis.circular.center = (graph->viewport.width/2.0
918 + graph->viewportxoff);
919 graph->grid.yaxis.circular.center = (graph->viewport.height/2.0
920 + graph->viewportyoff);
921 graph->grid.xaxis.circular.radius = graph->viewport.width / 2;
922
923
924 /* We don't handle cases where the range isn't centered about the
925 * X-axis. This is something that has to be fixed.
926 */
927 if (fabs(graph->datawindow.ymin) > fabs(graph->datawindow.ymax))
928 graph->datawindow.ymax = - graph->datawindow.ymin;
929 else
930 graph->datawindow.ymin = - graph->datawindow.ymax;
931
932
933 /* We have to make sure that the range is square. */
934 mx = graph->datawindow.xmax - graph->datawindow.xmin;
935 my = graph->datawindow.ymax - graph->datawindow.ymin;
936 if (mx > my) {
937 graph->datawindow.ymin -= (mx - my) / 2;
938 graph->datawindow.ymax += (mx - my) / 2;
939 }
940 else if (mx < my) {
941 graph->datawindow.xmin -= (my - mx) / 2;
942 graph->datawindow.xmax += (my - mx) / 2;
943 }
944
945 /* Figure out the minimum and maximum radii we're dealing with. */
946 mx = (graph->datawindow.xmin + graph->datawindow.xmax) / 2;
947 my = (graph->datawindow.ymin + graph->datawindow.ymax) / 2;
948 d = sqrt(mx * mx + my * my);
949 maxrad = d + (graph->datawindow.xmax - graph->datawindow.xmin) / 2;
950 minrad = d - (graph->datawindow.xmax - graph->datawindow.xmin) / 2;
951 }
952
953
954 static void
drawsmithgrid(graph)955 drawsmithgrid( graph )
956
957 GRAPH *graph;
958 {
959 double rad, tenpowmag, d, x;
960 double pixperunit;
961 int relcx, relcy;
962 int mt, mag, i, j, k, zheight;
963 int basemag;
964 char buf[64], plab[32], nlab[32];
965 bool centered = false;
966 int x1, y1, x2, y2;
967
968 DevSetColor(1);
969 DevSetLinestyle(0);
970
971 if ((graph->datawindow.xmin == - graph->datawindow.xmax) &&
972 (graph->datawindow.ymin == -
973 graph->datawindow.ymax) && (graph->datawindow.xmin ==
974 graph->datawindow.ymin))
975 centered = true;
976
977 mag = floor(mylog10(maxrad));
978 tenpowmag = pow(10.0, (double) mag);
979
980 pixperunit = 2 * graph->grid.xaxis.circular.radius / (maxrad - minrad);
981
982 relcx = - (graph->datawindow.xmin + graph->datawindow.xmax) / 2
983 * pixperunit;
984 relcy = - (graph->datawindow.ymin + graph->datawindow.ymax) / 2
985 * pixperunit;
986
987 /* Sweep the range from 10e-20 to 10e20. If any arcs fall into the
988 * picture, plot the arc set.
989 */
990 for (mag = -20; mag < 20; mag++) {
991
992 i = graph->grid.xaxis.circular.radius
993 * pow(10.0, (double) mag) / maxrad;
994 if (i > 10) {
995 j = 1;
996 break;
997 }
998 else if (i > 5) {
999 j = 2;
1000 break;
1001 }
1002 else if (i > 2) {
1003 j = 5;
1004 break;
1005 }
1006 }
1007 k = 0;
1008
1009 /* Now plot all the arc sets. Go as high as 5 times the radius that
1010 * will fit on the screen. The base magnitude is one more than
1011 * the least magnitude that will fit.
1012 */
1013 if (i > 20)
1014 basemag = mag;
1015 else
1016 basemag = mag + 1;
1017 mag = basemag;
1018 while (mag < 20) {
1019 /*
1020 i = j * gr_radius * pow(10.0, (double) mag) / (maxrad * 2);
1021 */
1022 i = j * pow(10.0, (double) mag) * pixperunit / 2;
1023 if (i / 5 > graph->grid.xaxis.circular.radius + pixperunit
1024 + ((relcx > 0) ? relcx : - relcx))
1025 break;
1026 x = j * pow(10.0, (double) (mag));
1027 if ( x == 0.0 ) {
1028 plab[0] = 0;
1029 nlab[0] = 0;
1030 }
1031 else {
1032 (void) sprintf(plab, "%lg", (2-x)/x);
1033 (void) sprintf(nlab, "-%lg",(1+x)/x);
1034 }
1035 arcset(i, k, graph,
1036 (int)(relcx + pixperunit), relcy,plab,nlab);
1037 if (i * 2.5 < graph->grid.xaxis.circular.radius
1038 + ((relcx > 0) ? relcx : - relcx))
1039 k = i;
1040 if (j == 5) {
1041 j = 1;
1042 mag++;
1043 }
1044 else if (j == 2)
1045 j = 5;
1046 else if (j == 1)
1047 j = 2;
1048 }
1049 if (mag == 20) {
1050 fprintf(cp_err, "smithgrid: Internal Error: screwed up\n");
1051 return;
1052 }
1053
1054 DevArc( graph->grid.xaxis.circular.center,
1055 graph->grid.yaxis.circular.center,
1056 graph->grid.xaxis.circular.radius, 0.0, 0.0);
1057 if ((relcx + pixperunit > - graph->grid.xaxis.circular.radius )
1058 && (relcx + pixperunit < graph->grid.xaxis.circular.radius )) {
1059 zheight = graph->grid.xaxis.circular.radius
1060 * sin(acos((double) (relcx + pixperunit)
1061 / graph->grid.xaxis.circular.radius ));
1062 if (zheight < 0)
1063 zheight = - zheight;
1064 x1 = graph->grid.xaxis.circular.center + relcx + pixperunit;
1065 y1 = graph->grid.yaxis.circular.center - zheight;
1066 x2 = graph->grid.xaxis.circular.center + relcx + pixperunit;
1067 y2 = graph->grid.yaxis.circular.center + zheight;
1068
1069 if (!clip_to_circle(&x1, &y1, &x2, &y2,
1070 graph->grid.xaxis.circular.center,
1071 graph->grid.yaxis.circular.center,
1072 graph->grid.xaxis.circular.radius))
1073 DevLine(x1, y1, x2, y2);
1074
1075 }
1076 if ((relcy > - graph->grid.xaxis.circular.radius )
1077 && (relcy < graph->grid.xaxis.circular.radius )) {
1078 zheight = graph->grid.xaxis.circular.radius
1079 * cos(asin((double) relcy / graph->grid.xaxis.circular.radius ));
1080 if (zheight < 0)
1081 zheight = - zheight;
1082 DevLine( graph->grid.xaxis.circular.center - zheight,
1083 graph->grid.yaxis.circular.center + relcy,
1084 graph->grid.xaxis.circular.center + zheight,
1085 graph->grid.yaxis.circular.center + relcy);
1086 }
1087
1088 /*
1089 (void) sprintf(buf, "e%d", basemag);
1090 DevText(buf,
1091 graph->grid.xaxis.circular.center + graph->grid.xaxis.circular.radius,
1092 graph->grid.yaxis.circular.center - graph->grid.xaxis.circular.radius
1093 );
1094
1095 */
1096 #ifdef notdef
1097 gi_text("0", gr_xcenter + gr_radius + gi_fntwidth, gr_ycenter -
1098 gi_fntheight / 2, 0, false);
1099 gi_text("o", gr_xcenter + gr_radius + gi_fntwidth * 2, gr_ycenter, 0,
1100 false);
1101 gi_text("90", gr_xcenter - gi_fntwidth, gr_ycenter + gr_radius +
1102 gi_fntheight / 2, 0, false);
1103 gi_text("o", gr_xcenter + gi_fntwidth, gr_ycenter + gr_radius +
1104 gi_fntheight, 0, false);
1105 gi_text("180", gr_xcenter - gr_radius - gi_fntwidth * 5, gr_ycenter
1106 - gi_fntheight / 2, 0, false);
1107 gi_text("o", gr_xcenter - gr_radius - gi_fntwidth * 2, gr_ycenter, 0,
1108 false);
1109 gi_text("-90", gr_xcenter - gi_fntwidth * 2, gr_ycenter - gr_radius -
1110 2 * gi_fntheight, 0, false);
1111 gi_text("o", gr_xcenter + gi_fntwidth, gr_ycenter - gr_radius -
1112 gi_fntheight - gi_fntheight / 2, 0, false);
1113 #endif
1114 DevUpdate();
1115 return;
1116 }
1117
1118
1119 /* Draw one arc set. The arcs should have radius rad. The outermost circle
1120 * is described by (centx, centy) and maxrad, and the distance from the
1121 * right side of the bounding circle to the logical center of the other
1122 * circles in pixels is xoffset (positive brings the negative plane into
1123 * the picture).
1124 * plab and nlab are the labels to put on the positive and negative X-arcs,
1125 * respectively. If the X-axis isn't on the screen, then we have to be
1126 * clever.
1127 */
1128
1129 static void
arcset(rad,prevrad,graph,xoffset,yoffset,plab,nlab)1130 arcset(rad, prevrad, graph, xoffset, yoffset, plab, nlab)
1131
1132 int rad, prevrad, xoffset, yoffset;
1133 GRAPH *graph;
1134 char *plab, *nlab;
1135 {
1136 double angle = atan(((double) prevrad) / rad);
1137 int x;
1138
1139 /* Let's be lazy and just draw everything -- we won't get called too
1140 * much and the circles get clipped anyway.
1141 */
1142 DevSetColor(5);
1143 cliparc(graph->grid.xaxis.circular.center + xoffset - rad,
1144 graph->grid.yaxis.circular.center + yoffset, rad,
1145 2 * angle, 2 * M_PI - 2 * angle,
1146 graph->grid.xaxis.circular.center,
1147 graph->grid.yaxis.circular.center,
1148 graph->grid.xaxis.circular.radius );
1149
1150 cliparc(graph->grid.xaxis.circular.center + xoffset + rad,
1151 graph->grid.yaxis.circular.center + yoffset, rad,
1152 M_PI + 2 * angle, M_PI - 2 * angle,
1153 graph->grid.xaxis.circular.center,
1154 graph->grid.yaxis.circular.center,
1155 graph->grid.xaxis.circular.radius );
1156
1157
1158 /* Draw the upper and lower circles. */
1159
1160 DevSetColor(3);
1161
1162 cliparc(graph->grid.xaxis.circular.center + xoffset,
1163 graph->grid.yaxis.circular.center + yoffset + rad, rad,
1164 M_PI * 1.5 + 2 * angle, M_PI * 1.5 - 2 * angle,
1165 graph->grid.xaxis.circular.center,
1166 graph->grid.yaxis.circular.center,
1167 graph->grid.xaxis.circular.radius );
1168
1169 cliparc(graph->grid.xaxis.circular.center + xoffset,
1170 graph->grid.yaxis.circular.center + yoffset - rad, rad,
1171 M_PI / 2 + 2 * angle, M_PI / 2 - 2 * angle,
1172 graph->grid.xaxis.circular.center,
1173 graph->grid.yaxis.circular.center,
1174 graph->grid.xaxis.circular.radius );
1175
1176
1177 /* Now toss the labels on. */
1178
1179 DevSetColor(1);
1180
1181 x = graph->grid.xaxis.circular.center + xoffset
1182 - 2 * rad - graph->fontwidth * strlen(plab) - 2;
1183 if ((x > graph->viewportxoff ) && (x < graph->viewportxoff +
1184 graph->viewport.width))
1185 DevText(plab, x,
1186 graph->grid.yaxis.circular.center - graph->fontheight - 1);
1187
1188 x = graph->grid.xaxis.circular.center + xoffset
1189 + 2 * rad - graph->fontwidth * strlen(nlab) - 2;
1190 if ((x > graph->viewportxoff )
1191 && (x < graph->viewportxoff + graph->viewport.width)
1192 && ( 2 * rad + xoffset < graph->grid.xaxis.circular.radius))
1193 DevText(nlab, x,
1194 graph->grid.yaxis.circular.center - graph->fontheight - 1);
1195
1196 return;
1197 }
1198
1199
1200 /* This routine draws an arc and clips it to a circle. It's hard to figure
1201 * out how it works without looking at the piece of scratch paper I have
1202 * in front of me, so let's hope it doesn't break.
1203 */
1204
1205 static void
cliparc(cx,cy,rad,start,end,clipx,clipy,cliprad)1206 cliparc(cx, cy, rad, start, end, clipx, clipy, cliprad)
1207
1208 int cx, cy, rad, clipx, clipy, cliprad;
1209 double start, end;
1210 {
1211 int x, y, tx, ty, dist;
1212 double alpha, theta, phi, a1, a2, d, l;
1213 bool in;
1214
1215 x = cx - clipx;
1216 y = cy - clipy;
1217 dist = sqrt((double) (x * x + y * y));
1218
1219 if (!rad || !cliprad)
1220 return;
1221 if (dist + rad < cliprad) {
1222 /* The arc is entirely in the boundary. */
1223 DevArc(cx, cy, rad, start, end);
1224 return;
1225 }
1226 else if ((dist - rad >= cliprad) || (rad - dist >= cliprad)) {
1227 /* The arc is outside of the boundary. */
1228 return;
1229 }
1230 /* Now let's figure out the angles at which the arc crosses the
1231 * circle. We know dist != 0.
1232 */
1233 if (x)
1234 phi = atan((double) y / x);
1235 else if (y > 0)
1236 phi = M_PI * 1.5;
1237 else
1238 phi = M_PI * .5;
1239 if (cx > clipx)
1240 theta = M_PI + phi;
1241 else
1242 theta = phi;
1243
1244 alpha = (double) (dist * dist + rad * rad - cliprad * cliprad) /
1245 (2 * dist * rad);
1246 if ( alpha < -1.0 )
1247 alpha = -1.0;
1248 else if (alpha > 1.0)
1249 alpha = 1.0;
1250 alpha = acos(alpha);
1251
1252 a1 = theta + alpha;
1253 a2 = theta - alpha;
1254 while (a1 < 0)
1255 a1 += M_PI * 2;
1256 while (a2 < 0)
1257 a2 += M_PI * 2;
1258 while (a1 >= M_PI * 2)
1259 a1 -= M_PI * 2;
1260 while (a2 >= M_PI * 2)
1261 a2 -= M_PI * 2;
1262
1263 tx = cos(start) * rad + x;
1264 ty = sin(start) * rad + y;
1265 d = sqrt((double) tx * tx + ty * ty);
1266 in = (d > cliprad) ? false : true;
1267
1268 /* Now begin with start. If the point is in, draw to either end, a1,
1269 * or a2, whichever comes first.
1270 */
1271 d = M_PI * 3;
1272 if ((end < d) && (end > start))
1273 d = end;
1274 if ((a1 < d) && (a1 > start))
1275 d = a1;
1276 if ((a2 < d) && (a2 > start))
1277 d = a2;
1278 if (d == M_PI * 3) {
1279 d = end;
1280 if (a1 < d)
1281 d = a1;
1282 if (a2 < d)
1283 d = a2;
1284 }
1285
1286 if (in)
1287 DevArc(cx, cy, rad, start, d);
1288 if (d == end)
1289 return;
1290 if (a1 != a2)
1291 in = in ? false : true;
1292
1293 /* Now go from here to the next point. */
1294 l = d;
1295 d = M_PI * 3;
1296 if ((end < d) && (end > l))
1297 d = end;
1298 if ((a1 < d) && (a1 > l))
1299 d = a1;
1300 if ((a2 < d) && (a2 > l))
1301 d = a2;
1302 if (d == M_PI * 3) {
1303 d = end;
1304 if (a1 < d)
1305 d = a1;
1306 if (a2 < d)
1307 d = a2;
1308 }
1309
1310 if (in)
1311 DevArc(cx, cy, rad, l, d);
1312 if (d == end)
1313 return;
1314 in = in ? false : true;
1315
1316 /* And from here to the end. */
1317 if (in)
1318 DevArc(cx, cy, rad, d, end);
1319 return;
1320 }
1321
1322
1323 #ifdef notdef
1324
1325 /* These routines are wrong on a lot of systems. */
1326
1327 double
floor(d)1328 floor(d)
1329 double d;
1330 {
1331 double fp, ip;
1332
1333 if (d >= 0.0) {
1334 fp = modf(d, &ip);
1335 return (ip);
1336 }
1337 else {
1338 /* This is the case that lattice C screws up */
1339 fp = modf(-d, &ip);
1340 if (fp != 0.0)
1341 return (-ip - 1);
1342 else
1343 return (-ip);
1344 }
1345 }
1346
1347 double
ceil(d)1348 ceil(d)
1349 double d;
1350 {
1351 return (- floor(-d));
1352 }
1353
1354 #endif
1355