lorenz.c

frasr200的win 版本源码(18.21),使用make文件,使用的vc版本较低,在我的环境下编译有问题! 很不错的分形程序代码!

/*
   This file contains two 3 dimensional orbit-type fractal
   generators - IFS and LORENZ3D, along with code to generate
   red/blue 3D images. Tim Wegner
*/

#include <stdio.h>
#include <stdlib.h>
#include <float.h>
#include <math.h>
#include <string.h>
#include "mpmath.h"
#include "fractint.h"
#include "fractype.h"
#include "prototyp.h"

#define RANDOM(x)  (rand()%(x))

/* BAD_PIXEL is used to cutoff orbits that are diverging. It might be better
to test the actual floating point prbit values, but this seems safe for now.
A higher value cannot be used - to test, turn off math coprocessor and
use +2.24 for type ICONS. If BAD_PIXEL is set to 20000, this will abort
Fractint with a math error. Note that this approach precludes zooming in very
far to an orbit type. */

#define BAD_PIXEL  10000L    /* pixels can't get this big */

struct affine
{
   /* weird order so a,b,e and c,d,f are vectors */
   double a;
   double b;
   double e;
   double c;
   double d;
   double f;
};
struct l_affine
{
   /* weird order so a,b,e and c,d,f are vectors */
   long a;
   long b;
   long e;
   long c;
   long d;
   long f;
};
struct long3dvtinf /* data used by 3d view transform subroutine */
{
   long ct;		/* iteration counter */
   long orbit[3];	/* interated function orbit value */
   long iview[3];	/* perspective viewer's coordinates */
   long viewvect[3];	/* orbit transformed for viewing */
   long viewvect1[3];	/* orbit transformed for viewing */
   long maxvals[3];
   long minvals[3];
   MATRIX doublemat;	/* transformation matrix */
   MATRIX doublemat1;	/* transformation matrix */
   long longmat[4][4];	/* long version of matrix */
   long longmat1[4][4]; /* long version of matrix */
   int row,col; 	/* results */
   int row1,col1;
   struct l_affine cvt;
};
struct float3dvtinf /* data used by 3d view transform subroutine */
{
   long ct;		/* iteration counter */
   double orbit[3];		   /* interated function orbit value */
   double viewvect[3];	      /* orbit transformed for viewing */
   double viewvect1[3];        /* orbit transformed for viewing */
   double maxvals[3];
   double minvals[3];
   MATRIX doublemat;	/* transformation matrix */
   MATRIX doublemat1;	/* transformation matrix */
   int row,col; 	/* results */
   int row1,col1;
   struct affine cvt;
};

/* Routines in this module	*/

static int  ifs3dlong(void);
static int  ifs3dfloat(void);
static double determinant(double mat[3][3]);
static int  solve3x3(double mat[3][3],double vec[3],double ans[3]);
extern int  setup_convert_to_screen(struct affine *);
static int  l_setup_convert_to_screen(struct l_affine *);
static void setupmatrix(MATRIX);
static int  long3dviewtransf(struct long3dvtinf *inf);
static int  float3dviewtransf(struct float3dvtinf *inf);
static FILE *open_orbitsave();
static void _fastcall plothist(int x, int y, int color);
extern char far insufficient_ifs_mem[];
extern int numaffine;
static int realtime;
extern int orbitsave;
static char orbitsavename[]    = {"orbits.raw"};
static char orbitsave_format[] = {"%g %g %g 15\n"};
extern int active_system;
extern int overflow;
extern int soundflag;
extern int basehertz;
extern int fractype;
extern int glassestype;
extern int whichimage;
extern int init3d[];
extern char floatflag;
extern VECTOR view;
extern int xxadjust, yyadjust;
extern int xxadjust1, yyadjust1;
extern int xshift,yshift;
extern int xshift1,yshift1;
extern int	debugflag;			/* for debugging purposes */
extern int	xdots, ydots;		/* coordinates of dots on the screen  */
extern int	maxit;				/* try this many iterations */
extern double param[];
extern double	xxmin,xxmax,yymin,yymax,xx3rd,yy3rd; /* selected screen corners  */
extern	int	diskvideo;			/* for disk-video klooges */
extern int	bitshift;			/* bit shift for fudge */
extern long	fudge;				/* fudge factor (2**n) */
extern int	colors; 			/* maximum colors available */
extern int display3d;
extern double dxsize,dysize,delxx,delxx2,delyy,delyy2;

