File: common\lorenz.c
1 /*
2 This file contains two 3 dimensional orbit-type fractal
3 generators - IFS and LORENZ3D, along with code to generate
4 red/blue 3D images. Tim Wegner
5 */
6
7 #include <string.h>
8 /* see Fractint.c for a description of the "include" hierarchy */
9 #include "port.h"
10 #include "prototyp.h"
11 #include "fractype.h"
12
13 /* orbitcalc is declared with no arguments so jump through hoops here */
14 #define LORBIT(x,y,z) \
15 (*(int(*)(long *,long *,long *))curfractalspecific->orbitcalc)(x,y,z)
16 #define FORBIT(x,y,z) \
17 (*(int(*)(double*,double*,double*))curfractalspecific->orbitcalc)(x,y,z)
18
19 #define RANDOM(x) (rand()%(x))
20 /* BAD_PIXEL is used to cutoff orbits that are diverging. It might be better
21 to test the actual floating point orbit values, but this seems safe for now.
22 A higher value cannot be used - to test, turn off math coprocessor and
23 use +2.24 for type ICONS. If BAD_PIXEL is set to 20000, this will abort
24 Fractint with a math error. Note that this approach precludes zooming in very
25 far to an orbit type. */
26
27 #define BAD_PIXEL 10000L /* pixels can't get this big */
28
29 struct l_affine
30 {
31 /* weird order so a,b,e and c,d,f are vectors */
32 long a;
33 long b;
34 long e;
35 long c;
36 long d;
37 long f;
38 };
39 struct long3dvtinf /* data used by 3d view transform subroutine */
40 {
41 long orbit[3]; /* interated function orbit value */
42 long iview[3]; /* perspective viewer's coordinates */
43 long viewvect[3]; /* orbit transformed for viewing */
44 long viewvect1[3]; /* orbit transformed for viewing */
45 long maxvals[3];
46 long minvals[3];
47 MATRIX doublemat; /* transformation matrix */
48 MATRIX doublemat1; /* transformation matrix */
49 long longmat[4][4]; /* long version of matrix */
50 long longmat1[4][4]; /* long version of matrix */
51 int row,col; /* results */
52 int row1,col1;
53 struct l_affine cvt;
54 };
55
56 struct float3dvtinf /* data used by 3d view transform subroutine */
57 {
58 double orbit[3]; /* interated function orbit value */
59 double viewvect[3]; /* orbit transformed for viewing */
60 double viewvect1[3]; /* orbit transformed for viewing */
61 double maxvals[3];
62 double minvals[3];
63 MATRIX doublemat; /* transformation matrix */
64 MATRIX doublemat1; /* transformation matrix */
65 int row,col; /* results */
66 int row1,col1;
67 struct affine cvt;
68 };
69
70 /* Routines in this module */
71
72 static int ifs2d(void);
73 static int ifs3d(void);
74 static int ifs3dlong(void);
75 static int ifs3dfloat(void);
76 static int l_setup_convert_to_screen(struct l_affine *);
77 static void setupmatrix(MATRIX);
78 static int long3dviewtransf(struct long3dvtinf *inf);
79 static int float3dviewtransf(struct float3dvtinf *inf);
80 static FILE *open_orbitsave(void);
81 static void _fastcall plothist(int x, int y, int color);
82 static int realtime;
83 static char orbitsavename[] = {"orbits.raw"};
84 static char orbitsave_format[] = {"%g %g %g 15\n"};
85
86 S32 maxct;
87 static int t;
88 static long l_dx,l_dy,l_dz,l_dt,l_a,l_b,l_c,l_d;
89 static long l_adt,l_bdt,l_cdt,l_xdt,l_ydt;
90 static long initorbitlong[3];
91
92 static double dx,dy,dz,dt,a,b,c,d;
93 static double adt,bdt,cdt,xdt,ydt,zdt;
94 static double initorbitfp[3];
95
96 /*
97 * The following declarations used for Inverse Julia. MVS
98 */
99
100 static FCODE NoQueue[] =
101 "Not enough memory: switching to random walk.\n";
102
103 static int mxhits;
104 int run_length;
105 /*
106 enum {breadth_first, depth_first, random_walk, random_run} major_method;
107 enum {left_first, right_first} minor_method;
108 */
109 enum Major major_method;
110 enum Minor minor_method;
111 struct affine cvt;
112 struct l_affine lcvt;
113
114 double Cx, Cy;
115 long CxLong, CyLong;
116
117 /*
118 * end of Inverse Julia declarations;
119 */
120
121 /* these are potential user parameters */
122 int connect = 1; /* flag to connect points with a line */
123 int euler = 0; /* use implicit euler approximation for dynamic system */
124 int waste = 100; /* waste this many points before plotting */
125 int projection = 2; /* projection plane - default is to plot x-y */
126
127 /******************************************************************/
128 /* zoom box conversion functions */
129 /******************************************************************/
130
131 /*
132 Conversion of complex plane to screen coordinates for rotating zoom box.
133 Assume there is an affine transformation mapping complex zoom parallelogram
134 to rectangular screen. We know this map must map parallelogram corners to
135 screen corners, so we have following equations:
136
137 a*xxmin+b*yymax+e == 0 (upper left)
138 c*xxmin+d*yymax+f == 0
139
140 a*xx3rd+b*yy3rd+e == 0 (lower left)
141 c*xx3rd+d*yy3rd+f == ydots-1
142
143 a*xxmax+b*yymin+e == xdots-1 (lower right)
144 c*xxmax+d*yymin+f == ydots-1
145
146 First we must solve for a,b,c,d,e,f - (which we do once per image),
147 then we just apply the transformation to each orbit value.
148 */
149
150 /*
151 Thanks to Sylvie Gallet for the following. The original code for
152 setup_convert_to_screen() solved for coefficients of the
153 complex-plane-to-screen transformation using a very straight-forward
154 application of determinants to solve a set of simulataneous
155 equations. The procedure was simple and general, but inefficient.
156 The inefficiecy wasn't hurting anything because the routine was called
157 only once per image, but it seemed positively sinful to use it
158 because the code that follows is SO much more compact, at the
159 expense of being less general. Here are Sylvie's notes. I have further
160 optimized the code a slight bit.
161 Tim Wegner
162 July, 1996
163 Sylvie's notes, slightly edited follow:
164
165 You don't need 3x3 determinants to solve these sets of equations because
166 the unknowns e and f have the same coefficient: 1.
167
168 First set of 3 equations:
169 a*xxmin+b*yymax+e == 0
170 a*xx3rd+b*yy3rd+e == 0
171 a*xxmax+b*yymin+e == xdots-1
172 To make things easy to read, I just replace xxmin, xxmax, xx3rd by x1,
173 x2, x3 (ditto for yy...) and xdots-1 by xd.
174
175 a*x1 + b*y2 + e == 0 (1)
176 a*x3 + b*y3 + e == 0 (2)
177 a*x2 + b*y1 + e == xd (3)
178
179 I subtract (1) to (2) and (3):
180 a*x1 + b*y2 + e == 0 (1)
181 a*(x3-x1) + b*(y3-y2) == 0 (2)-(1)
182 a*(x2-x1) + b*(y1-y2) == xd (3)-(1)
183
184 I just have to calculate a 2x2 determinant:
185 det == (x3-x1)*(y1-y2) - (y3-y2)*(x2-x1)
186
187 And the solution is:
188 a = -xd*(y3-y2)/det
189 b = xd*(x3-x1)/det
190 e = - a*x1 - b*y2
191
192 The same technique can be applied to the second set of equations:
193
194 c*xxmin+d*yymax+f == 0
195 c*xx3rd+d*yy3rd+f == ydots-1
196 c*xxmax+d*yymin+f == ydots-1
197
198 c*x1 + d*y2 + f == 0 (1)
199 c*x3 + d*y3 + f == yd (2)
200 c*x2 + d*y1 + f == yd (3)
201
202 c*x1 + d*y2 + f == 0 (1)
203 c*(x3-x2) + d*(y3-y1) == 0 (2)-(3)
204 c*(x2-x1) + d*(y1-y2) == yd (3)-(1)
205
206 det == (x3-x2)*(y1-y2) - (y3-y1)*(x2-x1)
207
208 c = -yd*(y3-y1)/det
209 d = yd*(x3-x2))det
210 f = - c*x1 - d*y2
211
212 - Sylvie
213 */
214
215 int setup_convert_to_screen(struct affine *scrn_cnvt)
216 {
217 double det, xd, yd;
218
219 if((det = (xx3rd-xxmin)*(yymin-yymax) + (yymax-yy3rd)*(xxmax-xxmin))==0)
220 return(-1);
221 xd = dxsize/det;
222 scrn_cnvt->a = xd*(yymax-yy3rd);
223 scrn_cnvt->b = xd*(xx3rd-xxmin);
224 scrn_cnvt->e = -scrn_cnvt->a*xxmin - scrn_cnvt->b*yymax;
225
226 if((det = (xx3rd-xxmax)*(yymin-yymax) + (yymin-yy3rd)*(xxmax-xxmin))==0)
227 return(-1);
228 yd = dysize/det;
229 scrn_cnvt->c = yd*(yymin-yy3rd);
230 scrn_cnvt->d = yd*(xx3rd-xxmax);
231 scrn_cnvt->f = -scrn_cnvt->c*xxmin - scrn_cnvt->d*yymax;
232 return(0);
233 }
234
235 static int l_setup_convert_to_screen(struct l_affine *l_cvt)
236 {
237 struct affine cvt;
238
239 /* MCP 7-7-91, This function should return a something! */
240 if(setup_convert_to_screen(&cvt))
241 return(-1);
242 l_cvt->a = (long)(cvt.a*fudge); l_cvt->b = (long)(cvt.b*fudge); l_cvt->c = (long)(cvt.c*fudge);
243 l_cvt->d = (long)(cvt.d*fudge); l_cvt->e = (long)(cvt.e*fudge); l_cvt->f = (long)(cvt.f*fudge);
244
245 /* MCP 7-7-91 */
246 return(0);
247 }
248
249 /******************************************************************/
250 /* setup functions - put in fractalspecific[fractype].per_image */
251 /******************************************************************/
252
253 static double orbit;
254 static long l_orbit;
255 static long l_sinx,l_cosx;
256
257 int orbit3dlongsetup()
258 {
259 maxct = 0L;
260 connect = 1;
261 waste = 100;
262 projection = 2;
263 if (fractype==LHENON || fractype==KAM || fractype==KAM3D ||
264 fractype==INVERSEJULIA)
265 connect=0;
266 if(fractype==LROSSLER)
267 waste = 500;
268 if(fractype==LLORENZ)
269 projection = 1;
270
271 initorbitlong[0] = fudge; /* initial conditions */
272 initorbitlong[1] = fudge;
273 initorbitlong[2] = fudge;
274
275 if(fractype==LHENON)
276 {
277 l_a = (long)(param[0]*fudge);
278 l_b = (long)(param[1]*fudge);
279 l_c = (long)(param[2]*fudge);
280 l_d = (long)(param[3]*fudge);
281 }
282 else if(fractype==KAM || fractype==KAM3D)
283 {
284 maxct = 1L;
285 a = param[0]; /* angle */
286 if(param[1] <= 0.0)
287 param[1] = .01;
288 l_b = (long)(param[1]*fudge); /* stepsize */
289 l_c = (long)(param[2]*fudge); /* stop */
290 t = (int)(l_d = (long)(param[3])); /* points per orbit */
291
292 l_sinx = (long)(sin(a)*fudge);
293 l_cosx = (long)(cos(a)*fudge);
294 l_orbit = 0;
295 initorbitlong[0] = initorbitlong[1] = initorbitlong[2] = 0;
296 }
297 else if (fractype == INVERSEJULIA)
298 {
299 LCMPLX Sqrt;
300
301 CxLong = (long)(param[0] * fudge);
302 CyLong = (long)(param[1] * fudge);
303
304 mxhits = (int) param[2];
305 run_length = (int) param[3];
306 if (mxhits <= 0)
307 mxhits = 1;
308 else if (mxhits >= colors)
309 mxhits = colors - 1;
310 param[2] = mxhits;
311
312 setup_convert_to_screen(&cvt);
313 /* Note: using bitshift of 21 for affine, 24 otherwise */
314
315 lcvt.a = (long)(cvt.a * (1L << 21));
316 lcvt.b = (long)(cvt.b * (1L << 21));
317 lcvt.c = (long)(cvt.c * (1L << 21));
318 lcvt.d = (long)(cvt.d * (1L << 21));
319 lcvt.e = (long)(cvt.e * (1L << 21));
320 lcvt.f = (long)(cvt.f * (1L << 21));
321
322 Sqrt = ComplexSqrtLong(fudge - 4 * CxLong, -4 * CyLong);
323
324 switch (major_method) {
325 case breadth_first:
326 if (Init_Queue((long)32*1024) == 0)
327 { /* can't get queue memory: fall back to random walk */
328 stopmsg(20, NoQueue);
329 major_method = random_walk;
330 goto lrwalk;
331 }
332 EnQueueLong((fudge + Sqrt.x) / 2, Sqrt.y / 2);
333 EnQueueLong((fudge - Sqrt.x) / 2, -Sqrt.y / 2);
334 break;
335 case depth_first:
336 if (Init_Queue((long)32*1024) == 0)
337 { /* can't get queue memory: fall back to random walk */
338 stopmsg(20, NoQueue);
339 major_method = random_walk;
340 goto lrwalk;
341 }
342 switch (minor_method) {
343 case left_first:
344 PushLong((fudge + Sqrt.x) / 2, Sqrt.y / 2);
345 PushLong((fudge - Sqrt.x) / 2, -Sqrt.y / 2);
346 break;
347 case right_first:
348 PushLong((fudge - Sqrt.x) / 2, -Sqrt.y / 2);
349 PushLong((fudge + Sqrt.x) / 2, Sqrt.y / 2);
350 break;
351 }
352 break;
353 case random_walk:
354 lrwalk:
355 lnew.x = initorbitlong[0] = fudge + Sqrt.x / 2;
356 lnew.y = initorbitlong[1] = Sqrt.y / 2;
357 break;
358 case random_run:
359 lnew.x = initorbitlong[0] = fudge + Sqrt.x / 2;
360 lnew.y = initorbitlong[1] = Sqrt.y / 2;
361 break;
362 }
363 }
364 else
365 {
366 l_dt = (long)(param[0]*fudge);
367 l_a = (long)(param[1]*fudge);
368 l_b = (long)(param[2]*fudge);
369 l_c = (long)(param[3]*fudge);
370 }
371
372 /* precalculations for speed */
373 l_adt = multiply(l_a,l_dt,bitshift);
374 l_bdt = multiply(l_b,l_dt,bitshift);
375 l_cdt = multiply(l_c,l_dt,bitshift);
376 return(1);
377 }
378
379 #define COSB dx
380 #define SINABC dy
381
382 int orbit3dfloatsetup()
383 {
384 maxct = 0L;
385 connect = 1;
386 waste = 100;
387 projection = 2;
388
389 if(fractype==FPHENON || fractype==FPPICKOVER || fractype==FPGINGERBREAD
390 || fractype == KAMFP || fractype == KAM3DFP
391 || fractype == FPHOPALONG || fractype == INVERSEJULIAFP)
392 connect=0;
393 if(fractype==FPLORENZ3D1 || fractype==FPLORENZ3D3 ||
394 fractype==FPLORENZ3D4)
395 waste = 750;
396 if(fractype==FPROSSLER)
397 waste = 500;
398 if(fractype==FPLORENZ)
399 projection = 1; /* plot x and z */
400
401 initorbitfp[0] = 1; /* initial conditions */
402 initorbitfp[1] = 1;
403 initorbitfp[2] = 1;
404 if(fractype==FPGINGERBREAD)
405 {
406 initorbitfp[0] = param[0]; /* initial conditions */
407 initorbitfp[1] = param[1];
408 }
409
410 if(fractype==ICON || fractype==ICON3D) /* DMF */
411 {
412 initorbitfp[0] = 0.01; /* initial conditions */
413 initorbitfp[1] = 0.003;
414 connect = 0;
415 waste = 2000;
416 }
417
418 if(fractype==LATOO) /* HB */
419 {
420 connect = 0;
421 }
422
423 if(fractype==FPHENON || fractype==FPPICKOVER)
424 {
425 a = param[0];
426 b = param[1];
427 c = param[2];
428 d = param[3];
429 }
430 else if(fractype==ICON || fractype==ICON3D) /* DMF */
431 {
432 initorbitfp[0] = 0.01; /* initial conditions */
433 initorbitfp[1] = 0.003;
434 connect = 0;
435 waste = 2000;
436 /* Initialize parameters */
437 a = param[0];
438 b = param[1];
439 c = param[2];
440 d = param[3];
441 }
442 else if(fractype==KAMFP || fractype==KAM3DFP)
443 {
444 maxct = 1L;
445 a = param[0]; /* angle */
446 if(param[1] <= 0.0)
447 param[1] = .01;
448 b = param[1]; /* stepsize */
449 c = param[2]; /* stop */
450 t = (int)(l_d = (long)(param[3])); /* points per orbit */
451 sinx = sin(a);
452 cosx = cos(a);
453 orbit = 0;
454 initorbitfp[0] = initorbitfp[1] = initorbitfp[2] = 0;
455 }
456 else if(fractype==FPHOPALONG || fractype==FPMARTIN || fractype==CHIP
457 || fractype==QUADRUPTWO || fractype==THREEPLY)
458 {
459 initorbitfp[0] = 0; /* initial conditions */
460 initorbitfp[1] = 0;
461 initorbitfp[2] = 0;
462 connect = 0;
463 a = param[0];
464 b = param[1];
465 c = param[2];
466 d = param[3];
467 if(fractype==THREEPLY)
468 {
469 COSB = cos(b);
470 SINABC = sin(a+b+c);
471 }
472 }
473 else if (fractype == INVERSEJULIAFP)
474 {
475 _CMPLX Sqrt;
476
477 Cx = param[0];
478 Cy = param[1];
479
480 mxhits = (int) param[2];
481 run_length = (int) param[3];
482 if (mxhits <= 0)
483 mxhits = 1;
484 else if (mxhits >= colors)
485 mxhits = colors - 1;
486 param[2] = mxhits;
487
488 setup_convert_to_screen(&cvt);
489
490 /* find fixed points: guaranteed to be in the set */
491 Sqrt = ComplexSqrtFloat(1 - 4 * Cx, -4 * Cy);
492 switch ((int) major_method) {
493 case breadth_first:
494 if (Init_Queue((long)32*1024) == 0)
495 { /* can't get queue memory: fall back to random walk */
496 stopmsg(20, NoQueue);
497 major_method = random_walk;
498 goto rwalk;
499 }
500 EnQueueFloat((float)((1 + Sqrt.x) / 2), (float)(Sqrt.y / 2));
501 EnQueueFloat((float)((1 - Sqrt.x) / 2), (float)(-Sqrt.y / 2));
502 break;
503 case depth_first: /* depth first (choose direction) */
504 if (Init_Queue((long)32*1024) == 0)
505 { /* can't get queue memory: fall back to random walk */
506 stopmsg(20, NoQueue);
507 major_method = random_walk;
508 goto rwalk;
509 }
510 switch (minor_method) {
511 case left_first:
512 PushFloat((float)((1 + Sqrt.x) / 2), (float)(Sqrt.y / 2));
513 PushFloat((float)((1 - Sqrt.x) / 2), (float)(-Sqrt.y / 2));
514 break;
515 case right_first:
516 PushFloat((float)((1 - Sqrt.x) / 2), (float)(-Sqrt.y / 2));
517 PushFloat((float)((1 + Sqrt.x) / 2), (float)(Sqrt.y / 2));
518 break;
519 }
520 break;
521 case random_walk:
522 rwalk:
523 new.x = initorbitfp[0] = 1 + Sqrt.x / 2;
524 new.y = initorbitfp[1] = Sqrt.y / 2;
525 break;
526 case random_run: /* random run, choose intervals */
527 major_method = random_run;
528 new.x = initorbitfp[0] = 1 + Sqrt.x / 2;
529 new.y = initorbitfp[1] = Sqrt.y / 2;
530 break;
531 }
532 }
533 else
534 {
535 dt = param[0];
536 a = param[1];
537 b = param[2];
538 c = param[3];
539
540 }
541
542 /* precalculations for speed */
543 adt = a*dt;
544 bdt = b*dt;
545 cdt = c*dt;
546
547 return(1);
548 }
549
550 /******************************************************************/
551 /* orbit functions - put in fractalspecific[fractype].orbitcalc */
552 /******************************************************************/
553
554 /* Julia sets by inverse iterations added by Juan J. Buhler 4/3/92 */
555 /* Integrated with Lorenz by Tim Wegner 7/20/92 */
556 /* Add Modified Inverse Iteration Method, 11/92 by Michael Snyder */
557
558 int
559 Minverse_julia_orbit()
560 {
561 static int random_dir = 0, random_len = 0;
562 int newrow, newcol;
563 int color, leftright;
564
565 /*
566 * First, compute new point
567 */
568 switch (major_method) {
569 case breadth_first:
570 if (QueueEmpty())
571 return -1;
572 new = DeQueueFloat();
573 break;
574 case depth_first:
575 if (QueueEmpty())
576 return -1;
577 new = PopFloat();
578 break;
579 case random_walk:
580 #if 0
587 #endif
588 break;
589 case random_run:
590 #if 0
612 #endif
613 break;
614 }
615
616 /*
617 * Next, find its pixel position
618 */
619 newcol = (int)(cvt.a * new.x + cvt.b * new.y + cvt.e);
620 newrow = (int)(cvt.c * new.x + cvt.d * new.y + cvt.f);
621
622 /*
623 * Now find the next point(s), and flip a coin to choose one.
624 */
625
626 new = ComplexSqrtFloat(new.x - Cx, new.y - Cy);
627 leftright = (RANDOM(2)) ? 1 : -1;
628
629 if (newcol < 1 || newcol >= xdots || newrow < 1 || newrow >= ydots)
630 {
631 /*
632 * MIIM must skip points that are off the screen boundary,
633 * since it cannot read their color.
634 */
635 switch (major_method) {
636 case breadth_first:
637 EnQueueFloat((float)(leftright * new.x), (float)(leftright * new.y));
638 return 1;
639 case depth_first:
640 PushFloat ((float)(leftright * new.x), (float)(leftright * new.y));
641 return 1;
642 case random_run:
643 case random_walk:
644 break;
645 }
646 }
647
648 /*
649 * Read the pixel's color:
650 * For MIIM, if color >= mxhits, discard the point
651 * else put the point's children onto the queue
652 */
653 color = getcolor(newcol, newrow);
654 switch (major_method) {
655 case breadth_first:
656 if (color < mxhits)
657 {
658 putcolor(newcol, newrow, color+1);
659 /* new = ComplexSqrtFloat(new.x - Cx, new.y - Cy); */
660 EnQueueFloat( (float)new.x, (float)new.y);
661 EnQueueFloat((float)-new.x, (float)-new.y);
662 }
663 break;
664 case depth_first:
665 if (color < mxhits)
666 {
667 putcolor(newcol, newrow, color+1);
668 /* new = ComplexSqrtFloat(new.x - Cx, new.y - Cy); */
669 if (minor_method == left_first)
670 {
671 if (QueueFullAlmost())
672 PushFloat((float)-new.x, (float)-new.y);
673 else
674 {
675 PushFloat( (float)new.x, (float)new.y);
676 PushFloat((float)-new.x, (float)-new.y);
677 }
678 }
679 else
680 {
681 if (QueueFullAlmost())
682 PushFloat( (float)new.x, (float)new.y);
683 else
684 {
685 PushFloat((float)-new.x, (float)-new.y);
686 PushFloat( (float)new.x, (float)new.y);
687 }
688 }
689 }
690 break;
691 case random_run:
692 if (random_len-- == 0)
693 {
694 random_len = RANDOM(run_length);
695 random_dir = RANDOM(3);
696 }
697 switch (random_dir) {
698 case 0: /* left */
699 break;
700 case 1: /* right */
701 new.x = -new.x;
702 new.y = -new.y;
703 break;
704 case 2: /* random direction */
705 new.x = leftright * new.x;
706 new.y = leftright * new.y;
707 break;
708 }
709 if (color < colors-1)
710 putcolor(newcol, newrow, color+1);
711 break;
712 case random_walk:
713 if (color < colors-1)
714 putcolor(newcol, newrow, color+1);
715 new.x = leftright * new.x;
716 new.y = leftright * new.y;
717 break;
718 }
719 return 1;
720
721 }
722
723 int
724 Linverse_julia_orbit()
725 {
726 static int random_dir = 0, random_len = 0;
727 int newrow, newcol;
728 int color;
729
730 /*
731 * First, compute new point
732 */
733 switch (major_method) {
734 case breadth_first:
735 if (QueueEmpty())
736 return -1;
737 lnew = DeQueueLong();
738 break;
739 case depth_first:
740 if (QueueEmpty())
741 return -1;
742 lnew = PopLong();
743 break;
744 case random_walk:
745 lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
746 if (RANDOM(2))
747 {
748 lnew.x = -lnew.x;
749 lnew.y = -lnew.y;
750 }
751 break;
752 case random_run:
753 lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
754 if (random_len == 0)
755 {
756 random_len = RANDOM(run_length);
757 random_dir = RANDOM(3);
758 }
759 switch (random_dir) {
760 case 0: /* left */
761 break;
762 case 1: /* right */
763 lnew.x = -lnew.x;
764 lnew.y = -lnew.y;
765 break;
766 case 2: /* random direction */
767 if (RANDOM(2))
768 {
769 lnew.x = -lnew.x;
770 lnew.y = -lnew.y;
771 }
772 break;
773 }
774 }
775
776 /*
777 * Next, find its pixel position
778 *
779 * Note: had to use a bitshift of 21 for this operation because
780 * otherwise the values of lcvt were truncated. Used bitshift
781 * of 24 otherwise, for increased precision.
782 */
783 newcol = (int)((multiply(lcvt.a, lnew.x >> (bitshift - 21), 21) +
784 multiply(lcvt.b, lnew.y >> (bitshift - 21), 21) + lcvt.e) >> 21);
785 newrow = (int)((multiply(lcvt.c, lnew.x >> (bitshift - 21), 21) +
786 multiply(lcvt.d, lnew.y >> (bitshift - 21), 21) + lcvt.f) >> 21);
787
788 if (newcol < 1 || newcol >= xdots || newrow < 1 || newrow >= ydots)
789 {
790 /*
791 * MIIM must skip points that are off the screen boundary,
792 * since it cannot read their color.
793 */
794 if (RANDOM(2))
795 color = 1;
796 else
797 color = -1;
798 switch (major_method) {
799 case breadth_first:
800 lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
801 EnQueueLong(color * lnew.x, color * lnew.y);
802 break;
803 case depth_first:
804 lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
805 PushLong(color * lnew.x, color * lnew.y);
806 break;
807 case random_run:
808 random_len--;
809 case random_walk:
810 break;
811 }
812 return 1;
813 }
814
815 /*
816 * Read the pixel's color:
817 * For MIIM, if color >= mxhits, discard the point
818 * else put the point's children onto the queue
819 */
820 color = getcolor(newcol, newrow);
821 switch (major_method) {
822 case breadth_first:
823 if (color < mxhits)
824 {
825 putcolor(newcol, newrow, color+1);
826 lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
827 EnQueueLong( lnew.x, lnew.y);
828 EnQueueLong(-lnew.x, -lnew.y);
829 }
830 break;
831 case depth_first:
832 if (color < mxhits)
833 {
834 putcolor(newcol, newrow, color+1);
835 lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
836 if (minor_method == left_first)
837 {
838 if (QueueFullAlmost())
839 PushLong(-lnew.x, -lnew.y);
840 else
841 {
842 PushLong( lnew.x, lnew.y);
843 PushLong(-lnew.x, -lnew.y);
844 }
845 }
846 else
847 {
848 if (QueueFullAlmost())
849 PushLong( lnew.x, lnew.y);
850 else
851 {
852 PushLong(-lnew.x, -lnew.y);
853 PushLong( lnew.x, lnew.y);
854 }
855 }
856 }
857 break;
858 case random_run:
859 random_len--;
860 /* fall through */
861 case random_walk:
862 if (color < colors-1)
863 putcolor(newcol, newrow, color+1);
864 break;
865 }
866 return 1;
867 }
868
869 #if 0
898 #endif
899
900 int lorenz3dlongorbit(long *l_x, long *l_y, long *l_z)
901 {
902 l_xdt = multiply(*l_x,l_dt,bitshift);
903 l_ydt = multiply(*l_y,l_dt,bitshift);
904 l_dx = -multiply(l_adt,*l_x,bitshift) + multiply(l_adt,*l_y,bitshift);
905 l_dy = multiply(l_bdt,*l_x,bitshift) -l_ydt -multiply(*l_z,l_xdt,bitshift);
906 l_dz = -multiply(l_cdt,*l_z,bitshift) + multiply(*l_x,l_ydt,bitshift);
907
908 *l_x += l_dx;
909 *l_y += l_dy;
910 *l_z += l_dz;
911 return(0);
912 }
913
914 int lorenz3d1floatorbit(double *x, double *y, double *z)
915 {
916 double norm;
917
918 xdt = (*x)*dt;
919 ydt = (*y)*dt;
920 zdt = (*z)*dt;
921
922 /* 1-lobe Lorenz */
923 norm = sqrt((*x)*(*x)+(*y)*(*y));
924 dx = (-adt-dt)*(*x) + (adt-bdt)*(*y) + (dt-adt)*norm + ydt*(*z);
925 dy = (bdt-adt)*(*x) - (adt+dt)*(*y) + (bdt+adt)*norm - xdt*(*z) -
926 norm*zdt;
927 dz = (ydt/2) - cdt*(*z);
928
929 *x += dx;
930 *y += dy;
931 *z += dz;
932 return(0);
933 }
934
935 int lorenz3dfloatorbit(double *x, double *y, double *z)
936 {
937 xdt = (*x)*dt;
938 ydt = (*y)*dt;
939 zdt = (*z)*dt;
940
941 /* 2-lobe Lorenz (the original) */
942 dx = -adt*(*x) + adt*(*y);
943 dy = bdt*(*x) - ydt - (*z)*xdt;
944 dz = -cdt*(*z) + (*x)*ydt;
945
946 *x += dx;
947 *y += dy;
948 *z += dz;
949 return(0);
950 }
951
952 int lorenz3d3floatorbit(double *x, double *y, double *z)
953 {
954 double norm;
955
956 xdt = (*x)*dt;
957 ydt = (*y)*dt;
958 zdt = (*z)*dt;
959
960 /* 3-lobe Lorenz */
961 norm = sqrt((*x)*(*x)+(*y)*(*y));
962 dx = (-(adt+dt)*(*x) + (adt-bdt+zdt)*(*y)) / 3 +
963 ((dt-adt)*((*x)*(*x)-(*y)*(*y)) +
964 2*(bdt+adt-zdt)*(*x)*(*y))/(3*norm);
965 dy = ((bdt-adt-zdt)*(*x) - (adt+dt)*(*y)) / 3 +
966 (2*(adt-dt)*(*x)*(*y) +
967 (bdt+adt-zdt)*((*x)*(*x)-(*y)*(*y)))/(3*norm);
968 dz = (3*xdt*(*x)*(*y)-ydt*(*y)*(*y))/2 - cdt*(*z);
969
970 *x += dx;
971 *y += dy;
972 *z += dz;
973 return(0);
974 }
975
976 int lorenz3d4floatorbit(double *x, double *y, double *z)
977 {
978 xdt = (*x)*dt;
979 ydt = (*y)*dt;
980 zdt = (*z)*dt;
981
982 /* 4-lobe Lorenz */
983 dx = (-adt*(*x)*(*x)*(*x) + (2*adt+bdt-zdt)*(*x)*(*x)*(*y) +
984 (adt-2*dt)*(*x)*(*y)*(*y) + (zdt-bdt)*(*y)*(*y)*(*y)) /
985 (2 * ((*x)*(*x)+(*y)*(*y)));
986 dy = ((bdt-zdt)*(*x)*(*x)*(*x) + (adt-2*dt)*(*x)*(*x)*(*y) +
987 (-2*adt-bdt+zdt)*(*x)*(*y)*(*y) - adt*(*y)*(*y)*(*y)) /
988 (2 * ((*x)*(*x)+(*y)*(*y)));
989 dz = (2*xdt*(*x)*(*x)*(*y) - 2*xdt*(*y)*(*y)*(*y) - cdt*(*z));
990
991 *x += dx;
992 *y += dy;
993 *z += dz;
994 return(0);
995 }
996
997 int henonfloatorbit(double *x, double *y, double *z)
998 {
999 double newx,newy;
1000 *z = *x; /* for warning only */
1001 newx = 1 + *y - a*(*x)*(*x);
1002 newy = b*(*x);
1003 *x = newx;
1004 *y = newy;
1005 return(0);
1006 }
1007
1008 int henonlongorbit(long *l_x, long *l_y, long *l_z)
1009 {
1010 long newx,newy;
1011 *l_z = *l_x; /* for warning only */
1012 newx = multiply(*l_x,*l_x,bitshift);
1013 newx = multiply(newx,l_a,bitshift);
1014 newx = fudge + *l_y - newx;
1015 newy = multiply(l_b,*l_x,bitshift);
1016 *l_x = newx;
1017 *l_y = newy;
1018 return(0);
1019 }
1020
1021 int rosslerfloatorbit(double *x, double *y, double *z)
1022 {
1023 xdt = (*x)*dt;
1024 ydt = (*y)*dt;
1025
1026 dx = -ydt - (*z)*dt;
1027 dy = xdt + (*y)*adt;
1028 dz = bdt + (*z)*xdt - (*z)*cdt;
1029
1030 *x += dx;
1031 *y += dy;
1032 *z += dz;
1033 return(0);
1034 }
1035
1036 int pickoverfloatorbit(double *x, double *y, double *z)
1037 {
1038 double newx,newy,newz;
1039 newx = sin(a*(*y)) - (*z)*cos(b*(*x));
1040 newy = (*z)*sin(c*(*x)) - cos(d*(*y));
1041 newz = sin(*x);
1042 *x = newx;
1043 *y = newy;
1044 *z = newz;
1045 return(0);
1046 }
1047 /* page 149 "Science of Fractal Images" */
1048 int gingerbreadfloatorbit(double *x, double *y, double *z)
1049 {
1050 double newx;
1051 *z = *x; /* for warning only */
1052 newx = 1 - (*y) + fabs(*x);
1053 *y = *x;
1054 *x = newx;
1055 return(0);
1056 }
1057
1058 int rosslerlongorbit(long *l_x, long *l_y, long *l_z)
1059 {
1060 l_xdt = multiply(*l_x,l_dt,bitshift);
1061 l_ydt = multiply(*l_y,l_dt,bitshift);
1062
1063 l_dx = -l_ydt - multiply(*l_z,l_dt,bitshift);
1064 l_dy = l_xdt + multiply(*l_y,l_adt,bitshift);
1065 l_dz = l_bdt + multiply(*l_z,l_xdt,bitshift)
1066 - multiply(*l_z,l_cdt,bitshift);
1067
1068 *l_x += l_dx;
1069 *l_y += l_dy;
1070 *l_z += l_dz;
1071
1072 return(0);
1073 }
1074
1075 /* OSTEP = Orbit Step (and inner orbit value) */
1076 /* NTURNS = Outside Orbit */
1077 /* TURN2 = Points per orbit */
1078 /* a = Angle */
1079
1080
1081 int kamtorusfloatorbit(double *r, double *s, double *z)
1082 {
1083 double srr;
1084 if(t++ >= l_d)
1085 {
1086 orbit += b;
1087 (*r) = (*s) = orbit/3;
1088 t = 0;
1089 *z = orbit;
1090 if(orbit > c)
1091 return(1);
1092 }
1093 srr = (*s)-(*r)*(*r);
1094 (*s)=(*r)*sinx+srr*cosx;
1095 (*r)=(*r)*cosx-srr*sinx;
1096 return(0);
1097 }
1098
1099 int kamtoruslongorbit(long *r, long *s, long *z)
1100 {
1101 long srr;
1102 if(t++ >= l_d)
1103 {
1104 l_orbit += l_b;
1105 (*r) = (*s) = l_orbit/3;
1106 t = 0;
1107 *z = l_orbit;
1108 if(l_orbit > l_c)
1109 return(1);
1110 }
1111 srr = (*s)-multiply((*r),(*r),bitshift);
1112 (*s)=multiply((*r),l_sinx,bitshift)+multiply(srr,l_cosx,bitshift);
1113 (*r)=multiply((*r),l_cosx,bitshift)-multiply(srr,l_sinx,bitshift);
1114 return(0);
1115 }
1116
1117 int hopalong2dfloatorbit(double *x, double *y, double *z)
1118 {
1119 double tmp;
1120 *z = *x; /* for warning only */
1121 tmp = *y - sign(*x)*sqrt(fabs(b*(*x)-c));
1122 *y = a - *x;
1123 *x = tmp;
1124 return(0);
1125 }
1126
1127 /* from Michael Peters and HOP */
1128 int chip2dfloatorbit(double *x, double *y, double *z)
1129 {
1130 double tmp;
1131 *z = *x; /* for warning only */
1132 tmp = *y - sign(*x) * cos(sqr(log(fabs(b*(*x)-c))))
1133 * atan(sqr(log(fabs(c*(*x)-b))));
1134 *y = a - *x;
1135 *x = tmp;
1136 return(0);
1137 }
1138
1139 /* from Michael Peters and HOP */
1140 int quadruptwo2dfloatorbit(double *x, double *y, double *z)
1141 {
1142 double tmp;
1143 *z = *x; /* for warning only */
1144 tmp = *y - sign(*x) * sin(log(fabs(b*(*x)-c)))
1145 * atan(sqr(log(fabs(c*(*x)-b))));
1146 *y = a - *x;
1147 *x = tmp;
1148 return(0);
1149 }
1150
1151 /* from Michael Peters and HOP */
1152 int threeply2dfloatorbit(double *x, double *y, double *z)
1153 {
1154 double tmp;
1155 *z = *x; /* for warning only */
1156 tmp = *y - sign(*x)*(fabs(sin(*x)*COSB+c-(*x)*SINABC));
1157 *y = a - *x;
1158 *x = tmp;
1159 return(0);
1160 }
1161
1162 int martin2dfloatorbit(double *x, double *y, double *z)
1163 {
1164 double tmp;
1165 *z = *x; /* for warning only */
1166 tmp = *y - sin(*x);
1167 *y = a - *x;
1168 *x = tmp;
1169 return(0);
1170 }
1171
1172 int mandelcloudfloat(double *x, double *y, double *z)
1173 {
1174 double newx,newy,x2,y2;
1175 #ifndef XFRACT
1176 newx = *z; /* for warning only */
1177 #endif
1178 x2 = (*x)*(*x);
1179 y2 = (*y)*(*y);
1180 if (x2+y2>2) return 1;
1181 newx = x2-y2+a;
1182 newy = 2*(*x)*(*y)+b;
1183 *x = newx;
1184 *y = newy;
1185 return(0);
1186 }
1187
1188 int dynamfloat(double *x, double *y, double *z)
1189 {
1190 _CMPLX cp,tmp;
1191 double newx,newy;
1192 #ifndef XFRACT
1193 newx = *z; /* for warning only */
1194 #endif
1195 cp.x = b* *x;
1196 cp.y = 0;
1197 CMPLXtrig0(cp, tmp);
1198 newy = *y + dt*sin(*x + a*tmp.x);
1199 if (euler) {
1200 *y = newy;
1201 }
1202
1203 cp.x = b* *y;
1204 cp.y = 0;
1205 CMPLXtrig0(cp, tmp);
1206 newx = *x - dt*sin(*y + a*tmp.x);
1207 *x = newx;
1208 *y = newy;
1209 return(0);
1210 }
1211
1212 /* dmf */
1213 #undef LAMBDA
1214 #define LAMBDA param[0]
1215 #define ALPHA param[1]
1216 #define BETA param[2]
1217 #define GAMMA param[3]
1218 #define OMEGA param[4]
1219 #define DEGREE param[5]
1220
1221 int iconfloatorbit(double *x, double *y, double *z)
1222 {
1223
1224 double oldx, oldy, zzbar, zreal, zimag, za, zb, zn, p;
1225 int i;
1226
1227 oldx = *x;
1228 oldy = *y;
1229
1230 zzbar = oldx * oldx + oldy * oldy;
1231 zreal = oldx;
1232 zimag = oldy;
1233
1234 for(i=1; i <= DEGREE-2; i++) {
1235 za = zreal * oldx - zimag * oldy;
1236 zb = zimag * oldx + zreal * oldy;
1237 zreal = za;
1238 zimag = zb;
1239 }
1240 zn = oldx * zreal - oldy * zimag;
1241 p = LAMBDA + ALPHA * zzbar + BETA * zn;
1242 *x = p * oldx + GAMMA * zreal - OMEGA * oldy;
1243 *y = p * oldy - GAMMA * zimag + OMEGA * oldx;
1244
1245 *z = zzbar;
1246 return(0);
1247 }
1248 #ifdef LAMBDA /* Tim says this will make me a "good citizen" */
1249 #undef LAMBDA /* Yeah, but you were bad, Dan - LAMBDA was already */
1250 #undef ALPHA /* defined! <grin!> TW */
1251 #undef BETA
1252 #undef GAMMA
1253 #endif
1254
1255 /* hb */
1256 #define PAR_A param[0]
1257 #define PAR_B param[1]
1258 #define PAR_C param[2]
1259 #define PAR_D param[3]
1260
1261 int latoofloatorbit(double *x, double *y, double *z)
1262 {
1263
1264 double xold, yold, tmp;
1265
1266 xold = *z; /* for warning only */
1267
1268 xold = *x;
1269 yold = *y;
1270
1271 /* *x = sin(yold * PAR_B) + PAR_C * sin(xold * PAR_B); */
1272 old.x = yold * PAR_B;
1273 old.y = 0; /* old = (y * B) + 0i (in the complex)*/
1274 CMPLXtrig0(old,new);
1275 tmp = (double) new.x;
1276 old.x = xold * PAR_B;
1277 old.y = 0; /* old = (x * B) + 0i */
1278 CMPLXtrig1(old,new);
1279 *x = PAR_C * new.x + tmp;
1280
1281 /* *y = sin(xold * PAR_A) + PAR_D * sin(yold * PAR_A); */
1282 old.x = xold * PAR_A;
1283 old.y = 0; /* old = (y * A) + 0i (in the complex)*/
1284 CMPLXtrig2(old,new);
1285 tmp = (double) new.x;
1286 old.x = yold * PAR_A;
1287 old.y = 0; /* old = (x * B) + 0i */
1288 CMPLXtrig3(old,new);
1289 *y = PAR_D * new.x + tmp;
1290
1291 return(0);
1292 }
1293
1294 #undef PAR_A
1295 #undef PAR_B
1296 #undef PAR_C
1297 #undef PAR_D
1298
1299 /**********************************************************************/
1300 /* Main fractal engines - put in fractalspecific[fractype].calctype */
1301 /**********************************************************************/
1302
1303 int inverse_julia_per_image()
1304 {
1305 int color = 0;
1306
1307 if (resuming) /* can't resume */
1308 return -1;
1309
1310 while (color >= 0) /* generate points */
1311 {
1312 if (check_key())
1313 {
1314 Free_Queue();
1315 return -1;
1316 }
1317 color = curfractalspecific->orbitcalc();
1318 old = new;
1319 }
1320 Free_Queue();
1321 return 0;
1322 }
1323
1324 int orbit2dfloat()
1325 {
1326 FILE *fp;
1327 double *soundvar;
1328 double x,y,z;
1329 int color,col,row;
1330 int count;
1331 int oldrow, oldcol;
1332 double *p0,*p1,*p2;
1333 struct affine cvt;
1334 int ret;
1335 soundvar = p0 = p1 = p2 = NULL;
1336
1337 fp = open_orbitsave();
1338 /* setup affine screen coord conversion */
1339 setup_convert_to_screen(&cvt);
1340
1341 /* set up projection scheme */
1342 if(projection==0)
1343 {
1344 p0 = &z;
1345 p1 = &x;
1346 p2 = &y;
1347 }
1348 else if(projection==1)
1349 {
1350 p0 = &x;
1351 p1 = &z;
1352 p2 = &y;
1353 }
1354 else if(projection==2)
1355 {
1356 p0 = &x;
1357 p1 = &y;
1358 p2 = &z;
1359 }
1360 if((soundflag&7) == 2)
1361 soundvar = &x;
1362 else if((soundflag&7) == 3)
1363 soundvar = &y;
1364 else if((soundflag&7) == 4)
1365 soundvar = &z;
1366
1367 if(inside > 0)
1368 color = inside;
1369 if(color >= colors)
1370 color = 1;
1371 oldcol = oldrow = -1;
1372 x = initorbitfp[0];
1373 y = initorbitfp[1];
1374 z = initorbitfp[2];
1375 coloriter = 0L;
1376 count = ret = 0;
1377 if(maxit > 0x1fffffL || maxct)
1378 maxct = 0x7fffffffL;
1379 else
1380 maxct = maxit*1024L;
1381
1382 if (resuming)
1383 {
1384 start_resume();
1385 get_resume(sizeof(count),&count,sizeof(color),&color,
1386 sizeof(oldrow),&oldrow,sizeof(oldcol),&oldcol,
1387 sizeof(x),&x,sizeof(y),&y,sizeof(z),&z,sizeof(t),&t,
1388 sizeof(orbit),&orbit,sizeof(coloriter),&coloriter,
1389 0);
1390 end_resume();
1391 }
1392
1393 while(coloriter++ <= maxct) /* loop until keypress or maxit */
1394 {
1395 if(keypressed())
1396 {
1397 mute();
1398 alloc_resume(100,1);
1399 put_resume(sizeof(count),&count,sizeof(color),&color,
1400 sizeof(oldrow),&oldrow,sizeof(oldcol),&oldcol,
1401 sizeof(x),&x,sizeof(y),&y,sizeof(z),&z,sizeof(t),&t,
1402 sizeof(orbit),&orbit,sizeof(coloriter),&coloriter,
1403 0);
1404 ret = -1;
1405 break;
1406 }
1407 if (++count > 1000)
1408 { /* time to switch colors? */
1409 count = 0;
1410 if (++color >= colors) /* another color to switch to? */
1411 color = 1; /* (don't use the background color) */
1412 }
1413
1414 col = (int)(cvt.a*x + cvt.b*y + cvt.e);
1415 row = (int)(cvt.c*x + cvt.d*y + cvt.f);
1416 if ( col >= 0 && col < xdots && row >= 0 && row < ydots )
1417 {
1418 if ((soundflag&7) > 1)
1419 w_snd((int)(*soundvar*100+basehertz));
1420 if((fractype!=ICON) && (fractype!=LATOO))
1421 {
1422 if(oldcol != -1 && connect)
1423 draw_line(col,row,oldcol,oldrow,color%colors);
1424 else
1425 (*plot)(col,row,color%colors);
1426 } else {
1427 /* should this be using plothist()? */
1428 color = getcolor(col,row)+1;
1429 if( color < colors ) /* color sticks on last value */
1430 (*plot)(col,row,color);
1431
1432 }
1433
1434 oldcol = col;
1435 oldrow = row;
1436 }
1437 else if((long)abs(row) + (long)abs(col) > BAD_PIXEL) /* sanity check */
1438 return(ret);
1439 else
1440 oldrow = oldcol = -1;
1441
1442 if(FORBIT(p0, p1, p2))
1443 break;
1444 if(fp)
1445 fprintf(fp,orbitsave_format,*p0,*p1,0.0);
1446 }
1447 if(fp)
1448 fclose(fp);
1449 return(ret);
1450 }
1451
1452 int orbit2dlong()
1453 {
1454 FILE *fp;
1455 long *soundvar;
1456 long x,y,z;
1457 int color,col,row;
1458 int count;
1459 int oldrow, oldcol;
1460 long *p0,*p1,*p2;
1461 struct l_affine cvt;
1462 int ret,start;
1463 start = 1;
1464 soundvar = p0 = p1 = p2 = NULL;
1465 fp = open_orbitsave();
1466
1467 /* setup affine screen coord conversion */
1468 l_setup_convert_to_screen(&cvt);
1469
1470 /* set up projection scheme */
1471 if(projection==0)
1472 {
1473 p0 = &z;
1474 p1 = &x;
1475 p2 = &y;
1476 }
1477 else if(projection==1)
1478 {
1479 p0 = &x;
1480 p1 = &z;
1481 p2 = &y;
1482 }
1483 else if(projection==2)
1484 {
1485 p0 = &x;
1486 p1 = &y;
1487 p2 = &z;
1488 }
1489 if((soundflag&7)==2)
1490 soundvar = &x;
1491 else if((soundflag&7)==3)
1492 soundvar = &y;
1493 else if((soundflag&7)==4)
1494 soundvar = &z;
1495 if(inside > 0)
1496 color = inside;
1497 if(color >= colors)
1498 color = 1;
1499 oldcol = oldrow = -1;
1500 x = initorbitlong[0];
1501 y = initorbitlong[1];
1502 z = initorbitlong[2];
1503 count = ret = 0;
1504 if(maxit > 0x1fffffL || maxct)
1505 maxct = 0x7fffffffL;
1506 else
1507 maxct = maxit*1024L;
1508 coloriter = 0L;
1509
1510 if (resuming)
1511 {
1512 start_resume();
1513 get_resume(sizeof(count),&count,sizeof(color),&color,
1514 sizeof(oldrow),&oldrow,sizeof(oldcol),&oldcol,
1515 sizeof(x),&x,sizeof(y),&y,sizeof(z),&z,sizeof(t),&t,
1516 sizeof(l_orbit),&l_orbit,sizeof(coloriter),&coloriter,
1517 0);
1518 end_resume();
1519 }
1520
1521 while(coloriter++ <= maxct) /* loop until keypress or maxit */
1522 {
1523 if(keypressed())
1524 {
1525 mute();
1526 alloc_resume(100,1);
1527 put_resume(sizeof(count),&count,sizeof(color),&color,
1528 sizeof(oldrow),&oldrow,sizeof(oldcol),&oldcol,
1529 sizeof(x),&x,sizeof(y),&y,sizeof(z),&z,sizeof(t),&t,
1530 sizeof(l_orbit),&l_orbit,sizeof(coloriter),&coloriter,
1531 0);
1532 ret = -1;
1533 break;
1534 }
1535 if (++count > 1000)
1536 { /* time to switch colors? */
1537 count = 0;
1538 if (++color >= colors) /* another color to switch to? */
1539 color = 1; /* (don't use the background color) */
1540 }
1541
1542 col = (int)((multiply(cvt.a,x,bitshift) + multiply(cvt.b,y,bitshift) + cvt.e) >> bitshift);
1543 row = (int)((multiply(cvt.c,x,bitshift) + multiply(cvt.d,y,bitshift) + cvt.f) >> bitshift);
1544 if(overflow)
1545 {
1546 overflow = 0;
1547 return(ret);
1548 }
1549 if ( col >= 0 && col < xdots && row >= 0 && row < ydots )
1550 {
1551 if ((soundflag&7) > 1)
1552 {
1553 double yy;
1554 yy = *soundvar;
1555 yy = yy/fudge;
1556 w_snd((int)(yy*100+basehertz));
1557 }
1558 if(oldcol != -1 && connect)
1559 draw_line(col,row,oldcol,oldrow,color%colors);
1560 else if(!start)
1561 (*plot)(col,row,color%colors);
1562 oldcol = col;
1563 oldrow = row;
1564 start = 0;
1565 }
1566 else if((long)abs(row) + (long)abs(col) > BAD_PIXEL) /* sanity check */
1567 return(ret);
1568 else
1569 oldrow = oldcol = -1;
1570
1571 /* Calculate the next point */
1572 if(LORBIT(p0, p1, p2))
1573 break;
1574 if(fp)
1575 fprintf(fp,orbitsave_format,(double)*p0/fudge,(double)*p1/fudge,0.0);
1576 }
1577 if(fp)
1578 fclose(fp);
1579 return(ret);
1580 }
1581
1582 static int orbit3dlongcalc(void)
1583 {
1584 FILE *fp;
1585 unsigned long count;
1586 int oldcol,oldrow;
1587 int oldcol1,oldrow1;
1588 struct long3dvtinf inf;
1589 int color;
1590 int ret;
1591
1592 /* setup affine screen coord conversion */
1593 l_setup_convert_to_screen(&inf.cvt);
1594
1595 oldcol1 = oldrow1 = oldcol = oldrow = -1;
1596 color = 2;
1597 if(color >= colors)
1598 color = 1;
1599
1600 inf.orbit[0] = initorbitlong[0];
1601 inf.orbit[1] = initorbitlong[1];
1602 inf.orbit[2] = initorbitlong[2];
1603
1604 if(diskvideo) /* this would KILL a disk drive! */
1605 notdiskmsg();
1606
1607 fp = open_orbitsave();
1608
1609 count = ret = 0;
1610 if(maxit > 0x1fffffL || maxct)
1611 maxct = 0x7fffffffL;
1612 else
1613 maxct = maxit*1024L;
1614 coloriter = 0L;
1615 while(coloriter++ <= maxct) /* loop until keypress or maxit */
1616 {
1617 /* calc goes here */
1618 if (++count > 1000)
1619 { /* time to switch colors? */
1620 count = 0;
1621 if (++color >= colors) /* another color to switch to? */
1622 color = 1; /* (don't use the background color) */
1623 }
1624 if(keypressed())
1625 {
1626 mute();
1627 ret = -1;
1628 break;
1629 }
1630
1631 LORBIT(&inf.orbit[0],&inf.orbit[1],&inf.orbit[2]);
1632 if(fp)
1633 fprintf(fp,orbitsave_format,(double)inf.orbit[0]/fudge,(double)inf.orbit[1]/fudge,(double)inf.orbit[2]/fudge);
1634 if (long3dviewtransf(&inf))
1635 {
1636 /* plot if inside window */
1637 if (inf.col >= 0)
1638 {
1639 if(realtime)
1640 whichimage=1;
1641 if ((soundflag&7) > 1)
1642 {
1643 double yy;
1644 yy = inf.viewvect[((soundflag&7) - 2)];
1645 yy = yy/fudge;
1646 w_snd((int)(yy*100+basehertz));
1647 }
1648 if(oldcol != -1 && connect)
1649 draw_line(inf.col,inf.row,oldcol,oldrow,color%colors);
1650 else
1651 (*plot)(inf.col,inf.row,color%colors);
1652 }
1653 else if (inf.col == -2)
1654 return(ret);
1655 oldcol = inf.col;
1656 oldrow = inf.row;
1657 if(realtime)
1658 {
1659 whichimage=2;
1660 /* plot if inside window */
1661 if (inf.col1 >= 0)
1662 {
1663 if(oldcol1 != -1 && connect)
1664 draw_line(inf.col1,inf.row1,oldcol1,oldrow1,color%colors);
1665 else
1666 (*plot)(inf.col1,inf.row1,color%colors);
1667 }
1668 else if (inf.col1 == -2)
1669 return(ret);
1670 oldcol1 = inf.col1;
1671 oldrow1 = inf.row1;
1672 }
1673 }
1674 }
1675 if(fp)
1676 fclose(fp);
1677 return(ret);
1678 }
1679
1680
1681 static int orbit3dfloatcalc(void)
1682 {
1683 FILE *fp;
1684 unsigned long count;
1685 int oldcol,oldrow;
1686 int oldcol1,oldrow1;
1687 int color;
1688 int ret;
1689 struct float3dvtinf inf;
1690
1691 /* setup affine screen coord conversion */
1692 setup_convert_to_screen(&inf.cvt);
1693
1694 oldcol = oldrow = -1;
1695 oldcol1 = oldrow1 = -1;
1696 color = 2;
1697 if(color >= colors)
1698 color = 1;
1699 inf.orbit[0] = initorbitfp[0];
1700 inf.orbit[1] = initorbitfp[1];
1701 inf.orbit[2] = initorbitfp[2];
1702
1703 if(diskvideo) /* this would KILL a disk drive! */
1704 notdiskmsg();
1705
1706 fp = open_orbitsave();
1707
1708 ret = 0;
1709 if(maxit > 0x1fffffL || maxct)
1710 maxct = 0x7fffffffL;
1711 else
1712 maxct = maxit*1024L;
1713 count = coloriter = 0L;
1714 while(coloriter++ <= maxct) /* loop until keypress or maxit */
1715 {
1716 /* calc goes here */
1717 if (++count > 1000)
1718 { /* time to switch colors? */
1719 count = 0;
1720 if (++color >= colors) /* another color to switch to? */
1721 color = 1; /* (don't use the background color) */
1722 }
1723
1724 if(keypressed())
1725 {
1726 mute();
1727 ret = -1;
1728 break;
1729 }
1730
1731 FORBIT(&inf.orbit[0],&inf.orbit[1],&inf.orbit[2]);
1732 if(fp)
1733 fprintf(fp,orbitsave_format,inf.orbit[0],inf.orbit[1],inf.orbit[2]);
1734 if (float3dviewtransf(&inf))
1735 {
1736 /* plot if inside window */
1737 if (inf.col >= 0)
1738 {
1739 if(realtime)
1740 whichimage=1;
1741 if ((soundflag&7) > 1) {
1742 w_snd((int)(inf.viewvect[((soundflag&7) - 2)]*100+basehertz));
1743 }
1744 if(oldcol != -1 && connect)
1745 draw_line(inf.col,inf.row,oldcol,oldrow,color%colors);
1746 else
1747 (*plot)(inf.col,inf.row,color%colors);
1748 }
1749 else if (inf.col == -2)
1750 return(ret);
1751 oldcol = inf.col;
1752 oldrow = inf.row;
1753 if(realtime)
1754 {
1755 whichimage=2;
1756 /* plot if inside window */
1757 if (inf.col1 >= 0)
1758 {
1759 if(oldcol1 != -1 && connect)
1760 draw_line(inf.col1,inf.row1,oldcol1,oldrow1,color%colors);
1761 else
1762 (*plot)(inf.col1,inf.row1,color%colors);
1763 }
1764 else if (inf.col1 == -2)
1765 return(ret);
1766 oldcol1 = inf.col1;
1767 oldrow1 = inf.row1;
1768 }
1769 }
1770 }
1771 if(fp)
1772 fclose(fp);
1773 return(ret);
1774 }
1775
1776 int dynam2dfloatsetup()
1777 {
1778 connect = 0;
1779 euler = 0;
1780 d = param[0]; /* number of intervals */
1781 if (d<0) {
1782 d = -d;
1783 connect = 1;
1784 }
1785 else if (d==0) {
1786 d = 1;
1787 }
1788 if (fractype==DYNAMICFP) {
1789 a = param[2]; /* parameter */
1790 b = param[3]; /* parameter */
1791 dt = param[1]; /* step size */
1792 if (dt<0) {
1793 dt = -dt;
1794 euler = 1;
1795 }
1796 if (dt==0) dt = 0.01;
1797 }
1798 if (outside == SUM) {
1799 plot = plothist;
1800 }
1801 return(1);
1802 }
1803
1804 /*
1805 * This is the routine called to perform a time-discrete dynamical
1806 * system image.
1807 * The starting positions are taken by stepping across the image in steps
1808 * of parameter1 pixels. maxit differential equation steps are taken, with
1809 * a step size of parameter2.
1810 */
1811 int dynam2dfloat()
1812 {
1813 FILE *fp;
1814 double *soundvar = NULL;
1815 double x,y,z;
1816 int color,col,row;
1817 long count;
1818 int oldrow, oldcol;
1819 double *p0,*p1;
1820 struct affine cvt;
1821 int ret;
1822 int xstep, ystep; /* The starting position step number */
1823 double xpixel, ypixel; /* Our pixel position on the screen */
1824
1825 fp = open_orbitsave();
1826 /* setup affine screen coord conversion */
1827 setup_convert_to_screen(&cvt);
1828
1829 p0 = &x;
1830 p1 = &y;
1831
1832
1833 if((soundflag&7)==2)
1834 soundvar = &x;
1835 else if((soundflag&7)==3)
1836 soundvar = &y;
1837 else if((soundflag&7)==4)
1838 soundvar = &z;
1839
1840 count = 0;
1841 if(inside > 0)
1842 color = inside;
1843 if(color >= colors)
1844 color = 1;
1845 oldcol = oldrow = -1;
1846
1847 xstep = -1;
1848 ystep = 0;
1849
1850 if (resuming) {
1851 start_resume();
1852 get_resume(sizeof(count),&count, sizeof(color),&color,
1853 sizeof(oldrow),&oldrow, sizeof(oldcol),&oldcol,
1854 sizeof(x),&x, sizeof(y), &y, sizeof(xstep), &xstep,
1855 sizeof(ystep), &ystep, 0);
1856 end_resume();
1857 }
1858
1859 ret = 0;
1860 for(;;)
1861 {
1862 if(keypressed())
1863 {
1864 mute();
1865 alloc_resume(100,1);
1866 put_resume(sizeof(count),&count, sizeof(color),&color,
1867 sizeof(oldrow),&oldrow, sizeof(oldcol),&oldcol,
1868 sizeof(x),&x, sizeof(y), &y, sizeof(xstep), &xstep,
1869 sizeof(ystep), &ystep, 0);
1870 ret = -1;
1871 break;
1872 }
1873
1874 xstep ++;
1875 if (xstep>=d) {
1876 xstep = 0;
1877 ystep ++;
1878 if (ystep>d) {
1879 mute();
1880 ret = -1;
1881 break;
1882 }
1883 }
1884
1885 xpixel = dxsize*(xstep+.5)/d;
1886 ypixel = dysize*(ystep+.5)/d;
1887 x = (double)((xxmin+delxx*xpixel) + (delxx2*ypixel));
1888 y = (double)((yymax-delyy*ypixel) + (-delyy2*xpixel));
1889 if (fractype==MANDELCLOUD) {
1890 a = x;
1891 b = y;
1892 }
1893 oldcol = -1;
1894
1895 if (++color >= colors) /* another color to switch to? */
1896 color = 1; /* (don't use the background color) */
1897
1898 for (count=0;count<maxit;count++) {
1899
1900 if (count % 2048L == 0)
1901 if (keypressed())
1902 break;
1903
1904 col = (int)(cvt.a*x + cvt.b*y + cvt.e);
1905 row = (int)(cvt.c*x + cvt.d*y + cvt.f);
1906 if ( col >= 0 && col < xdots && row >= 0 && row < ydots )
1907 {
1908 if ((soundflag&7) > 1)
1909 w_snd((int)(*soundvar*100+basehertz));
1910
1911 if (count>=orbit_delay) {
1912 if(oldcol != -1 && connect)
1913 draw_line(col,row,oldcol,oldrow,color%colors);
1914 else if(count > 0 || fractype != MANDELCLOUD)
1915 (*plot)(col,row,color%colors);
1916 }
1917 oldcol = col;
1918 oldrow = row;
1919 }
1920 else if((long)abs(row) + (long)abs(col) > BAD_PIXEL) /* sanity check */
1921 return(ret);
1922 else
1923 oldrow = oldcol = -1;
1924
1925 if(FORBIT(p0, p1, NULL))
1926 break;
1927 if(fp)
1928 fprintf(fp,orbitsave_format,*p0,*p1,0.0);
1929 }
1930 }
1931 if(fp)
1932 fclose(fp);
1933 return(ret);
1934 }
1935
1936 int keep_scrn_coords = 0;
1937 int set_orbit_corners = 0;
1938 long orbit_interval;
1939 double oxmin, oymin, oxmax, oymax, ox3rd, oy3rd;
1940 struct affine o_cvt;
1941 static int o_color;
1942
1943 int setup_orbits_to_screen(struct affine *scrn_cnvt)
1944 {
1945 double det, xd, yd;
1946
1947 if((det = (ox3rd-oxmin)*(oymin-oymax) + (oymax-oy3rd)*(oxmax-oxmin))==0)
1948 return(-1);
1949 xd = dxsize/det;
1950 scrn_cnvt->a = xd*(oymax-oy3rd);
1951 scrn_cnvt->b = xd*(ox3rd-oxmin);
1952 scrn_cnvt->e = -scrn_cnvt->a*oxmin - scrn_cnvt->b*oymax;
1953
1954 if((det = (ox3rd-oxmax)*(oymin-oymax) + (oymin-oy3rd)*(oxmax-oxmin))==0)
1955 return(-1);
1956 yd = dysize/det;
1957 scrn_cnvt->c = yd*(oymin-oy3rd);
1958 scrn_cnvt->d = yd*(ox3rd-oxmax);
1959 scrn_cnvt->f = -scrn_cnvt->c*oxmin - scrn_cnvt->d*oymax;
1960 return(0);
1961 }
1962
1963 int plotorbits2dsetup(void)
1964 {
1965
1966 #ifndef XFRACT
1967 if (curfractalspecific->isinteger != 0) {
1968 int tofloat;
1969 if ((tofloat = curfractalspecific->tofloat) == NOFRACTAL)
1970 return(-1);
1971 floatflag = usr_floatflag = 1; /* force floating point */
1972 curfractalspecific = &fractalspecific[tofloat];
1973 fractype = tofloat;
1974 }
1975 #endif
1976
1977 PER_IMAGE();
1978
1979 /* setup affine screen coord conversion */
1980 if (keep_scrn_coords) {
1981 if(setup_orbits_to_screen(&o_cvt))
1982 return(-1);
1983 } else {
1984 if(setup_convert_to_screen(&o_cvt))
1985 return(-1);
1986 }
1987 /* set so truncation to int rounds to nearest */
1988 o_cvt.e += 0.5;
1989 o_cvt.f += 0.5;
1990
1991 if (orbit_delay >= maxit) /* make sure we get an image */
1992 orbit_delay = (int)(maxit - 1);
1993
1994 o_color = 1;
1995
1996 if (outside == SUM) {
1997 plot = plothist;
1998 }
1999 return(1);
2000 }
2001
2002 int plotorbits2dfloat(void)
2003 {
2004 double *soundvar = NULL;
2005 double x,y,z;
2006 int col,row;
2007 long count;
2008
2009 if(keypressed())
2010 {
2011 mute();
2012 alloc_resume(100,1);
2013 put_resume(sizeof(o_color),&o_color, 0);
2014 return(-1);
2015 }
2016
2017 #if 0
2023 #endif
2024
2025 if((soundflag&7)==2)
2026 soundvar = &x;
2027 else if((soundflag&7)==3)
2028 soundvar = &y;
2029 else if((soundflag&7)==4)
2030 soundvar = &z;
2031
2032 if (resuming) {
2033 start_resume();
2034 get_resume(sizeof(o_color),&o_color, 0);
2035 end_resume();
2036 }
2037
2038 if(inside > 0)
2039 o_color = inside;
2040 else { /* inside <= 0 */
2041 o_color++;
2042 if (o_color >= colors) /* another color to switch to? */
2043 o_color = 1; /* (don't use the background color) */
2044 }
2045
2046 PER_PIXEL(); /* initialize the calculations */
2047
2048 for (count = 0; count < maxit; count++) {
2049
2050 if (ORBITCALC() == 1 && periodicitycheck)
2051 continue; /* bailed out, don't plot */
2052
2053 if (count < orbit_delay || count%orbit_interval)
2054 continue; /* don't plot it */
2055
2056 /* else count >= orbit_delay and we want to plot it */
2057 col = (int)(o_cvt.a*new.x + o_cvt.b*new.y + o_cvt.e);
2058 row = (int)(o_cvt.c*new.x + o_cvt.d*new.y + o_cvt.f);
2059 #ifdef XFRACT
2061 #else
2062 /* don't know why the next line is necessary, the one above should work */
2063 if ( col > 0 && col < xdots && row > 0 && row < ydots )
2064 #endif
2065 { /* plot if on the screen */
2066 if ((soundflag&7) > 1)
2067 w_snd((int)(*soundvar*100+basehertz));
2068
2069 (*plot)(col,row,o_color%colors);
2070 }
2071 else
2072 { /* off screen, don't continue unless periodicity=0 */
2073 if (periodicitycheck)
2074 return(0); /* skip to next pixel */
2075 }
2076 }
2077 return(0);
2078 }
2079
2080 /* this function's only purpose is to manage funnyglasses related */
2081 /* stuff so the code is not duplicated for ifs3d() and lorenz3d() */
2082 int funny_glasses_call(int (*calc)(void))
2083 {
2084 int status;
2085 status = 0;
2086 if(glassestype)
2087 whichimage = 1;
2088 else
2089 whichimage = 0;
2090 plot_setup();
2091 plot = standardplot;
2092 status = calc();
2093 if(realtime && glassestype < 3)
2094 {
2095 realtime = 0;
2096 goto done;
2097 }
2098 if(glassestype && status == 0 && display3d)
2099 {
2100 if(glassestype==3) { /* photographer's mode */
2101 if(active_system == 0) { /* dos version */
2102 int i;
2103 static FCODE firstready[]={"\
2104 First image (left eye) is ready. Hit any key to see it,\n\
2105 then hit <s> to save, hit any other key to create second image."};
2106 stopmsg(16,firstready);
2107 while ((i = getakey()) == 's' || i == 'S') {
2108 diskisactive = 1;
2109 savetodisk(savename);
2110 diskisactive = 0;
2111 }
2112 /* is there a better way to clear the screen in graphics mode? */
2113 setvideomode(videoentry.videomodeax,
2114 videoentry.videomodebx,
2115 videoentry.videomodecx,
2116 videoentry.videomodedx);
2117 }
2118 else { /* Windows version */
2119 static FCODE firstready2[]={"First (Left Eye) image is complete"};
2120 stopmsg(0,firstready2);
2121 clear_screen(0);
2122 }
2123 }
2124 whichimage = 2;
2125 if(curfractalspecific->flags & INFCALC)
2126 curfractalspecific->per_image(); /* reset for 2nd image */
2127 plot_setup();
2128 plot = standardplot;
2129 /* is there a better way to clear the graphics screen ? */
2130 if((status = calc()) != 0)
2131 goto done;
2132 if(glassestype==3) /* photographer's mode */
2133 if(active_system == 0) { /* dos version */
2134 static FCODE secondready[]={"Second image (right eye) is ready"};
2135 stopmsg(16,secondready);
2136 }
2137 }
2138 done:
2139 if(glassestype == 4 && sxdots >= 2*xdots)
2140 {
2141 /* turn off view windows so will save properly */
2142 sxoffs = syoffs = 0;
2143 xdots = sxdots;
2144 ydots = sydots;
2145 viewwindow = 0;
2146 }
2147 return(status);
2148 }
2149
2150 /* double version - mainly for testing */
2151 static int ifs3dfloat(void)
2152 {
2153 int color_method;
2154 FILE *fp;
2155 int color;
2156
2157 double newx,newy,newz,r,sum;
2158
2159 int k;
2160 int ret;
2161
2162 struct float3dvtinf inf;
2163
2164 float far *ffptr;
2165
2166 /* setup affine screen coord conversion */
2167 setup_convert_to_screen(&inf.cvt);
2168 srand(1);
2169 color_method = (int)param[0];
2170 if(diskvideo) /* this would KILL a disk drive! */
2171 notdiskmsg();
2172
2173 inf.orbit[0] = 0;
2174 inf.orbit[1] = 0;
2175 inf.orbit[2] = 0;
2176
2177 fp = open_orbitsave();
2178
2179 ret = 0;
2180 if(maxit > 0x1fffffL)
2181 maxct = 0x7fffffffL;
2182 else
2183 maxct = maxit*1024;
2184 coloriter = 0L;
2185 while(coloriter++ <= maxct) /* loop until keypress or maxit */
2186 {
2187 if( keypressed() ) /* keypress bails out */
2188 {
2189 ret = -1;
2190 break;
2191 }
2192 r = rand(); /* generate a random number between 0 and 1 */
2193 r /= RAND_MAX;
2194
2195 /* pick which iterated function to execute, weighted by probability */
2196 sum = ifs_defn[12]; /* [0][12] */
2197 k = 0;
2198 while ( sum < r && ++k < numaffine*IFS3DPARM)
2199 {
2200 sum += ifs_defn[k*IFS3DPARM+12];
2201 if (ifs_defn[(k+1)*IFS3DPARM+12] == 0) break; /* for safety */
2202 }
2203
2204 /* calculate image of last point under selected iterated function */
2205 ffptr = ifs_defn + k*IFS3DPARM; /* point to first parm in row */
2206 newx = *ffptr * inf.orbit[0] +
2207 *(ffptr+1) * inf.orbit[1] +
2208 *(ffptr+2) * inf.orbit[2] + *(ffptr+9);
2209 newy = *(ffptr+3) * inf.orbit[0] +
2210 *(ffptr+4) * inf.orbit[1] +
2211 *(ffptr+5) * inf.orbit[2] + *(ffptr+10);
2212 newz = *(ffptr+6) * inf.orbit[0] +
2213 *(ffptr+7) * inf.orbit[1] +
2214 *(ffptr+8) * inf.orbit[2] + *(ffptr+11);
2215
2216 inf.orbit[0] = newx;
2217 inf.orbit[1] = newy;
2218 inf.orbit[2] = newz;
2219 if(fp)
2220 fprintf(fp,orbitsave_format,newx,newy,newz);
2221 if (float3dviewtransf(&inf))
2222 {
2223 /* plot if inside window */
2224 if (inf.col >= 0)
2225 {
2226 if(realtime)
2227 whichimage=1;
2228 if(color_method)
2229 color = (k%colors)+1;
2230 else
2231 color = getcolor(inf.col,inf.row)+1;
2232 if( color < colors ) /* color sticks on last value */
2233 (*plot)(inf.col,inf.row,color);
2234 }
2235 else if (inf.col == -2)
2236 return(ret);
2237 if(realtime)
2238 {
2239 whichimage=2;
2240 /* plot if inside window */
2241 if (inf.col1 >= 0)
2242 {
2243 if(color_method)
2244 color = (k%colors)+1;
2245 else
2246 color = getcolor(inf.col1,inf.row1)+1;
2247 if( color < colors ) /* color sticks on last value */
2248 (*plot)(inf.col1,inf.row1,color);
2249 }
2250 else if (inf.col1 == -2)
2251 return(ret);
2252 }
2253 }
2254 } /* end while */
2255 if(fp)
2256 fclose(fp);
2257 return(ret);
2258 }
2259
2260 int ifs() /* front-end for ifs2d and ifs3d */
2261 {
2262 if (ifs_defn == NULL && ifsload() < 0)
2263 return(-1);
2264 if(diskvideo) /* this would KILL a disk drive! */
2265 notdiskmsg();
2266 return((ifs_type == 0) ? ifs2d() : ifs3d());
2267 }
2268
2269
2270 /* IFS logic shamelessly converted to integer math */
2271 static int ifs2d(void)
2272 {
2273 int color_method;
2274 FILE *fp;
2275 int col;
2276 int row;
2277 int color;
2278 int ret;
2279 long far *localifs;
2280 long far *lfptr;
2281 long x,y,newx,newy,r,sum, tempr;
2282
2283 int i,j,k;
2284 struct l_affine cvt;
2285 /* setup affine screen coord conversion */
2286 l_setup_convert_to_screen(&cvt);
2287
2288 srand(1);
2289 color_method = (int)param[0];
2290 if((localifs=(long far *)farmemalloc((long)numaffine*IFSPARM*sizeof(long)))==NULL)
2291 {
2292 stopmsg(0,insufficient_ifs_mem);
2293 return(-1);
2294 }
2295
2296 for (i = 0; i < numaffine; i++) /* fill in the local IFS array */
2297 for (j = 0; j < IFSPARM; j++)
2298 localifs[i*IFSPARM+j] = (long)(ifs_defn[i*IFSPARM+j] * fudge);
2299
2300 tempr = fudge / 32767; /* find the proper rand() fudge */
2301
2302 fp = open_orbitsave();
2303
2304 x = y = 0;
2305 ret = 0;
2306 if(maxit > 0x1fffffL)
2307 maxct = 0x7fffffffL;
2308 else
2309 maxct = maxit*1024L;
2310 coloriter = 0L;
2311 while(coloriter++ <= maxct) /* loop until keypress or maxit */
2312 {
2313 if( keypressed() ) /* keypress bails out */
2314 {
2315 ret = -1;
2316 break;
2317 }
2318 r = rand15(); /* generate fudged random number between 0 and 1 */
2319 r *= tempr;
2320
2321 /* pick which iterated function to execute, weighted by probability */
2322 sum = localifs[6]; /* [0][6] */
2323 k = 0;
2324 while ( sum < r && k < numaffine-1) /* fixed bug of error if sum < 1 */
2325 sum += localifs[++k*IFSPARM+6];
2326 /* calculate image of last point under selected iterated function */
2327 lfptr = localifs + k*IFSPARM; /* point to first parm in row */
2328 newx = multiply(lfptr[0],x,bitshift) +
2329 multiply(lfptr[1],y,bitshift) + lfptr[4];
2330 newy = multiply(lfptr[2],x,bitshift) +
2331 multiply(lfptr[3],y,bitshift) + lfptr[5];
2332 x = newx;
2333 y = newy;
2334 if(fp)
2335 fprintf(fp,orbitsave_format,(double)newx/fudge,(double)newy/fudge,0.0);
2336
2337 /* plot if inside window */
2338 col = (int)((multiply(cvt.a,x,bitshift) + multiply(cvt.b,y,bitshift) + cvt.e) >> bitshift);
2339 row = (int)((multiply(cvt.c,x,bitshift) + multiply(cvt.d,y,bitshift) + cvt.f) >> bitshift);
2340 if ( col >= 0 && col < xdots && row >= 0 && row < ydots )
2341 {
2342 /* color is count of hits on this pixel */
2343 if(color_method)
2344 color = (k%colors)+1;
2345 else
2346 color = getcolor(col,row)+1;
2347 if( color < colors ) /* color sticks on last value */
2348 (*plot)(col,row,color);
2349 }
2350 else if((long)abs(row) + (long)abs(col) > BAD_PIXEL) /* sanity check */
2351 return(ret);
2352 }
2353 if(fp)
2354 fclose(fp);
2355 farmemfree(localifs);
2356 return(ret);
2357 }
2358
2359 static int ifs3dlong(void)
2360 {
2361 int color_method;
2362 FILE *fp;
2363 int color;
2364 int ret;
2365
2366 long far *localifs;
2367 long far *lfptr;
2368 long newx,newy,newz,r,sum, tempr;
2369
2370 int i,j,k;
2371
2372 struct long3dvtinf inf;
2373 srand(1);
2374 color_method = (int)param[0];
2375 if((localifs=(long far *)farmemalloc((long)numaffine*IFS3DPARM*sizeof(long)))==NULL)
2376 {
2377 stopmsg(0,insufficient_ifs_mem);
2378 return(-1);
2379 }
2380
2381 /* setup affine screen coord conversion */
2382 l_setup_convert_to_screen(&inf.cvt);
2383
2384 for (i = 0; i < numaffine; i++) /* fill in the local IFS array */
2385 for (j = 0; j < IFS3DPARM; j++)
2386 localifs[i*IFS3DPARM+j] = (long)(ifs_defn[i*IFS3DPARM+j] * fudge);
2387
2388 tempr = fudge / 32767; /* find the proper rand() fudge */
2389
2390 inf.orbit[0] = 0;
2391 inf.orbit[1] = 0;
2392 inf.orbit[2] = 0;
2393
2394 fp = open_orbitsave();
2395
2396 ret = 0;
2397 if(maxit > 0x1fffffL)
2398 maxct = 0x7fffffffL;
2399 else
2400 maxct = maxit*1024L;
2401 coloriter = 0L;
2402 while(coloriter++ <= maxct) /* loop until keypress or maxit */
2403 {
2404 if( keypressed() ) /* keypress bails out */
2405 {
2406 ret = -1;
2407 break;
2408 }
2409 r = rand15(); /* generate fudged random number between 0 and 1 */
2410 r *= tempr;
2411
2412 /* pick which iterated function to execute, weighted by probability */
2413 sum = localifs[12]; /* [0][12] */
2414 k = 0;
2415 while ( sum < r && ++k < numaffine*IFS3DPARM)
2416 {
2417 sum += localifs[k*IFS3DPARM+12];
2418 if (ifs_defn[(k+1)*IFS3DPARM+12] == 0) break; /* for safety */
2419 }
2420
2421 /* calculate image of last point under selected iterated function */
2422 lfptr = localifs + k*IFS3DPARM; /* point to first parm in row */
2423
2424 /* calculate image of last point under selected iterated function */
2425 newx = multiply(lfptr[0], inf.orbit[0], bitshift) +
2426 multiply(lfptr[1], inf.orbit[1], bitshift) +
2427 multiply(lfptr[2], inf.orbit[2], bitshift) + lfptr[9];
2428 newy = multiply(lfptr[3], inf.orbit[0], bitshift) +
2429 multiply(lfptr[4], inf.orbit[1], bitshift) +
2430 multiply(lfptr[5], inf.orbit[2], bitshift) + lfptr[10];
2431 newz = multiply(lfptr[6], inf.orbit[0], bitshift) +
2432 multiply(lfptr[7], inf.orbit[1], bitshift) +
2433 multiply(lfptr[8], inf.orbit[2], bitshift) + lfptr[11];
2434
2435 inf.orbit[0] = newx;
2436 inf.orbit[1] = newy;
2437 inf.orbit[2] = newz;
2438 if(fp)
2439 fprintf(fp,orbitsave_format,(double)newx/fudge,(double)newy/fudge,(double)newz/fudge);
2440
2441 if (long3dviewtransf(&inf))
2442 {
2443 if((long)abs(inf.row) + (long)abs(inf.col) > BAD_PIXEL) /* sanity check */
2444 return(ret);
2445 /* plot if inside window */
2446 if (inf.col >= 0)
2447 {
2448 if(realtime)
2449 whichimage=1;
2450 if(color_method)
2451 color = (k%colors)+1;
2452 else
2453 color = getcolor(inf.col,inf.row)+1;
2454 if( color < colors ) /* color sticks on last value */
2455 (*plot)(inf.col,inf.row,color);
2456 }
2457 if(realtime)
2458 {
2459 whichimage=2;
2460 /* plot if inside window */
2461 if (inf.col1 >= 0)
2462 {
2463 if(color_method)
2464 color = (k%colors)+1;
2465 else
2466 color = getcolor(inf.col1,inf.row1)+1;
2467 if( color < colors ) /* color sticks on last value */
2468 (*plot)(inf.col1,inf.row1,color);
2469 }
2470 }
2471 }
2472 }
2473 if(fp)
2474 fclose(fp);
2475 farmemfree(localifs);
2476 return(ret);
2477 }
2478
2479 static void setupmatrix(MATRIX doublemat)
2480 {
2481 /* build transformation matrix */
2482 identity (doublemat);
2483
2484 /* apply rotations - uses the same rotation variables as line3d.c */
2485 xrot ((double)XROT / 57.29577,doublemat);
2486 yrot ((double)YROT / 57.29577,doublemat);
2487 zrot ((double)ZROT / 57.29577,doublemat);
2488
2489 /* apply scale */
2490 /* scale((double)XSCALE/100.0,(double)YSCALE/100.0,(double)ROUGH/100.0,doublemat);*/
2491
2492 }
2493
2494 int orbit3dfloat()
2495 {
2496 display3d = -1;
2497 if(0 < glassestype && glassestype < 3)
2498 realtime = 1;
2499 else
2500 realtime = 0;
2501 return(funny_glasses_call(orbit3dfloatcalc));
2502 }
2503
2504 int orbit3dlong()
2505 {
2506 display3d = -1;
2507 if(0 < glassestype && glassestype < 3)
2508 realtime = 1;
2509 else
2510 realtime = 0;
2511 return(funny_glasses_call(orbit3dlongcalc));
2512 }
2513
2514 static int ifs3d(void)
2515 {
2516 display3d = -1;
2517
2518 if(0 < glassestype && glassestype < 3)
2519 realtime = 1;
2520 else
2521 realtime = 0;
2522 if(floatflag)
2523 return(funny_glasses_call(ifs3dfloat)); /* double version of ifs3d */
2524 else
2525 return(funny_glasses_call(ifs3dlong)); /* long version of ifs3d */
2526 }
2527
2528
2529
2530 static int long3dviewtransf(struct long3dvtinf *inf)
2531 {
2532 int i,j;
2533 double tmpx, tmpy, tmpz;
2534 long tmp;
2535
2536 if (coloriter == 1) /* initialize on first call */
2537 {
2538 for(i=0;i<3;i++)
2539 {
2540 inf->minvals[i] = 1L << 30;
2541 inf->maxvals[i] = -inf->minvals[i];
2542 }
2543 setupmatrix(inf->doublemat);
2544 if(realtime)
2545 setupmatrix(inf->doublemat1);
2546 /* copy xform matrix to long for for fixed point math */
2547 for (i = 0; i < 4; i++)
2548 for (j = 0; j < 4; j++)
2549 {
2550 inf->longmat[i][j] = (long)(inf->doublemat[i][j] * fudge);
2551 if(realtime)
2552 inf->longmat1[i][j] = (long)(inf->doublemat1[i][j] * fudge);
2553 }
2554 }
2555
2556 /* 3D VIEWING TRANSFORM */
2557 longvmult(inf->orbit,inf->longmat,inf->viewvect,bitshift);
2558 if(realtime)
2559 longvmult(inf->orbit,inf->longmat1,inf->viewvect1,bitshift);
2560
2561 if(coloriter <= waste) /* waste this many points to find minz and maxz */
2562 {
2563 /* find minz and maxz */
2564 for(i=0;i<3;i++)
2565 if ((tmp = inf->viewvect[i]) < inf->minvals[i])
2566 inf->minvals[i] = tmp;
2567 else if (tmp > inf->maxvals[i])
2568 inf->maxvals[i] = tmp;
2569
2570 if(coloriter == waste) /* time to work it out */
2571 {
2572 inf->iview[0] = inf->iview[1] = 0L; /* center viewer on origin */
2573
2574 /* z value of user's eye - should be more negative than extreme
2575 negative part of image */
2576 inf->iview[2] = (long)((inf->minvals[2]-inf->maxvals[2])*(double)ZVIEWER/100.0);
2577
2578 /* center image on origin */
2579 tmpx = (-inf->minvals[0]-inf->maxvals[0])/(2.0*fudge); /* center x */
2580 tmpy = (-inf->minvals[1]-inf->maxvals[1])/(2.0*fudge); /* center y */
2581
2582 /* apply perspective shift */
2583 tmpx += ((double)xshift*(xxmax-xxmin))/(xdots);
2584 tmpy += ((double)yshift*(yymax-yymin))/(ydots);
2585 tmpz = -((double)inf->maxvals[2]) / fudge;
2586 trans(tmpx,tmpy,tmpz,inf->doublemat);
2587
2588 if(realtime)
2589 {
2590 /* center image on origin */
2591 tmpx = (-inf->minvals[0]-inf->maxvals[0])/(2.0*fudge); /* center x */
2592 tmpy = (-inf->minvals[1]-inf->maxvals[1])/(2.0*fudge); /* center y */
2593
2594 tmpx += ((double)xshift1*(xxmax-xxmin))/(xdots);
2595 tmpy += ((double)yshift1*(yymax-yymin))/(ydots);
2596 tmpz = -((double)inf->maxvals[2]) / fudge;
2597 trans(tmpx,tmpy,tmpz,inf->doublemat1);
2598 }
2599 for(i=0;i<3;i++)
2600 view[i] = (double)inf->iview[i] / fudge;
2601
2602 /* copy xform matrix to long for for fixed point math */
2603 for (i = 0; i < 4; i++)
2604 for (j = 0; j < 4; j++)
2605 {
2606 inf->longmat[i][j] = (long)(inf->doublemat[i][j] * fudge);
2607 if(realtime)
2608 inf->longmat1[i][j] = (long)(inf->doublemat1[i][j] * fudge);
2609 }
2610 }
2611 return(0);
2612 }
2613
2614 /* apply perspective if requested */
2615 if(ZVIEWER)
2616 {
2617 if(debugflag==22 || ZVIEWER < 100) /* use float for small persp */
2618 {
2619 /* use float perspective calc */
2620 VECTOR tmpv;
2621 for(i=0;i<3;i++)
2622 tmpv[i] = (double)inf->viewvect[i] / fudge;
2623 perspective(tmpv);
2624 for(i=0;i<3;i++)
2625 inf->viewvect[i] = (long)(tmpv[i]*fudge);
2626 if(realtime)
2627 {
2628 for(i=0;i<3;i++)
2629 tmpv[i] = (double)inf->viewvect1[i] / fudge;
2630 perspective(tmpv);
2631 for(i=0;i<3;i++)
2632 inf->viewvect1[i] = (long)(tmpv[i]*fudge);
2633 }
2634 }
2635 else
2636 {
2637 longpersp(inf->viewvect,inf->iview,bitshift);
2638 if(realtime)
2639 longpersp(inf->viewvect1,inf->iview,bitshift);
2640 }
2641 }
2642
2643 /* work out the screen positions */
2644 inf->row = (int)(((multiply(inf->cvt.c,inf->viewvect[0],bitshift) +
2645 multiply(inf->cvt.d,inf->viewvect[1],bitshift) + inf->cvt.f)
2646 >> bitshift)
2647 + yyadjust);
2648 inf->col = (int)(((multiply(inf->cvt.a,inf->viewvect[0],bitshift) +
2649 multiply(inf->cvt.b,inf->viewvect[1],bitshift) + inf->cvt.e)
2650 >> bitshift)
2651 + xxadjust);
2652 if (inf->col < 0 || inf->col >= xdots || inf->row < 0 || inf->row >= ydots)
2653 {
2654 if((long)abs(inf->col)+(long)abs(inf->row) > BAD_PIXEL)
2655 inf->col= inf->row = -2;
2656 else
2657 inf->col= inf->row = -1;
2658 }
2659 if(realtime)
2660 {
2661 inf->row1 = (int)(((multiply(inf->cvt.c,inf->viewvect1[0],bitshift) +
2662 multiply(inf->cvt.d,inf->viewvect1[1],bitshift) +
2663 inf->cvt.f) >> bitshift)
2664 + yyadjust1);
2665 inf->col1 = (int)(((multiply(inf->cvt.a,inf->viewvect1[0],bitshift) +
2666 multiply(inf->cvt.b,inf->viewvect1[1],bitshift) +
2667 inf->cvt.e) >> bitshift)
2668 + xxadjust1);
2669 if (inf->col1 < 0 || inf->col1 >= xdots || inf->row1 < 0 || inf->row1 >= ydots)
2670 {
2671 if((long)abs(inf->col1)+(long)abs(inf->row1) > BAD_PIXEL)
2672 inf->col1= inf->row1 = -2;
2673 else
2674 inf->col1= inf->row1 = -1;
2675 }
2676 }
2677 return(1);
2678 }
2679
2680 static int float3dviewtransf(struct float3dvtinf *inf)
2681 {
2682 int i;
2683 double tmpx, tmpy, tmpz;
2684 double tmp;
2685
2686 if (coloriter == 1) /* initialize on first call */
2687 {
2688 for(i=0;i<3;i++)
2689 {
2690 inf->minvals[i] = 100000.0; /* impossible value */
2691 inf->maxvals[i] = -100000.0;
2692 }
2693 setupmatrix(inf->doublemat);
2694 if(realtime)
2695 setupmatrix(inf->doublemat1);
2696 }
2697
2698 /* 3D VIEWING TRANSFORM */
2699 vmult(inf->orbit,inf->doublemat,inf->viewvect );
2700 if(realtime)
2701 vmult(inf->orbit,inf->doublemat1,inf->viewvect1);
2702
2703 if(coloriter <= waste) /* waste this many points to find minz and maxz */
2704 {
2705 /* find minz and maxz */
2706 for(i=0;i<3;i++)
2707 if ((tmp = inf->viewvect[i]) < inf->minvals[i])
2708 inf->minvals[i] = tmp;
2709 else if (tmp > inf->maxvals[i])
2710 inf->maxvals[i] = tmp;
2711 if(coloriter == waste) /* time to work it out */
2712 {
2713 view[0] = view[1] = 0; /* center on origin */
2714 /* z value of user's eye - should be more negative than extreme
2715 negative part of image */
2716 view[2] = (inf->minvals[2]-inf->maxvals[2])*(double)ZVIEWER/100.0;
2717
2718 /* center image on origin */
2719 tmpx = (-inf->minvals[0]-inf->maxvals[0])/(2.0); /* center x */
2720 tmpy = (-inf->minvals[1]-inf->maxvals[1])/(2.0); /* center y */
2721
2722 /* apply perspective shift */
2723 tmpx += ((double)xshift*(xxmax-xxmin))/(xdots);
2724 tmpy += ((double)yshift*(yymax-yymin))/(ydots);
2725 tmpz = -(inf->maxvals[2]);
2726 trans(tmpx,tmpy,tmpz,inf->doublemat);
2727
2728 if(realtime)
2729 {
2730 /* center image on origin */
2731 tmpx = (-inf->minvals[0]-inf->maxvals[0])/(2.0); /* center x */
2732 tmpy = (-inf->minvals[1]-inf->maxvals[1])/(2.0); /* center y */
2733
2734 tmpx += ((double)xshift1*(xxmax-xxmin))/(xdots);
2735 tmpy += ((double)yshift1*(yymax-yymin))/(ydots);
2736 tmpz = -(inf->maxvals[2]);
2737 trans(tmpx,tmpy,tmpz,inf->doublemat1);
2738 }
2739 }
2740 return(0);
2741 }
2742
2743 /* apply perspective if requested */
2744 if(ZVIEWER)
2745 {
2746 perspective(inf->viewvect);
2747 if(realtime)
2748 perspective(inf->viewvect1);
2749 }
2750 inf->row = (int)(inf->cvt.c*inf->viewvect[0] + inf->cvt.d*inf->viewvect[1]
2751 + inf->cvt.f + yyadjust);
2752 inf->col = (int)(inf->cvt.a*inf->viewvect[0] + inf->cvt.b*inf->viewvect[1]
2753 + inf->cvt.e + xxadjust);
2754 if (inf->col < 0 || inf->col >= xdots || inf->row < 0 || inf->row >= ydots)
2755 {
2756 if((long)abs(inf->col)+(long)abs(inf->row) > BAD_PIXEL)
2757 inf->col= inf->row = -2;
2758 else
2759 inf->col= inf->row = -1;
2760 }
2761 if(realtime)
2762 {
2763 inf->row1 = (int)(inf->cvt.c*inf->viewvect1[0] + inf->cvt.d*inf->viewvect1[1]
2764 + inf->cvt.f + yyadjust1);
2765 inf->col1 = (int)(inf->cvt.a*inf->viewvect1[0] + inf->cvt.b*inf->viewvect1[1]
2766 + inf->cvt.e + xxadjust1);
2767 if (inf->col1 < 0 || inf->col1 >= xdots || inf->row1 < 0 || inf->row1 >= ydots)
2768 {
2769 if((long)abs(inf->col1)+(long)abs(inf->row1) > BAD_PIXEL)
2770 inf->col1= inf->row1 = -2;
2771 else
2772 inf->col1= inf->row1 = -1;
2773 }
2774 }
2775 return(1);
2776 }
2777
2778 static FILE *open_orbitsave(void)
2779 {
2780 FILE *fp;
2781 if ((orbitsave&1) && (fp = fopen(orbitsavename,"w")) != NULL)
2782 {
2783 fprintf(fp,"pointlist x y z color\n");
2784 return fp;
2785 }
2786 return NULL;
2787 }
2788
2789 /* Plot a histogram by incrementing the pixel each time it it touched */
2790 static void _fastcall plothist(int x, int y, int color)
2791 {
2792 color = getcolor(x,y)+1;
2793 if (color >= colors)
2794 color = 1;
2795 putcolor(x,y,color);
2796 }
2797
2798