Annotation of OpenXM_contrib2/asir2000/plot/calc.c, Revision 1.15
1.2 noro 1: /*
2: * Copyright (c) 1994-2000 FUJITSU LABORATORIES LIMITED
3: * All rights reserved.
4: *
5: * FUJITSU LABORATORIES LIMITED ("FLL") hereby grants you a limited,
6: * non-exclusive and royalty-free license to use, copy, modify and
7: * redistribute, solely for non-commercial and non-profit purposes, the
8: * computer program, "Risa/Asir" ("SOFTWARE"), subject to the terms and
9: * conditions of this Agreement. For the avoidance of doubt, you acquire
10: * only a limited right to use the SOFTWARE hereunder, and FLL or any
11: * third party developer retains all rights, including but not limited to
12: * copyrights, in and to the SOFTWARE.
13: *
14: * (1) FLL does not grant you a license in any way for commercial
15: * purposes. You may use the SOFTWARE only for non-commercial and
16: * non-profit purposes only, such as academic, research and internal
17: * business use.
18: * (2) The SOFTWARE is protected by the Copyright Law of Japan and
19: * international copyright treaties. If you make copies of the SOFTWARE,
20: * with or without modification, as permitted hereunder, you shall affix
21: * to all such copies of the SOFTWARE the above copyright notice.
22: * (3) An explicit reference to this SOFTWARE and its copyright owner
23: * shall be made on your publication or presentation in any form of the
24: * results obtained by use of the SOFTWARE.
25: * (4) In the event that you modify the SOFTWARE, you shall notify FLL by
1.3 noro 26: * e-mail at risa-admin@sec.flab.fujitsu.co.jp of the detailed specification
1.2 noro 27: * for such modification or the source code of the modified part of the
28: * SOFTWARE.
29: *
30: * THE SOFTWARE IS PROVIDED AS IS WITHOUT ANY WARRANTY OF ANY KIND. FLL
31: * MAKES ABSOLUTELY NO WARRANTIES, EXPRESSED, IMPLIED OR STATUTORY, AND
32: * EXPRESSLY DISCLAIMS ANY IMPLIED WARRANTY OF MERCHANTABILITY, FITNESS
33: * FOR A PARTICULAR PURPOSE OR NONINFRINGEMENT OF THIRD PARTIES'
34: * RIGHTS. NO FLL DEALER, AGENT, EMPLOYEES IS AUTHORIZED TO MAKE ANY
35: * MODIFICATIONS, EXTENSIONS, OR ADDITIONS TO THIS WARRANTY.
36: * UNDER NO CIRCUMSTANCES AND UNDER NO LEGAL THEORY, TORT, CONTRACT,
37: * OR OTHERWISE, SHALL FLL BE LIABLE TO YOU OR ANY OTHER PERSON FOR ANY
38: * DIRECT, INDIRECT, SPECIAL, INCIDENTAL, PUNITIVE OR CONSEQUENTIAL
39: * DAMAGES OF ANY CHARACTER, INCLUDING, WITHOUT LIMITATION, DAMAGES
40: * ARISING OUT OF OR RELATING TO THE SOFTWARE OR THIS AGREEMENT, DAMAGES
41: * FOR LOSS OF GOODWILL, WORK STOPPAGE, OR LOSS OF DATA, OR FOR ANY
42: * DAMAGES, EVEN IF FLL SHALL HAVE BEEN INFORMED OF THE POSSIBILITY OF
43: * SUCH DAMAGES, OR FOR ANY CLAIM BY ANY OTHER PARTY. EVEN IF A PART
44: * OF THE SOFTWARE HAS BEEN DEVELOPED BY A THIRD PARTY, THE THIRD PARTY
45: * DEVELOPER SHALL HAVE NO LIABILITY IN CONNECTION WITH THE USE,
46: * PERFORMANCE OR NON-PERFORMANCE OF THE SOFTWARE.
47: *
1.15 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/plot/calc.c,v 1.14 2018/03/29 01:32:55 noro Exp $
1.2 noro 49: */
1.1 noro 50: #include "ca.h"
1.5 noro 51: #include "parse.h"
1.1 noro 52: #include "ifplot.h"
53: #include <math.h>
1.7 ohara 54: #if defined(PARI)
1.1 noro 55: #include "genpari.h"
56: #endif
57:
1.6 noro 58: #ifndef MAXSHORT
59: #define MAXSHORT ((short)0x7fff)
60: #endif
61:
1.15 ! noro 62: void todouble(Obj,Obj *);
! 63: void Psetprec();
! 64:
1.9 saito 65: void calc(double **tab,struct canvas *can,int nox){
1.14 noro 66: //memory_plot,IFPLOTD,INEQND,INEQNANDD,INEQNORD
67: //INEQNXORD,conplotmainD
68: double x,y,xstep,ystep;
69: int ix,iy;
70: Real r,rx,ry;
71: Obj fr,g,t,s;
1.1 noro 72:
1.14 noro 73: if(!nox)initmarker(can,"Evaluating...");
1.12 noro 74: todouble((Obj)can->formula,&fr);
1.14 noro 75: xstep=(can->xmax-can->xmin)/can->width;
76: ystep=(can->ymax-can->ymin)/can->height;
77: MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
78: BDY(rx)=can->xmin;
79: substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&t); devalr(CO,t,&g);
80: BDY(ry)=can->ymin;
81: substr(CO,0,g,can->vy,can->ymin?(Obj)ry:0,&t); devalr(CO,t,&s);
82: can->vmax=can->vmin=ToReal(s);
83: for(ix=0,x=can->xmin; ix<can->width; ix++,x+=xstep){
84: BDY(rx)=x; substr(CO,0,fr,can->vx,x?(Obj)rx:0,&t);
85: devalr(CO,t,&g);
86: if(!nox)marker(can,DIR_X,ix);
87: for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
88: BDY(ry)=y;
89: substr(CO,0,g,can->vy,y?(Obj)ry:0,&t);
90: devalr(CO,t,&s);
91: tab[ix][iy]=ToReal(s);
92: if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
93: if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
94: }
95: }
1.9 saito 96: }
97:
98: void calcq(double **tab,struct canvas *can,int nox){
1.14 noro 99: //IFPLOTQ,INEQNQ,INEQNANDQ,INEQNORQ,INEQNXORQ
100: //plotoverD
101: Q dx,dy,xstep,ystep,q1,w,h,c;
102: P g,g1,f1,f2,x,y;
103: int ix,iy;
104: Obj fr,gm,t,s;
105: Real r,rx,ry;
1.9 saito 106:
1.12 noro 107: todouble((Obj)can->formula,&fr);
1.14 noro 108: MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
109: BDY(rx)=can->xmin;
110: substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&t); devalr(CO,t,&gm);
111: BDY(ry)=can->ymin;
112: substr(CO,0,gm,can->vy,can->ymin?(Obj)ry:0,&t); devalr(CO,t,&s);
113: can->vmax=can->vmin=ToReal(s);
114:
115: subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
116: subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
117: MKV(can->vx,x); mulp(CO,(P)xstep,x,(P *)&t);
118: addp(CO,(P)can->qxmin,(P)t,(P *)&s); substp(CO,can->formula,can->vx,(P)s,&f1);
119: MKV(can->vy,y); mulp(CO,(P)ystep,y,(P *)&t);
120: addp(CO,(P)can->qymin,(P)t,(P *)&s); substp(CO,f1,can->vy,(P)s,&f2);
121: ptozp(f2,1,&c,&g);
122: if(!nox) initmarker(can,"Evaluating...");
123: for(iy=0;iy<can->height;iy++){
124: marker(can,DIR_Y,iy);
125: STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,(P *)&t); ptozp((P)t,1,&c,&g1);
126: for(ix=0;ix<can->width;ix++){
127: STOQ(ix,q1);substp(CO,g1,can->vx,(P)q1,(P *)&t);
128: devalr(CO,t,&s);
129: tab[ix][iy]=ToReal(s);
130: if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
131: if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
132: }
133: }
1.9 saito 134: }
135:
136: void calcb(double **tab,struct canvas *can,int nox){
1.14 noro 137: //IFPLOTB,INEQNB,INEQNANDB,INEQNORB,INEQNXORB
138: Q dx,dy,xstep,ystep,q1,w,h,c;
139: P g,g1,f1,f2,x,y,t,s;
140: int ix,iy,*a,*pa;
141: VECT ss;
142: Obj fr,gm,tm,sm;
143: Real r,rx,ry;
1.9 saito 144:
1.12 noro 145: todouble((Obj)can->formula,&fr);
1.14 noro 146: MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
147: BDY(rx)=can->xmin;
148: substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&tm); devalr(CO,tm,&gm);
149: BDY(ry)=can->ymin;
150: substr(CO,0,gm,can->vy,can->ymin?(Obj)ry:0,&tm); devalr(CO,tm,&sm);
151: can->vmax=can->vmin=ToReal(sm);
152:
153: for(iy=0;iy<can->height;iy++)for(ix=0;ix<can->width;ix++)tab[ix][iy]=1.0;
154: subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
155: subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
156: MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
157: addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
158: MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
159: addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
160: ptozp(f2,1,&c,&g);
161: a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
162: for(iy=0;iy<can->height;iy++)for(ix=0;ix<can->width;ix++)tab[ix][iy]=1.0;
163: subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
164: subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
165: MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
166: addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
167: MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
168: addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
169: ptozp(f2,1,&c,&g);
170: a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
171: for(ix=0;ix<can->width;ix++){
172: STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
173: if(!g1)for(iy=0;iy<can->height;iy++)tab[ix][iy]=0.0;
174: else if(!NUM(g1)){
175: sturmseq(CO,g1,&ss);
176: seproot(ss,0,can->width,a);
177: for(iy=0,pa=a;iy<can->height;iy++,pa++){
178: if(*pa<0||(*(pa+1)>=0&&(*pa>*(pa+1))))tab[ix][iy]=0.0;
179: else {
180: STOQ(iy,q1);substp(CO,g1,can->vy,(P)q1,&t);
181: devalr(CO,(Obj)t,(Obj *)&s);
182: tab[ix][iy]=ToReal(s);
183: if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
184: if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
185: }
186: }
187: }
188: }
189: for(iy=0;iy<can->height;iy++){
190: STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
191: if(!g1) for(ix=0;ix<can->width;ix++)tab[ix][iy]=0.0;
192: else if(!NUM(g1)){
193: sturmseq(CO,g1,&ss);
194: seproot(ss,0,can->height,a);
195: for(ix=0,pa=a;ix<can->width;ix++,pa++){
196: if(tab[ix][iy]!=0.0){
197: if(*pa<0||(*(pa+1)>=0&&(*pa>*(pa+1))))tab[ix][iy]=0.0;
198: else {
199: STOQ(ix,q1);substp(CO,g1,can->vx,(P)q1,&t);
200: devalr(CO,(Obj)t,(Obj *)&s);
201: tab[ix][iy]=ToReal(s);
202: if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
203: if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
204: }
205: }
206: }
207: }
208: }
1.1 noro 209: }
210:
1.9 saito 211: double usubstrp(P p,double r){
1.14 noro 212: DCP dc;
213: int d;
214: double t,pwrreal0();
215:
216: if(!p) t=0.0;
217: else if(NUM(p))t=BDY((Real)p);
218: else {
219: dc=DC(p); t=BDY((Real)COEF(dc));
220: for(d=QTOS(DEG(dc)),dc=NEXT(dc);dc;d=QTOS(DEG(dc)),dc=NEXT(dc)){
221: t=t*pwrreal0(r,(d-QTOS(DEG(dc))))+BDY((Real)COEF(dc));
222: }
223: if(d)t*=pwrreal0(r,d);
224: }
225: return t;
1.1 noro 226: }
227:
1.9 saito 228: void qcalc(char **tab,struct canvas *can){
1.14 noro 229: //qifplotmain(Old type)
230: Q dx,dy,w,h,xstep,ystep,c,q1;
231: P g,g1,f1,f2,x,y,t,s;
232: int ix,iy;
233: int *a,*pa;
234: VECT ss;
235:
236: subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
237: subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
238: MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
239: addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
240: MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
241: addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
242: ptozp(f2,1,&c,&g);
243: a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
244: initmarker(can,"Horizontal scan...");
245: for( ix=0; ix < can->width; ix++ ){
246: marker(can,DIR_X,ix);
247: STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
248: if( !g1 )
249: for(iy=0; iy < can->height; iy++ )
250: tab[ix][iy]=1;
251: else if( !NUM(g1) ){
252: sturmseq(CO,g1,&ss); seproot(ss,0,can->height,a);
253: for(iy=0, pa=a; iy < can->height; iy++, pa++ )
254: if( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
255: tab[ix][iy]=1;
256: }
257: }
258: initmarker(can,"Vertical scan...");
259: for( iy=0; iy < can->height; iy++ ){
260: marker(can,DIR_Y,iy);
261: STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
262: if( !g1 )
263: for(ix=0; ix < can->width; ix++ )
264: tab[ix][iy]=1;
265: else if( !NUM(g1) ){
266: sturmseq(CO,g1,&ss); seproot(ss,0,can->width,a);
267: for(ix=0, pa=a; ix < can->width; ix++, pa++ )
268: if( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
269: tab[ix][iy]=1;
270: }
271: }
1.1 noro 272: }
273:
1.9 saito 274: void sturmseq(VL vl,P p,VECT *rp){
1.14 noro 275: P g1,g2,q,r,s,*t;
276: V v;
277: VECT ret;
278: int i,j;
279: Q a,b,c,d,h,l,m,x;
280:
281: v=VR(p);t=(P *)ALLOCA((deg(v,p)+1)*sizeof(P));
282: g1=t[0]=p;diffp(vl,p,v,(P *)&a);ptozp((P)a,1,&c,&g2);t[1]=g2;
283: for(i=1,h=ONE,x=ONE;;){
284: if(NUM(g2)) break;
285: subq(DEG(DC(g1)),DEG(DC(g2)),&d);
286: l=(Q)LC(g2);
287: if(SGN(l)<0){
288: chsgnq(l,&a);l=a;
289: }
290: addq(d,ONE,&a);pwrq(l,a,&b);mulp(vl,(P)b,g1,(P *)&a);
291: divsrp(vl,(P)a,g2,&q,&r);
292: if(!r) break;
293: chsgnp(r,&s);r=s;i++;
294: if(NUM(r)){
295: t[i]=r;break;
296: }
297: pwrq(h,d,&m);g1=g2;
298: mulq(m,x,&a);divsp(vl,r,(P)a,&g2);t[i]=g2;
299: x=(Q)LC(g1);
300: if(SGN(x)<0){
301: chsgnq(x,&a);x=a;
302: }
303: pwrq(x,d,&a);mulq(a,h,&b);divq(b,m,&h);
304: }
305: MKVECT(ret,i+1);
306: for(j=0;j<=i;j++)
307: ret->body[j]=(pointer)t[j];
308: *rp=ret;
1.1 noro 309: }
310:
1.9 saito 311: void seproot(VECT s,int min,int max,int *ar){
1.14 noro 312: P f,*ss;
313: Q q,t;
314: int i,j,k;
315:
316: ss=(P *)s->body;f=ss[0];
317: for(i=min;i<=max;i++){
318: STOQ(i,q);usubstqp(f,q,&t);
319: if(!t)ar[i]=-1;
320: else {
321: ar[i]=numch(s,q,t);break;
322: }
323: }
324: if(i>max) return;
325: for(j=max;j>= min;j--){
326: STOQ(j,q); usubstqp(f,q,&t);
327: if(!t)ar[j]=-1;
328: else {
329: if(i!=j)ar[j]=numch(s,q,t);
330: break;
331: }
332: }
333: if(j<=i+1) return;
334: if(ar[i]==ar[j]){
335: for(k=i+1;k<j;k++)ar[k]=ar[i];
336: return;
337: }
338: k=(i+j)/2;
339: seproot(s,i,k,ar);
340: seproot(s,k,j,ar);
1.1 noro 341: }
342:
1.9 saito 343: int numch(VECT s,Q n,Q a0){
1.14 noro 344: int len,i,c;
345: Q a;
346: P *ss;
347:
348: len=s->len;ss=(P *)s->body;
349: for(i=1,c=0;i<len;i++){
350: usubstqp(ss[i],n,&a);
351: if(a){
352: if((SGN(a)>0 && SGN(a0)<0)||(SGN(a)<0&&SGN(a0)>0))c++;
353: a0=a;
354: }
355: }
356: return c;
1.1 noro 357: }
358:
1.9 saito 359: void usubstqp(P p,Q r,Q *v){
1.14 noro 360: Q d,d1,a,b,t;
361: DCP dc;
1.1 noro 362:
1.14 noro 363: if(!p)
364: *v=0;
365: else if(NUM(p))*v=(Q)p;
366: else {
367: dc=DC(p);t=(Q)COEF(dc);
368: for(d=DEG(dc),dc=NEXT(dc);dc;d=DEG(dc),dc=NEXT(dc)){
369: subq(d,DEG(dc),&d1);pwrq(r,d1,&a);
370: mulq(t,a,&b);addq(b,(Q)COEF(dc),&t);
371: }
372: if(d){
373: pwrq(r,d,&a);mulq(t,a,&b);t=b;
374: }
375: *v=t;
376: }
1.1 noro 377: }
378:
1.13 noro 379: Num tobf(Num,int);
380:
381: void plotcalcbf(struct canvas *can){
382: Obj fr,s,t;
383: Num xmin,xmax,ymin,ymax,xstep;
384: Num u,v,ha,dx,dy,x;
385: Num *tab;
386: Real r;
387: Q w,h1;
388: int ix;
389: POINT *pa;
390: double rr;
391: Q prec;
392: NODE arg;
393:
394: STOQ(can->prec,prec); arg = mknode(1,prec); Psetprec(arg,&t);
395: evalr(CO,(Obj)can->formula,can->prec,&fr);
396: MKReal(can->xmin,r); xmin = tobf((Num)r,can->prec);
397: MKReal(can->xmax,r); xmax = tobf((Num)r,can->prec);
398: MKReal(can->ymin,r); ymin = tobf((Num)r,can->prec);
399: MKReal(can->ymax,r); ymax = tobf((Num)r,can->prec);
400: STOQ(can->width,w);
401: subbf(xmax,xmin,&dx); divbf(dx,(Num)w,&xstep);
1.14 noro 402: tab=(Num *)MALLOC(can->width*sizeof(Num));
403: for(ix=0,x=xmin;ix<can->width;ix++){
404: substr(CO,0,fr,can->vx,(Obj)x,(Obj *)&s);
405: evalr(CO,(Obj)s,can->prec,&t);
406: if(t&&(OID(t)!=O_N))
407: error("plotcalcbf : invalid evaluation");
408: tab[ix]=(Num)t;
1.13 noro 409: addbf(x,xstep,&u); x = u;
1.14 noro 410: }
411: if(!cmpbf(ymax,ymin)){
412: for(ymax=ymin=tab[0],ix=1;ix<can->width;ix++){
413: if(cmpbf(tab[ix],ymax)>0)ymax=tab[ix];
414: if(cmpbf(tab[ix],ymin)<0)ymin=tab[ix];
415: }
416: can->ymax=ToReal(ymax);can->ymin=ToReal(ymin);
417: }
1.13 noro 418: subbf(ymax,ymin,&dy);
1.14 noro 419: can->pa=(struct pa *)MALLOC(sizeof(struct pa));
420: can->pa[0].length=can->width;
421: can->pa[0].pos=pa=(POINT *)MALLOC(can->width*sizeof(POINT));
1.13 noro 422: STOQ(can->height-1,h1);
1.14 noro 423: for(ix=0;ix<can->width;ix++){
424: XC(pa[ix])=ix;
1.13 noro 425: subbf(ymax,tab[ix],&u); divbf(u,dy,&v); mulbf(v,(Num)h1,&u);
426: rr = ToReal(u);
1.14 noro 427: if(rr>MAXSHORT)YC(pa[ix])=MAXSHORT;
428: else if(rr<-MAXSHORT)YC(pa[ix])=-MAXSHORT;
429: else YC(pa[ix])=(long)rr;
430: }
1.13 noro 431: }
432:
1.9 saito 433: void plotcalc(struct canvas *can){
1.14 noro 434: //plot,memory_plot,plotover,plot_resize
435: double x,xmin,xstep,ymax,ymin,dy,*tab,usubstrp();
436: int ix,w,h;
437: Real r,rx;
438: Obj fr,t,s;
439: POINT *pa;
1.1 noro 440:
1.13 noro 441: if ( can->prec ) {
442: plotcalcbf(can);
443: return;
444: }
1.12 noro 445: todouble((Obj)can->formula,&fr);
1.14 noro 446: w=can->width;h=can->height;
447: xmin=can->xmin;xstep=(can->xmax-can->xmin)/w;
448: tab=(double *)ALLOCA(w*sizeof(double));
449: MKReal(1,rx); // dummy real number
450: for(ix=0,x=xmin;ix<w;ix++,x+=xstep){
451: // full substitution
452: BDY(rx)=x;
453: substr(CO,0,fr,can->vx,x?(Obj)rx:0,&s);
454: devalr(CO,(Obj)s,&t);
455: if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
456: error("plotcalc : invalid evaluation");
457: tab[ix]=ToReal((Num)t);
458: }
459: if(can->ymax==can->ymin){
460: for(ymax=ymin=tab[0],ix=1;ix<w;ix++){
461: if(tab[ix]>ymax)ymax=tab[ix];
462: if(tab[ix]<ymin)ymin=tab[ix];
463: }
464: can->ymax=ymax;can->ymin=ymin;
465: } else {
466: ymax=can->ymax;ymin=can->ymin;
467: }
468: dy=ymax-ymin;
469: can->pa=(struct pa *)MALLOC(sizeof(struct pa));
470: can->pa[0].length=w;
471: can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
472: for(ix=0;ix<w;ix++){
473: double t;
474: XC(pa[ix])=ix;
475: t=(h-1)*(ymax-tab[ix])/dy;
476: if(t>MAXSHORT)YC(pa[ix])=MAXSHORT;
477: else if(t<-MAXSHORT)YC(pa[ix])=-MAXSHORT;
478: else YC(pa[ix])=(long)t;
479: }
1.6 noro 480: }
481:
1.10 saito 482: void polarcalc(struct canvas *can){
1.14 noro 483: double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep,tr,p,*tabx,*taby;
484: double usubstrp();
485: int i,nstep,w,h;
486: POINT *pa;
487: Real r;
488: Obj fr,t,s;
1.11 saito 489:
1.12 noro 490: todouble((Obj)can->formula,&fr);
1.14 noro 491: w=can->width; h=can->height; nstep=can->nzstep;
492: pmin=can->zmin; pstep=(can->zmax-can->zmin)/nstep;
493: tabx=(double *)ALLOCA(nstep*sizeof(double));
494: taby=(double *)ALLOCA(nstep*sizeof(double));
495: MKReal(1,r); // dummy real number
496:
497: for(i=0,p=pmin;i<nstep;i++,p+= pstep){
498: // full substitution
499: BDY(r)=p;
500: substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
501: devalr(CO,(Obj)s,&t);
502: if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
503: error("polarcalc : invalid evaluation");
504: tr=ToReal((Num)t);
505: tabx[i]=tr*cos(p);
506: taby[i]=tr*sin(p);
507: }
508: xmax=xmin=tabx[0];
509: ymax=ymin=taby[0];
510: for(i=1;i<nstep;i++){
511: if(tabx[i]>xmax)xmax=tabx[i];
512: if(tabx[i]<xmin)xmin=tabx[i];
513: if(taby[i]>ymax)ymax=taby[i];
514: if(taby[i]<ymin)ymin=taby[i];
515: }
516: can->xmax=xmax;can->xmin=xmin;
517: can->ymax=ymax;can->ymin=ymin;
518: dx=xmax-xmin;
519: dy=ymax-ymin;
520: can->pa=(struct pa *)MALLOC(sizeof(struct pa));
521: can->pa[0].length=nstep;
522: can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
523: for(i=0;i<nstep;i++){
524: XC(pa[i])=(w-1)*(tabx[i]-xmin)/dx;
525: YC(pa[i])=(h-1)*(ymax-taby[i])/dy;
526: }
1.11 saito 527: }
528:
529: void polarcalcNG(struct canvas *can){
1.14 noro 530: //polarplotNG
531: double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep,tr,p, *tabx,*taby;
532: double usubstrp();
533: int i,ix,iy,nstep,w,h;
534: POINT *pa;
535: Real r;
536: Obj fr,t,s;
1.6 noro 537:
1.12 noro 538: todouble((Obj)can->formula,&fr);
1.14 noro 539: w=can->width; h=can->height; nstep=can->nzstep;
540: pmin=can->zmin; pstep=(can->zmax-can->zmin)/nstep;
541: tabx=(double *)ALLOCA(nstep*sizeof(double));
542: taby=(double *)ALLOCA(nstep*sizeof(double));
543: MKReal(1,r); // dummy real number
544:
545: for(i=0,p=pmin;i<nstep;i++,p+= pstep){
546: // full substitution
547: BDY(r)=p;
548: substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
549: devalr(CO,(Obj)s,&t);
550: if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
551: error("polarcalc : invalid evaluation");
552: tr=ToReal((Num)t);
553: tabx[i]=tr*cos(p);
554: taby[i]=tr*sin(p);
555: if(i==0){
556: xmax=xmin=tabx[0];
557: ymax=ymin=taby[0];
558: } else {
559: if(tabx[i]>xmax)xmax=tabx[i];
560: if(tabx[i]<xmin)xmin=tabx[i];
561: if(taby[i]>ymax)ymax=taby[i];
562: if(taby[i]<ymin)ymin=taby[i];
563: }
564: }
565: can->xmax=xmax;can->xmin=xmin;
566: can->ymax=ymax;can->ymin=ymin;
567: dx=xmax-xmin;
568: dy=ymax-ymin;
569: can->pa=(struct pa *)MALLOC(sizeof(struct pa));
570: can->pa[0].length=nstep;
571: can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
572: for(i=0;i<nstep;i++){
573: XC(pa[i])=(w-1)*(tabx[i]-xmin)/dx;
574: YC(pa[i])=(h-1)*(ymax-taby[i])/dy;
575: }
1.1 noro 576: }
1.8 saito 577:
1.9 saito 578: /*
579: void ineqncalc(double **tab,struct canvas *can,int nox){
1.14 noro 580: double x,y,xmin,ymin,xstep,ystep;
581: int ix,iy,w,h;
582: Real r,rx,ry;
583: Obj fr,g,t,s;
584: V vx,vy;
1.8 saito 585:
1.14 noro 586: if(!nox) initmarker(can,"Evaluating...");
1.12 noro 587: todouble((Obj)can->formula,&fr);
1.14 noro 588: vx=can->vx;vy=can->vy;
589: w=can->width;h=can->height;
590: xmin=can->xmin;xstep=(can->xmax-can->xmin)/w;
591: ymin=can->ymin;ystep=(can->ymin-can->ymin)/h;
592: MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
593:
594: for(ix=0,x=xmin;ix<=w;ix++,x+=xstep){
595: BDY(rx)=x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
596: devalr(CO,t,&g);
597: if(!nox) marker(can,DIR_X,ix);
598: for(iy=0,y=ymin;iy<=h;iy++,y+=ystep){
599: BDY(ry)=y;
600: substr(CO,0,g,vy,y?(Obj)ry:0,&t);
601: devalr(CO,t,&s);
602: tab[ix][iy]=ToReal(s);
603: }
604: }
1.8 saito 605: }
1.9 saito 606: */
1.8 saito 607:
1.9 saito 608: #if defined(INTERVAL)
609: void itvcalc(double **mask, struct canvas *can, int nox){
1.14 noro 610: //ITVIFPLOT
611: double x,y,xstep,ystep,dx,dy,wx,wy;
612: int idv,ix,iy,idx,idy;
613: Itv ity,itx,ddx,ddy;
614: Real r,rx,ry,rx1,ry1,rdx,rdy,rdx1,rdy1;
615: V vx,vy;
616: Obj fr,g,t,s;
1.8 saito 617:
1.14 noro 618: idv=can->division;
1.12 noro 619: todouble((Obj)can->formula,&fr);
1.14 noro 620: vx=can->vx; vy=can->vy;
621: xstep=(can->xmax-can->xmin)/can->width;
622: ystep=(can->ymax-can->ymin)/can->height;
623: if(idv!=0){
624: wx=xstep/can->division;
625: wy=ystep/can->division;
626: }
627: MKReal(can->ymin,ry1);
628: for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
629: ry=ry1;
630: MKReal(y+ystep,ry1);
631: istoitv((Num)(ry1),(Num)ry,&ity);
632: substr(CO,0,(Obj)fr,vy,(Obj)ity,&t);
633: MKReal(can->xmin,rx1);
634: for(ix=0,x=can->xmin; ix<can->width; ix++,x+=xstep){
635: rx=rx1;
636: MKReal(x+xstep,rx1);
637: istoitv((Num)(rx1),(Num)rx,&itx);
638: substr(CO,0,(Obj)fr,vx,(Obj)itx,&t);
639: MKReal(can->ymin,ry1);
640: for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
641: ry=ry1;
642: MKReal(y+ystep,ry1);
643: istoitv((Num)ry,(Num)ry1,&ity);
644: substr(CO,0,(Obj)t,vy,(Obj)ity,&g);
645: if(compnum(0,0,(Num)g))mask[ix][iy]=-1;
646: else {
647: mask[ix][iy]=0;
1.9 saito 648: /*
1.14 noro 649: if(idv==0) mask[ix][iy]=0;
650: else {
651: MKReal(y,rdy1);
652: for(idy=0,dy=y;idy<idv;dy+=wy,idy++){
653: rdy=rdy1;
654: MKReal(dy+wy,rdy1);
655: istoitv((Num)rdy,(Num)rdy1,&ddy);
656: substr(CO,0,(Obj)fr,vy,(Obj)ddy,&t);
657: MKReal(x,rdx1);
658: for(idx=0,dx=x;idx<idx;dx+=wx,idx++){
659: rdx=rdx1;
660: MKReal(dx+wx,rdx1);
661: istoitv((Num)rdx,(Num)rdx1,&ddx);
662: substr(CO,0,(Obj)t,vx,(Obj)ddx,&g);
663: if(!compnum(0,0,(Num)g)){
664: mask[ix][iy]=0;
665: break;
666: }
667: }
668: if(mask[ix][iy]==0)break;
669: }
670: }
1.9 saito 671: */
1.14 noro 672: }
673: }
674: }
675: }
1.8 saito 676: }
677: #endif
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to calc.c CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib2 / asir2000 / plot