extern float far *ifs_defn;
extern int ifs_type;

static int t;
static long l_dx,l_dy,l_dz,l_dt,l_a,l_b,l_c,l_d;
static long l_adt,l_bdt,l_cdt,l_xdt,l_ydt;
static long l_initx,l_inity,l_initz;
static long initorbitlong[3];

static double dx,dy,dz,dt,a,b,c,d;
static double adt,bdt,cdt,xdt,ydt,zdt;
static double initx,inity,initz;
static double initorbit[3];
extern int inside;
extern int outside;
extern int orbit_delay;

extern void (_fastcall * standardplot)(int,int,int);
extern int  calc_status, resuming;
extern int  diskisactive;
extern char savename[];

/*
 * The following declarations used for Inverse Julia.  MVS
 */

extern int    Init_Queue      (unsigned long);
extern void   Free_Queue      (void);
extern void   ClearQueue      (void);
extern int    QueueEmpty      (void);
extern int    QueueFull       (void);
extern int    QueueFullAlmost (void);
extern int    PushLong        (long ,  long);
extern int    PushFloat       (float,  float);
extern int    EnQueueLong     (long ,  long);
extern int    EnQueueFloat    (float,  float);
extern LCMPLX PopLong         (void);
extern _CMPLX PopFloat        (void);
extern LCMPLX DeQueueLong     (void);
extern _CMPLX DeQueueFloat    (void);
extern LCMPLX ComplexSqrtLong (long ,  long);
extern _CMPLX ComplexSqrtFloat(double, double);
static char far NoQueue[] =
  "Not enough memory: switching to random walk.\n";

static int    maxhits;
int    run_length;
enum   {breadth_first, depth_first, random_walk, random_run} major_method;
enum   {left_first, right_first}                             minor_method;
struct affine cvt;
struct l_affine lcvt;

extern _CMPLX  old,  new;
extern LCMPLX lold, lnew;

double Cx, Cy;
long   CxLong, CyLong;

/*
 * end of Inverse Julia declarations;
 */

/* these are potential user parameters */
int connect = 1;    /* flag to connect points with a line */
int euler = 0;	    /* use implicit euler approximation for dynamic system */
int waste = 100;    /* waste this many points before plotting */
int projection = 2; /* projection plane - default is to plot x-y */

/******************************************************************/
/*		   zoom box conversion functions		  */
/******************************************************************/

static double determinant(mat) /* determinant of 3x3 matrix */
double mat[3][3];
{
   /* calculate determinant of 3x3 matrix */
   return(mat[0][0]*mat[1][1]*mat[2][2] +
	  mat[0][2]*mat[1][0]*mat[2][1] +
	  mat[0][1]*mat[1][2]*mat[2][0] -
	  mat[2][0]*mat[1][1]*mat[0][2] -
	  mat[1][0]*mat[0][1]*mat[2][2] -
	  mat[0][0]*mat[1][2]*mat[2][1]);

}

static int solve3x3(mat,vec,ans) /* solve 3x3 inhomogeneous linear equations */
double mat[3][3], vec[3], ans[3];
{
   /* solve 3x3 linear equation [mat][ans] = [vec] */
   double denom;
   double tmp[3][3];
   int i;
   denom = determinant(mat);
   if(fabs(denom) < DBL_EPSILON) /* test if can solve */
     return(-1);
   memcpy(tmp,mat,sizeof(double)*9);
   for(i=0;i<3;i++)
   {
      tmp[0][i] = vec[0];
      tmp[1][i] = vec[1];
      tmp[2][i] = vec[2];
      ans[i]  =  determinant(tmp)/denom;
      tmp[0][i] = mat[0][i];
      tmp[1][i] = mat[1][i];
      tmp[2][i] = mat[2][i];
    }
    return(0);
}


/* Conversion of complex plane to screen coordinates for rotating zoom box.
   Assume there is an affine transformation mapping complex zoom parallelogram
   to rectangular screen. We know this map must map parallelogram corners to
   screen corners, so we have following equations:

      a*xxmin+b*yymax+e == 0	    (upper left)
      c*xxmin+d*yymax+f == 0

      a*xx3rd+b*yy3rd+e == 0	    (lower left)
      c*xx3rd+d*yy3rd+f == ydots-1

      a*xxmax+b*yymin+e == xdots-1  (lower right)
      c*xxmax+d*yymin+f == ydots-1

      First we must solve for a,b,c,d,e,f - (which we do once per image),
      then we just apply the transformation to each orbit value.
*/

