Annotation of OpenXM_contrib2/asir2000/plot/calc.c, Revision 1.13
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.13 ! noro 48: * $OpenXM: OpenXM_contrib2/asir2000/plot/calc.c,v 1.12 2017/08/31 04:21:48 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.9 saito 62: void calc(double **tab,struct canvas *can,int nox){
1.10 saito 63: //memory_plot,IFPLOTD,INEQND,INEQNANDD,INEQNORD
64: //INEQNXORD,conplotmainD
1.9 saito 65: double x,y,xstep,ystep;
1.1 noro 66: int ix,iy;
67: Real r,rx,ry;
1.9 saito 68: Obj fr,g,t,s;
1.1 noro 69:
1.9 saito 70: if(!nox)initmarker(can,"Evaluating...");
1.12 noro 71: todouble((Obj)can->formula,&fr);
1.9 saito 72: xstep=(can->xmax-can->xmin)/can->width;
73: ystep=(can->ymax-can->ymin)/can->height;
74: MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
75: BDY(rx)=can->xmin;
76: substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&t); devalr(CO,t,&g);
77: BDY(ry)=can->ymin;
78: substr(CO,0,g,can->vy,can->ymin?(Obj)ry:0,&t); devalr(CO,t,&s);
79: can->vmax=can->vmin=ToReal(s);
80: for(ix=0,x=can->xmin; ix<can->width; ix++,x+=xstep){
81: BDY(rx)=x; substr(CO,0,fr,can->vx,x?(Obj)rx:0,&t);
1.5 noro 82: devalr(CO,t,&g);
1.9 saito 83: if(!nox)marker(can,DIR_X,ix);
84: for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
85: BDY(ry)=y;
86: substr(CO,0,g,can->vy,y?(Obj)ry:0,&t);
87: devalr(CO,t,&s);
88: tab[ix][iy]=ToReal(s);
89: if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
90: if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
91: }
92: }
93: }
94:
95: void calcq(double **tab,struct canvas *can,int nox){
1.10 saito 96: //IFPLOTQ,INEQNQ,INEQNANDQ,INEQNORQ,INEQNXORQ
1.9 saito 97: //plotoverD
98: Q dx,dy,xstep,ystep,q1,w,h,c;
99: P g,g1,f1,f2,x,y;
100: int ix,iy;
101: Obj fr,gm,t,s;
102: Real r,rx,ry;
103:
1.12 noro 104: todouble((Obj)can->formula,&fr);
1.9 saito 105: MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
106: BDY(rx)=can->xmin;
107: substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&t); devalr(CO,t,&gm);
108: BDY(ry)=can->ymin;
109: substr(CO,0,gm,can->vy,can->ymin?(Obj)ry:0,&t); devalr(CO,t,&s);
110: can->vmax=can->vmin=ToReal(s);
111:
112: subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
113: subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
114: MKV(can->vx,x); mulp(CO,(P)xstep,x,(P *)&t);
115: addp(CO,(P)can->qxmin,(P)t,(P *)&s); substp(CO,can->formula,can->vx,(P)s,&f1);
116: MKV(can->vy,y); mulp(CO,(P)ystep,y,(P *)&t);
117: addp(CO,(P)can->qymin,(P)t,(P *)&s); substp(CO,f1,can->vy,(P)s,&f2);
118: ptozp(f2,1,&c,&g);
119: if(!nox) initmarker(can,"Evaluating...");
120: for(iy=0;iy<can->height;iy++){
121: marker(can,DIR_Y,iy);
122: STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,(P *)&t); ptozp((P)t,1,&c,&g1);
123: for(ix=0;ix<can->width;ix++){
124: STOQ(ix,q1);substp(CO,g1,can->vx,(P)q1,(P *)&t);
1.5 noro 125: devalr(CO,t,&s);
1.9 saito 126: tab[ix][iy]=ToReal(s);
127: if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
128: if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
129: }
130: }
131: }
132:
133: void calcb(double **tab,struct canvas *can,int nox){
1.10 saito 134: //IFPLOTB,INEQNB,INEQNANDB,INEQNORB,INEQNXORB
1.9 saito 135: Q dx,dy,xstep,ystep,q1,w,h,c;
136: P g,g1,f1,f2,x,y,t,s;
137: int ix,iy,*a,*pa;
138: VECT ss;
139: Obj fr,gm,tm,sm;
140: Real r,rx,ry;
141:
1.12 noro 142: todouble((Obj)can->formula,&fr);
1.9 saito 143: MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
144: BDY(rx)=can->xmin;
145: substr(CO,0,fr,can->vx,can->xmin?(Obj)rx:0,&tm); devalr(CO,tm,&gm);
146: BDY(ry)=can->ymin;
147: substr(CO,0,gm,can->vy,can->ymin?(Obj)ry:0,&tm); devalr(CO,tm,&sm);
148: can->vmax=can->vmin=ToReal(sm);
149:
150: for(iy=0;iy<can->height;iy++)for(ix=0;ix<can->width;ix++)tab[ix][iy]=1.0;
151: subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
152: subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
153: MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
154: addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
155: MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
156: addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
157: ptozp(f2,1,&c,&g);
158: a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
159: for(iy=0;iy<can->height;iy++)for(ix=0;ix<can->width;ix++)tab[ix][iy]=1.0;
160: subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
161: subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
162: MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
163: addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
164: MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
165: addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
166: ptozp(f2,1,&c,&g);
167: a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
168: for(ix=0;ix<can->width;ix++){
169: STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
170: if(!g1)for(iy=0;iy<can->height;iy++)tab[ix][iy]=0.0;
171: else if(!NUM(g1)){
172: sturmseq(CO,g1,&ss);
173: seproot(ss,0,can->width,a);
174: for(iy=0,pa=a;iy<can->height;iy++,pa++){
175: if(*pa<0||(*(pa+1)>=0&&(*pa>*(pa+1))))tab[ix][iy]=0.0;
176: else {
177: STOQ(iy,q1);substp(CO,g1,can->vy,(P)q1,&t);
178: devalr(CO,(Obj)t,(Obj *)&s);
179: tab[ix][iy]=ToReal(s);
180: if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
181: if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
182: }
183: }
184: }
185: }
186: for(iy=0;iy<can->height;iy++){
187: STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
188: if(!g1) for(ix=0;ix<can->width;ix++)tab[ix][iy]=0.0;
189: else if(!NUM(g1)){
190: sturmseq(CO,g1,&ss);
191: seproot(ss,0,can->height,a);
192: for(ix=0,pa=a;ix<can->width;ix++,pa++){
193: if(tab[ix][iy]!=0.0){
194: if(*pa<0||(*(pa+1)>=0&&(*pa>*(pa+1))))tab[ix][iy]=0.0;
195: else {
196: STOQ(ix,q1);substp(CO,g1,can->vx,(P)q1,&t);
197: devalr(CO,(Obj)t,(Obj *)&s);
198: tab[ix][iy]=ToReal(s);
199: if(can->vmax<tab[ix][iy])can->vmax=tab[ix][iy];
200: if(can->vmin>tab[ix][iy])can->vmin=tab[ix][iy];
201: }
202: }
203: }
1.1 noro 204: }
205: }
206: }
207:
1.9 saito 208: double usubstrp(P p,double r){
1.1 noro 209: DCP dc;
210: int d;
1.9 saito 211: double t,pwrreal0();
1.1 noro 212:
1.9 saito 213: if(!p) t=0.0;
214: else if(NUM(p))t=BDY((Real)p);
1.1 noro 215: else {
1.9 saito 216: dc=DC(p); t=BDY((Real)COEF(dc));
217: for(d=QTOS(DEG(dc)),dc=NEXT(dc);dc;d=QTOS(DEG(dc)),dc=NEXT(dc)){
218: t=t*pwrreal0(r,(d-QTOS(DEG(dc))))+BDY((Real)COEF(dc));
1.1 noro 219: }
1.9 saito 220: if(d)t*=pwrreal0(r,d);
1.1 noro 221: }
222: return t;
223: }
224:
1.9 saito 225: void qcalc(char **tab,struct canvas *can){
226: //qifplotmain(Old type)
1.5 noro 227: Q dx,dy,w,h,xstep,ystep,c,q1;
1.1 noro 228: P g,g1,f1,f2,x,y,t,s;
229: int ix,iy;
230: int *a,*pa;
231: VECT ss;
232:
233: subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
234: subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
235: MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
236: addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
237: MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
238: addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
239: ptozp(f2,1,&c,&g);
1.9 saito 240: a=(int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
1.1 noro 241: initmarker(can,"Horizontal scan...");
1.9 saito 242: for( ix=0; ix < can->width; ix++ ){
1.1 noro 243: marker(can,DIR_X,ix);
244: STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
1.9 saito 245: if( !g1 )
246: for(iy=0; iy < can->height; iy++ )
247: tab[ix][iy]=1;
248: else if( !NUM(g1) ){
1.5 noro 249: sturmseq(CO,g1,&ss); seproot(ss,0,can->height,a);
1.9 saito 250: for(iy=0, pa=a; iy < can->height; iy++, pa++ )
251: if( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
252: tab[ix][iy]=1;
1.1 noro 253: }
254: }
255: initmarker(can,"Vertical scan...");
1.9 saito 256: for( iy=0; iy < can->height; iy++ ){
1.1 noro 257: marker(can,DIR_Y,iy);
258: STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
1.9 saito 259: if( !g1 )
260: for(ix=0; ix < can->width; ix++ )
261: tab[ix][iy]=1;
262: else if( !NUM(g1) ){
1.5 noro 263: sturmseq(CO,g1,&ss); seproot(ss,0,can->width,a);
1.9 saito 264: for(ix=0, pa=a; ix < can->width; ix++, pa++ )
265: if( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
266: tab[ix][iy]=1;
1.1 noro 267: }
268: }
269: }
270:
1.9 saito 271: void sturmseq(VL vl,P p,VECT *rp){
272: P g1,g2,q,r,s,*t;
1.1 noro 273: V v;
274: VECT ret;
275: int i,j;
276: Q a,b,c,d,h,l,m,x;
277:
1.9 saito 278: v=VR(p);t=(P *)ALLOCA((deg(v,p)+1)*sizeof(P));
279: g1=t[0]=p;diffp(vl,p,v,(P *)&a);ptozp((P)a,1,&c,&g2);t[1]=g2;
280: for(i=1,h=ONE,x=ONE;;){
281: if(NUM(g2)) break;
1.1 noro 282: subq(DEG(DC(g1)),DEG(DC(g2)),&d);
1.9 saito 283: l=(Q)LC(g2);
284: if(SGN(l)<0){
285: chsgnq(l,&a);l=a;
1.1 noro 286: }
1.9 saito 287: addq(d,ONE,&a);pwrq(l,a,&b);mulp(vl,(P)b,g1,(P *)&a);
1.1 noro 288: divsrp(vl,(P)a,g2,&q,&r);
1.9 saito 289: if(!r) break;
290: chsgnp(r,&s);r=s;i++;
291: if(NUM(r)){
292: t[i]=r;break;
293: }
294: pwrq(h,d,&m);g1=g2;
295: mulq(m,x,&a);divsp(vl,r,(P)a,&g2);t[i]=g2;
296: x=(Q)LC(g1);
297: if(SGN(x)<0){
298: chsgnq(x,&a);x=a;
1.1 noro 299: }
1.9 saito 300: pwrq(x,d,&a);mulq(a,h,&b);divq(b,m,&h);
1.1 noro 301: }
302: MKVECT(ret,i+1);
1.9 saito 303: for(j=0;j<=i;j++)
304: ret->body[j]=(pointer)t[j];
305: *rp=ret;
1.1 noro 306: }
307:
1.9 saito 308: void seproot(VECT s,int min,int max,int *ar){
309: P f,*ss;
1.1 noro 310: Q q,t;
311: int i,j,k;
312:
1.9 saito 313: ss=(P *)s->body;f=ss[0];
314: for(i=min;i<=max;i++){
315: STOQ(i,q);usubstqp(f,q,&t);
316: if(!t)ar[i]=-1;
1.1 noro 317: else {
1.9 saito 318: ar[i]=numch(s,q,t);break;
1.1 noro 319: }
320: }
1.9 saito 321: if(i>max) return;
322: for(j=max;j>= min;j--){
1.1 noro 323: STOQ(j,q); usubstqp(f,q,&t);
1.9 saito 324: if(!t)ar[j]=-1;
1.1 noro 325: else {
1.9 saito 326: if(i!=j)ar[j]=numch(s,q,t);
1.1 noro 327: break;
328: }
329: }
1.9 saito 330: if(j<=i+1) return;
331: if(ar[i]==ar[j]){
332: for(k=i+1;k<j;k++)ar[k]=ar[i];
1.1 noro 333: return;
334: }
1.9 saito 335: k=(i+j)/2;
1.1 noro 336: seproot(s,i,k,ar);
337: seproot(s,k,j,ar);
338: }
339:
1.9 saito 340: int numch(VECT s,Q n,Q a0){
1.1 noro 341: int len,i,c;
342: Q a;
343: P *ss;
344:
1.9 saito 345: len=s->len;ss=(P *)s->body;
346: for(i=1,c=0;i<len;i++){
1.1 noro 347: usubstqp(ss[i],n,&a);
1.9 saito 348: if(a){
349: if((SGN(a)>0 && SGN(a0)<0)||(SGN(a)<0&&SGN(a0)>0))c++;
350: a0=a;
1.1 noro 351: }
352: }
353: return c;
354: }
355:
1.9 saito 356: void usubstqp(P p,Q r,Q *v){
1.1 noro 357: Q d,d1,a,b,t;
358: DCP dc;
359:
1.9 saito 360: if(!p)
361: *v=0;
362: else if(NUM(p))*v=(Q)p;
1.1 noro 363: else {
1.9 saito 364: dc=DC(p);t=(Q)COEF(dc);
365: for(d=DEG(dc),dc=NEXT(dc);dc;d=DEG(dc),dc=NEXT(dc)){
366: subq(d,DEG(dc),&d1);pwrq(r,d1,&a);
367: mulq(t,a,&b);addq(b,(Q)COEF(dc),&t);
1.1 noro 368: }
1.9 saito 369: if(d){
370: pwrq(r,d,&a);mulq(t,a,&b);t=b;
1.1 noro 371: }
1.9 saito 372: *v=t;
1.1 noro 373: }
374: }
375:
1.13 ! noro 376: Num tobf(Num,int);
! 377:
! 378: void plotcalcbf(struct canvas *can){
! 379: Obj fr,s,t;
! 380: Num xmin,xmax,ymin,ymax,xstep;
! 381: Num u,v,ha,dx,dy,x;
! 382: Num *tab;
! 383: Real r;
! 384: Q w,h1;
! 385: int ix;
! 386: POINT *pa;
! 387: double rr;
! 388: Q prec;
! 389: NODE arg;
! 390:
! 391: STOQ(can->prec,prec); arg = mknode(1,prec); Psetprec(arg,&t);
! 392: evalr(CO,(Obj)can->formula,can->prec,&fr);
! 393: MKReal(can->xmin,r); xmin = tobf((Num)r,can->prec);
! 394: MKReal(can->xmax,r); xmax = tobf((Num)r,can->prec);
! 395: MKReal(can->ymin,r); ymin = tobf((Num)r,can->prec);
! 396: MKReal(can->ymax,r); ymax = tobf((Num)r,can->prec);
! 397: STOQ(can->width,w);
! 398: subbf(xmax,xmin,&dx); divbf(dx,(Num)w,&xstep);
! 399: tab=(Num *)MALLOC(can->width*sizeof(Num));
! 400: for(ix=0,x=xmin;ix<can->width;ix++){
! 401: substr(CO,0,fr,can->vx,(Obj)x,(Obj *)&s);
! 402: evalr(CO,(Obj)s,can->prec,&t);
! 403: if(t&&(OID(t)!=O_N))
! 404: error("plotcalcbf : invalid evaluation");
! 405: tab[ix]=(Num)t;
! 406: addbf(x,xstep,&u); x = u;
! 407: }
! 408: if(!cmpbf(ymax,ymin)){
! 409: for(ymax=ymin=tab[0],ix=1;ix<can->width;ix++){
! 410: if(cmpbf(tab[ix],ymax)>0)ymax=tab[ix];
! 411: if(cmpbf(tab[ix],ymin)<0)ymin=tab[ix];
! 412: }
! 413: can->ymax=ToReal(ymax);can->ymin=ToReal(ymin);
! 414: }
! 415: subbf(ymax,ymin,&dy);
! 416: can->pa=(struct pa *)MALLOC(sizeof(struct pa));
! 417: can->pa[0].length=can->width;
! 418: can->pa[0].pos=pa=(POINT *)MALLOC(can->width*sizeof(POINT));
! 419: STOQ(can->height-1,h1);
! 420: for(ix=0;ix<can->width;ix++){
! 421: XC(pa[ix])=ix;
! 422: subbf(ymax,tab[ix],&u); divbf(u,dy,&v); mulbf(v,(Num)h1,&u);
! 423: rr = ToReal(u);
! 424: if(rr>MAXSHORT)YC(pa[ix])=MAXSHORT;
! 425: else if(rr<-MAXSHORT)YC(pa[ix])=-MAXSHORT;
! 426: else YC(pa[ix])=(long)rr;
! 427: }
! 428: }
! 429:
1.9 saito 430: void plotcalc(struct canvas *can){
431: //plot,memory_plot,plotover,plot_resize
432: double x,xmin,xstep,ymax,ymin,dy,*tab,usubstrp();
433: int ix,w,h;
434: Real r,rx;
435: Obj fr,t,s;
1.1 noro 436: POINT *pa;
437:
1.13 ! noro 438: if ( can->prec ) {
! 439: plotcalcbf(can);
! 440: return;
! 441: }
1.12 noro 442: todouble((Obj)can->formula,&fr);
1.9 saito 443: w=can->width;h=can->height;
444: xmin=can->xmin;xstep=(can->xmax-can->xmin)/w;
445: tab=(double *)ALLOCA(w*sizeof(double));
446: MKReal(1,rx); // dummy real number
447: for(ix=0,x=xmin;ix<w;ix++,x+=xstep){
448: // full substitution
449: BDY(rx)=x;
1.5 noro 450: substr(CO,0,fr,can->vx,x?(Obj)rx:0,&s);
451: devalr(CO,(Obj)s,&t);
1.9 saito 452: if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
1.1 noro 453: error("plotcalc : invalid evaluation");
1.9 saito 454: tab[ix]=ToReal((Num)t);
1.1 noro 455: }
1.9 saito 456: if(can->ymax==can->ymin){
457: for(ymax=ymin=tab[0],ix=1;ix<w;ix++){
458: if(tab[ix]>ymax)ymax=tab[ix];
459: if(tab[ix]<ymin)ymin=tab[ix];
1.1 noro 460: }
1.9 saito 461: can->ymax=ymax;can->ymin=ymin;
1.1 noro 462: } else {
1.9 saito 463: ymax=can->ymax;ymin=can->ymin;
1.1 noro 464: }
1.9 saito 465: dy=ymax-ymin;
466: can->pa=(struct pa *)MALLOC(sizeof(struct pa));
467: can->pa[0].length=w;
468: can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
469: for(ix=0;ix<w;ix++){
1.1 noro 470: double t;
1.9 saito 471: XC(pa[ix])=ix;
472: t=(h-1)*(ymax-tab[ix])/dy;
473: if(t>MAXSHORT)YC(pa[ix])=MAXSHORT;
474: else if(t<-MAXSHORT)YC(pa[ix])=-MAXSHORT;
475: else YC(pa[ix])=(long)t;
1.6 noro 476: }
477: }
478:
1.10 saito 479: void polarcalc(struct canvas *can){
1.11 saito 480: double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep,tr,p,*tabx,*taby;
481: double usubstrp();
482: int i,nstep,w,h;
483: POINT *pa;
484: Real r;
485: Obj fr,t,s;
486:
1.12 noro 487: todouble((Obj)can->formula,&fr);
1.11 saito 488: w=can->width; h=can->height; nstep=can->nzstep;
489: pmin=can->zmin; pstep=(can->zmax-can->zmin)/nstep;
490: tabx=(double *)ALLOCA(nstep*sizeof(double));
491: taby=(double *)ALLOCA(nstep*sizeof(double));
492: MKReal(1,r); // dummy real number
493:
494: for(i=0,p=pmin;i<nstep;i++,p+= pstep){
495: // full substitution
496: BDY(r)=p;
497: substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
498: devalr(CO,(Obj)s,&t);
499: if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
500: error("polarcalc : invalid evaluation");
501: tr=ToReal((Num)t);
502: tabx[i]=tr*cos(p);
503: taby[i]=tr*sin(p);
504: }
505: xmax=xmin=tabx[0];
506: ymax=ymin=taby[0];
507: for(i=1;i<nstep;i++){
508: if(tabx[i]>xmax)xmax=tabx[i];
509: if(tabx[i]<xmin)xmin=tabx[i];
510: if(taby[i]>ymax)ymax=taby[i];
511: if(taby[i]<ymin)ymin=taby[i];
512: }
513: can->xmax=xmax;can->xmin=xmin;
514: can->ymax=ymax;can->ymin=ymin;
515: dx=xmax-xmin;
516: dy=ymax-ymin;
517: can->pa=(struct pa *)MALLOC(sizeof(struct pa));
518: can->pa[0].length=nstep;
519: can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
520: for(i=0;i<nstep;i++){
521: XC(pa[i])=(w-1)*(tabx[i]-xmin)/dx;
522: YC(pa[i])=(h-1)*(ymax-taby[i])/dy;
523: }
524: }
525:
526: void polarcalcNG(struct canvas *can){
1.9 saito 527: //polarplotNG
528: double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep,tr,p, *tabx,*taby;
1.6 noro 529: double usubstrp();
1.10 saito 530: int i,ix,iy,nstep,w,h;
1.6 noro 531: POINT *pa;
532: Real r;
533: Obj fr,t,s;
534:
1.12 noro 535: todouble((Obj)can->formula,&fr);
1.9 saito 536: w=can->width; h=can->height; nstep=can->nzstep;
537: pmin=can->zmin; pstep=(can->zmax-can->zmin)/nstep;
538: tabx=(double *)ALLOCA(nstep*sizeof(double));
539: taby=(double *)ALLOCA(nstep*sizeof(double));
540: MKReal(1,r); // dummy real number
1.10 saito 541:
1.9 saito 542: for(i=0,p=pmin;i<nstep;i++,p+= pstep){
543: // full substitution
544: BDY(r)=p;
1.6 noro 545: substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
546: devalr(CO,(Obj)s,&t);
1.9 saito 547: if(t&&(OID(t)!=O_N||NID((Num)t)!=N_R))
1.10 saito 548: error("polarcalc : invalid evaluation");
1.9 saito 549: tr=ToReal((Num)t);
550: tabx[i]=tr*cos(p);
551: taby[i]=tr*sin(p);
1.10 saito 552: if(i==0){
553: xmax=xmin=tabx[0];
554: ymax=ymin=taby[0];
555: } else {
556: if(tabx[i]>xmax)xmax=tabx[i];
557: if(tabx[i]<xmin)xmin=tabx[i];
558: if(taby[i]>ymax)ymax=taby[i];
559: if(taby[i]<ymin)ymin=taby[i];
560: }
1.9 saito 561: }
562: can->xmax=xmax;can->xmin=xmin;
563: can->ymax=ymax;can->ymin=ymin;
564: dx=xmax-xmin;
565: dy=ymax-ymin;
566: can->pa=(struct pa *)MALLOC(sizeof(struct pa));
567: can->pa[0].length=nstep;
568: can->pa[0].pos=pa=(POINT *)MALLOC(w*sizeof(POINT));
569: for(i=0;i<nstep;i++){
570: XC(pa[i])=(w-1)*(tabx[i]-xmin)/dx;
571: YC(pa[i])=(h-1)*(ymax-taby[i])/dy;
1.1 noro 572: }
573: }
1.8 saito 574:
1.9 saito 575: /*
576: void ineqncalc(double **tab,struct canvas *can,int nox){
577: double x,y,xmin,ymin,xstep,ystep;
578: int ix,iy,w,h;
579: Real r,rx,ry;
580: Obj fr,g,t,s;
581: V vx,vy;
1.8 saito 582:
1.9 saito 583: if(!nox) initmarker(can,"Evaluating...");
1.12 noro 584: todouble((Obj)can->formula,&fr);
1.9 saito 585: vx=can->vx;vy=can->vy;
586: w=can->width;h=can->height;
587: xmin=can->xmin;xstep=(can->xmax-can->xmin)/w;
588: ymin=can->ymin;ystep=(can->ymin-can->ymin)/h;
589: MKReal(1.0,rx); MKReal(1.0,ry); // dummy real
1.8 saito 590:
1.9 saito 591: for(ix=0,x=xmin;ix<=w;ix++,x+=xstep){
592: BDY(rx)=x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
1.8 saito 593: devalr(CO,t,&g);
1.9 saito 594: if(!nox) marker(can,DIR_X,ix);
595: for(iy=0,y=ymin;iy<=h;iy++,y+=ystep){
596: BDY(ry)=y;
597: substr(CO,0,g,vy,y?(Obj)ry:0,&t);
1.8 saito 598: devalr(CO,t,&s);
1.9 saito 599: tab[ix][iy]=ToReal(s);
1.8 saito 600: }
601: }
602: }
1.9 saito 603: */
1.8 saito 604:
1.9 saito 605: #if defined(INTERVAL)
606: void itvcalc(double **mask, struct canvas *can, int nox){
1.10 saito 607: //ITVIFPLOT
1.9 saito 608: double x,y,xstep,ystep,dx,dy,wx,wy;
609: int idv,ix,iy,idx,idy;
610: Itv ity,itx,ddx,ddy;
611: Real r,rx,ry,rx1,ry1,rdx,rdy,rdx1,rdy1;
612: V vx,vy;
613: Obj fr,g,t,s;
1.8 saito 614:
1.9 saito 615: idv=can->division;
1.12 noro 616: todouble((Obj)can->formula,&fr);
1.9 saito 617: vx=can->vx; vy=can->vy;
618: xstep=(can->xmax-can->xmin)/can->width;
619: ystep=(can->ymax-can->ymin)/can->height;
620: if(idv!=0){
621: wx=xstep/can->division;
622: wy=ystep/can->division;
623: }
624: MKReal(can->ymin,ry1);
625: for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
626: ry=ry1;
627: MKReal(y+ystep,ry1);
1.8 saito 628: istoitv((Num)(ry1),(Num)ry,&ity);
629: substr(CO,0,(Obj)fr,vy,(Obj)ity,&t);
1.9 saito 630: MKReal(can->xmin,rx1);
631: for(ix=0,x=can->xmin; ix<can->width; ix++,x+=xstep){
632: rx=rx1;
1.8 saito 633: MKReal(x+xstep,rx1);
1.9 saito 634: istoitv((Num)(rx1),(Num)rx,&itx);
635: substr(CO,0,(Obj)fr,vx,(Obj)itx,&t);
636: MKReal(can->ymin,ry1);
637: for(iy=0,y=can->ymin; iy<can->height; iy++,y+=ystep){
638: ry=ry1;
639: MKReal(y+ystep,ry1);
640: istoitv((Num)ry,(Num)ry1,&ity);
641: substr(CO,0,(Obj)t,vy,(Obj)ity,&g);
642: if(compnum(0,0,(Num)g))mask[ix][iy]=-1;
643: else {
644: mask[ix][iy]=0;
645: /*
646: if(idv==0) mask[ix][iy]=0;
647: else {
648: MKReal(y,rdy1);
649: for(idy=0,dy=y;idy<idv;dy+=wy,idy++){
650: rdy=rdy1;
651: MKReal(dy+wy,rdy1);
652: istoitv((Num)rdy,(Num)rdy1,&ddy);
653: substr(CO,0,(Obj)fr,vy,(Obj)ddy,&t);
654: MKReal(x,rdx1);
655: for(idx=0,dx=x;idx<idx;dx+=wx,idx++){
656: rdx=rdx1;
657: MKReal(dx+wx,rdx1);
658: istoitv((Num)rdx,(Num)rdx1,&ddx);
659: substr(CO,0,(Obj)t,vx,(Obj)ddx,&g);
660: if(!compnum(0,0,(Num)g)){
661: mask[ix][iy]=0;
662: break;
663: }
664: }
665: if(mask[ix][iy]==0)break;
666: }
667: }
668: */
669: }
670: }
1.8 saito 671: }
672: }
673: }
674: #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