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