int setup_convert_to_screen(struct affine *scrn_cnvt)
{
   /* we do this twice - rather than having six equations with six unknowns,
      everything partitions to two sets of three equations with three
      unknowns. Nice, because who wants to calculate a 6x6 determinant??
    */
   double mat[3][3];
   double vec[3];
   /*
      first these equations - solve for a,b,e
      a*xxmin+b*yymax+e == 0	    (upper left)
      a*xx3rd+b*yy3rd+e == 0	    (lower left)
      a*xxmax+b*yymin+e == xdots-1  (lower right)
   */
   mat[0][0] = xxmin;
   mat[0][1] = yymax;
   mat[0][2] = 1.0;
   mat[1][0] = xx3rd;
   mat[1][1] = yy3rd;
   mat[1][2] = 1.0;
   mat[2][0] = xxmax;
   mat[2][1] = yymin;
   mat[2][2] = 1.0;
   vec[0]    = 0.0;
   vec[1]    = 0.0;
   vec[2]    = (double)(xdots-1);

   if(solve3x3(mat,vec, &(scrn_cnvt->a)))
      return(-1);
   /*
      now solve these:
      c*xxmin+d*yymax+f == 0
      c*xx3rd+d*yy3rd+f == ydots-1
      c*xxmax+d*yymin+f == ydots-1
      (mat[][] has not changed - only vec[])
   */
   vec[0]    = 0.0;
   vec[1]    = (double)(ydots-1);
   vec[2]    = (double)(ydots-1);

   if(solve3x3(mat,vec, &scrn_cnvt->c))
      return(-1);
   return(0);
}

static int l_setup_convert_to_screen(struct l_affine *l_cvt)
{
   struct affine cvt;

   /* MCP 7-7-91, This function should return a something! */
   if(setup_convert_to_screen(&cvt))
      return(-1);
   l_cvt->a = cvt.a*fudge; l_cvt->b = cvt.b*fudge; l_cvt->c = cvt.c*fudge;
   l_cvt->d = cvt.d*fudge; l_cvt->e = cvt.e*fudge; l_cvt->f = cvt.f*fudge;

   /* MCP 7-7-91 */
   return(0);
}


/******************************************************************/
/*   setup functions - put in fractalspecific[fractype].per_image */
/******************************************************************/

static double orbit;
static long   l_orbit;

extern double sinx,cosx;
static long l_sinx,l_cosx;

