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