Annotation of OpenXM_contrib2/asir2000/plot/calc.c, Revision 1.8
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.8 ! saito 48: * $OpenXM: OpenXM_contrib2/asir2000/plot/calc.c,v 1.7 2003/02/14 22:29:19 ohara 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.5 noro 62: void calc(double **tab,struct canvas *can,int nox)
1.1 noro 63: {
64: double x,y,xmin,ymin,xstep,ystep;
65: int ix,iy;
66: Real r,rx,ry;
67: Obj fr,g;
68: int w,h;
69: V vx,vy;
70: Obj t,s;
71:
1.4 noro 72: if ( !nox ) initmarker(can,"Evaluating...");
1.1 noro 73: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
74: vx = can->vx;
75: vy = can->vy;
76: w = can->width; h = can->height;
77: xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;
78: ymin = can->ymin; ystep = (can->ymax-can->ymin)/h;
79: MKReal(1.0,rx); MKReal(1.0,ry); /* dummy real */
80: for( ix = 0, x = xmin; ix < w ; ix++, x += xstep ) {
1.5 noro 81: BDY(rx) = x; substr(CO,0,fr,vx,x?(Obj)rx:0,&t);
82: devalr(CO,t,&g);
1.4 noro 83: if ( !nox ) marker(can,DIR_X,ix);
1.1 noro 84: for( iy = 0, y = ymin; iy < h ; iy++, y += ystep ) {
85: BDY(ry) = y;
1.5 noro 86: substr(CO,0,g,vy,y?(Obj)ry:0,&t);
87: devalr(CO,t,&s);
1.1 noro 88: tab[ix][iy] = ToReal(s);
89: }
90: }
91: }
92:
1.5 noro 93: double usubstrp(P p,double r)
1.1 noro 94: {
95: double t;
96: DCP dc;
97: int d;
98: double pwrreal0();
99:
100: if ( !p )
101: t = 0.0;
102: else if ( NUM(p) )
103: t = BDY((Real)p);
104: else {
105: dc = DC(p); t = BDY((Real)COEF(dc));
106: for ( d = QTOS(DEG(dc)), dc = NEXT(dc); dc;
107: d = QTOS(DEG(dc)), dc = NEXT(dc) ) {
108: t = t*pwrreal0(r,(d-QTOS(DEG(dc))))+BDY((Real)COEF(dc));
109: }
110: if ( d )
111: t *= pwrreal0(r,d);
112: }
113: return t;
114: }
115:
1.5 noro 116: void qcalc(char **tab,struct canvas *can)
1.1 noro 117: {
1.5 noro 118: Q dx,dy,w,h,xstep,ystep,c,q1;
1.1 noro 119: P g,g1,f1,f2,x,y,t,s;
120: int ix,iy;
121: int *a,*pa;
122: VECT ss;
123:
124: subq(can->qxmax,can->qxmin,&dx); STOQ(can->width,w); divq(dx,w,&xstep);
125: subq(can->qymax,can->qymin,&dy); STOQ(can->height,h); divq(dy,h,&ystep);
126: MKV(can->vx,x); mulp(CO,(P)xstep,x,&t);
127: addp(CO,(P)can->qxmin,t,&s); substp(CO,can->formula,can->vx,s,&f1);
128: MKV(can->vy,y); mulp(CO,(P)ystep,y,&t);
129: addp(CO,(P)can->qymin,t,&s); substp(CO,f1,can->vy,s,&f2);
130: ptozp(f2,1,&c,&g);
131: a = (int *)ALLOCA((MAX(can->width,can->height)+1)*sizeof(int));
132: initmarker(can,"Horizontal scan...");
133: for( ix = 0; ix < can->width; ix++ ) {
134: marker(can,DIR_X,ix);
135: STOQ(ix,q1); substp(CO,g,can->vx,(P)q1,&t); ptozp(t,1,&c,&g1);
136: if ( !g1 )
137: for ( iy = 0; iy < can->height; iy++ )
138: tab[ix][iy] = 1;
139: else if ( !NUM(g1) ) {
1.5 noro 140: sturmseq(CO,g1,&ss); seproot(ss,0,can->height,a);
1.1 noro 141: for ( iy = 0, pa = a; iy < can->height; iy++, pa++ )
142: if ( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
143: tab[ix][iy] = 1;
144: }
145: }
146: initmarker(can,"Vertical scan...");
147: for( iy = 0; iy < can->height; iy++ ) {
148: marker(can,DIR_Y,iy);
149: STOQ(iy,q1); substp(CO,g,can->vy,(P)q1,&t); ptozp(t,1,&c,&g1);
150: if ( !g1 )
151: for ( ix = 0; ix < can->width; ix++ )
152: tab[ix][iy] = 1;
153: else if ( !NUM(g1) ) {
1.5 noro 154: sturmseq(CO,g1,&ss); seproot(ss,0,can->width,a);
1.1 noro 155: for ( ix = 0, pa = a; ix < can->width; ix++, pa++ )
156: if ( *pa < 0 || (*(pa+1) >= 0 && (*pa > *(pa+1))) )
157: tab[ix][iy] = 1;
158: }
159: }
160: }
161:
1.5 noro 162: void sturmseq(VL vl,P p,VECT *rp)
1.1 noro 163: {
164: P g1,g2,q,r,s;
165: P *t;
166: V v;
167: VECT ret;
168: int i,j;
169: Q a,b,c,d,h,l,m,x;
170:
171: v = VR(p); t = (P *)ALLOCA((deg(v,p)+1)*sizeof(P));
172: g1 = t[0] = p; diffp(vl,p,v,(P *)&a); ptozp((P)a,1,&c,&g2); t[1] = g2;
173: for ( i = 1, h = ONE, x = ONE; ; ) {
174: if ( NUM(g2) )
175: break;
176: subq(DEG(DC(g1)),DEG(DC(g2)),&d);
177: l = (Q)LC(g2);
178: if ( SGN(l) < 0 ) {
179: chsgnq(l,&a); l = a;
180: }
181: addq(d,ONE,&a); pwrq(l,a,&b); mulp(vl,(P)b,g1,(P *)&a);
182: divsrp(vl,(P)a,g2,&q,&r);
183: if ( !r )
184: break;
185: chsgnp(r,&s); r = s; i++;
186: if ( NUM(r) ) {
187: t[i] = r; break;
188: }
189: pwrq(h,d,&m); g1 = g2;
190: mulq(m,x,&a); divsp(vl,r,(P)a,&g2); t[i] = g2;
191: x = (Q)LC(g1);
192: if ( SGN(x) < 0 ) {
193: chsgnq(x,&a); x = a;
194: }
195: pwrq(x,d,&a); mulq(a,h,&b); divq(b,m,&h);
196: }
197: MKVECT(ret,i+1);
198: for ( j = 0; j <= i; j++ )
199: ret->body[j] = (pointer)t[j];
200: *rp = ret;
201: }
202:
1.5 noro 203: void seproot(VECT s,int min,int max,int *ar)
1.1 noro 204: {
205: P f;
206: P *ss;
207: Q q,t;
208: int i,j,k;
209:
210: ss = (P *)s->body; f = ss[0];
211: for ( i = min; i <= max; i++ ) {
212: STOQ(i,q); usubstqp(f,q,&t);
213: if ( !t )
214: ar[i] = -1;
215: else {
216: ar[i] = numch(s,q,t); break;
217: }
218: }
219: if ( i > max )
220: return;
221: for ( j = max; j >= min; j-- ) {
222: STOQ(j,q); usubstqp(f,q,&t);
223: if ( !t )
224: ar[j] = -1;
225: else {
226: if ( i != j )
227: ar[j] = numch(s,q,t);
228: break;
229: }
230: }
231: if ( j <= i+1 )
232: return;
233: if ( ar[i] == ar[j] ) {
234: for ( k = i+1; k < j; k++ )
235: ar[k] = ar[i];
236: return;
237: }
238: k = (i+j)/2;
239: seproot(s,i,k,ar);
240: seproot(s,k,j,ar);
241: }
242:
1.5 noro 243: int numch(VECT s,Q n,Q a0)
1.1 noro 244: {
245: int len,i,c;
246: Q a;
247: P *ss;
248:
249: len = s->len; ss = (P *)s->body;
250: for ( i = 1, c = 0; i < len; i++ ) {
251: usubstqp(ss[i],n,&a);
252: if ( a ) {
253: if ( (SGN(a)>0 && SGN(a0)<0) || (SGN(a)<0 && SGN(a0)>0) )
254: c++;
255: a0 = a;
256: }
257: }
258: return c;
259: }
260:
1.5 noro 261: void usubstqp(P p,Q r,Q *v)
1.1 noro 262: {
263: Q d,d1,a,b,t;
264: DCP dc;
265:
266: if ( !p )
267: *v = 0;
268: else if ( NUM(p) )
269: *v = (Q)p;
270: else {
271: dc = DC(p); t = (Q)COEF(dc);
272: for ( d = DEG(dc), dc = NEXT(dc); dc;
273: d = DEG(dc), dc = NEXT(dc) ) {
274: subq(d,DEG(dc),&d1); pwrq(r,d1,&a);
275: mulq(t,a,&b); addq(b,(Q)COEF(dc),&t);
276: }
277: if ( d ) {
278: pwrq(r,d,&a); mulq(t,a,&b); t = b;
279: }
280: *v = t;
281: }
282: }
283:
1.5 noro 284: void plotcalc(struct canvas *can)
1.1 noro 285: {
286: double x,xmin,xstep,ymax,ymin,dy;
287: int ix;
288: Real r;
289: Obj fr;
290: double usubstrp();
291: int w,h;
292: double *tab;
293: POINT *pa;
294: Real rx;
295: Obj t,s;
296:
297: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
298: w = can->width; h = can->height;
299: xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;
300: tab = (double *)ALLOCA(w*sizeof(double));
301: MKReal(1,rx); /* dummy real number */
302: for( ix = 0, x = xmin; ix < w ; ix++, x += xstep ) {
303: /* full substitution */
304: BDY(rx) = x;
1.5 noro 305: substr(CO,0,fr,can->vx,x?(Obj)rx:0,&s);
306: devalr(CO,(Obj)s,&t);
1.1 noro 307: if ( t && (OID(t)!=O_N || NID((Num)t)!=N_R) )
308: error("plotcalc : invalid evaluation");
309: tab[ix] = ToReal((Num)t);
310: }
311: if ( can->ymax == can->ymin ) {
312: for ( ymax = ymin = tab[0], ix = 1; ix < w; ix++ ) {
313: if ( tab[ix] > ymax )
314: ymax = tab[ix];
315: if ( tab[ix] < ymin )
316: ymin = tab[ix];
317: }
318: can->ymax = ymax; can->ymin = ymin;
319: } else {
320: ymax = can->ymax; ymin = can->ymin;
321: }
322: dy = ymax-ymin;
323: can->pa = (struct pa *)MALLOC(sizeof(struct pa));
324: can->pa[0].length = w;
325: can->pa[0].pos = pa = (POINT *)MALLOC(w*sizeof(POINT));
326: for ( ix = 0; ix < w; ix++ ) {
327: double t;
328:
329: XC(pa[ix]) = ix;
330: t = (h - 1)*(ymax - tab[ix])/dy;
331: if ( t > MAXSHORT )
332: YC(pa[ix]) = MAXSHORT;
333: else if ( t < -MAXSHORT )
334: YC(pa[ix]) = -MAXSHORT;
335: else
1.5 noro 336: YC(pa[ix]) = (long)t;
1.6 noro 337: }
338: }
339:
340: void polarplotcalc(struct canvas *can)
341: {
342: double xmax,xmin,ymax,ymin,dx,dy,pmin,pstep;
343: int i,nstep;
344: double usubstrp();
345: int w,h;
346: double tr,p;
347: double *tabx,*taby;
348: POINT *pa;
349: Real r;
350: Obj fr,t,s;
351:
352: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
353: w = can->width; h = can->height; nstep = can->nzstep;
354: pmin = can->zmin; pstep = (can->zmax-can->zmin)/nstep;
355: tabx = (double *)ALLOCA(nstep*sizeof(double));
356: taby = (double *)ALLOCA(nstep*sizeof(double));
357: MKReal(1,r); /* dummy real number */
358: for( i = 0, p = pmin; i < nstep ; i++, p += pstep ) {
359: /* full substitution */
360: BDY(r) = p;
361: substr(CO,0,fr,can->vx,p?(Obj)r:0,&s);
362: devalr(CO,(Obj)s,&t);
363: if ( t && (OID(t)!=O_N || NID((Num)t)!=N_R) )
364: error("polarplotcalc : invalid evaluation");
365: tr = ToReal((Num)t);
366: tabx[i] = tr*cos(p);
367: taby[i] = tr*sin(p);
368: }
369: xmax = xmin = tabx[0];
370: ymax = ymin = taby[0];
371: for ( i = 1; i < nstep; i++ ) {
372: if ( tabx[i] > xmax ) xmax = tabx[i];
373: if ( tabx[i] < xmin ) xmin = tabx[i];
374: if ( taby[i] > ymax ) ymax = taby[i];
375: if ( taby[i] < ymin ) ymin = taby[i];
376: }
377: can->xmax = xmax; can->xmin = xmin;
378: can->ymax = ymax; can->ymin = ymin;
379: dx = xmax-xmin;
380: dy = ymax-ymin;
381: can->pa = (struct pa *)MALLOC(sizeof(struct pa));
382: can->pa[0].length = nstep;
383: can->pa[0].pos = pa = (POINT *)MALLOC(w*sizeof(POINT));
384: for ( i = 0; i < nstep; i++ ) {
385: XC(pa[i]) = (w-1)*(tabx[i]-xmin)/dx;
386: YC(pa[i]) = (h-1)*(ymax-taby[i])/dy;
1.1 noro 387: }
388: }
1.8 ! saito 389:
! 390: #if defined(INTERVAL)
! 391: #define NORMAL 1
! 392: #define TRANSFER 1
! 393: #define RECURSION 1
! 394: #define RECTRANS 1
! 395:
! 396: void ineqncalc(double **tab,struct canvas *can,int nox)
! 397: {
! 398: double x, y, xmin, ymax, xstep, ystep;
! 399: int ix, iy;
! 400: Real r, rx, ry;
! 401: Obj fr, g, t, s;
! 402: int w, h;
! 403: V vx, vy;
! 404:
! 405: if ( !nox ) initmarker(can,"Evaluating...");
! 406: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
! 407: vx = can->vx; vy = can->vy;
! 408: w = can->width; h = can->height;
! 409: xmin = can->xmin; xstep = (can->xmax-can->xmin)/w;
! 410: ymax = can->ymax; ystep = (can->ymax-can->ymin)/h;
! 411: MKReal(1.0,rx); MKReal(1.0,ry); /* dummy real */
! 412:
! 413: for( iy = 0, y = ymax; iy <= h ; iy++, y -= ystep ) {
! 414: BDY(ry) = y; substr(CO,0,fr,vy,y?(Obj)ry:0,&t);
! 415: devalr(CO,t,&g);
! 416: if ( !nox ) marker(can,DIR_Y,iy);
! 417: for( ix = 0, x = xmin; ix <= w ; ix++, x += xstep ) {
! 418: BDY(rx) = x;
! 419: substr(CO,0,g,vx,x?(Obj)rx:0,&t);
! 420: devalr(CO,t,&s);
! 421: tab[iy][ix] = ToReal(s);
! 422: }
! 423: }
! 424: }
! 425:
! 426: void reccalc(P fr, V vx, V vy, int ixlw, int ixhg, int iylw, int iyhg,
! 427: double *xval, double *yval, int w, int **mask, int itvsize)
! 428: {
! 429: int i, j;
! 430: int ixmd, iymd;
! 431: double inf, sup, dt;
! 432: Real rxwid, rywid, rxlw, rylw, rzero;
! 433: Itv itx, ity;
! 434: Obj fr1;
! 435: P py, px, tmp;
! 436:
! 437: double **tbl , xs, ys;
! 438: int ywidth, xwidth;
! 439: Real rx, ry;
! 440: Obj g, t, s;
! 441:
! 442: if ( (ixlw > ixhg) || (iylw > iyhg) ) return;
! 443: NEWItvP(itx); NEWItvP(ity);
! 444: MKReal(xval[ixhg + 1] - xval[ixlw], rxwid);
! 445: MKReal(yval[iyhg + 1] - yval[iylw], rywid);
! 446: MKReal(0.0,rzero);
! 447: istoitv((Num)rzero, (Num)rxwid, &itx);
! 448: istoitv((Num)rzero, (Num)rywid, &ity);
! 449:
! 450: MKReal(xval[ixlw],rxlw);MKReal(yval[iylw],rylw);
! 451: MKV(vx,tmp); addp(CO,(P)tmp,(P)rxlw,&px);
! 452: MKV(vy,tmp); addp(CO,(P)tmp,(P)rylw,&py);
! 453: substp(CO,(P)fr,(V)vx,(P)px,(P*)&fr1);
! 454: substp(CO,(P)fr1,(V)vy,(P)py,(P*)&fr1);
! 455: substr(CO,0,(Obj)fr1,vx,(Obj)itx,&fr1);
! 456: substr(CO,0,(Obj)fr1,vy,(Obj)ity,&fr1);
! 457: Num2double((Num)fr1,&inf,&sup);
! 458:
! 459: if(inf > sup){dt=inf;inf=sup;sup=dt;}
! 460: if ( inf <= 0.0 && sup >= 0.0 ) {
! 461: if ( (ixhg - ixlw <= itvsize) && (iyhg - iylw <= itvsize) ){
! 462: ywidth = iyhg - iylw + 2; xwidth = ixhg - ixlw + 2;
! 463: MKReal(1.0, rx); MKReal(1.0, ry);
! 464: tbl = (double **)ALLOCA((ywidth)*sizeof(double *));
! 465: for (i=0;i<ywidth;i++)
! 466: tbl[i] = (double *)ALLOCA((xwidth)*sizeof(double));
! 467: for (i = 0; i < ywidth; i++){
! 468: BDY(ry) = yval[iylw+i];
! 469: substr(CO,0,(Obj)fr,vy,(Obj)ry,&t);
! 470: devalr(CO, t, &g);
! 471: for (j = 0; j < xwidth; j++){
! 472: BDY(rx) = xval[ixlw+j];
! 473: substr(CO,0,g,vx,(Obj)rx,&t);
! 474: devalr(CO, t, &s);
! 475: tbl[i][j] = ToReal(s);
! 476: }
! 477: }
! 478: for ( i = 0; i < ywidth - 1; i++){
! 479: for ( j = 0; j < xwidth - 1; j++){
! 480: if ( tbl[i][j] >= 0.0 ){
! 481: if ( (tbl[i+1][j] <= 0.0) ||
! 482: (tbl[i+1][j+1] <= 0.0) ||
! 483: (tbl[i][j+1] <= 0.0) ) mask[w-(i+iylw)][j+ixlw] = 0;
! 484: else mask[w-(i+iylw)][j+ixlw] = 1;
! 485: } else {
! 486: if ( (tbl[i+1][j] >= 0.0) ||
! 487: (tbl[i+1][j+1] >= 0.0) ||
! 488: (tbl[i][j+1] >= 0.0) ) mask[w-(i+iylw)][j+ixlw] = 0;
! 489: else mask[w-(i+iylw)][j+ixlw] = 1;
! 490: }
! 491: }
! 492: }
! 493: } else {
! 494: ixmd = (ixhg + ixlw)/2; iymd = (iyhg + iylw)/2;
! 495: reccalc((P)fr,vx,vy,ixlw,ixmd,iylw,iymd,xval,yval,w,mask,itvsize);
! 496: reccalc((P)fr,vx,vy,ixmd+1,ixhg,iylw,iymd,xval,yval,w,mask,itvsize);
! 497: reccalc((P)fr,vx,vy,ixlw,ixmd,iymd+1,iyhg,xval,yval,w,mask,itvsize);
! 498: reccalc((P)fr,vx,vy,ixmd+1,ixhg,iymd+1,iyhg,xval,yval,w,mask,itvsize);
! 499: }
! 500: } else {
! 501: for (i=w-iyhg; i> w-iylw; i--)
! 502: for (j=ixlw; j <= ixhg; j++) mask[i][j] = -1;
! 503: }
! 504: }
! 505:
! 506: void itvcalc(int **mask, struct canvas *can, int nox, int itvsize)
! 507: {
! 508: double xstep, ystep, xv, yv, *xval, *yval;
! 509: int imx, imy, i, w, h;
! 510: Real r;
! 511: Obj fr;
! 512: V vx, vy;
! 513:
! 514: if ( !nox ) initmarker(can,"Evaluating...");
! 515: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
! 516: vx = can->vx; vy = can->vy;
! 517: w = can->width; h = can->height;
! 518: xstep = (can->xmax - can->xmin)/w;
! 519: ystep = (can->ymax - can->ymin)/h;
! 520:
! 521: xval = (double *)ALLOCA((w+1)*sizeof(double));
! 522: yval = (double *)ALLOCA((h+1)*sizeof(double));
! 523: for (i=0, xv=can->xmin; i<= w; i++, xv += xstep) xval[i] = xv;
! 524: for (i=0, yv=can->ymin; i<= h; i++, yv += ystep) yval[i] = yv;
! 525: imx = w/2; imy = h/2;
! 526:
! 527: reccalc((P)fr,vx,vy,0,imx,0,imy,xval,yval,w-1,mask,itvsize);
! 528: reccalc((P)fr,vx,vy,imx+1,w-1,0,imy,xval,yval,w-1,mask,itvsize);
! 529: reccalc((P)fr,vx,vy,0,imx,imy+1,h-1,xval,yval,w-1,mask,itvsize);
! 530: reccalc((P)fr,vx,vy,imx+1,w-1,imy+1,h-1,xval,yval,w-1,mask,itvsize);
! 531: }
! 532:
! 533: #if NORMAL
! 534: void itvcalc1(int **mask, struct canvas *can, int nox)
! 535: {
! 536: double x, y, xstep, ystep;
! 537: int ix, iy, w, h;
! 538: Itv ity, itx;
! 539: Real r, rx, ry, rx1, ry1;
! 540: Obj fr, g, t;
! 541: V vx, vy;
! 542:
! 543: if ( !nox ) initmarker(can,"Evaluating...");
! 544: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
! 545: vx = can->vx; vy = can->vy;
! 546: w = can->width; h = can->height;
! 547: xstep = (can->xmax - can->xmin)/w;
! 548: ystep = (can->ymax - can->ymin)/h;
! 549:
! 550: for( iy = 0, y = can->ymax; iy < h ; iy++, y -= ystep ) {
! 551: MKReal(y, ry);
! 552: MKReal(y-ystep,ry1);
! 553: istoitv((Num)(ry1),(Num)ry,&ity);
! 554: substr(CO,0,(Obj)fr,vy,(Obj)ity,&t);
! 555: for( ix = 0, x = can->xmin; ix < w ; ix++, x += xstep ) {
! 556: MKReal(x,rx);
! 557: MKReal(x+xstep,rx1);
! 558: istoitv((Num)rx,(Num)rx1,&itx);
! 559: substr(CO,0,(Obj)t,vx,(Obj)itx,&g);
! 560: if (compnum(0,0,(Num)g)) mask[iy][ix] = -1;
! 561: else mask[iy][ix] = 0;
! 562: }
! 563: }
! 564: }
! 565: #endif
! 566:
! 567: #if TRANSFER
! 568: void itvcalc2(int **mask, struct canvas *can, int nox)
! 569: {
! 570: double x, y, xstep, ystep;
! 571: int ix, iy, w, h;
! 572: Itv ity, itx;
! 573: Real r, rx, ry, rzero;
! 574: Obj fr, g, t, s;
! 575: V vx, vy;
! 576: P mp, fr2;
! 577: Real qx,qy;
! 578:
! 579: if ( !nox ) initmarker(can,"Evaluating...");
! 580: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
! 581: vx = can->vx; vy = can->vy;
! 582: w = can->width; h = can->height;
! 583: xstep = (can->xmax - can->xmin)/w;
! 584: ystep = (can->ymax - can->ymin)/h;
! 585:
! 586: MKReal(0.0, rzero);
! 587: MKReal(ystep,ry);
! 588: istoitv((Num)rzero,(Num)ry,&ity);
! 589: MKReal(xstep,rx);
! 590: istoitv((Num)rzero,(Num)rx,&itx);
! 591: for( iy = 0, y = can->ymin; iy < h ; iy++, y += ystep ) {
! 592: MKReal(y,qy);
! 593: MKV(vy,mp);subp(CO,mp,(P)qy,&mp);
! 594: substp(CO,(P)fr,(V)vy,(P)mp,&fr2);
! 595: substr(CO,0,(Obj)fr2,vy,(Obj)ity,&t);
! 596: for( ix = 0, x = can->xmin; ix < w ; ix++, x += xstep ) {
! 597: MKReal(x,qx);
! 598: MKV(vx,mp);addp(CO,mp,(P)qx,&mp);
! 599: substp(CO,(P)t,(V)vx,(P)mp,(P*)&s);
! 600: substr(CO,0,(Obj)s,vx,(Obj)itx,&g);
! 601: if (compnum(0,0,(Num)g)) mask[iy][ix] = -1;
! 602: else mask[iy][ix] = 0;
! 603: }
! 604: }
! 605: }
! 606: #endif
! 607:
! 608: #if RECURSION
! 609: void reccalc3(P fr, V vx, V vy, int ixlw, int ixhg, int iylw, int iyhg,
! 610: double *xval, double *yval,int w, int **mask, int itvsize)
! 611: {
! 612: int i, j, ixmd, iymd;
! 613: double inf, sup, dt;
! 614: Real rxlw, rxhg, rylw, ryhg;
! 615: Itv itx, ity;
! 616: Obj fr1;
! 617:
! 618: NEWItvP(itx); NEWItvP(ity);
! 619:
! 620: if ( (ixlw > ixhg) || (iylw > iyhg) ) return;
! 621:
! 622: MKReal(xval[ixlw], rxlw); MKReal(xval[ixhg + 1], rxhg);
! 623: MKReal(yval[iylw], rylw); MKReal(yval[iyhg + 1], ryhg);
! 624: istoitv((Num)rxlw, (Num)rxhg, &itx);
! 625: istoitv((Num)rylw, (Num)ryhg, &ity);
! 626: substr(CO, 0, (Obj)fr, vy, (Obj)ity, &fr1);
! 627: substr(CO, 0, (Obj)fr1, vx, (Obj)itx, &fr1);
! 628: Num2double((Num)fr1, &inf, &sup);
! 629:
! 630: if( itvsize <= 0 ) itvsize = 1;
! 631: if(inf > sup){dt = inf; inf = sup; sup = dt;}
! 632: if ( inf <= 0.0 && sup >= 0.0 ) {
! 633: if ( (ixhg - ixlw < itvsize) || (iyhg - iylw < itvsize) ){
! 634: for(i = iylw; i <=iyhg; i++)
! 635: for(j = ixlw; j <=ixhg; j++) mask[w-i][j] = 0;
! 636: } else {
! 637: ixmd = (ixhg + ixlw)/2; iymd = (iyhg + iylw)/2;
! 638: reccalc3((P)fr,vx,vy,ixlw,ixmd,iylw,iymd,xval,yval,w,mask,itvsize);
! 639: reccalc3((P)fr,vx,vy,ixmd+1,ixhg,iylw,iymd,xval,yval,w,mask,itvsize);
! 640: reccalc3((P)fr,vx,vy,ixlw,ixmd,iymd+1,iyhg,xval,yval,w,mask,itvsize);
! 641: reccalc3((P)fr,vx,vy,ixmd+1,ixhg,iymd+1,iyhg,xval,yval,w,mask,itvsize);
! 642: }
! 643: } else {
! 644: for (i = iylw; i<= iyhg; i++)
! 645: for (j = ixlw; j <= ixhg; j++) mask[w-i][j] = -1;
! 646: }
! 647: }
! 648:
! 649: void itvcalc3(int **mask, struct canvas *can, int nox, int itvsize)
! 650: {
! 651: double xstep, ystep, xv, yv, *xval, *yval;
! 652: int i, imx, imy, w, h;
! 653: Real r;
! 654: Obj fr;
! 655: V vx, vy;
! 656:
! 657: if ( !nox ) initmarker(can,"Evaluating...");
! 658: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
! 659: vx = can->vx; vy = can->vy;
! 660: w = can->width; h = can->height;
! 661: xstep = (can->xmax - can->xmin)/w;
! 662: ystep = (can->ymax - can->ymin)/h;
! 663:
! 664: xval = (double *)ALLOCA((w+1)*sizeof(double));
! 665: yval = (double *)ALLOCA((h+1)*sizeof(double));
! 666: for (i = 0, xv = can->xmin; i <= w; i++, xv += xstep) xval[i] = xv;
! 667: for (i = 0, yv = can->ymin; i <= h; i++, yv += ystep) yval[i] = yv;
! 668: imx = w/2; imy = h/2;
! 669:
! 670: reccalc3((P)fr,vx,vy,0,imx,0,imy,xval,yval,w-1,mask,itvsize);
! 671: reccalc3((P)fr,vx,vy,imx+1,w-1,0,imy,xval,yval,w-1,mask,itvsize);
! 672: reccalc3((P)fr,vx,vy,0,imx,imy+1,h-1,xval,yval,w-1,mask,itvsize);
! 673: reccalc3((P)fr,vx,vy,imx+1,w-1,imy+1,h-1,xval,yval,w-1,mask,itvsize);
! 674: }
! 675: #endif
! 676:
! 677: #if RECTRANS
! 678: void reccalc4(P fr, V vx, V vy, int ixlw, int ixhg, int iylw, int iyhg,
! 679: double *xval, double *yval, int w, int **mask, int itvsize)
! 680: {
! 681: int i, j;
! 682: int ixmd, iymd;
! 683: double inf, sup, dt;
! 684: Real rxlw, rylw, rzero, rxwid, rywid;
! 685: Itv itx, ity;
! 686: Obj fr1;
! 687: P py, px, tmp;
! 688:
! 689: if ( (ixlw > ixhg) || (iylw > iyhg)) return;
! 690: NEWItvP(itx); NEWItvP(ity);
! 691:
! 692: MKReal(0.0,rzero);
! 693: MKReal(xval[ixhg + 1] - xval[ixlw], rxwid);
! 694: MKReal(yval[iyhg + 1] - yval[iylw], rywid);
! 695: istoitv((Num)rzero, (Num)rxwid, &itx);
! 696: istoitv((Num)rzero, (Num)rywid, &ity);
! 697:
! 698: MKReal(xval[ixlw], rxlw); MKReal(yval[iylw], rylw);
! 699: MKV(vx,tmp);addp(CO, (P)tmp, (P)rxlw, &px);
! 700: MKV(vy,tmp);addp(CO, (P)tmp, (P)rylw, &py);
! 701: substp(CO, (P)fr, (V)vx, (P)px, (P*)&fr1);
! 702: substp(CO, (P)fr1, (V)vy, (P)py, (P*)&fr1);
! 703: substr(CO, 0, (Obj)fr1, vx, (Obj)itx, &fr1);
! 704: substr(CO, 0, (Obj)fr1, vy, (Obj)ity, &fr1);
! 705: Num2double((Num)fr1, &inf, &sup);
! 706:
! 707: if(inf > sup){dt=inf;inf=sup;sup=dt;}
! 708: if ( inf <= 0.0 && sup >= 0.0 ) {
! 709: if ( (ixhg - ixlw <= itvsize) && (iyhg - iylw <= itvsize) ){
! 710: for(i = iylw; i <= iyhg; i++)
! 711: for(j = ixlw; j <= ixhg; j++) mask[w-i][j]=0;
! 712: } else {
! 713: ixmd = (ixhg + ixlw)/2; iymd = (iyhg + iylw)/2;
! 714:
! 715: reccalc4((P)fr,vx,vy,ixlw,ixmd,iylw,iymd,xval,yval,w,mask,itvsize);
! 716: reccalc4((P)fr,vx,vy,ixmd+1,ixhg,iylw,iymd,xval,yval,w,mask,itvsize);
! 717: reccalc4((P)fr,vx,vy,ixlw,ixmd,iymd+1,iyhg,xval,yval,w,mask,itvsize);
! 718: reccalc4((P)fr,vx,vy,ixmd+1,ixhg,iymd+1,iyhg,xval,yval,w,mask,itvsize);
! 719: }
! 720: } else {
! 721: for (i = iylw; i<= iyhg; i++)
! 722: for (j = ixlw; j <= ixhg; j++) mask[w-i][j] = -1;
! 723: }
! 724: }
! 725:
! 726: void itvcalc4(int **mask, struct canvas *can, int nox, int itvsize)
! 727: {
! 728: double xstep, ystep, xv, yv, *xval, *yval;
! 729: int i, imx, imy, w, h;
! 730: Real r;
! 731: Obj fr;
! 732: V vx, vy;
! 733:
! 734: if ( !nox ) initmarker(can,"Evaluating...");
! 735: MKReal(1.0,r); mulr(CO,(Obj)can->formula,(Obj)r,&fr);
! 736: vx = can->vx; vy = can->vy;
! 737: w = can->width; h = can->height;
! 738: xstep = (can->xmax - can->xmin)/w;
! 739: ystep = (can->ymax - can->ymin)/h;
! 740:
! 741: xval = (double *)ALLOCA((w+1)*sizeof(double));
! 742: yval = (double *)ALLOCA((h+1)*sizeof(double));
! 743: for (i=0, xv=can->xmin; i <= w; i++, xv += xstep) xval[i] = xv;
! 744: for (i=0, yv=can->ymin; i <= h; i++, yv += ystep) yval[i] = yv;
! 745: imx = w/2; imy = h/2;
! 746:
! 747: reccalc4((P)fr,vx,vy,0,imx,0,imy,xval,yval,w-1,mask,itvsize);
! 748: reccalc4((P)fr,vx,vy,imx+1,w-1,0,imy,xval,yval,w-1,mask,itvsize);
! 749: reccalc4((P)fr,vx,vy,0,imx,imy+1,h-1,xval,yval,w-1,mask,itvsize);
! 750: reccalc4((P)fr,vx,vy,imx+1,w-1,imy+1,h-1,xval,yval,w-1,mask,itvsize);
! 751: }
! 752: #endif
! 753: #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