Annotation of OpenXM_contrib/pari-2.2/src/graph/plotport.c, Revision 1.2
1.2 ! noro 1: /* $Id: plotport.c,v 1.25 2002/06/09 18:49:11 karim Exp $
1.1 noro 2:
3: Copyright (C) 2000 The PARI group.
4:
5: This file is part of the PARI/GP package.
6:
7: PARI/GP is free software; you can redistribute it and/or modify it under the
8: terms of the GNU General Public License as published by the Free Software
9: Foundation. It is distributed in the hope that it will be useful, but WITHOUT
10: ANY WARRANTY WHATSOEVER.
11:
12: Check the License for details. You should have received a copy of it, along
13: with the package; see the file 'COPYING'. If not, write to the Free Software
14: Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */
15:
16: /*******************************************************************/
17: /* */
18: /* PLOT ROUTINES */
19: /* */
20: /*******************************************************************/
21: #include "pari.h"
22: #include "rect.h"
23:
24: extern void push_val(entree *ep, GEN a);
25: extern void pop_val(entree *ep);
26: extern void postdraw0(long *w, long *x, long *y, long lw);
27: extern void postdraw00(long *w, long *x, long *y, long lw, long scale);
28: extern void rectdraw0(long *w, long *x, long *y, long lw, long do_free);
1.2 ! noro 29: static void PARI_get_psplot(void);
1.1 noro 30:
31: static long current_color[NUMRECT];
32: PariRect **rectgraph = NULL;
33: PARI_plot pari_plot, pari_psplot;
1.2 ! noro 34: PARI_plot *pari_plot_engine = &pari_plot;
! 35: PARI_plot pari_X11plot; /* Used if BOTH_GNUPLOT_AND_X11 */
1.1 noro 36: long rectpoint_itype = 0;
37: long rectline_itype = 0;
38:
39: #define STRINGRECT (NUMRECT-2)
40: #define DRAWRECT (NUMRECT-1)
41:
42: #define PLOTH_NUMPOINTS 1000
43: #define PARAM_NUMPOINTS 1500
44: #define RECUR_NUMPOINTS 8
45:
46: #define RECUR_MAXDEPTH 10
47: #define RECUR_PREC 0.001
48: #define PARAMR_MAXDEPTH 10
49:
50: /********************************************************************/
51: /** **/
52: /** LOW-RES PLOT **/
53: /** **/
54: /********************************************************************/
55: #define ISCR 64
56: #define JSCR 22
57: #define BLANK ' '
58: #define ZERO1 ','
59: #define ZERO2 '-'
60: #define ZERO3 '`'
61: #define PICTZERO(j) ((j) % 3 ? ((j) % 3 == 2 ? ZERO3 : ZERO2) : ZERO1)
62: #define YY '|'
63: #define XX_UPPER '\''
64: #define XX_LOWER '.'
65: #define FF1 '_'
66: #define FF2 'x'
67: #define FF3 '"'
68: #define PICT(j) ((j) % 3 ? ((j) % 3 == 2 ? FF3 : FF2) : FF1)
69:
70: static char *
71: dsprintf9(double d, char *buf)
72: {
73: int i = 10;
74:
75: while (--i >= 0) {
76: sprintf(buf, "%9.*g", i, d);
77: if (strlen(buf) <= 9) return buf;
78: }
79: return buf; /* Should not happen? */
80: }
81:
82: typedef unsigned char screen[ISCR+1][JSCR+1];
83:
84: static void
85: fill_gap(screen scr, long i, int jnew, int jpre)
86: {
87: int mid, i_up, i_lo, up, lo;
88:
89: if (jpre < jnew - 2) {
90: up = jnew - 1; i_up = i;
91: lo = jpre + 1; i_lo = i - 1;
92: } else if (jnew < jpre - 2) {
93: up = jpre - 1; i_up = i - 1;
94: lo = jnew + 1; i_lo = i;
95: } else return; /* if gap < 2, leave it as it is. */
96:
97: mid = (jpre+jnew)/2;
98: if (mid>JSCR) mid=JSCR; else if (mid<0) mid=0;
99: if (lo<0) lo=0;
100: if (lo<=JSCR) while (lo <= mid) scr[i_lo][lo++] = ':';
101: if (up>JSCR) up=JSCR;
102: if (up>=0) while (up > mid) scr[i_up][up--] = ':';
103: }
104:
105: #define QUARK ((char*)NULL) /* Used as a special-case */
106: static GEN quark_gen;
107:
108: #define READ_EXPR(s) ((s)==QUARK? quark_gen : lisexpr(s))
109:
110: void
111: plot(entree *ep, GEN a, GEN b, char *ch,GEN ysmlu,GEN ybigu, long prec)
112: {
1.2 ! noro 113: long jz, j, i, sig;
! 114: gpmem_t av = avma, av2, limite;
1.1 noro 115: int jnew, jpre = 0; /* for lint */
116: GEN p1,p2,ysml,ybig,x,diff,dyj,dx,y[ISCR+1];
117: screen scr;
118: char buf[80], z;
119:
120: sig=gcmp(b,a); if (!sig) return;
121: if (sig<0) { x=a; a=b; b=x; }
122: x = cgetr(prec); gaffect(a,x); push_val(ep, x);
123: for (i=1; i<=ISCR; i++) y[i]=cgetr(3);
124: p1 = gdivgs(gsub(b,a), ISCR-1);
125: dx = cgetr(prec); gaffect(p1, dx);
126: ysml=gzero; ybig=gzero;
127: for (j=1; j<=JSCR; j++) scr[1][j]=scr[ISCR][j]=YY;
128: for (i=2; i<ISCR; i++)
129: {
130: scr[i][1] = XX_LOWER;
131: scr[i][JSCR]=XX_UPPER;
132: for (j=2; j<JSCR; j++) scr[i][j]=BLANK;
133: }
134: av2=avma; limite=stack_lim(av2,1);
135: for (i=1; i<=ISCR; i++)
136: {
137: gaffect(READ_EXPR(ch),y[i]);
138: if (gcmp(y[i],ysml)<0) ysml=y[i];
139: if (gcmp(y[i],ybig)>0) ybig=y[i];
140: x = addrr(x,dx);
141: if (low_stack(limite, stack_lim(av2,1)))
142: {
1.2 ! noro 143: gpmem_t tetpil=avma;
1.1 noro 144: if (DEBUGMEM>1) err(warnmem,"plot");
145: x = gerepile(av2,tetpil,rcopy(x));
146: }
147: ep->value = (void*)x;
148: }
149: if (ysmlu) ysml=ysmlu;
150: if (ybigu) ybig=ybigu;
151: avma=av2; diff=gsub(ybig,ysml);
152: if (gcmp0(diff)) { ybig=gaddsg(1,ybig); diff=gun; }
153: dyj = gdivsg((JSCR-1)*3+2,diff);
154: jz = 3-gtolong(gmul(ysml,dyj));
155: av2=avma; z = PICTZERO(jz); jz = jz/3;
156: for (i=1; i<=ISCR; i++)
157: {
158: if (0<=jz && jz<=JSCR) scr[i][jz]=z;
159: j = 3+gtolong(gmul(gsub(y[i],ysml),dyj));
160: jnew = j/3;
161: if (i > 1) fill_gap(scr, i, jnew, jpre);
162: if (0<=jnew && jnew<=JSCR) scr[i][jnew] = PICT(j);
163: avma = av2;
164: jpre = jnew;
165: }
166: p1=cgetr(3); gaffect(ybig,p1); pariputc('\n');
167: pariputsf("%s ", dsprintf9(rtodbl(p1), buf));
168: for (i=1; i<=ISCR; i++) pariputc(scr[i][JSCR]);
169: pariputc('\n');
170: for (j=(JSCR-1); j>=2; j--)
171: {
172: pariputs(" ");
173: for (i=1; i<=ISCR; i++) pariputc(scr[i][j]);
174: pariputc('\n');
175: }
176: p1=cgetr(3); gaffect(ysml,p1);
177: pariputsf("%s ", dsprintf9(rtodbl(p1), buf));
178: for (i=1; i<=ISCR; i++) pariputc(scr[i][1]);
179: pariputc('\n');
180: p1=cgetr(3); gaffect(a,p1);
181: p2=cgetr(3); gaffect(b,p2);
182: pariputsf("%10s%-9.7g%*.7g\n"," ",rtodbl(p1),ISCR-9,rtodbl(p2));
183: pop_val(ep); avma=av;
184: }
185:
186: /********************************************************************/
187: /** **/
188: /** RECTPLOT FUNCTIONS **/
189: /** **/
190: /********************************************************************/
191: void
1.2 ! noro 192: init_graph(void)
1.1 noro 193: {
194: int n;
195:
196: rectgraph = (PariRect**) gpmalloc(sizeof(PariRect*)*NUMRECT);
197: for (n=0; n<NUMRECT; n++)
198: {
199: PariRect *e = (PariRect*) gpmalloc(sizeof(PariRect));
200:
201: e->head = e->tail = NULL;
202: e->sizex = e->sizey = 0;
203: current_color[n] = DEFAULT_COLOR;
204: rectgraph[n] = e;
205: }
206: }
207:
208: void
1.2 ! noro 209: free_graph(void)
1.1 noro 210: {
211: int i;
212:
1.2 ! noro 213: if (!rectgraph)
! 214: return;
1.1 noro 215: for (i=0; i<NUMRECT; i++)
216: {
217: PariRect *e=rectgraph[i];
218:
219: if (RHead(e)) killrect(i);
220: free((void *)e);
221: }
222: free((void *)rectgraph);
1.2 ! noro 223: rectgraph = 0;
1.1 noro 224: }
225:
226: static PariRect *
227: check_rect(long ne)
228: {
229: if (!GOODRECT(ne))
1.2 ! noro 230: err(talker,
! 231: "incorrect rectwindow number in graphic function (%ld not in [0, %ld])",
! 232: ne, NUMRECT-1);
1.1 noro 233: return rectgraph[ne];
234: }
235:
236: static PariRect *
237: check_rect_init(long ne)
238: {
239: PariRect *e = check_rect(ne);
1.2 ! noro 240: if (!RHead(e)) err(talker,"you must initialize the rectwindow first");
1.1 noro 241: return e;
242: }
243:
244: void
245: initrect_gen(long ne, GEN x, GEN y, long flag)
246: {
247: long xi, yi;
248: if (flag) {
249: double xd = gtodouble(x), yd = gtodouble(y);
250:
251: PARI_get_plot(0);
252: xi = w_width - 1; yi = w_height - 1;
253: if (xd)
254: xi = (long)(xd*xi + 0.5);
255: if (yd)
256: yi = (long)(yd*yi + 0.5);
257: } else {
258: xi = itos(x); yi = itos(y);
259: if (!xi || !yi)
260: PARI_get_plot(0);
261: if (!xi)
262: xi = w_width - 1;
263: if (!yi)
264: yi = w_height - 1;
265: }
266: initrect(ne, xi, yi);
267: }
268:
269: void
270: initrect(long ne, long x, long y)
271: {
272: PariRect *e;
273: RectObj *z;
274:
275: if (x<=1 || y<=1) err(talker,"incorrect dimensions in initrect");
276: e = check_rect(ne); if (RHead(e)) killrect(ne);
277:
278: z = (RectObj*) gpmalloc(sizeof(RectObj));
279: RoNext(z) = NULL;
280: RoType(z) = ROt_NULL;
281: RHead(e)=RTail(e)=z;
282: RXsize(e)=x; RYsize(e)=y;
283: RXcursor(e)=0; RYcursor(e)=0;
284: RXscale(e)=1; RXshift(e)=0;
285: RYscale(e)=1; RYshift(e)=0;
286: RHasGraph(e) = 0;
287: }
288:
289: GEN
290: rectcursor(long ne)
291: {
292: PariRect *e = check_rect_init(ne);
293: GEN z=cgetg(3,t_VEC);
294:
295: z[1] = lstoi((long)RXcursor(e)); z[2] = lstoi((long)RYcursor(e));
296: return z;
297: }
298:
299: static void
300: rectscale0(long ne, double x1, double x2, double y1, double y2)
301: {
302: PariRect *e = check_rect_init(ne);
303: double p2,p3;
304:
305: p2 = RXshift(e) + RXscale(e) * RXcursor(e);
306: p3 = RYshift(e) + RYscale(e) * RYcursor(e);
307: RXscale(e) = RXsize(e)/(x2-x1); RXshift(e) = -x1*RXscale(e);
308: RYscale(e) = RYsize(e)/(y1-y2); RYshift(e) = -y2*RYscale(e);
309: RXcursor(e) = (p2 - RXshift(e)) / RXscale(e);
310: RYcursor(e) = (p3 - RYshift(e)) / RYscale(e);
311: }
312:
313: void
314: rectscale(long ne, GEN x1, GEN x2, GEN y1, GEN y2)
315: {
316: rectscale0(ne, gtodouble(x1), gtodouble(x2), gtodouble(y1), gtodouble(y2));
317: }
318:
319: static void
320: rectmove0(long ne, double x, double y, long relative)
321: {
322: PariRect *e = check_rect_init(ne);
323: RectObj *z = (RectObj*) gpmalloc(sizeof(RectObj1P));
324:
325: if (relative)
326: { RXcursor(e) += x; RYcursor(e) += y; }
327: else
328: { RXcursor(e) = x; RYcursor(e) = y; }
329: RoNext(z) = 0; RoType(z) = ROt_MV;
330: RoMVx(z) = RXcursor(e) * RXscale(e) + RXshift(e);
331: RoMVy(z) = RYcursor(e) * RYscale(e) + RYshift(e);
332: if (!RHead(e)) RHead(e)=RTail(e)=z;
333: else { RoNext(RTail(e))=z; RTail(e)=z; }
334: }
335:
336: void
337: rectmove(long ne, GEN x, GEN y)
338: {
339: rectmove0(ne,gtodouble(x),gtodouble(y),0);
340: }
341:
342: void
343: rectrmove(long ne, GEN x, GEN y)
344: {
345: rectmove0(ne,gtodouble(x),gtodouble(y),1);
346: }
347:
348: void
349: rectpoint0(long ne, double x, double y,long relative) /* code = ROt_MV/ROt_PT */
350: {
351: PariRect *e = check_rect_init(ne);
352: RectObj *z = (RectObj*) gpmalloc(sizeof(RectObj1P));
353:
354: if (relative)
355: { RXcursor(e) += x; RYcursor(e) += y; }
356: else
357: { RXcursor(e) = x; RYcursor(e) = y; }
358: RoNext(z)=0;
359: RoPTx(z) = RXcursor(e)*RXscale(e) + RXshift(e);
360: RoPTy(z) = RYcursor(e)*RYscale(e) + RYshift(e);
361: RoType(z) = ( DTOL(RoPTx(z)) < 0
362: || DTOL(RoPTy(z)) < 0 || DTOL(RoPTx(z)) > RXsize(e)
363: || DTOL(RoPTy(z)) > RYsize(e) ) ? ROt_MV : ROt_PT;
364: if (!RHead(e)) RHead(e)=RTail(e)=z;
365: else { RoNext(RTail(e))=z; RTail(e)=z; }
366: RoCol(z)=current_color[ne];
367: }
368:
369: void
370: rectpoint(long ne, GEN x, GEN y)
371: {
372: rectpoint0(ne,gtodouble(x),gtodouble(y),0);
373: }
374:
375: void
376: rectrpoint(long ne, GEN x, GEN y)
377: {
378: rectpoint0(ne,gtodouble(x),gtodouble(y),1);
379: }
380:
381: void
382: rectcolor(long ne, long color)
383: {
384: check_rect(ne);
385: if (!GOODCOLOR(color)) err(talker,"This is not a valid color");
386: current_color[ne]=color;
387: }
388:
389: void
390: rectline0(long ne, double gx2, double gy2, long relative) /* code = ROt_MV/ROt_LN */
391: {
392: double dx,dy,dxy,xmin,xmax,ymin,ymax,x1,y1,x2,y2;
393: PariRect *e = check_rect_init(ne);
394: RectObj *z = (RectObj*) gpmalloc(sizeof(RectObj2P));
395: const double c = 1 + 1e-10;
396:
397: x1 = RXcursor(e)*RXscale(e) + RXshift(e);
398: y1 = RYcursor(e)*RYscale(e) + RYshift(e);
399: if (relative)
400: { RXcursor(e)+=gx2; RYcursor(e)+=gy2; }
401: else
402: { RXcursor(e)=gx2; RYcursor(e)=gy2; }
403: x2 = RXcursor(e)*RXscale(e) + RXshift(e);
404: y2 = RYcursor(e)*RYscale(e) + RYshift(e);
405: xmin = max(min(x1,x2),0); xmax = min(max(x1,x2),RXsize(e));
406: ymin = max(min(y1,y2),0); ymax = min(max(y1,y2),RYsize(e));
407: dxy = x1*y2 - y1*x2; dx = x2-x1; dy = y2-y1;
408: if (dy)
409: {
410: if (dx*dy<0)
411: { xmin = max(xmin,(dxy+RYsize(e)*dx)/dy); xmax=min(xmax,dxy/dy); }
412: else
413: { xmin=max(xmin,dxy/dy); xmax=min(xmax,(dxy+RYsize(e)*dx)/dy); }
414: }
415: if (dx)
416: {
417: if (dx*dy<0)
418: { ymin=max(ymin,(RXsize(e)*dy-dxy)/dx); ymax=min(ymax,-dxy/dx); }
419: else
420: { ymin=max(ymin,-dxy/dx); ymax=min(ymax,(RXsize(e)*dy-dxy)/dx); }
421: }
422: RoNext(z)=0;
423: RoLNx1(z) = xmin; RoLNx2(z) = xmax;
424: if (dx*dy<0) { RoLNy1(z) = ymax; RoLNy2(z) = ymin; }
425: else { RoLNy1(z) = ymin; RoLNy2(z) = ymax; }
426: RoType(z) = (xmin>xmax*c || ymin>ymax*c) ? ROt_MV : ROt_LN;
427: if (!RHead(e)) RHead(e)=RTail(e)=z;
428: else { RoNext(RTail(e))=z; RTail(e)=z; }
429: RoCol(z)=current_color[ne];
430: }
431:
432: /* Given coordinates of ends of a line, and labels l1 l2 attached to the
433: ends, plot ticks where the label coordinate takes "round" values */
434:
435: static void
436: rectticks(PARI_plot *WW, long ne,
437: double dx1, double dy1, double dx2, double dy2,
438: double l1, double l2, long flags)
439: {
440: long dx,dy,dxy,dxy1,x1,y1,x2,y2,nticks,n,n1,dn;
441: double minstep, maxstep, step, l_min, l_max, minl, maxl, dl, dtx, dty, x, y;
442: double ddx, ddy;
443: const double mult[3] = { 2./1., 5./2., 10./5. };
444: PariRect *e = check_rect_init(ne);
445: int do_double = !(flags & TICKS_NODOUBLE);
446:
447: x1 = DTOL(dx1*RXscale(e) + RXshift(e));
448: y1 = DTOL(dy1*RYscale(e) + RYshift(e));
449: x2 = DTOL(dx2*RXscale(e) + RXshift(e));
450: y2 = DTOL(dy2*RYscale(e) + RYshift(e));
451: dx = x2 - x1;
452: dy = y2 - y1;
453: if (dx < 0) dx = -dx;
454: if (dy < 0) dy = -dy;
455: dxy1 = max(dx, dy);
456: if (WW) {
457: dx /= WW->hunit;
458: dy /= WW->vunit;
459: }
460: else {
461: PARI_get_plot(0);
462: dx /= h_unit;
463: dy /= v_unit;
464: }
465: dxy = (long)sqrt(dx*dx + dy*dy);
466: nticks = (long) ((dxy + 2.5)/4);
467: if (!nticks) return;
468:
469: /* Now we want to find nticks (or less) "round" numbers between l1 and l2.
470: For our purpose round numbers have "last significant" digit either
471: *) any;
472: *) even;
473: *) divisible by 5.
474: We need to choose which alternative is better.
475: */
476: if (l1 < l2)
477: l_min = l1, l_max = l2;
478: else
479: l_min = l2, l_max = l1;
480: minstep = (l_max - l_min)/(nticks + 1);
481: maxstep = 2.5*(l_max - l_min);
482: step = exp(log(10) * floor(log10(minstep)));
483: if (!(flags & TICKS_ENDSTOO)) {
484: double d = 2*(l_max - l_min)/dxy1; /* Two pixels off */
485:
486: l_min += d;
487: l_max -= d;
488: }
489: for (n = 0; ; n++) {
490: if (step >= maxstep) return;
491:
492: if (step >= minstep) {
493: minl = ceil(l_min/step);
494: maxl = floor(l_max/step);
495: if (minl <= maxl && maxl - minl + 1 <= nticks) {
496: nticks = (long) (maxl - minl + 1);
497: l_min = minl * step;
498: l_max = maxl * step; break;
499: }
500: }
501: step *= mult[ n % 3 ];
502: }
503: /* Where to position doubleticks, variants:
504: small: each 5, double: each 10 (n===2 mod 3)
505: small: each 2, double: each 10 (n===1 mod 3)
506: small: each 1, double: each 5 */
507: dn = (n % 3 == 2)? 2: 5;
508: n1 = ((long)minl) % dn; /* unused if do_double = FALSE */
509:
510: /* now l_min and l_max keep min/max values of l with ticks, and nticks is
511: the number of ticks to draw. */
512: if (nticks == 1) ddx = ddy = 0; /* unused: for lint */
513: else {
514: dl = (l_max - l_min)/(nticks - 1);
515: ddx = (dx2 - dx1) * dl / (l2 - l1);
516: ddy = (dy2 - dy1) * dl / (l2 - l1);
517: }
518: x = dx1 + (dx2 - dx1) * (l_min - l1) / (l2 - l1);
519: y = dy1 + (dy2 - dy1) * (l_min - l1) / (l2 - l1);
520: /* assume h_unit and v_unit form a square. For clockwise ticks: */
521: dtx = h_unit * dy/dxy * (y2 > y1 ? 1 : -1); /* y-coord runs down */
522: dty = v_unit * dx/dxy * (x2 > x1 ? 1 : -1);
523: for (n = 0; n < nticks; n++) {
524: RectObj *z = (RectObj*) gpmalloc(sizeof(RectObj2P));
525: double lunit = h_unit > 1 ? 1.5 : 2;
526: double l = (do_double && (n + n1) % dn == 0) ? lunit: 1;
527:
528: RoNext(z) = 0;
529: RoLNx1(z) = RoLNx2(z) = x*RXscale(e) + RXshift(e);
530: RoLNy1(z) = RoLNy2(z) = y*RYscale(e) + RYshift(e);
531:
532: if (flags & TICKS_CLOCKW) {
533: RoLNx1(z) += dtx*l;
534: RoLNy1(z) -= dty*l; /* y-coord runs down */
535: }
536: if (flags & TICKS_ACLOCKW) {
537: RoLNx2(z) -= dtx*l;
538: RoLNy2(z) += dty*l; /* y-coord runs down */
539: }
540: RoType(z) = ROt_LN;
541:
542: if (!RHead(e))
543: RHead(e)=RTail(e)=z;
544: else {
545: RoNext(RTail(e))=z; RTail(e)=z;
546: }
547: RoCol(z)=current_color[ne];
548: x += ddx;
549: y += ddy;
550: }
551: }
552:
553: void
554: rectline(long ne, GEN gx2, GEN gy2)
555: {
556: rectline0(ne, gtodouble(gx2), gtodouble(gy2),0);
557: }
558:
559: void
560: rectrline(long ne, GEN gx2, GEN gy2)
561: {
562: rectline0(ne, gtodouble(gx2), gtodouble(gy2),1);
563: }
564:
565: void
566: rectbox0(long ne, double gx2, double gy2, long relative)
567: {
568: double x1,y1,x2,y2,xmin,ymin,xmax,ymax;
569: double xx,yy;
570: PariRect *e = check_rect_init(ne);
571: RectObj *z = (RectObj*) gpmalloc(sizeof(RectObj2P));
572:
573: x1 = RXcursor(e)*RXscale(e) + RXshift(e);
574: y1 = RYcursor(e)*RYscale(e) + RYshift(e);
575: if (relative)
576: { xx = RXcursor(e)+gx2; yy = RYcursor(e)+gy2; }
577: else
578: { xx = gx2; yy = gy2; }
579: x2 = xx*RXscale(e) + RXshift(e);
580: y2 = yy*RYscale(e) + RYshift(e);
581: xmin = max(min(x1,x2),0); xmax = min(max(x1,x2),RXsize(e));
582: ymin = max(min(y1,y2),0); ymax = min(max(y1,y2),RYsize(e));
583:
584: RoNext(z)=0; RoType(z) = ROt_BX;
585: RoBXx1(z) = xmin; RoBXy1(z) = ymin;
586: RoBXx2(z) = xmax; RoBXy2(z) = ymax;
587: if (!RHead(e)) RHead(e)=RTail(e)=z;
588: else { RoNext(RTail(e))=z; RTail(e)=z; }
589: RoCol(z)=current_color[ne];
590: }
591:
592: void
593: rectbox(long ne, GEN gx2, GEN gy2)
594: {
595: rectbox0(ne, gtodouble(gx2), gtodouble(gy2), 0);
596: }
597:
598: void
599: rectrbox(long ne, GEN gx2, GEN gy2)
600: {
601: rectbox0(ne, gtodouble(gx2), gtodouble(gy2), 1);
602: }
603:
604: void
605: killrect(long ne)
606: {
607: RectObj *p1,*p2;
608: PariRect *e = check_rect_init(ne);
609:
610: current_color[ne]=DEFAULT_COLOR;
611: p1=RHead(e);
612: RHead(e) = RTail(e) = NULL;
613: RXsize(e) = RYsize(e) = 0;
614: RXcursor(e) = RYcursor(e) = 0;
615: RXscale(e) = RYscale(e) = 1;
616: RXshift(e) = RYshift(e) = 0;
617: while (p1)
618: {
619: if (RoType(p1)==ROt_MP || RoType(p1)==ROt_ML)
620: {
621: free(RoMPxs(p1)); free(RoMPys(p1));
622: }
623: if (RoType(p1)==ROt_ST) free(RoSTs(p1));
624: p2=RoNext(p1); free(p1); p1=p2;
625: }
626: }
627:
628: void
629: rectpoints0(long ne, double *listx, double *listy, long lx) /* code = ROt_MP */
630: {
631: double *ptx, *pty, x, y;
632: long i, cp=0;
633: PariRect *e = check_rect_init(ne);
634: RectObj *z = (RectObj*) gpmalloc(sizeof(RectObjMP));
635:
636: ptx=(double*) gpmalloc(lx*sizeof(double));
637: pty=(double*) gpmalloc(lx*sizeof(double));
638: for (i=0; i<lx; i++)
639: {
640: x = RXscale(e)*listx[i] + RXshift(e);
641: y = RYscale(e)*listy[i] + RYshift(e);
642: if ((x>=0)&&(y>=0)&&(x<=RXsize(e))&&(y<=RYsize(e)))
643: {
644: ptx[cp]=x; pty[cp]=y; cp++;
645: }
646: }
647: RoNext(z)=0; RoType(z) = ROt_MP;
648: RoMPcnt(z) = cp; RoMPxs(z) = ptx; RoMPys(z) = pty;
649: if (!RHead(e)) RHead(e)=RTail(e)=z;
650: else { RoNext(RTail(e))=z; RTail(e)=z; }
651: RoCol(z)=current_color[ne];
652: }
653:
654: void
655: rectpoints(long ne, GEN listx, GEN listy)
656: {
657: long i,lx, tx=typ(listx), ty=typ(listy);
658: double *px,*py;
659:
660: if (!is_matvec_t(tx) || !is_matvec_t(ty))
661: {
662: rectpoint(ne, listx, listy); return;
663: }
664: if (tx == t_MAT || ty == t_MAT) err(ploter4);
665: lx=lg(listx); if (lg(listy)!=lx) err(ploter5);
666: lx--; if (!lx) return;
667:
668: px = (double*) gpmalloc(lx*sizeof(double));
669: py = (double*) gpmalloc(lx*sizeof(double));
670: for (i=0; i<lx; i++)
671: {
672: px[i]=gtodouble((GEN)listx[i+1]); py[i]=gtodouble((GEN)listy[i+1]);
673: }
674: rectpoints0(ne,px,py,lx);
675: free(px); free(py);
676: }
677:
678: void
679: rectlines0(long ne, double *x, double *y, long lx, long flag) /* code = ROt_ML */
680: {
681: long i,I;
682: double *ptx,*pty;
683: PariRect *e = check_rect_init(ne);
684: RectObj *z = (RectObj*) gpmalloc(sizeof(RectObj2P));
685:
686: I = flag ? lx+1 : lx;
687: ptx = (double*) gpmalloc(I*sizeof(double));
688: pty = (double*) gpmalloc(I*sizeof(double));
689: for (i=0; i<lx; i++)
690: {
691: ptx[i] = RXscale(e)*x[i] + RXshift(e);
692: pty[i] = RYscale(e)*y[i] + RYshift(e);
693: }
694: if (flag)
695: {
696: ptx[i] = RXscale(e)*x[0] + RXshift(e);
697: pty[i] = RYscale(e)*y[0] + RYshift(e);
698: }
699: RoNext(z)=0; RoType(z)=ROt_ML;
700: RoMLcnt(z)=lx; RoMLxs(z)=ptx; RoMLys(z)=pty;
701: if (!RHead(e)) RHead(e)=RTail(e)=z;
702: else { RoNext(RTail(e))=z; RTail(e)=z; }
703: RoCol(z) = current_color[ne];
704: }
705:
706: void
707: rectlines(long ne, GEN listx, GEN listy, long flag)
708: {
709: long tx=typ(listx), ty=typ(listy), lx=lg(listx), i;
710: double *x, *y;
711:
712: if (!is_matvec_t(tx) || !is_matvec_t(ty))
713: {
714: rectline(ne, listx, listy); return;
715: }
716: if (tx == t_MAT || ty == t_MAT) err(ploter4);
717: if (lg(listy)!=lx) err(ploter5);
718: lx--; if (!lx) return;
719:
720: x = (double*) gpmalloc(lx*sizeof(double));
721: y = (double*) gpmalloc(lx*sizeof(double));
722: for (i=0; i<lx; i++)
723: {
724: x[i] = gtodouble((GEN)listx[i+1]); y[i] = gtodouble((GEN)listy[i+1]);
725: }
726: rectlines0(ne,x,y,lx,flag);
727: free(x); free(y);
728: }
729:
730: static void
731: put_string(long win, long x, long y, char *str, long dir)
732: {
733: rectmove0(win,(double)x,(double)y,0); rectstring3(win,str,dir);
734: }
735:
736: void
737: rectstring(long ne, char *str)
738: {
739: rectstring3(ne,str,RoSTdirLEFT);
740: }
741:
742: /* Allocate memory, then put string */
743: void
744: rectstring3(long ne, char *str, long dir) /* code = ROt_ST */
745: {
746: PariRect *e = check_rect_init(ne);
747: RectObj *z = (RectObj*) gpmalloc(sizeof(RectObjST));
748: long l = strlen(str);
749: char *s = (char *) gpmalloc(l+1);
750:
751: strcpy(s,str);
752: RoNext(z)=0; RoType(z) = ROt_ST;
753: RoSTl(z) = l; RoSTs(z) = s;
754: RoSTx(z) = RXscale(e)*RXcursor(e)+RXshift(e);
755: RoSTy(z) = RYscale(e)*RYcursor(e)+RYshift(e);
756: RoSTdir(z) = dir;
757: if (!RHead(e)) RHead(e)=RTail(e)=z;
758: else { RoNext(RTail(e))=z; RTail(e)=z; }
759: RoCol(z)=current_color[ne];
760: }
761:
762: void
763: rectpointtype(long ne, long type) /* code = ROt_PTT */
764: {
765: if (ne == -1) {
766: rectpoint_itype = type;
767: } else {
768: PariRect *e = check_rect_init(ne);
769: RectObj *z = (RectObj*) gpmalloc(sizeof(RectObjPN));
770:
771: RoNext(z) = 0; RoType(z) = ROt_PTT;
772: RoPTTpen(z) = type;
773: if (!RHead(e)) RHead(e)=RTail(e)=z;
774: else { RoNext(RTail(e))=z; RTail(e)=z; }
775: }
776: }
777:
778: void
779: rectpointsize(long ne, GEN size) /* code = ROt_PTS */
780: {
781: if (ne == -1) {
782: set_pointsize(gtodouble(size)); /* Immediate set */
783: } else {
784: PariRect *e = check_rect_init(ne);
785: RectObj *z = (RectObj*) gpmalloc(sizeof(RectObjPS));
786:
787: RoNext(z) = 0; RoType(z) = ROt_PTS;
788: RoPTSsize(z) = gtodouble(size);
789: if (!RHead(e)) RHead(e)=RTail(e)=z;
790: else { RoNext(RTail(e))=z; RTail(e)=z; }
791: }
792: }
793:
794: void
795: rectlinetype(long ne, long type)
796: {
797: if (ne == -1) {
798: rectline_itype = type;
799: } else {
800: PariRect *e = check_rect_init(ne);
801: RectObj *z = (RectObj*) gpmalloc(sizeof(RectObjPN));
802:
803: RoNext(z) = 0; RoType(z) = ROt_LNT;
804: RoLNTpen(z) = type;
805: if (!RHead(e)) RHead(e)=RTail(e)=z;
806: else { RoNext(RTail(e))=z; RTail(e)=z; }
807: }
808: }
809:
810: void
811: rectcopy_gen(long source, long dest, GEN xoff, GEN yoff, long flag)
812: {
813: long xi, yi;
814: if (flag & RECT_CP_RELATIVE) {
815: double xd = gtodouble(xoff), yd = gtodouble(yoff);
816:
817: PARI_get_plot(0);
818: xi = w_width - 1; yi = w_height - 1;
819: xi = (long)(xd*xi + 0.5);
820: yi = (long)(yd*yi + 0.5);
821: } else {
822: xi = itos(xoff); yi = itos(yoff);
823: }
824: if (flag & ~RECT_CP_RELATIVE) {
825: PariRect *s = check_rect_init(source), *d = check_rect_init(dest);
826:
827: switch (flag & ~RECT_CP_RELATIVE) {
828: case RECT_CP_NW:
829: break;
830: case RECT_CP_SW:
831: yi = RYsize(d) - RYsize(s) - yi;
832: break;
833: case RECT_CP_SE:
834: yi = RYsize(d) - RYsize(s) - yi;
835: /* FALL THROUGH */
836: case RECT_CP_NE:
837: xi = RXsize(d) - RXsize(s) - xi;
838: break;
839: }
840: }
841: rectcopy(source, dest, xi, yi);
842: }
843:
844: void
845: rectcopy(long source, long dest, long xoff, long yoff)
846: {
847: PariRect *s = check_rect_init(source), *d = check_rect_init(dest);
848: RectObj *p1 = RHead(s);
849: RectObj *tail = RTail(d);
850: RectObj *next;
851: int i;
852:
853: while (p1)
854: {
855: switch(RoType(p1))
856: {
857: case ROt_PT:
858: next = (RectObj*) gpmalloc(sizeof(RectObj1P));
859: memcpy(next,p1,sizeof(RectObj1P));
860: RoPTx(next) += xoff; RoPTy(next) += yoff;
861: RoNext(tail) = next; tail = next;
862: break;
863: case ROt_LN: case ROt_BX:
864: next = (RectObj*) gpmalloc(sizeof(RectObj2P));
865: memcpy(next,p1,sizeof(RectObj2P));
866: RoLNx1(next) += xoff; RoLNy1(next) += yoff;
867: RoLNx2(next) += xoff; RoLNy2(next) += yoff;
868: RoNext(tail) = next; tail = next;
869: break;
870: case ROt_MP: case ROt_ML:
871: next = (RectObj*) gpmalloc(sizeof(RectObjMP));
872: memcpy(next,p1,sizeof(RectObjMP));
873: RoMPxs(next) = (double*) gpmalloc(sizeof(double)*RoMPcnt(next));
874: RoMPys(next) = (double*) gpmalloc(sizeof(double)*RoMPcnt(next));
875: memcpy(RoMPxs(next),RoMPxs(p1),sizeof(double)*RoMPcnt(next));
876: memcpy(RoMPys(next),RoMPys(p1),sizeof(double)*RoMPcnt(next));
877: for (i=0; i<RoMPcnt(next); i++)
878: {
879: RoMPxs(next)[i] += xoff; RoMPys(next)[i] += yoff;
880: }
881: RoNext(tail) = next; tail = next;
882: break;
883: case ROt_ST:
884: next = (RectObj*) gpmalloc(sizeof(RectObjST));
885: memcpy(next,p1,sizeof(RectObjMP));
886: RoSTs(next) = (char*) gpmalloc(RoSTl(p1)+1);
887: memcpy(RoSTs(next),RoSTs(p1),RoSTl(p1)+1);
888: RoSTx(next) += xoff; RoSTy(next) += yoff;
889: RoNext(tail) = next; tail = next;
890: break;
891: case ROt_PTT: case ROt_LNT: case ROt_PTS:
892: next = (RectObj*) gpmalloc(sizeof(RectObjPN));
893: memcpy(next,p1,sizeof(RectObjPN));
894: RoNext(tail) = next; tail = next;
895: break;
896: }
897: p1=RoNext(p1);
898: }
899: RoNext(tail) = NULL; RTail(d) = tail;
900: }
901:
902: #define CLIPLINE_NONEMPTY 1
903: #define CLIPLINE_CLIP_1 2
904: #define CLIPLINE_CLIP_2 4
905:
906: /* A simpler way is to clip by 4 half-planes */
907: static int
908: clipline(double xmin, double xmax, double ymin, double ymax,
909: double *x1p, double *y1p, double *x2p, double *y2p)
910: {
911: int xy_exch = 0, rc = CLIPLINE_NONEMPTY;
912: double t, sl;
913: double xi, xmn, xmx;
914: double yi, ymn, ymx;
915: int x1_is_ymn, x1_is_xmn;
916: double x1 = *x1p, x2 = *x2p, y1 = *y1p, y2 = *y2p;
917:
918: if ((x1 < xmin && x2 < xmin) || (x1 > xmax && x2 > xmax))
919: return 0;
920: if (fabs(x1 - x2) < fabs(y1 - y2)) { /* Exchange x and y */
921: xy_exch = 1;
922: t = xmin; xmin = ymin; ymin = t;
923: t = xmax; xmax = ymax; ymax = t;
924: t = x1; x1 = y1; y1 = t;
925: t = x2; x2 = y2; y2 = t;
926: }
927:
928: /* Build y as a function of x */
929: xi = x1, yi = y1;
930: sl = ( (x1==x2) ? 0 : (y2 - yi)/(x2 - xi) );
931:
932: xmn = x1, xmx = x2;
933: if (x1 > x2) {
934: x1_is_xmn = 0;
935: xmn = x2, xmx = x1;
936: } else
937: x1_is_xmn = 1;
938:
939: if (xmn < xmin) {
940: xmn = xmin;
941: rc |= (x1_is_xmn ? CLIPLINE_CLIP_1 : CLIPLINE_CLIP_2);
942: }
943:
944: if (xmx > xmax) {
945: xmx = xmax;
946: rc |= (x1_is_xmn ? CLIPLINE_CLIP_2 : CLIPLINE_CLIP_1);
947: }
948:
949: if (xmn > xmx)
950: return 0;
951:
952: ymn = yi + (xmn - xi)*sl;
953: ymx = yi + (xmx - xi)*sl;
954:
955: if (sl < 0)
956: t = ymn, ymn = ymx, ymx = t;
957: if (ymn > ymax || ymx < ymin)
958: return 0;
959:
960: if (rc & CLIPLINE_CLIP_1)
961: x1 = (x1_is_xmn ? xmn : xmx);
962: if (rc & CLIPLINE_CLIP_2)
963: x2 = (x1_is_xmn ? xmx : xmn);
964:
965: /* Now we know there is an intersection, need to move x1 and x2 */
966: x1_is_ymn = ((sl >= 0) == (x1 < x2));
967: if (ymn < ymin) {
968: double x = (ymin - yi)/sl + xi; /* slope != 0 ! */
969:
970: if (x1_is_ymn)
971: x1 = x, rc |= CLIPLINE_CLIP_1;
972: else
973: x2 = x, rc |= CLIPLINE_CLIP_2;
974: }
975: if (ymx > ymax) {
976: double x = (ymax - yi)/sl + xi; /* slope != 0 ! */
977:
978: if (x1_is_ymn)
979: x2 = x, rc |= CLIPLINE_CLIP_2;
980: else
981: x1 = x, rc |= CLIPLINE_CLIP_1;
982: }
983: if (rc & CLIPLINE_CLIP_1)
984: y1 = yi + (x1 - xi)*sl;
985: if (rc & CLIPLINE_CLIP_2)
986: y2 = yi + (x2 - xi)*sl;
987: if (xy_exch) /* Exchange x and y */
988: *x1p = y1, *x2p = y2, *y1p = x1, *y2p = x2;
989: else
990: *x1p = x1, *x2p = x2, *y1p = y1, *y2p = y2;
991: return rc;
992: }
993:
994: void
995: rectclip(long rect)
996: {
997: PariRect *s = check_rect_init(rect);
998: RectObj *p1 = RHead(s);
999: RectObj **prevp = &RHead(s);
1000: RectObj *next;
1001: double xmin = 0;
1002: double xmax = RXsize(s);
1003: double ymin = 0;
1004: double ymax = RYsize(s);
1005:
1006: while (p1) {
1007: int did_clip = 0;
1008:
1009: next = RoNext(p1);
1010: switch(RoType(p1)) {
1011: case ROt_PT:
1012: if ( DTOL(RoPTx(p1)) < xmin || DTOL(RoPTx(p1)) > xmax
1013: || DTOL(RoPTy(p1)) < ymin || DTOL(RoPTy(p1)) > ymax) {
1014: remove:
1015: *prevp = next;
1016: free(p1);
1017: break;
1018: }
1019: goto do_next;
1020: case ROt_BX:
1021: if (RoLNx1(p1) < xmin)
1022: RoLNx1(p1) = xmin, did_clip = 1;
1023: if (RoLNx2(p1) < xmin)
1024: RoLNx2(p1) = xmin, did_clip = 1;
1025: if (RoLNy1(p1) < ymin)
1026: RoLNy1(p1) = ymin, did_clip = 1;
1027: if (RoLNy2(p1) < ymin)
1028: RoLNy2(p1) = ymin, did_clip = 1;
1029: if (RoLNx1(p1) > xmax)
1030: RoLNx1(p1) = xmax, did_clip = 1;
1031: if (RoLNx2(p1) > xmax)
1032: RoLNx2(p1) = xmax, did_clip = 1;
1033: if (RoLNy1(p1) > ymax)
1034: RoLNy1(p1) = ymax, did_clip = 1;
1035: if (RoLNy2(p1) > ymax)
1036: RoLNy2(p1) = ymax, did_clip = 1;
1037: /* Remove zero-size clipped boxes */
1038: if ( did_clip
1039: && RoLNx1(p1) == RoLNx2(p1) && RoLNy1(p1) == RoLNy2(p1) )
1040: goto remove;
1041: goto do_next;
1042: case ROt_LN:
1043: if (!clipline(xmin, xmax, ymin, ymax,
1044: &RoLNx1(p1), &RoLNy1(p1), &RoLNx2(p1), &RoLNy2(p1)))
1045: goto remove;
1046: goto do_next;
1047: case ROt_MP:
1048: {
1049: int c = RoMPcnt(p1);
1050: int f = 0, t = 0;
1051:
1052: while (f < c) {
1053: if ( DTOL(RoMPxs(p1)[f]) >= xmin && DTOL(RoMPxs(p1)[f]) <= xmax
1054: && DTOL(RoMPys(p1)[f]) >= ymin
1055: && DTOL(RoMPys(p1)[f]) <= ymax) {
1056: if (t != f) {
1057: RoMPxs(p1)[t] = RoMPxs(p1)[f];
1058: RoMPys(p1)[t] = RoMPys(p1)[f];
1059: }
1060: t++;
1061: }
1062: f++;
1063: }
1064: if (t == 0)
1065: goto remove;
1066: RoMPcnt(p1) = t;
1067: goto do_next;
1068: }
1069: case ROt_ML:
1070: {
1071: /* Hard case. Here we need to break a multiline into
1072: several pieces if some part is clipped. */
1073: int c = RoMPcnt(p1) - 1;
1074: int f = 0, t = 0, had_lines = 0, had_hole = 0, rc;
1075: double ox = RoMLxs(p1)[0], oy = RoMLys(p1)[0], oxn, oyn;
1076:
1077: while (f < c) {
1078: /* Endpoint of this segment is the startpoint of the
1079: next one, so we need to preserve it if it is clipped. */
1080: oxn = RoMLxs(p1)[f + 1], oyn = RoMLys(p1)[f + 1];
1081: rc = clipline(xmin, xmax, ymin, ymax,
1082: /* &RoMLxs(p1)[f], &RoMLys(p1)[f], */
1083: &ox, &oy,
1084: &RoMLxs(p1)[f+1], &RoMLys(p1)[f+1]);
1085: RoMLxs(p1)[f] = ox; RoMLys(p1)[f] = oy;
1086: ox = oxn; oy = oyn;
1087: if (!rc) {
1088: if (had_lines)
1089: had_hole = 1;
1090: f++;
1091: continue;
1092: }
1093:
1094: if ( !had_lines || (!(rc & CLIPLINE_CLIP_1) && !had_hole) ) {
1095: /* Continuous */
1096: had_lines = 1;
1097: if (t != f) {
1098: if (t == 0) {
1099: RoMPxs(p1)[t] = RoMPxs(p1)[f];
1100: RoMPys(p1)[t] = RoMPys(p1)[f];
1101: }
1102: RoMPxs(p1)[t+1] = RoMPxs(p1)[f+1];
1103: RoMPys(p1)[t+1] = RoMPys(p1)[f+1];
1104: }
1105: t++;
1106: f++;
1107: if ( rc & CLIPLINE_CLIP_2)
1108: had_hole = 1, RoMLcnt(p1) = t + 1;
1109: continue;
1110: }
1111: /* Is not continuous, automatically p1 is not free()ed. */
1112: RoMLcnt(p1) = t + 1;
1113: if ( rc & CLIPLINE_CLIP_2) { /* Needs separate entry */
1114: RectObj *n = (RectObj*) gpmalloc(sizeof(RectObj2P));
1115:
1116: RoType(n) = ROt_LN;
1117: RoCol(n) = RoCol(p1);
1118: RoLNx1(n) = RoMLxs(p1)[f]; RoLNy1(n) = RoMLys(p1)[f];
1119: RoLNx2(n) = RoMLxs(p1)[f+1]; RoLNy2(n) = RoMLys(p1)[f+1];
1120: RoNext(n) = next;
1121: RoNext(p1) = n;
1122: /* Restore the unclipped value: */
1123: RoMLxs(p1)[f + 1] = oxn; RoMLys(p1)[f + 1] = oyn;
1124: f++;
1125: prevp = &RoNext(n);
1126: }
1127: if (f + 1 < c) { /* Are other lines */
1128: RectObj *n = (RectObj*) gpmalloc(sizeof(RectObjMP));
1129: RoType(n) = ROt_ML;
1130: RoCol(n) = RoCol(p1);
1131: RoMLcnt(n) = c - f;
1132: RoMLxs(n) = (double*) gpmalloc(sizeof(double)*(c - f));
1133: RoMLys(n) = (double*) gpmalloc(sizeof(double)*(c - f));
1134: memcpy(RoMPxs(n),RoMPxs(p1) + f, sizeof(double)*(c - f));
1135: memcpy(RoMPys(n),RoMPys(p1) + f, sizeof(double)*(c - f));
1136: RoMPxs(n)[0] = oxn; RoMPys(n)[0] = oyn;
1137: RoNext(n) = next;
1138: RoNext(p1) = n;
1139: next = n;
1140: }
1141: goto do_next;
1142: }
1143: if (t == 0)
1144: goto remove;
1145: goto do_next;
1146: }
1147: #if 0
1148: case ROt_ST:
1149: next = (RectObj*) gpmalloc(sizeof(RectObjMP));
1150: memcpy(next,p1,sizeof(RectObjMP));
1151: RoSTs(next) = (char*) gpmalloc(RoSTl(p1)+1);
1152: memcpy(RoSTs(next),RoSTs(p1),RoSTl(p1)+1);
1153: RoSTx(next) += xoff; RoSTy(next) += yoff;
1154: RoNext(tail) = next; tail = next;
1155: break;
1156: #endif
1157: default: {
1158: do_next:
1159: prevp = &RoNext(p1);
1160: break;
1161: }
1162: }
1163: p1 = next;
1164: }
1165: }
1166:
1167: /********************************************************************/
1168: /** **/
1169: /** HI-RES PLOT **/
1170: /** **/
1171: /********************************************************************/
1172:
1173: static void
1174: Appendx(dblPointList *f, dblPointList *l,double x)
1175: {
1176: (l->d)[l->nb++]=x;
1177: if (x < f->xsml) f->xsml=x;
1178: else if (x > f->xbig) f->xbig=x;
1179: }
1180:
1181: static void
1182: Appendy(dblPointList *f, dblPointList *l,double y)
1183: {
1184: (l->d)[l->nb++]=y;
1185: if (y < f->ysml) f->ysml=y;
1186: else if (y > f->ybig) f->ybig=y;
1187: }
1188:
1189: /* Convert data from GEN to double before we call rectplothrawin */
1190: static dblPointList*
1191: gtodblList(GEN data, long flags)
1192: {
1193: dblPointList *l;
1194: double xsml,xbig,ysml,ybig;
1195: long tx=typ(data), ty, nl=lg(data)-1, lx1,lx;
1196: register long i,j,u,v;
1197: long param=(flags & PLOT_PARAMETRIC);
1198: GEN x,y;
1199:
1200: if (! is_vec_t(tx)) err(talker,"not a vector in gtodblList");
1201: if (!nl) return NULL;
1.2 ! noro 1202: lx1 = lg(data[1]);
1.1 noro 1203:
1.2 ! noro 1204: if (nl == 1) err(talker,"single vector in gtodblList");
1.1 noro 1205: /* Allocate memory, then convert coord. to double */
1.2 ! noro 1206: l = (dblPointList*) gpmalloc(nl*sizeof(dblPointList));
1.1 noro 1207: for (i=0; i<nl-1; i+=2)
1208: {
1.2 ! noro 1209: u = i+1;
! 1210: x = (GEN)data[u]; tx = typ(x);
! 1211: y = (GEN)data[u+1]; ty = typ(y);
1.1 noro 1212: if (!is_vec_t(tx) || !is_vec_t(ty)) err(ploter4);
1.2 ! noro 1213: lx = lg(x); if (lg(y) != lx) err(ploter5);
1.1 noro 1214: if (!param && lx != lx1) err(ploter5);
1.2 ! noro 1215:
1.1 noro 1216: lx--;
1.2 ! noro 1217: l[i].d = (double*) gpmalloc(lx*sizeof(double));
! 1218: l[u].d = (double*) gpmalloc(lx*sizeof(double));
! 1219: for (j=0; j<lx; j=v)
1.1 noro 1220: {
1.2 ! noro 1221: v = j+1;
! 1222: l[i].d[j] = gtodouble((GEN)x[v]);
! 1223: l[u].d[j] = gtodouble((GEN)y[v]);
1.1 noro 1224: }
1.2 ! noro 1225: l[i].nb = l[u].nb = lx;
1.1 noro 1226: }
1227:
1228: /* Now compute extremas */
1229: if (param)
1230: {
1.2 ! noro 1231: l[0].nb = nl/2;
! 1232: for (i=0; i < l[0].nb; i+=2)
! 1233: if (l[i+1].nb) break;
! 1234: if (i >= l[0].nb) { free(l); return NULL; }
! 1235: xsml = xbig = l[i ].d[0];
! 1236: ysml = ybig = l[i+1].d[0];
! 1237:
! 1238: for (i=0; i < l[0].nb; i+=2)
1.1 noro 1239: {
1.2 ! noro 1240: u = i+1;
! 1241: for (j=0; j < l[u].nb; j++)
1.1 noro 1242: {
1.2 ! noro 1243: if (l[i].d[j] < xsml) xsml = l[i].d[j];
! 1244: else if (l[i].d[j] > xbig) xbig = l[i].d[j];
1.1 noro 1245:
1.2 ! noro 1246: if (l[u].d[j] < ysml) ysml = l[u].d[j];
! 1247: else if (l[u].d[j] > ybig) ybig = l[u].d[j];
1.1 noro 1248: }
1249: }
1250: }
1251: else
1252: {
1.2 ! noro 1253: if (!l[0].nb) { free(l); return NULL; }
! 1254: l[0].nb = nl-1;
! 1255:
! 1256: xsml = xbig = l[0].d[0];
! 1257: ysml = ybig = l[1].d[0];
! 1258:
! 1259: for (j=0; j < l[1].nb; j++)
1.1 noro 1260: {
1.2 ! noro 1261: if (l[0].d[j] < xsml) xsml = l[0].d[j];
! 1262: else if (l[0].d[j] > xbig) xbig = l[0].d[j];
1.1 noro 1263: }
1.2 ! noro 1264: for (i=1; i <= l[0].nb; i++)
! 1265: for (j=0; j < l[i].nb; j++)
1.1 noro 1266: {
1.2 ! noro 1267: if (l[i].d[j] < ysml) ysml = l[i].d[j];
! 1268: else if (l[i].d[j] > ybig) ybig = l[i].d[j];
1.1 noro 1269: }
1270: }
1.2 ! noro 1271: l[0].xsml = xsml; l[0].xbig = xbig;
! 1272: l[0].ysml = ysml; l[0].ybig = ybig;
1.1 noro 1273: return l;
1274: }
1275:
1276: static void
1277: single_recursion(dblPointList *pl,char *ch,entree *ep,GEN xleft,GEN yleft,
1278: GEN xright,GEN yright,long depth)
1279: {
1280: GEN xx,yy;
1.2 ! noro 1281: gpmem_t av=avma;
1.1 noro 1282: double dy=pl[0].ybig - pl[0].ysml;
1283:
1284: if (depth==RECUR_MAXDEPTH) return;
1285:
1286: yy=cgetr(3); xx=gmul2n(gadd(xleft,xright),-1);
1287: ep->value = (void*) xx; gaffect(READ_EXPR(ch), yy);
1288:
1289: if (dy)
1290: {
1291: if (fabs(rtodbl(yleft)+rtodbl(yright)-2*rtodbl(yy))/dy < RECUR_PREC)
1292: return;
1293: }
1294: single_recursion(pl,ch,ep, xleft,yleft, xx,yy, depth+1);
1295:
1296: Appendx(&pl[0],&pl[0],rtodbl(xx));
1297: Appendy(&pl[0],&pl[1],rtodbl(yy));
1298:
1299: single_recursion(pl,ch,ep, xx,yy, xright,yright, depth+1);
1300:
1301: avma=av;
1302: }
1303:
1304: static void
1305: param_recursion(dblPointList *pl,char *ch,entree *ep, GEN tleft,GEN xleft,
1306: GEN yleft, GEN tright,GEN xright,GEN yright, long depth)
1307: {
1308: GEN tt,xx,yy, p1;
1.2 ! noro 1309: gpmem_t av=avma;
1.1 noro 1310: double dy=pl[0].ybig - pl[0].ysml;
1311: double dx=pl[0].xbig - pl[0].xsml;
1312:
1313: if (depth==PARAMR_MAXDEPTH) return;
1314:
1315: xx=cgetr(3); yy=cgetr(3); tt=gmul2n(gadd(tleft,tright),-1);
1316: ep->value = (void*)tt; p1=READ_EXPR(ch);
1317: gaffect((GEN)p1[1],xx); gaffect((GEN)p1[2],yy);
1318:
1319: if (dx && dy)
1320: {
1321: if (fabs(rtodbl(xleft)+rtodbl(xright)-2*rtodbl(xx))/dx < RECUR_PREC &&
1322: fabs(rtodbl(yleft)+rtodbl(yright)-2*rtodbl(yy))/dy < RECUR_PREC)
1323: return;
1324: }
1325: param_recursion(pl,ch,ep, tleft,xleft,yleft, tt,xx,yy, depth+1);
1326:
1327: Appendx(&pl[0],&pl[0],rtodbl(xx));
1328: Appendy(&pl[0],&pl[1],rtodbl(yy));
1329:
1330: param_recursion(pl,ch,ep, tt,xx,yy, tright,xright,yright, depth+1);
1331:
1332: avma=av;
1333: }
1334:
1335: /*
1336: * Pure graphing. If testpoints is 0, it is set to the default.
1337: * Returns a dblPointList of (absolute) coordinates.
1338: */
1339: static dblPointList *
1340: rectplothin(entree *ep, GEN a, GEN b, char *ch, long prec, ulong flags,
1341: long testpoints)
1342: {
1343: long single_c;
1344: long param=flags & PLOT_PARAMETRIC;
1345: long recur=flags & PLOT_RECURSIVE;
1346: GEN p1,dx,x,xleft,xright,yleft,yright,tleft,tright;
1347: dblPointList *pl;
1.2 ! noro 1348: long tx, i, j, sig, nc, nl, nbpoints;
! 1349: gpmem_t av = avma, av2;
1.1 noro 1350: double xsml,xbig,ysml,ybig,fx,fy;
1351:
1352: if (!testpoints)
1353: {
1354: if (recur)
1355: testpoints = RECUR_NUMPOINTS;
1356: else
1357: testpoints = param? PARAM_NUMPOINTS : PLOTH_NUMPOINTS;
1358: }
1359: if (recur)
1360: nbpoints = testpoints << (param? PARAMR_MAXDEPTH : RECUR_MAXDEPTH);
1361: else
1362: nbpoints = testpoints;
1363:
1364: sig=gcmp(b,a); if (!sig) return 0;
1365: if (sig<0) { x=a; a=b; b=x; }
1366: dx=gdivgs(gsub(b,a),testpoints-1);
1367:
1368: x = cgetr(prec); gaffect(a,x); push_val(ep, x);
1369: av2=avma; p1=READ_EXPR(ch); tx=typ(p1);
1370: if (!is_matvec_t(tx))
1371: {
1372: xsml = gtodouble(a);
1373: xbig = gtodouble(b);
1374: ysml = ybig = gtodouble(p1);
1375: nc=1; nl=2; /* nc = nb of curves; nl = nb of coord. lists */
1376: if (param)
1377: err(warner,"flag PLOT_PARAMETRIC ignored");
1378: single_c=1; param=0;
1379: }
1380: else
1381: {
1382: if (tx != t_VEC) err(talker,"not a row vector in ploth");
1383: nl=lg(p1)-1; if (!nl) { avma=av; return 0; }
1384: single_c=0;
1385: if (param) nc=nl/2; else { nc=nl; nl++; }
1386: if (recur && nc > 1) err(talker,"multi-curves cannot be plot recursively");
1387:
1388: if (param)
1389: {
1390: xbig=xsml=gtodouble((GEN)p1[1]);
1391: ybig=ysml=gtodouble((GEN)p1[2]);
1392: for (i=3; i<=nl; i++)
1393: {
1394: fx=gtodouble((GEN)p1[i]); i++;
1395: fy=gtodouble((GEN)p1[i]);
1396: if (xsml>fx) xsml=fx; else if (xbig<fx) xbig=fx;
1397: if (ysml>fy) ysml=fy; else if (ybig<fy) ybig=fy;
1398: }
1399: }
1400: else
1401: {
1402: xsml=gtodouble(a); xbig=gtodouble(b);
1403: ysml=gtodouble(vecmin(p1)); ybig=gtodouble(vecmax(p1));
1404: }
1405: }
1406:
1407: pl=(dblPointList*) gpmalloc(nl*sizeof(dblPointList));
1408: for (i = 0; i < nl; i++)
1409: {
1410: pl[i].d = (double*) gpmalloc((nbpoints+1)*sizeof(dblPointList));
1411: pl[i].nb=0;
1412: }
1413: pl[0].xsml=xsml; pl[0].ysml=ysml; pl[0].xbig=xbig; pl[0].ybig=ybig;
1414:
1415: if (recur) /* recursive plot */
1416: {
1417: xleft=cgetr(3); xright=cgetr(3); yleft=cgetr(3); yright=cgetr(3);
1418: if (param)
1419: {
1420: tleft=cgetr(prec); tright=cgetr(prec);
1421: av2=avma;
1422: gaffect(a,tleft); ep->value = (void*)tleft; p1=READ_EXPR(ch);
1423: gaffect((GEN)p1[1],xleft); gaffect((GEN)p1[2],yleft);
1424: for (i=0; i<testpoints-1; i++)
1425: {
1426: if (i) {
1427: gaffect(tright,tleft); gaffect(xright,xleft); gaffect(yright,yleft);
1428: }
1429: gaddz(tleft,dx,tright); ep->value = (void*)tright;
1430: p1 = READ_EXPR(ch);
1431: if (lg(p1) != 3) err(talker,"inconsistent data in rectplothin");
1432: gaffect((GEN)p1[1],xright); gaffect((GEN)p1[2],yright);
1433:
1434: Appendx(&pl[0],&pl[0],rtodbl(xleft));
1435: Appendy(&pl[0],&pl[1],rtodbl(yleft));
1436:
1437: param_recursion(pl,ch,ep, tleft,xleft,yleft, tright,xright,yright, 0);
1438: avma=av2;
1439: }
1440: Appendx(&pl[0],&pl[0],rtodbl(xright));
1441: Appendy(&pl[0],&pl[1],rtodbl(yright));
1442: }
1443: else /* single_c */
1444: {
1445: av2=avma;
1446: gaffect(a,xleft); ep->value = (void*) xleft;
1447: gaffect(READ_EXPR(ch),yleft);
1448: for (i=0; i<testpoints-1; i++)
1449: {
1450: gaddz(xleft,dx,xright); ep->value = (void*) xright;
1451: gaffect(READ_EXPR(ch),yright);
1452:
1453: Appendx(&pl[0],&pl[0],rtodbl(xleft));
1454: Appendy(&pl[0],&pl[1],rtodbl(yleft));
1455:
1456: single_recursion(pl,ch,ep,xleft,yleft,xright,yright,0);
1457: avma=av2;
1458: gaffect(xright,xleft); gaffect(yright,yleft);
1459: }
1460: Appendx(&pl[0],&pl[0],rtodbl(xright));
1461: Appendy(&pl[0],&pl[1],rtodbl(yright));
1462: }
1463: }
1464: else /* non-recursive plot */
1465: {
1466: if (single_c)
1467: for (i=0; i<testpoints; i++)
1468: {
1469: p1 = READ_EXPR(ch);
1470: pl[0].d[i]=gtodouble(x);
1471: Appendy(&pl[0],&pl[1],gtodouble(p1));
1472: gaddz(x,dx,x); avma=av2;
1473: }
1474: else if (param)
1475: {
1476: long k;
1477: double z;
1478:
1479: for (i=0; i<testpoints; i++)
1480: {
1481: p1 = READ_EXPR(ch);
1482: if (lg(p1) != nl+1) err(talker,"inconsistent data in rectplothin");
1483: for (j=0; j<nl; j=k)
1484: {
1485: k=j+1; z=gtodouble((GEN)p1[k]);
1486: Appendx(&pl[0], &pl[j],z);
1487:
1488: j=k; k++; z=gtodouble((GEN)p1[k]);
1489: Appendy(&pl[0], &pl[j],z);
1490: }
1491: gaddz(x,dx,x); avma=av2;
1492: }
1493: }
1494: else /* plothmult */
1495: for (i=0; i<testpoints; i++)
1496: {
1497: p1 = READ_EXPR(ch);
1498: if (lg(p1) != nl) err(talker,"inconsistent data in rectplothin");
1499: pl[0].d[i]=gtodouble(x);
1500: for (j=1; j<nl; j++) { Appendy(&pl[0],&pl[j],gtodouble((GEN)p1[j])); }
1501: gaddz(x,dx,x); avma=av2;
1502: }
1503: }
1504: pl[0].nb=nc; pop_val(ep); avma = av;
1505: return pl;
1506: }
1507:
1508: extern GEN polint_i(GEN xa, GEN ya, GEN x, long n, GEN *ptdy);
1509:
1510: /* Uses highlevel plotting functions to implement splines as
1511: a low-level plotting function.
1512: Most probably one does not need to make varn==0 into pure variable,
1513: since plotting functions should take care of this. */
1514: static void
1515: rectsplines(long ne, double *x, double *y, long lx, long flag)
1516: {
1.2 ! noro 1517: long i, j;
! 1518: gpmem_t oldavma = avma;
1.1 noro 1519: GEN tas, xa = cgetg(lx+1, t_VEC), ya = cgetg(lx+1, t_VEC);
1520: entree *var0 = varentries[ordvar[0]];
1521:
1522: if (lx < 4)
1523: err(talker, "Too few points (%ld) for spline plot", lx);
1524: for (i = 1; i <= lx; i++) {
1525: xa[i] = (long) dbltor(x[i-1]);
1526: ya[i] = (long) dbltor(y[i-1]);
1527: }
1528: if (flag & PLOT_PARAMETRIC) {
1529: tas = new_chunk(4);
1530: for (j = 1; j <= 4; j++) tas[j-1] = lstoi(j);
1531: quark_gen = cgetg(2 + 1, t_VEC);
1532: }
1533: else tas = NULL; /* for lint */
1534: for (i = 0; i <= lx - 4; i++) {
1.2 ! noro 1535: gpmem_t oavma = avma;
1.1 noro 1536:
1537: xa++; ya++;
1538: if (flag & PLOT_PARAMETRIC) {
1539: quark_gen[1] = (long)polint_i(tas, xa, polx[0], 4, NULL);
1540: quark_gen[2] = (long)polint_i(tas, ya, polx[0], 4, NULL);
1541: } else {
1542: quark_gen = polint_i(xa, ya, polx[0], 4, NULL);
1543: tas = xa;
1544: }
1545: rectploth(ne, var0,
1546: (GEN)(i==0 ? tas[0] : tas[1]),
1547: (GEN)(i==lx-4 ? tas[3] : tas[2]),
1548: QUARK,
1549: DEFAULTPREC, /* XXXX precision */
1550: PLOT_RECURSIVE
1551: | PLOT_NO_RESCALE
1552: | PLOT_NO_FRAME
1553: | PLOT_NO_AXE_Y
1554: | PLOT_NO_AXE_X
1555: | (flag & PLOT_PARAMETRIC),
1556: 2); /* Start with 3 points */
1557: avma = oavma;
1558: }
1559: avma = oldavma;
1560: }
1561:
1562: /*
1563: * Plot a dblPointList. Complete with axes, bounding box, etc.
1564: * We use two drawing rectangles: one for strings, another
1565: * for graphs.
1566: *
1567: * data is an array of structs. Its meaning depends on flags :
1568: *
1569: * + data[0] contains global extremas, the number of curves to plot
1570: * (data[0].nb) and a list of doubles (first set of x-coordinates).
1571: *
1572: * + data[i].nb (i>0) contains the number of points in the list
1573: * data[i].d (hopefully, data[2i].nb=data[2i+1].nb when i>0...)
1574: *
1575: * + If flags contain PLOT_PARAMETRIC, the array length should be
1576: * even, and successive pairs (data[2i].d, data[2i+1].d) represent
1577: * curves to plot.
1578: *
1579: * + If there is no such flag, the first element is an array with
1580: * x-coordinates and the following ones contain y-coordinates.
1581: *
1582: * Additional flags: PLOT_NO_AXE_X, PLOT_NO_AXE_Y, PLOT_NO_FRAME.
1583: */
1584:
1585: static GEN
1586: rectplothrawin(long stringrect, long drawrect, dblPointList *data,
1587: long flags, PARI_plot *WW)
1588: {
1589: GEN res;
1590: dblPointList y,x;
1591: double xsml,xbig,ysml,ybig,tmp;
1.2 ! noro 1592: long ltype;
! 1593: gpmem_t ltop=avma;
1.1 noro 1594: long i,nc,nbpoints, w[2], wx[2], wy[2];
1595:
1596: w[0]=stringrect; w[1]=drawrect;
1597: if (!data) return cgetg(1,t_VEC);
1598: x = data[0]; nc = x.nb;
1599: xsml = x.xsml; xbig = x.xbig;
1600: ysml = x.ysml; ybig = x.ybig;
1601: if (xbig-xsml < 1.e-9)
1602: {
1603: tmp=fabs(xsml)/10; if (!tmp) tmp=0.1;
1604: xbig+=tmp; xsml-=tmp;
1605: }
1606: if (ybig-ysml < 1.e-9)
1607: {
1608: tmp=fabs(ysml)/10; if (!tmp) tmp=0.1;
1609: ybig+=tmp; ysml-=tmp;
1610: }
1611:
1612: if (WW)
1613: { /* no rectwindow supplied ==> ps or screen output */
1614: char c1[16],c2[16],c3[16],c4[16];
1615: PARI_plot W = *WW;
1616: long lm = W.fwidth*10; /* left margin */
1617: long rm = W.hunit-1; /* right margin */
1618: long tm = W.vunit-1; /* top margin */
1619: long bm = W.vunit+W.fheight-1; /* bottom margin */
1620: long is = W.width - (lm+rm);
1621: long js = W.height - (tm+bm);
1622:
1623: wx[0]=wy[0]=0; wx[1]=lm; wy[1]=tm;
1624: /* Window size (W.width x W.height) is given in pixels, and
1625: * correct pixels are 0..w_width-1.
1626: * On the other hand, rect functions work with windows whose pixel
1627: * range is [0,width]. */
1628: initrect(stringrect, W.width-1, W.height-1);
1629: if (drawrect != stringrect) initrect(drawrect, is-1, js-1);
1630:
1631: /* draw labels on stringrect */
1632: sprintf(c1,"%.5g",ybig); sprintf(c2,"%.5g",ysml);
1633: sprintf(c3,"%.5g",xsml); sprintf(c4,"%.5g",xbig);
1634:
1635: rectlinetype(stringrect,-2); /* Frame */
1636: current_color[stringrect]=BLACK;
1637: put_string( stringrect, lm, 0, c1,
1638: RoSTdirRIGHT | RoSTdirHGAP | RoSTdirTOP);
1639: put_string(stringrect, lm, W.height - bm, c2,
1640: RoSTdirRIGHT | RoSTdirHGAP | RoSTdirVGAP);
1641: put_string(stringrect, lm, W.height - bm, c3,
1642: RoSTdirLEFT | RoSTdirTOP);
1643: put_string(stringrect, W.width - rm - 1, W.height - bm, c4,
1644: RoSTdirRIGHT | RoSTdirTOP);
1645: }
1646: RHasGraph(check_rect(drawrect)) = 1;
1647:
1648: if (!(flags & PLOT_NO_RESCALE))
1649: rectscale0(drawrect, xsml, xbig, ysml, ybig);
1650:
1651: if (!(flags & PLOT_NO_FRAME))
1652: {
1653: int do_double = (flags & PLOT_NODOUBLETICK) ? TICKS_NODOUBLE : 0;
1654:
1655: rectlinetype(drawrect, -2); /* Frame. */
1656: current_color[drawrect]=BLACK;
1657: rectmove0(drawrect,xsml,ysml,0);
1658: rectbox0(drawrect,xbig,ybig,0);
1659: if (!(flags & PLOT_NO_TICK_X)) {
1660: rectticks(WW, drawrect, xsml, ysml, xbig, ysml, xsml, xbig,
1661: TICKS_CLOCKW | do_double);
1662: rectticks(WW, drawrect, xbig, ybig, xsml, ybig, xbig, xsml,
1663: TICKS_CLOCKW | do_double);
1664: }
1665: if (!(flags & PLOT_NO_TICK_Y)) {
1666: rectticks(WW, drawrect, xbig, ysml, xbig, ybig, ysml, ybig,
1667: TICKS_CLOCKW | do_double);
1668: rectticks(WW, drawrect, xsml, ybig, xsml, ysml, ybig, ysml,
1669: TICKS_CLOCKW | do_double);
1670: }
1671: }
1672:
1673: if (!(flags & PLOT_NO_AXE_Y) && (xsml<=0 && xbig >=0))
1674: {
1675: rectlinetype(drawrect, -1); /* Axes. */
1676: current_color[drawrect]=BLUE;
1677: rectmove0(drawrect,0.0,ysml,0);
1678: rectline0(drawrect,0.0,ybig,0);
1679: }
1680:
1681: if (!(flags & PLOT_NO_AXE_X) && (ysml<=0 && ybig >=0))
1682: {
1683: rectlinetype(drawrect, -1); /* Axes. */
1684: current_color[drawrect]=BLUE;
1685: rectmove0(drawrect,xsml,0.0,0);
1686: rectline0(drawrect,xbig,0.0,0);
1687: }
1688:
1689: i = (flags & PLOT_PARAMETRIC)? 0: 1;
1690: for (ltype = 0; ltype < nc; ltype++)
1691: {
1692: current_color[drawrect] = ltype&1 ? SIENNA: RED;
1693: if (flags & PLOT_PARAMETRIC) x = data[i++];
1694:
1695: y=data[i++]; nbpoints=y.nb;
1696: if ((flags & PLOT_POINTS_LINES) || (flags & PLOT_POINTS)) {
1697: rectlinetype(drawrect, rectpoint_itype + ltype); /* Graphs. */
1698: rectpointtype(drawrect, rectpoint_itype + ltype); /* Graphs. */
1699: rectpoints0(drawrect,x.d,y.d,nbpoints);
1700: }
1701: if ((flags & PLOT_POINTS_LINES) || !(flags & PLOT_POINTS)) {
1702: if (flags & PLOT_SPLINES) {
1703: /* rectsplines will call us back with ltype == 0 */
1704: int old = rectline_itype;
1705:
1706: rectline_itype = rectline_itype + ltype;
1707: rectsplines(drawrect,x.d,y.d,nbpoints,flags);
1708: rectline_itype = old;
1709: } else {
1710: rectlinetype(drawrect, rectline_itype + ltype); /* Graphs. */
1711: rectlines0(drawrect,x.d,y.d,nbpoints,0);
1712: }
1713: }
1714: }
1715: for (i--; i>=0; i--) free(data[i].d);
1716: free(data);
1717:
1718: if (WW)
1719: {
1720: if (flags & PLOT_POSTSCRIPT)
1721: postdraw0(w,wx,wy,2);
1722: else
1723: rectdraw0(w,wx,wy,2, 0);
1724:
1725: killrect(drawrect); if (stringrect != drawrect) killrect(stringrect);
1726: }
1727:
1728: avma=ltop;
1729: res = cgetg(5,t_VEC);
1730: res[1]=(long)dbltor(xsml); res[2]=(long)dbltor(xbig);
1731: res[3]=(long)dbltor(ysml); res[4]=(long)dbltor(ybig);
1732: return res;
1733: }
1734:
1735: /*************************************************************************/
1736: /* */
1737: /* HI-RES FUNCTIONS */
1738: /* */
1739: /*************************************************************************/
1740:
1741: GEN
1742: rectploth(long drawrect,entree *ep,GEN a,GEN b,char *ch,
1743: long prec,ulong flags,long testpoints)
1744: {
1745: dblPointList *pl=rectplothin(ep, a,b, ch, prec, flags,testpoints);
1746: return rectplothrawin(0,drawrect, pl, flags,NULL);
1747: }
1748:
1749: GEN
1750: rectplothraw(long drawrect, GEN data, long flags)
1751: {
1752: dblPointList *pl=gtodblList(data,flags);
1753: return rectplothrawin(0,drawrect,pl,flags,NULL);
1754: }
1755:
1756: static PARI_plot*
1757: init_output(long flags)
1758: {
1759: if (flags & PLOT_POSTSCRIPT)
1760: { PARI_get_psplot(); return &pari_psplot; }
1761: else
1762: { PARI_get_plot(0); return &pari_plot; }
1763: }
1764:
1765: static GEN
1766: ploth0(long stringrect,long drawrect,entree *ep,GEN a,GEN b,char *ch,
1767: long prec,ulong flags,long testpoints)
1768: {
1769: PARI_plot *output = init_output(flags);
1770: dblPointList *pl=rectplothin(ep, a,b, ch, prec, flags,testpoints);
1771: return rectplothrawin(stringrect,drawrect, pl, flags, output);
1772: }
1773:
1774: static GEN
1775: plothraw0(long stringrect, long drawrect, GEN listx, GEN listy, long flags)
1776: {
1777: PARI_plot *output = init_output(flags);
1.2 ! noro 1778: long data[] = {evaltyp(t_VEC) | _evallg(3), 0, 0};
1.1 noro 1779: dblPointList *pl;
1780:
1781: data[1] = (long) listx;
1782: data[2] = (long) listy;
1783: pl=gtodblList(data,PLOT_PARAMETRIC);
1784: if (!pl) return cgetg(1,t_VEC);
1785: return rectplothrawin(stringrect,drawrect,pl,flags | PLOT_PARAMETRIC,output);
1786: }
1787:
1788: GEN
1789: plothraw(GEN listx, GEN listy, long flags)
1790: {
1791: flags = (flags > 1 ? flags : (flags ? 0: PLOT_POINTS));
1792: return plothraw0(STRINGRECT, DRAWRECT, listx, listy, flags);
1793: }
1794:
1795: GEN
1796: ploth(entree *ep, GEN a, GEN b, char *ch, long prec,long flags,long numpoints)
1797: {
1798: return ploth0(STRINGRECT, DRAWRECT, ep,a,b,ch,prec,flags,numpoints);
1799: }
1800:
1801: GEN
1802: ploth2(entree *ep, GEN a, GEN b, char *ch, long prec)
1803: {
1804: return ploth0(STRINGRECT, DRAWRECT, ep,a,b,ch,prec,PLOT_PARAMETRIC,0);
1805: }
1806:
1807: GEN
1808: plothmult(entree *ep, GEN a, GEN b, char *ch, long prec)
1809: {
1810: return ploth0(STRINGRECT, DRAWRECT, ep,a,b,ch,prec,0,0);
1811: }
1812:
1813: GEN
1814: postplothraw(GEN listx, GEN listy, long flags)
1815: {
1816: flags=(flags)? 0: PLOT_POINTS;
1817: return plothraw0(STRINGRECT, DRAWRECT, listx, listy, flags|PLOT_POSTSCRIPT);
1818: }
1819:
1820: GEN
1821: postploth(entree *ep, GEN a, GEN b, char *ch, long prec,long flags,
1822: long numpoints)
1823: {
1824: return ploth0(STRINGRECT,DRAWRECT,ep,a,b,ch,prec,flags|PLOT_POSTSCRIPT,
1825: numpoints);
1826: }
1827:
1828: GEN
1829: postploth2(entree *ep, GEN a, GEN b, char *ch, long prec,
1830: long numpoints)
1831: {
1832: return ploth0(STRINGRECT,DRAWRECT,ep,a,b,ch,prec,
1833: PLOT_PARAMETRIC|PLOT_POSTSCRIPT,numpoints);
1834: }
1835:
1836: GEN
1.2 ! noro 1837: plothsizes(void)
1.1 noro 1838: {
1839: return plothsizes_flag(0);
1840: }
1841:
1842: GEN
1843: plothsizes_flag(long flag)
1844: {
1845: GEN vect = cgetg(1+6,t_VEC);
1846:
1847: PARI_get_plot(0);
1848: vect[1] = lstoi(w_width);
1849: vect[2] = lstoi(w_height);
1850: if (flag) {
1851: vect[3] = (long)dbltor(h_unit*1.0/w_width);
1852: vect[4] = (long)dbltor(v_unit*1.0/w_height);
1853: vect[5] = (long)dbltor(f_width*1.0/w_width);
1854: vect[6] = (long)dbltor(f_height*1.0/w_height);
1855: } else {
1856: vect[3] = lstoi(h_unit);
1857: vect[4] = lstoi(v_unit);
1858: vect[5] = lstoi(f_width);
1859: vect[6] = lstoi(f_height);
1860: }
1861: return vect;
1862: }
1863:
1864: /*************************************************************************/
1865: /* */
1866: /* POSTSCRIPT OUTPUT */
1867: /* */
1868: /*************************************************************************/
1869:
1870: typedef struct spoint {
1871: int x,y; } SPoint;
1872: typedef struct ssegment {
1873: int x1,y1,x2,y2; } SSegment;
1874: typedef struct srectangle {
1875: int x,y,width,height; } SRectangle;
1876:
1877: static void ps_point(FILE *psfile, int x, int y);
1878: static void ps_line(FILE *psfile, int x1, int y1, int x2, int y2);
1879: static void ps_rect(FILE *psfile, int x1, int y1, int x2, int y2);
1880: static void ps_string(FILE *psfile, int x, int y, char *c, long dir);
1881:
1882: #undef ISCR
1883: #undef JSCR
1884: #define ISCR 1120 /* 1400 en haute resolution */
1885: #define JSCR 800 /* 1120 en haute resolution */
1886:
1887: static void
1.2 ! noro 1888: PARI_get_psplot(void)
1.1 noro 1889: {
1890: pari_psplot.height = JSCR - 60;
1891: pari_psplot.width = ISCR - 40;
1892: pari_psplot.fheight = 15;
1893: pari_psplot.fwidth = 6;
1894: pari_psplot.hunit = 5;
1895: pari_psplot.vunit = 5;
1896: }
1897:
1898: static void
1899: gendraw(GEN list, long ps, long flag)
1900: {
1901: long i,n,ne,*w,*x,*y;
1902: GEN x0,y0,win;
1903:
1904: if (typ(list) != t_VEC) err(talker,"not a vector in rectdraw");
1905: n = lg(list)-1;
1906: if (n%3) err(talker,"incorrect number of components in rectdraw");
1907: n = n/3; if (!n) return;
1908: w = (long*)gpmalloc(n*sizeof(long));
1909: x = (long*)gpmalloc(n*sizeof(long));
1910: y = (long*)gpmalloc(n*sizeof(long));
1911: if (flag)
1912: PARI_get_plot(0);
1913: for (i=0; i<n; i++)
1914: {
1915: win=(GEN)list[3*i+1]; x0=(GEN)list[3*i+2]; y0=(GEN)list[3*i+3];
1916: if (typ(win)!=t_INT || (!flag && (typ(x0)!=t_INT || typ(y0)!= t_INT)))
1917: err(talker, "not an integer type in rectdraw");
1918: if (flag) {
1919: double xd = gtodouble(x0), yd = gtodouble(y0);
1920: long xi, yi;
1921:
1922: xi = w_width - 1; yi = w_height - 1;
1923: xi = (long)(xd*xi + 0.5);
1924: yi = (long)(yd*yi + 0.5);
1925: x[i] = xi; y[i] = yi;
1926: } else {
1927: x[i]=itos(x0); y[i]=itos(y0);
1928: }
1929: ne=itos(win); check_rect(ne);
1930: w[i]=ne;
1931: }
1932: if (ps) postdraw00(w,x,y,n,flag); else rectdraw0(w,x,y,n, 1);
1933: free(x); free(y); free(w);
1934: }
1935:
1936: void
1937: postdraw(GEN list) { gendraw(list, 1, 0); }
1938:
1939: void
1940: rectdraw(GEN list) { gendraw(list, 0, 0); }
1941:
1942: void
1943: postdraw_flag(GEN list, long flag) { gendraw(list, 1, flag); }
1944:
1945: void
1946: rectdraw_flag(GEN list, long flag) { gendraw(list, 0, flag); }
1947:
1948: static char*
1949: zmalloc(size_t x)
1950: {
1951: return x? gpmalloc(x): NULL;
1952: }
1953:
1954: void
1955: postdraw0(long *w, long *x, long *y, long lw)
1956: {
1957: postdraw00(w, x, y, lw, 0);
1958: }
1959:
1960: void
1961: postdraw00(long *w, long *x, long *y, long lw, long scale)
1962: {
1963: double *ptx,*pty;
1964: long *numpoints,*numtexts,*xtexts,*ytexts,*dirtexts;
1965: RectObj *p1;
1966: PariRect *e;
1967: long i,j,x0,y0;
1968: long a,b,c,d,nd[ROt_MAX+1];
1969: char **texts;
1970: FILE *psfile;
1971: double xscale = 0.65, yscale = 0.65;
1972: long fontsize = 16, xtick = 5, ytick = 5;
1973:
1974: SPoint *points, **lines, *SLine;
1975: SSegment *seg;
1976: SRectangle *rect, SRec;
1977:
1978: if (scale) {
1979: double postxsize, postysize;
1980:
1981: PARI_get_psplot();
1982: postxsize = pari_psplot.width;
1983: postysize = pari_psplot.height;
1984: PARI_get_plot(0);
1985: xscale *= postxsize * 1.0/w_width;
1986: fontsize = (long) (fontsize/(postxsize * 1.0/w_width));
1987: yscale *= postysize * 1.0/w_height;
1988: xtick = h_unit; ytick = v_unit;
1989: }
1990: psfile = fopen(current_psfile, "a");
1991: if (!psfile)
1992: err(openfiler,"postscript",current_psfile);
1993:
1994: for (i=0; i<=ROt_MAX; i++) nd[i]=0;
1995:
1996: for (i=0; i<lw; i++)
1997: {
1998: e = check_rect_init(w[i]); p1=RHead(e);
1999: while (p1)
2000: {
2001: if (RoType(p1) != ROt_MP) nd[RoType(p1)]++;
2002: else nd[ROt_PT]+=RoMPcnt(p1);
2003: p1=RoNext(p1);
2004: }
2005: }
2006: points=(SPoint*) zmalloc(nd[ROt_PT]*sizeof(SPoint));
2007: seg=(SSegment*) zmalloc(nd[ROt_LN]*sizeof(SSegment));
2008: rect=(SRectangle*) zmalloc(nd[ROt_BX]*sizeof(SRectangle));
2009: lines=(SPoint**) zmalloc(nd[ROt_ML]*sizeof(SPoint*));
2010: numpoints=(long*) zmalloc(nd[ROt_ML]*sizeof(long));
2011: texts=(char**) zmalloc(nd[ROt_ST]*sizeof(char*));
2012: numtexts=(long*) zmalloc(nd[ROt_ST]*sizeof(long));
2013: xtexts = (long*) zmalloc(nd[ROt_ST]*sizeof(long));
2014: ytexts = (long*) zmalloc(nd[ROt_ST]*sizeof(long));
2015: dirtexts = (long*) zmalloc(nd[ROt_ST]*sizeof(long));
2016: for (i=0; i<=ROt_MAX; i++) nd[i]=0;
2017:
2018: for (i=0; i<lw; i++)
2019: {
2020: e=rectgraph[w[i]]; p1=RHead(e); x0=x[i]; y0=y[i];
2021: while (p1)
2022: {
2023: switch(RoType(p1))
2024: {
2025: case ROt_PT:
2026: points[nd[ROt_PT]].x = DTOL(RoPTx(p1)+x0);
2027: points[nd[ROt_PT]].y = DTOL(RoPTy(p1)+y0);
2028: nd[ROt_PT]++; break;
2029: case ROt_LN:
2030: seg[nd[ROt_LN]].x1 = DTOL(RoLNx1(p1)+x0);
2031: seg[nd[ROt_LN]].y1 = DTOL(RoLNy1(p1)+y0);
2032: seg[nd[ROt_LN]].x2 = DTOL(RoLNx2(p1)+x0);
2033: seg[nd[ROt_LN]].y2 = DTOL(RoLNy2(p1)+y0);
2034: nd[ROt_LN]++; break;
2035: case ROt_BX:
2036: a=rect[nd[ROt_BX]].x = DTOL(RoBXx1(p1)+x0);
2037: b=rect[nd[ROt_BX]].y = DTOL(RoBXy1(p1)+y0);
2038: rect[nd[ROt_BX]].width = DTOL(RoBXx2(p1)+x0-a);
2039: rect[nd[ROt_BX]].height = DTOL(RoBXy2(p1)+y0-b);
2040: nd[ROt_BX]++; break;
2041: case ROt_MP:
2042: ptx = RoMPxs(p1); pty = RoMPys(p1);
2043: for (j=0; j<RoMPcnt(p1); j++)
2044: {
2045: points[nd[ROt_PT]+j].x = DTOL(ptx[j]+x0);
2046: points[nd[ROt_PT]+j].y = DTOL(pty[j]+y0);
2047: }
2048: nd[ROt_PT]+=RoMPcnt(p1); break;
2049: case ROt_ML:
2050: ptx=RoMLxs(p1); pty=RoMLys(p1);
2051: numpoints[nd[ROt_ML]]=RoMLcnt(p1);
2052: lines[nd[ROt_ML]]=(SPoint*)gpmalloc(RoMLcnt(p1)*sizeof(SPoint));
2053: for (j=0; j<RoMLcnt(p1); j++)
2054: {
2055: lines[nd[ROt_ML]][j].x = DTOL(ptx[j]+x0);
2056: lines[nd[ROt_ML]][j].y = DTOL(pty[j]+y0);
2057: }
2058: nd[ROt_ML]++; break;
2059: case ROt_ST:
2060: texts[nd[ROt_ST]] = RoSTs(p1);
2061: numtexts[nd[ROt_ST]] = RoSTl(p1);
2062: xtexts[nd[ROt_ST]] = DTOL(RoSTx(p1)+x0);
2063: ytexts[nd[ROt_ST]] = DTOL(RoSTy(p1)+y0);
2064: dirtexts[nd[ROt_ST]] = RoSTdir(p1);
2065: nd[ROt_ST]++; break;
2066: default: break;
2067: }
2068: p1=RoNext(p1);
2069: }
2070: }
2071: /* Definitions taken from post terminal of Gnuplot. */
2072: fprintf(psfile,"%%!\n50 50 translate\n/Times-Roman findfont %ld scalefont setfont\n%g %g scale\n", fontsize, yscale,xscale);
2073: fprintf(psfile,"/Lshow { moveto 90 rotate show -90 rotate } def\n/Rshow { 3 -1 roll dup 4 1 roll stringwidth pop sub Lshow } def\n/Cshow { 3 -1 roll dup 4 1 roll stringwidth pop 2 div sub Lshow } def\n");
2074: fprintf(psfile,"/Xgap %ld def\n/Ygap %ld def\n", xtick, ytick);
2075: fprintf(psfile,"/Bbox { gsave newpath 0 0 moveto true charpath pathbbox grestore } def\n");
2076: fprintf(psfile,"/Height { Bbox 4 1 roll pop pop pop } def\n");
2077: fprintf(psfile,"/TopAt { 3 -1 roll dup 4 1 roll Height 3 -1 roll add exch } def\n");
2078: fprintf(psfile,"/VCenter { 3 -1 roll dup 4 1 roll Height 2 div 3 -1 roll add exch } def\n");
2079: fprintf(psfile,"/Tgap { exch Ygap add exch } def\n");
2080: fprintf(psfile,"/Bgap { exch Ygap sub exch } def\n");
2081: fprintf(psfile,"/Lgap { Xgap add } def\n");
2082: fprintf(psfile,"/Rgap { Xgap sub } def\n");
2083:
2084: for (i=0; i<nd[ROt_PT]; i++)
2085: ps_point(psfile,points[i].x,points[i].y);
2086: for (i=0; i<nd[ROt_LN]; i++)
2087: ps_line(psfile,seg[i].x1,seg[i].y1,seg[i].x2,seg[i].y2);
2088: for (i=0; i<nd[ROt_BX]; i++)
2089: {
2090: SRec=rect[i]; a=SRec.x; b=SRec.y; c=a+SRec.width;
2091: d=b+SRec.height; ps_rect(psfile,a,b,c,d);
2092: }
2093: for (i=0; i<nd[ROt_ML]; i++)
2094: {
2095: SLine=lines[i];
2096: for (j=0; j<numpoints[i]; j++)
2097: {
2098: if (!j) fprintf(psfile,"%d %d moveto\n",SLine[0].y,SLine[0].x);
2099: else fprintf(psfile,"%d %d lineto\n",SLine[j].y,SLine[j].x);
2100: }
2101: }
2102: for (i=0; i<nd[ROt_ST]; i++)
2103: ps_string(psfile,xtexts[i],ytexts[i],texts[i], dirtexts[i]);
2104: fprintf(psfile,"stroke showpage\n"); fclose(psfile);
2105: #define xfree(pointer) if (pointer) free(pointer)
2106: xfree(points); xfree(seg); xfree(rect); xfree(numpoints);
2107: for (i=0; i<nd[ROt_ML]; i++) xfree(lines[i]);
2108: xfree(lines); xfree(texts); xfree(numtexts); xfree(xtexts); xfree(ytexts);
2109: #undef xfree
2110: }
2111:
2112: static void
2113: ps_point(FILE *psfile, int x, int y)
2114: {
2115: fprintf(psfile,"%d %d moveto\n0 2 rlineto 2 0 rlineto 0 -2 rlineto closepath fill\n",y,x);
2116: }
2117:
2118: static void
2119: ps_line(FILE *psfile, int x1, int y1, int x2, int y2)
2120: {
2121: fprintf(psfile,"%d %d moveto\n%d %d lineto\n",y1,x1,y2,x2);
2122: }
2123:
2124: static void
2125: ps_rect(FILE *psfile, int x1, int y1, int x2, int y2)
2126: {
2127: fprintf(psfile,"%d %d moveto\n%d %d lineto\n%d %d lineto\n%d %d lineto\nclosepath\n",y1,x1,y1,x2,y2,x2,y2,x1);
2128: }
2129:
2130: static void
2131: ps_string(FILE *psfile, int x, int y, char *s, long dir)
2132: {
2133: if (strpbrk(s, "(\\)")) {
2134: fprintf(psfile,"(");
2135: while (*s) {
2136: if ( *s=='(' || *s==')' || *s=='\\' )
2137: fputc('\\', psfile);
2138: fputc(*s, psfile);
2139: s++;
2140: }
2141: } else
2142: fprintf(psfile,"(%s", s);
2143: fprintf(psfile,") %d %d %s%s%s%sshow\n",
2144: y, x,
2145: ((dir & RoSTdirVGAP)
2146: ? ((dir & RoSTdirVPOS_mask) == RoSTdirTOP ? "Tgap " : "Bgap ")
2147: : ""),
2148: ((dir & RoSTdirHGAP)
2149: ? ((dir & RoSTdirHPOS_mask) == RoSTdirRIGHT ? "Rgap " : "Lgap ")
2150: : ""),
2151: ((dir & RoSTdirVPOS_mask) == RoSTdirBOTTOM ? ""
2152: : (dir & RoSTdirVPOS_mask) == RoSTdirTOP ? "TopAt " : "VCenter "),
2153: ((dir & RoSTdirHPOS_mask) == RoSTdirLEFT ? "L"
2154: : ((dir & RoSTdirHPOS_mask) == RoSTdirRIGHT ? "R" : "C")));
2155: }
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to plotport.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / pari-2.2 / src / graph