int orbit3dlongsetup()
{
   connect = 1;
   waste = 100;
   projection = 2;
   if (fractype==LHENON || fractype==KAM || fractype==KAM3D || 
       fractype==INVERSEJULIA)
      connect=0;
   if(fractype==LROSSLER)
      waste = 500;
   if(fractype==LLORENZ)
      projection = 1;

   initorbitlong[0] = fudge;  /* initial conditions */
   initorbitlong[1] = fudge;
   initorbitlong[2] = fudge;

   if(fractype==LHENON)
   {
      l_a =  param[0]*fudge;
      l_b =  param[1]*fudge;
      l_c =  param[2]*fudge;
      l_d =  param[3]*fudge;
   }
   else if(fractype==KAM || fractype==KAM3D)
   {
      a   = param[0];		/* angle */
      if(param[1] <= 0.0)
	 param[1] = .01;
      l_b =  param[1]*fudge;	/* stepsize */
      l_c =  param[2]*fudge;	/* stop */
      t = l_d =  param[3];     /* points per orbit */

      l_sinx = sin(a)*fudge;
      l_cosx = cos(a)*fudge;
      l_orbit = 0;
      initorbitlong[0] = initorbitlong[1] = initorbitlong[2] = 0;
   }
   else if (fractype == INVERSEJULIA)
   {
      LCMPLX Sqrt;

      CxLong = param[0] * fudge;
      CyLong = param[1] * fudge;

      maxhits    = (int) param[2];
      run_length = (int) param[3];
      if (maxhits == 0)
	  maxhits = 1;
      else if (maxhits >= colors)
	  maxhits = colors - 1;

      setup_convert_to_screen(&cvt);
      /* Note: using bitshift of 21 for affine, 24 otherwise */

      lcvt.a = cvt.a * (1L << 21);
      lcvt.b = cvt.b * (1L << 21);
      lcvt.c = cvt.c * (1L << 21);
      lcvt.d = cvt.d * (1L << 21);
      lcvt.e = cvt.e * (1L << 21);
      lcvt.f = cvt.f * (1L << 21);

      Sqrt = ComplexSqrtLong(fudge - 4 * CxLong, -4 * CyLong);

      switch (major_method) {
	 case breadth_first:
	    if (Init_Queue((long)32*1024) == 0)
	    { /* can't get queue memory: fall back to random walk */
	       stopmsg(20, NoQueue);
	       major_method = random_walk;
	       goto lrwalk;
	    }
	    EnQueueLong((fudge + Sqrt.x) / 2,  Sqrt.y / 2);
	    EnQueueLong((fudge - Sqrt.x) / 2, -Sqrt.y / 2);
	    break;
	 case depth_first:
	    if (Init_Queue((long)32*1024) == 0)
	    { /* can't get queue memory: fall back to random walk */
	       stopmsg(20, NoQueue);
	       major_method = random_walk;
	       goto lrwalk;
	    }
	    switch (minor_method) {
	       case left_first:
		  PushLong((fudge + Sqrt.x) / 2,  Sqrt.y / 2);
		  PushLong((fudge - Sqrt.x) / 2, -Sqrt.y / 2);
		  break;
	       case right_first:
		  PushLong((fudge - Sqrt.x) / 2, -Sqrt.y / 2);
		  PushLong((fudge + Sqrt.x) / 2,  Sqrt.y / 2);
		  break;
	    }
	    break;
	 case random_walk:
lrwalk:
	    lnew.x = initorbitlong[0] = fudge + Sqrt.x / 2;
	    lnew.y = initorbitlong[1] =         Sqrt.y / 2;
	    break;
	 case random_run:
	    lnew.x = initorbitlong[0] = fudge + Sqrt.x / 2;
	    lnew.y = initorbitlong[1] =         Sqrt.y / 2;
	    break;
      }
   }
   else
   {
      l_dt = param[0]*fudge;
      l_a =  param[1]*fudge;
      l_b =  param[2]*fudge;
      l_c =  param[3]*fudge;
   }

   /* precalculations for speed */
   l_adt = multiply(l_a,l_dt,bitshift);
   l_bdt = multiply(l_b,l_dt,bitshift);
   l_cdt = multiply(l_c,l_dt,bitshift);
   return(1);
}

int orbit3dfloatsetup()
{
   connect = 1;
   waste = 100;
   projection = 2;

   if(fractype==FPHENON || fractype==FPPICKOVER || fractype==FPGINGERBREAD
	    || fractype == KAMFP || fractype == KAM3DFP
	    || fractype == FPHOPALONG || fractype == INVERSEJULIAFP)
      connect=0;
   if(fractype==FPLORENZ3D1 || fractype==FPLORENZ3D3 ||
      fractype==FPLORENZ3D4)
      waste = 750;
   if(fractype==FPROSSLER)
      waste = 500;
   if(fractype==FPLORENZ)
      projection = 1; /* plot x and z */

   initorbit[0] = 1;  /* initial conditions */
   initorbit[1] = 1;
   initorbit[2] = 1;
   if(fractype==FPGINGERBREAD)
   {
      initorbit[0] = param[0];	/* initial conditions */
      initorbit[1] = param[1];
   }

   if(fractype==ICON || fractype==ICON3D)        /* DMF */
   {
      initorbit[0] = 0.01;  /* initial conditions */
      initorbit[1] = 0.003;
      connect = 0;
      waste = 2000;
   }

   if(fractype==FPHENON || fractype==FPPICKOVER)
   {
      a =  param[0];
      b =  param[1];
      c =  param[2];
      d =  param[3];
   }
   else if(fractype==ICON || fractype==ICON3D)        /* DMF */
   {
      initorbit[0] = 0.01;  /* initial conditions */
      initorbit[1] = 0.003;
      connect = 0;
      waste = 2000;
      /* Initialize parameters */
      a  =   param[0];
      b  =   param[1];
      c  =   param[2];
      d  =   param[3];
   }
   else if(fractype==KAMFP || fractype==KAM3DFP)