1/* Last non-groff version: hgraph.c  1.14 (Berkeley) 84/11/27
2 *
3 * This file contains the graphics routines for converting gremlin pictures
4 * to troff input.
5 */
6
7#include "lib.h"
8
9#include "gprint.h"
10
11#define MAXVECT	40
12#define MAXPOINTS	200
13#define LINELENGTH	1
14#define PointsPerInterval 64
15#define pi		3.14159265358979324
16#define twopi		(2.0 * pi)
17#define len(a, b)	groff_hypot((double)(b.x-a.x), (double)(b.y-a.y))
18
19
20extern int dotshifter;		/* for the length of dotted curves */
21
22extern int style[];		/* line and character styles */
23extern double thick[];
24extern char *tfont[];
25extern int tsize[];
26extern int stipple_index[];	/* stipple font index for stipples 0 - 16 */
27extern char *stipple;		/* stipple type (cf or ug) */
28
29
30extern double troffscale;	/* imports from main.c */
31extern double linethickness;
32extern int linmod;
33extern int lastx;
34extern int lasty;
35extern int lastyline;
36extern int ytop;
37extern int ybottom;
38extern int xleft;
39extern int xright;
40extern enum E {
41  OUTLINE, FILL, BOTH
42} polyfill;
43
44extern double adj1;
45extern double adj2;
46extern double adj3;
47extern double adj4;
48extern int res;
49
50void HGSetFont(int font, int size);
51void HGPutText(int justify, POINT pnt, register char *string);
52void HGSetBrush(int mode);
53void tmove2(int px, int py);
54void doarc(POINT cp, POINT sp, int angle);
55void tmove(POINT * ptr);
56void cr();
57void drawwig(POINT * ptr, int type);
58void HGtline(int x1, int y1);
59void deltax(double x);
60void deltay(double y);
61void HGArc(register int cx, register int cy, int px, int py, int angle);
62void picurve(register int *x, register int *y, int npts);
63void HGCurve(int *x, int *y, int numpoints);
64void Paramaterize(int x[], int y[], double h[], int n);
65void PeriodicSpline(double h[], int z[],
66		    double dz[], double d2z[], double d3z[],
67		    int npoints);
68void NaturalEndSpline(double h[], int z[],
69		      double dz[], double d2z[], double d3z[],
70		      int npoints);
71
72
73
74/*----------------------------------------------------------------------------*
75 | Routine:	HGPrintElt (element_pointer, baseline)
76 |
77 | Results:	Examines a picture element and calls the appropriate
78 |		routine(s) to print them according to their type.  After the
79 |		picture is drawn, current position is (lastx, lasty).
80 *----------------------------------------------------------------------------*/
81
82void
83HGPrintElt(ELT *element,
84	   int /* baseline */)
85{
86  register POINT *p1;
87  register POINT *p2;
88  register int length;
89  register int graylevel;
90
91  if (!DBNullelt(element) && !Nullpoint((p1 = element->ptlist))) {
92    /* p1 always has first point */
93    if (TEXT(element->type)) {
94      HGSetFont(element->brushf, element->size);
95      switch (element->size) {
96      case 1:
97	p1->y += adj1;
98	break;
99      case 2:
100	p1->y += adj2;
101	break;
102      case 3:
103	p1->y += adj3;
104	break;
105      case 4:
106	p1->y += adj4;
107	break;
108      default:
109	break;
110      }
111      HGPutText(element->type, *p1, element->textpt);
112    } else {
113      if (element->brushf)		/* if there is a brush, the */
114	HGSetBrush(element->brushf);	/* graphics need it set     */
115
116      switch (element->type) {
117
118      case ARC:
119	p2 = PTNextPoint(p1);
120	tmove(p2);
121	doarc(*p1, *p2, element->size);
122	cr();
123	break;
124
125      case CURVE:
126	length = 0;	/* keep track of line length */
127	drawwig(p1, CURVE);
128	cr();
129	break;
130
131      case BSPLINE:
132	length = 0;	/* keep track of line length */
133	drawwig(p1, BSPLINE);
134	cr();
135	break;
136
137      case VECTOR:
138	length = 0;		/* keep track of line length so */
139	tmove(p1);		/* single lines don't get long  */
140	while (!Nullpoint((p1 = PTNextPoint(p1)))) {
141	  HGtline((int) (p1->x * troffscale),
142		  (int) (p1->y * troffscale));
143	  if (length++ > LINELENGTH) {
144	    length = 0;
145	    printf("\\\n");
146	  }
147	}			/* end while */
148	cr();
149	break;
150
151      case POLYGON:
152	{
153	  /* brushf = style of outline; size = color of fill:
154	   * on first pass (polyfill=FILL), do the interior using 'P'
155	   *    unless size=0
156	   * on second pass (polyfill=OUTLINE), do the outline using a series
157	   *    of vectors. It might make more sense to use \D'p ...',
158	   *    but there is no uniform way to specify a 'fill character'
159	   *    that prints as 'no fill' on all output devices (and
160	   *    stipple fonts).
161	   * If polyfill=BOTH, just use the \D'p ...' command.
162	   */
163	  double firstx = p1->x;
164	  double firsty = p1->y;
165
166	  length = 0;		/* keep track of line length so */
167				/* single lines don't get long  */
168
169	  if (polyfill == FILL || polyfill == BOTH) {
170	    /* do the interior */
171	    char command = (polyfill == BOTH && element->brushf) ? 'p' : 'P';
172
173	    /* include outline, if there is one and */
174	    /* the -p flag was set                  */
175
176	    /* switch based on what gremlin gives */
177	    switch (element->size) {
178	    case 1:
179	      graylevel = 1;
180	      break;
181	    case 3:
182	      graylevel = 2;
183	      break;
184	    case 12:
185	      graylevel = 3;
186	      break;
187	    case 14:
188	      graylevel = 4;
189	      break;
190	    case 16:
191	      graylevel = 5;
192	      break;
193	    case 19:
194	      graylevel = 6;
195	      break;
196	    case 21:
197	      graylevel = 7;
198	      break;
199	    case 23:
200	      graylevel = 8;
201	      break;
202	    default:		/* who's giving something else? */
203	      graylevel = NSTIPPLES;
204	      break;
205	    }
206	    /* int graylevel = element->size; */
207
208	    if (graylevel < 0)
209	      break;
210	    if (graylevel > NSTIPPLES)
211	      graylevel = NSTIPPLES;
212	    printf("\\D'Fg %.3f'",
213		   double(1000 - stipple_index[graylevel]) / 1000.0);
214	    cr();
215	    tmove(p1);
216	    printf("\\D'%c", command);
217
218	    while (!Nullpoint((PTNextPoint(p1)))) {
219	      p1 = PTNextPoint(p1);
220	      deltax((double) p1->x);
221	      deltay((double) p1->y);
222	      if (length++ > LINELENGTH) {
223		length = 0;
224		printf("\\\n");
225	      }
226	    } /* end while */
227
228	    /* close polygon if not done so by user */
229	    if ((firstx != p1->x) || (firsty != p1->y)) {
230	      deltax((double) firstx);
231	      deltay((double) firsty);
232	    }
233	    putchar('\'');
234	    cr();
235	    break;
236	  }
237	  /* else polyfill == OUTLINE; only draw the outline */
238	  if (!(element->brushf))
239	    break;
240	  length = 0;		/* keep track of line length */
241	  tmove(p1);
242
243	  while (!Nullpoint((PTNextPoint(p1)))) {
244	    p1 = PTNextPoint(p1);
245	    HGtline((int) (p1->x * troffscale),
246		    (int) (p1->y * troffscale));
247	    if (length++ > LINELENGTH) {
248	      length = 0;
249	      printf("\\\n");
250	    }
251	  }			/* end while */
252
253	  /* close polygon if not done so by user */
254	  if ((firstx != p1->x) || (firsty != p1->y)) {
255	    HGtline((int) (firstx * troffscale),
256		    (int) (firsty * troffscale));
257	  }
258	  cr();
259	  break;
260	}			/* end case POLYGON */
261      }				/* end switch */
262    }				/* end else Text */
263  }				/* end if */
264}				/* end PrintElt */
265
266
267/*----------------------------------------------------------------------------*
268 | Routine:	HGPutText (justification, position_point, string)
269 |
270 | Results:	Given the justification, a point to position with, and a
271 |		string to put, HGPutText first sends the string into a
272 |		diversion, moves to the positioning point, then outputs
273 |		local vertical and horizontal motions as needed to justify
274 |		the text.  After all motions are done, the diversion is
275 |		printed out.
276 *----------------------------------------------------------------------------*/
277
278void
279HGPutText(int justify,
280	  POINT pnt,
281	  register char *string)
282{
283  int savelasty = lasty;	/* vertical motion for text is to be */
284				/* ignored.  Save current y here     */
285
286  printf(".nr g8 \\n(.d\n");	/* save current vertical position. */
287  printf(".ds g9 \"");		/* define string containing the text. */
288  while (*string) {		/* put out the string */
289    if (*string == '\\' &&
290	*(string + 1) == '\\') {	/* one character at a */
291      printf("\\\\\\");			/* time replacing //  */
292      string++;				/* by //// to prevent */
293    }					/* interpretation at  */
294    printf("%c", *(string++));		/* printout time      */
295  }
296  printf("\n");
297
298  tmove(&pnt);			/* move to positioning point */
299
300  switch (justify) {
301    /* local vertical motions                                            */
302    /* (the numbers here are used to be somewhat compatible with gprint) */
303  case CENTLEFT:
304  case CENTCENT:
305  case CENTRIGHT:
306    printf("\\v'0.85n'");	/* down half */
307    break;
308
309  case TOPLEFT:
310  case TOPCENT:
311  case TOPRIGHT:
312    printf("\\v'1.7n'");	/* down whole */
313  }
314
315  switch (justify) {
316    /* local horizontal motions */
317  case BOTCENT:
318  case CENTCENT:
319  case TOPCENT:
320    printf("\\h'-\\w'\\*(g9'u/2u'");	/* back half */
321    break;
322
323  case BOTRIGHT:
324  case CENTRIGHT:
325  case TOPRIGHT:
326    printf("\\h'-\\w'\\*(g9'u'");	/* back whole */
327  }
328
329  printf("\\&\\*(g9\n");	/* now print the text. */
330  printf(".sp |\\n(g8u\n");	/* restore vertical position */
331  lasty = savelasty;		/* vertical position restored to where it */
332  lastx = xleft;		/* was before text, also horizontal is at */
333				/* left                                   */
334}				/* end HGPutText */
335
336
337/*----------------------------------------------------------------------------*
338 | Routine:	doarc (center_point, start_point, angle)
339 |
340 | Results:	Produces either drawarc command or a drawcircle command
341 |		depending on the angle needed to draw through.
342 *----------------------------------------------------------------------------*/
343
344void
345doarc(POINT cp,
346      POINT sp,
347      int angle)
348{
349  if (angle)			/* arc with angle */
350    HGArc((int) (cp.x * troffscale), (int) (cp.y * troffscale),
351	  (int) (sp.x * troffscale), (int) (sp.y * troffscale), angle);
352  else				/* a full circle (angle == 0) */
353    HGArc((int) (cp.x * troffscale), (int) (cp.y * troffscale),
354	  (int) (sp.x * troffscale), (int) (sp.y * troffscale), 0);
355}
356
357
358/*----------------------------------------------------------------------------*
359 | Routine:	HGSetFont (font_number, Point_size)
360 |
361 | Results:	ALWAYS outputs a .ft and .ps directive to troff.  This is
362 |		done because someone may change stuff inside a text string.
363 |		Changes thickness back to default thickness.  Default
364 |		thickness depends on font and pointsize.
365 *----------------------------------------------------------------------------*/
366
367void
368HGSetFont(int font,
369	  int size)
370{
371  printf(".ft %s\n"
372	 ".ps %d\n", tfont[font - 1], tsize[size - 1]);
373  linethickness = DEFTHICK;
374}
375
376
377/*----------------------------------------------------------------------------*
378 | Routine:	HGSetBrush (line_mode)
379 |
380 | Results:	Generates the troff commands to set up the line width and
381 |		style of subsequent lines.  Does nothing if no change is
382 |              needed.
383 |
384 | Side Efct:	Sets `linmode' and `linethicknes'.
385 *----------------------------------------------------------------------------*/
386
387void
388HGSetBrush(int mode)
389{
390  register int printed = 0;
391
392  if (linmod != style[--mode]) {
393    /* Groff doesn't understand \Ds, so we take it out */
394    /* printf ("\\D's %du'", linmod = style[mode]); */
395    linmod = style[mode];
396    printed = 1;
397  }
398  if (linethickness != thick[mode]) {
399    linethickness = thick[mode];
400    printf("\\h'-%.2fp'\\D't %.2fp'", linethickness, linethickness);
401    printed = 1;
402  }
403  if (printed)
404    cr();
405}
406
407
408/*----------------------------------------------------------------------------*
409 | Routine:	deltax (x_destination)
410 |
411 | Results:	Scales and outputs a number for delta x (with a leading
412 |		space) given `lastx' and x_destination.
413 |
414 | Side Efct:	Resets `lastx' to x_destination.
415 *----------------------------------------------------------------------------*/
416
417void
418deltax(double x)
419{
420  register int ix = (int) (x * troffscale);
421
422  printf(" %du", ix - lastx);
423  lastx = ix;
424}
425
426
427/*----------------------------------------------------------------------------*
428 | Routine:	deltay (y_destination)
429 |
430 | Results:	Scales and outputs a number for delta y (with a leading
431 |		space) given `lastyline' and y_destination.
432 |
433 | Side Efct:	Resets `lastyline' to y_destination.  Since `line' vertical
434 |		motions don't affect `page' ones, `lasty' isn't updated.
435 *----------------------------------------------------------------------------*/
436
437void
438deltay(double y)
439{
440  register int iy = (int) (y * troffscale);
441
442  printf(" %du", iy - lastyline);
443  lastyline = iy;
444}
445
446
447/*----------------------------------------------------------------------------*
448 | Routine:	tmove2 (px, py)
449 |
450 | Results:	Produces horizontal and vertical moves for troff given the
451 |		pair of points to move to and knowing the current position.
452 |		Also puts out a horizontal move to start the line.  This is
453 |		a variation without the .sp command.
454 *----------------------------------------------------------------------------*/
455
456void
457tmove2(int px,
458       int py)
459{
460  register int dx;
461  register int dy;
462
463  if ((dy = py - lasty)) {
464    printf("\\v'%du'", dy);
465  }
466  lastyline = lasty = py;	/* lasty is always set to current */
467  if ((dx = px - lastx)) {
468    printf("\\h'%du'", dx);
469    lastx = px;
470  }
471}
472
473
474/*----------------------------------------------------------------------------*
475 | Routine:	tmove (point_pointer)
476 |
477 | Results:	Produces horizontal and vertical moves for troff given the
478 |		pointer of a point to move to and knowing the current
479 |		position.  Also puts out a horizontal move to start the
480 |		line.
481 *----------------------------------------------------------------------------*/
482
483void
484tmove(POINT * ptr)
485{
486  register int ix = (int) (ptr->x * troffscale);
487  register int iy = (int) (ptr->y * troffscale);
488  register int dx;
489  register int dy;
490
491  if ((dy = iy - lasty)) {
492    printf(".sp %du\n", dy);
493  }
494  lastyline = lasty = iy;	/* lasty is always set to current */
495  if ((dx = ix - lastx)) {
496    printf("\\h'%du'", dx);
497    lastx = ix;
498  }
499}
500
501
502/*----------------------------------------------------------------------------*
503 | Routine:	cr ( )
504 |
505 | Results:	Ends off an input line.  `.sp -1' is also added to counteract
506 |		the vertical move done at the end of text lines.
507 |
508 | Side Efct:	Sets `lastx' to `xleft' for troff's return to left margin.
509 *----------------------------------------------------------------------------*/
510
511void
512cr()
513{
514  printf("\n.sp -1\n");
515  lastx = xleft;
516}
517
518
519/*----------------------------------------------------------------------------*
520 | Routine:	line ( )
521 |
522 | Results:	Draws a single solid line to (x,y).
523 *----------------------------------------------------------------------------*/
524
525void
526line(int px,
527     int py)
528{
529  printf("\\D'l");
530  printf(" %du", px - lastx);
531  printf(" %du'", py - lastyline);
532  lastx = px;
533  lastyline = lasty = py;
534}
535
536
537/*----------------------------------------------------------------------------
538 | Routine:	drawwig (ptr, type)
539 |
540 | Results:	The point sequence found in the structure pointed by ptr is
541 |		placed in integer arrays for further manipulation by the
542 |		existing routing.  With the corresponding type parameter,
543 |		either picurve or HGCurve are called.
544 *----------------------------------------------------------------------------*/
545
546void
547drawwig(POINT * ptr,
548	int type)
549{
550  register int npts;			/* point list index */
551  int x[MAXPOINTS], y[MAXPOINTS];	/* point list */
552
553  for (npts = 1; !Nullpoint(ptr); ptr = PTNextPoint(ptr), npts++) {
554    x[npts] = (int) (ptr->x * troffscale);
555    y[npts] = (int) (ptr->y * troffscale);
556  }
557  if (--npts) {
558    if (type == CURVE) /* Use the 2 different types of curves */
559      HGCurve(&x[0], &y[0], npts);
560    else
561      picurve(&x[0], &y[0], npts);
562  }
563}
564
565
566/*----------------------------------------------------------------------------
567 | Routine:	HGArc (xcenter, ycenter, xstart, ystart, angle)
568 |
569 | Results:	This routine plots an arc centered about (cx, cy) counter
570 |		clockwise starting from the point (px, py) through `angle'
571 |		degrees.  If angle is 0, a full circle is drawn.  It does so
572 |		by creating a draw-path around the arc whose density of
573 |		points depends on the size of the arc.
574 *----------------------------------------------------------------------------*/
575
576void
577HGArc(register int cx,
578      register int cy,
579      int px,
580      int py,
581      int angle)
582{
583  double xs, ys, resolution, fullcircle;
584  int m;
585  register int mask;
586  register int extent;
587  register int nx;
588  register int ny;
589  register int length;
590  register double epsilon;
591
592  xs = px - cx;
593  ys = py - cy;
594
595  length = 0;
596
597  resolution = (1.0 + groff_hypot(xs, ys) / res) * PointsPerInterval;
598  /* mask = (1 << (int) log10(resolution + 1.0)) - 1; */
599  (void) frexp(resolution, &m);		/* A bit more elegant than log10 */
600  for (mask = 1; mask < m; mask = mask << 1);
601  mask -= 1;
602
603  epsilon = 1.0 / resolution;
604  fullcircle = (2.0 * pi) * resolution;
605  if (angle == 0)
606    extent = (int) fullcircle;
607  else
608    extent = (int) (angle * fullcircle / 360.0);
609
610  HGtline(px, py);
611  while (--extent >= 0) {
612    xs += epsilon * ys;
613    nx = cx + (int) (xs + 0.5);
614    ys -= epsilon * xs;
615    ny = cy + (int) (ys + 0.5);
616    if (!(extent & mask)) {
617      HGtline(nx, ny);		/* put out a point on circle */
618      if (length++ > LINELENGTH) {
619	length = 0;
620	printf("\\\n");
621      }
622    }
623  }				/* end for */
624}				/* end HGArc */
625
626
627/*----------------------------------------------------------------------------
628 | Routine:	picurve (xpoints, ypoints, num_of_points)
629 |
630 | Results:	Draws a curve delimited by (not through) the line segments
631 |		traced by (xpoints, ypoints) point list.  This is the `Pic'
632 |		style curve.
633 *----------------------------------------------------------------------------*/
634
635void
636picurve(register int *x,
637	register int *y,
638	int npts)
639{
640  register int nseg;		/* effective resolution for each curve */
641  register int xp;		/* current point (and temporary) */
642  register int yp;
643  int pxp, pyp;			/* previous point (to make lines from) */
644  int i;			/* inner curve segment traverser */
645  int length = 0;
646  double w;			/* position factor */
647  double t1, t2, t3;		/* calculation temps */
648
649  if (x[1] == x[npts] && y[1] == y[npts]) {
650    x[0] = x[npts - 1];		/* if the lines' ends meet, make */
651    y[0] = y[npts - 1];		/* sure the curve meets          */
652    x[npts + 1] = x[2];
653    y[npts + 1] = y[2];
654  } else {			/* otherwise, make the ends of the  */
655    x[0] = x[1];		/* curve touch the ending points of */
656    y[0] = y[1];		/* the line segments                */
657    x[npts + 1] = x[npts];
658    y[npts + 1] = y[npts];
659  }
660
661  pxp = (x[0] + x[1]) / 2;	/* make the last point pointers       */
662  pyp = (y[0] + y[1]) / 2;	/* point to the start of the 1st line */
663  tmove2(pxp, pyp);
664
665  for (; npts--; x++, y++) {	/* traverse the line segments */
666    xp = x[0] - x[1];
667    yp = y[0] - y[1];
668    nseg = (int) groff_hypot((double) xp, (double) yp);
669    xp = x[1] - x[2];
670    yp = y[1] - y[2];
671				/* `nseg' is the number of line    */
672				/* segments that will be drawn for */
673				/* each curve segment.             */
674    nseg = (int) ((double) (nseg + (int) groff_hypot((double) xp, (double) yp)) /
675		  res * PointsPerInterval);
676
677    for (i = 1; i < nseg; i++) {
678      w = (double) i / (double) nseg;
679      t1 = w * w;
680      t3 = t1 + 1.0 - (w + w);
681      t2 = 2.0 - (t3 + t1);
682      xp = (((int) (t1 * x[2] + t2 * x[1] + t3 * x[0])) + 1) / 2;
683      yp = (((int) (t1 * y[2] + t2 * y[1] + t3 * y[0])) + 1) / 2;
684
685      HGtline(xp, yp);
686      if (length++ > LINELENGTH) {
687	length = 0;
688	printf("\\\n");
689      }
690    }
691  }
692}
693
694
695/*----------------------------------------------------------------------------
696 | Routine:	HGCurve(xpoints, ypoints, num_points)
697 |
698 | Results:	This routine generates a smooth curve through a set of
699 |		points.  The method used is the parametric spline curve on
700 |		unit knot mesh described in `Spline Curve Techniques' by
701 |		Patrick Baudelaire, Robert Flegal, and Robert Sproull --
702 |		Xerox Parc.
703 *----------------------------------------------------------------------------*/
704
705void
706HGCurve(int *x,
707	int *y,
708	int numpoints)
709{
710  double h[MAXPOINTS], dx[MAXPOINTS], dy[MAXPOINTS];
711  double d2x[MAXPOINTS], d2y[MAXPOINTS], d3x[MAXPOINTS], d3y[MAXPOINTS];
712  double t, t2, t3;
713  register int j;
714  register int k;
715  register int nx;
716  register int ny;
717  int lx, ly;
718  int length = 0;
719
720  lx = x[1];
721  ly = y[1];
722  tmove2(lx, ly);
723
724  /*
725   * Solve for derivatives of the curve at each point separately for x and y
726   * (parametric).
727   */
728  Paramaterize(x, y, h, numpoints);
729
730  /* closed curve */
731  if ((x[1] == x[numpoints]) && (y[1] == y[numpoints])) {
732    PeriodicSpline(h, x, dx, d2x, d3x, numpoints);
733    PeriodicSpline(h, y, dy, d2y, d3y, numpoints);
734  } else {
735    NaturalEndSpline(h, x, dx, d2x, d3x, numpoints);
736    NaturalEndSpline(h, y, dy, d2y, d3y, numpoints);
737  }
738
739  /*
740   * generate the curve using the above information and PointsPerInterval
741   * vectors between each specified knot.
742   */
743
744  for (j = 1; j < numpoints; ++j) {
745    if ((x[j] == x[j + 1]) && (y[j] == y[j + 1]))
746      continue;
747    for (k = 0; k <= PointsPerInterval; ++k) {
748      t = (double) k *h[j] / (double) PointsPerInterval;
749      t2 = t * t;
750      t3 = t * t * t;
751      nx = x[j] + (int) (t * dx[j] + t2 * d2x[j] / 2 + t3 * d3x[j] / 6);
752      ny = y[j] + (int) (t * dy[j] + t2 * d2y[j] / 2 + t3 * d3y[j] / 6);
753      HGtline(nx, ny);
754      if (length++ > LINELENGTH) {
755	length = 0;
756	printf("\\\n");
757      }
758    }				/* end for k */
759  }				/* end for j */
760}				/* end HGCurve */
761
762
763/*----------------------------------------------------------------------------
764 | Routine:	Paramaterize (xpoints, ypoints, hparams, num_points)
765 |
766 | Results:	This routine calculates parameteric values for use in
767 |		calculating curves.  The parametric values are returned
768 |		in the array h.  The values are an approximation of
769 |		cumulative arc lengths of the curve (uses cord length).
770 |		For additional information, see paper cited below.
771 *----------------------------------------------------------------------------*/
772
773void
774Paramaterize(int x[],
775	     int y[],
776	     double h[],
777	     int n)
778{
779  register int dx;
780  register int dy;
781  register int i;
782  register int j;
783  double u[MAXPOINTS];
784
785  for (i = 1; i <= n; ++i) {
786    u[i] = 0;
787    for (j = 1; j < i; j++) {
788      dx = x[j + 1] - x[j];
789      dy = y[j + 1] - y[j];
790      /* Here was overflowing, so I changed it.       */
791      /* u[i] += sqrt ((double) (dx * dx + dy * dy)); */
792      u[i] += groff_hypot((double) dx, (double) dy);
793    }
794  }
795  for (i = 1; i < n; ++i)
796    h[i] = u[i + 1] - u[i];
797}				/* end Paramaterize */
798
799
800/*----------------------------------------------------------------------------
801 | Routine:	PeriodicSpline (h, z, dz, d2z, d3z, npoints)
802 |
803 | Results:	This routine solves for the cubic polynomial to fit a spline
804 |		curve to the the points specified by the list of values.
805 |		The Curve generated is periodic.  The algorithms for this
806 |		curve are from the `Spline Curve Techniques' paper cited
807 |		above.
808 *----------------------------------------------------------------------------*/
809
810void
811PeriodicSpline(double h[],	/* paramaterization  */
812	       int z[],		/* point list */
813	       double dz[],	/* to return the 1st derivative */
814	       double d2z[],	/* 2nd derivative */
815	       double d3z[],	/* 3rd derivative */
816	       int npoints)	/* number of valid points */
817{
818  double d[MAXPOINTS];
819  double deltaz[MAXPOINTS], a[MAXPOINTS], b[MAXPOINTS];
820  double c[MAXPOINTS], r[MAXPOINTS], s[MAXPOINTS];
821  int i;
822
823  /* step 1 */
824  for (i = 1; i < npoints; ++i) {
825    deltaz[i] = h[i] ? ((double) (z[i + 1] - z[i])) / h[i] : 0;
826  }
827  h[0] = h[npoints - 1];
828  deltaz[0] = deltaz[npoints - 1];
829
830  /* step 2 */
831  for (i = 1; i < npoints - 1; ++i) {
832    d[i] = deltaz[i + 1] - deltaz[i];
833  }
834  d[0] = deltaz[1] - deltaz[0];
835
836  /* step 3a */
837  a[1] = 2 * (h[0] + h[1]);
838  b[1] = d[0];
839  c[1] = h[0];
840  for (i = 2; i < npoints - 1; ++i) {
841    a[i] = 2 * (h[i - 1] + h[i]) -
842	   pow((double) h[i - 1], (double) 2.0) / a[i - 1];
843    b[i] = d[i - 1] - h[i - 1] * b[i - 1] / a[i - 1];
844    c[i] = -h[i - 1] * c[i - 1] / a[i - 1];
845  }
846
847  /* step 3b */
848  r[npoints - 1] = 1;
849  s[npoints - 1] = 0;
850  for (i = npoints - 2; i > 0; --i) {
851    r[i] = -(h[i] * r[i + 1] + c[i]) / a[i];
852    s[i] = (6 * b[i] - h[i] * s[i + 1]) / a[i];
853  }
854
855  /* step 4 */
856  d2z[npoints - 1] = (6 * d[npoints - 2] - h[0] * s[1]
857		      - h[npoints - 1] * s[npoints - 2])
858		     / (h[0] * r[1] + h[npoints - 1] * r[npoints - 2]
859		      + 2 * (h[npoints - 2] + h[0]));
860  for (i = 1; i < npoints - 1; ++i) {
861    d2z[i] = r[i] * d2z[npoints - 1] + s[i];
862  }
863  d2z[npoints] = d2z[1];
864
865  /* step 5 */
866  for (i = 1; i < npoints; ++i) {
867    dz[i] = deltaz[i] - h[i] * (2 * d2z[i] + d2z[i + 1]) / 6;
868    d3z[i] = h[i] ? (d2z[i + 1] - d2z[i]) / h[i] : 0;
869  }
870}				/* end PeriodicSpline */
871
872
873/*----------------------------------------------------------------------------
874 | Routine:	NaturalEndSpline (h, z, dz, d2z, d3z, npoints)
875 |
876 | Results:	This routine solves for the cubic polynomial to fit a spline
877 |		curve the the points specified by the list of values.  The
878 |		alogrithms for this curve are from the `Spline Curve
879 |		Techniques' paper cited above.
880 *----------------------------------------------------------------------------*/
881
882void
883NaturalEndSpline(double h[],	/* parameterization */
884		 int z[],	/* Point list */
885		 double dz[],	/* to return the 1st derivative */
886		 double d2z[],	/* 2nd derivative */
887		 double d3z[],	/* 3rd derivative */
888		 int npoints)	/* number of valid points */
889{
890  double d[MAXPOINTS];
891  double deltaz[MAXPOINTS], a[MAXPOINTS], b[MAXPOINTS];
892  int i;
893
894  /* step 1 */
895  for (i = 1; i < npoints; ++i) {
896    deltaz[i] = h[i] ? ((double) (z[i + 1] - z[i])) / h[i] : 0;
897  }
898  deltaz[0] = deltaz[npoints - 1];
899
900  /* step 2 */
901  for (i = 1; i < npoints - 1; ++i) {
902    d[i] = deltaz[i + 1] - deltaz[i];
903  }
904  d[0] = deltaz[1] - deltaz[0];
905
906  /* step 3 */
907  a[0] = 2 * (h[2] + h[1]);
908  b[0] = d[1];
909  for (i = 1; i < npoints - 2; ++i) {
910    a[i] = 2 * (h[i + 1] + h[i + 2]) -
911	    pow((double) h[i + 1], (double) 2.0) / a[i - 1];
912    b[i] = d[i + 1] - h[i + 1] * b[i - 1] / a[i - 1];
913  }
914
915  /* step 4 */
916  d2z[npoints] = d2z[1] = 0;
917  for (i = npoints - 1; i > 1; --i) {
918    d2z[i] = (6 * b[i - 2] - h[i] * d2z[i + 1]) / a[i - 2];
919  }
920
921  /* step 5 */
922  for (i = 1; i < npoints; ++i) {
923    dz[i] = deltaz[i] - h[i] * (2 * d2z[i] + d2z[i + 1]) / 6;
924    d3z[i] = h[i] ? (d2z[i + 1] - d2z[i]) / h[i] : 0;
925  }
926}				/* end NaturalEndSpline */
927
928
929/*----------------------------------------------------------------------------*
930 | Routine:	change (x_position, y_position, visible_flag)
931 |
932 | Results:	As HGtline passes from the invisible to visible (or vice
933 |		versa) portion of a line, change is called to either draw
934 |		the line, or initialize the beginning of the next one.
935 |		Change calls line to draw segments if visible_flag is set
936 |		(which means we're leaving a visible area).
937 *----------------------------------------------------------------------------*/
938
939void
940change(register int x,
941       register int y,
942       register int vis)
943{
944  static int length = 0;
945
946  if (vis) {			/* leaving a visible area, draw it. */
947    line(x, y);
948    if (length++ > LINELENGTH) {
949      length = 0;
950      printf("\\\n");
951    }
952  } else {			/* otherwise, we're entering one, remember */
953				/* beginning                               */
954    tmove2(x, y);
955  }
956}
957
958
959/*----------------------------------------------------------------------------
960 | Routine:	HGtline (xstart, ystart, xend, yend)
961 |
962 | Results:	Draws a line from current position to (x1,y1) using line(x1,
963 |		y1) to place individual segments of dotted or dashed lines.
964 *----------------------------------------------------------------------------*/
965
966void
967HGtline(int x_1,
968	int y_1)
969{
970  register int x_0 = lastx;
971  register int y_0 = lasty;
972  register int dx;
973  register int dy;
974  register int oldcoord;
975  register int res1;
976  register int visible;
977  register int res2;
978  register int xinc;
979  register int yinc;
980  register int dotcounter;
981
982  if (linmod == SOLID) {
983    line(x_1, y_1);
984    return;
985  }
986
987  /* for handling different resolutions */
988  dotcounter = linmod << dotshifter;
989
990  xinc = 1;
991  yinc = 1;
992  if ((dx = x_1 - x_0) < 0) {
993    xinc = -xinc;
994    dx = -dx;
995  }
996  if ((dy = y_1 - y_0) < 0) {
997    yinc = -yinc;
998    dy = -dy;
999  }
1000  res1 = 0;
1001  res2 = 0;
1002  visible = 0;
1003  if (dx >= dy) {
1004    oldcoord = y_0;
1005    while (x_0 != x_1) {
1006      if ((x_0 & dotcounter) && !visible) {
1007	change(x_0, y_0, 0);
1008	visible = 1;
1009      } else if (visible && !(x_0 & dotcounter)) {
1010	change(x_0 - xinc, oldcoord, 1);
1011	visible = 0;
1012      }
1013      if (res1 > res2) {
1014	oldcoord = y_0;
1015	res2 += dx - res1;
1016	res1 = 0;
1017	y_0 += yinc;
1018      }
1019      res1 += dy;
1020      x_0 += xinc;
1021    }
1022  } else {
1023    oldcoord = x_0;
1024    while (y_0 != y_1) {
1025      if ((y_0 & dotcounter) && !visible) {
1026	change(x_0, y_0, 0);
1027	visible = 1;
1028      } else if (visible && !(y_0 & dotcounter)) {
1029	change(oldcoord, y_0 - yinc, 1);
1030	visible = 0;
1031      }
1032      if (res1 > res2) {
1033	oldcoord = x_0;
1034	res2 += dy - res1;
1035	res1 = 0;
1036	x_0 += xinc;
1037      }
1038      res1 += dx;
1039      y_0 += yinc;
1040    }
1041  }
1042  if (visible)
1043    change(x_1, y_1, 1);
1044  else
1045    change(x_1, y_1, 0);
1046}
1047
1048/* EOF */
1049