lorenz.c

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

   {
      a = param[0];	      /* angle */
      if(param[1] <= 0.0)
	 param[1] = .01;
      b =  param[1];	/* stepsize */
      c =  param[2];	/* stop */
      t = l_d =  param[3];     /* points per orbit */
      sinx = sin(a);
      cosx = cos(a);
      orbit = 0;
      initorbit[0] = initorbit[1] = initorbit[2] = 0;
   }
   else if(fractype==FPHOPALONG || fractype==FPMARTIN)
   {
      initorbit[0] = 0;  /* initial conditions */
      initorbit[1] = 0;
      initorbit[2] = 0;
      connect = 0;
      a =  param[0];
      b =  param[1];
      c =  param[2];
      d =  param[3];
   }
   else if (fractype == INVERSEJULIAFP)
   {
      _CMPLX Sqrt;

      Cx = param[0];
      Cy = param[1];

      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);

      /* find fixed points: guaranteed to be in the set */
      Sqrt = ComplexSqrtFloat(1 - 4 * Cx, -4 * Cy);
      switch ((int) 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 rwalk;
	    }
	    EnQueueFloat((1 + Sqrt.x) / 2,  Sqrt.y / 2);
	    EnQueueFloat((1 - Sqrt.x) / 2, -Sqrt.y / 2);
	    break;
	 case depth_first:			/* depth first (choose direction) */
	    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 rwalk;
	    }
	    switch (minor_method) {
	       case left_first:
		  PushFloat((1 + Sqrt.x) / 2,  Sqrt.y / 2);
		  PushFloat((1 - Sqrt.x) / 2, -Sqrt.y / 2);
		  break;
	       case right_first:
		  PushFloat((1 - Sqrt.x) / 2, -Sqrt.y / 2);
		  PushFloat((1 + Sqrt.x) / 2,  Sqrt.y / 2);
		  break;
	    }
	    break;
	 case random_walk:
rwalk:
	    new.x = initorbit[0] = 1 + Sqrt.x / 2;
	    new.y = initorbit[1] = Sqrt.y / 2;
	    break;
	 case random_run:	/* random run, choose intervals */
	    major_method = random_run;
	    new.x = initorbit[0] = 1 + Sqrt.x / 2;
	    new.y = initorbit[1] = Sqrt.y / 2;
	    break;
      }
   }
   else
   {
      dt = param[0];
      a =  param[1];
      b =  param[2];
      c =  param[3];

   }

   /* precalculations for speed */
   adt = a*dt;
   bdt = b*dt;
   cdt = c*dt;

   return(1);
}

/******************************************************************/
/*   orbit functions - put in fractalspecific[fractype].orbitcalc */
/******************************************************************/

/* Julia sets by inverse iterations added by Juan J. Buhler 4/3/92 */
/* Integrated with Lorenz by Tim Wegner 7/20/92 */
/* Add Modified Inverse Iteration Method, 11/92 by Michael Snyder  */

int
Minverse_julia_orbit()
{
   static int   random_dir = 0, random_len = 0;
   int    newrow, newcol;
   int    color,  leftright;

   /*
    * First, compute new point
    */
   switch (major_method) {
      case breadth_first:
	 if (QueueEmpty())
	    return -1;
	 new = DeQueueFloat();
	 break;
      case depth_first:
	 if (QueueEmpty())
	    return -1;
	 new = PopFloat();
	 break;
      case random_walk:
#if 0
	 new = ComplexSqrtFloat(new.x - Cx, new.y - Cy);
	 if (RANDOM(2))
	 {
	    new.x = -new.x;
	    new.y = -new.y;
	 }
#endif
	 break;
      case random_run:
#if 0
	 new = ComplexSqrtFloat(new.x - Cx, new.y - Cy);
	 if (random_len == 0)
	 {
	    random_len = RANDOM(run_length);
	    random_dir = RANDOM(3);
	 }
	 switch (random_dir) {
	    case 0:	/* left */
	       break;
	    case 1:	/* right */
	       new.x = -new.x;
	       new.y = -new.y;
	       break;
	    case 2:	/* random direction */
	       if (RANDOM(2))
	       {
		  new.x = -new.x;
		  new.y = -new.y;
	       }
	       break;
	 }
#endif
	 break;
   }

   /*
    * Next, find its pixel position
    */
   newcol = cvt.a * new.x + cvt.b * new.y + cvt.e;
   newrow = cvt.c * new.x + cvt.d * new.y + cvt.f;

   /*
    * Now find the next point(s), and flip a coin to choose one.
    */

   new       = ComplexSqrtFloat(new.x - Cx, new.y - Cy);
   leftright = (RANDOM(2)) ? 1 : -1;

   if (newcol < 1 || newcol >= xdots || newrow < 1 || newrow >= ydots)
   {
      /*
       * MIIM must skip points that are off the screen boundary,
       * since it cannot read their color.
       */
      switch (major_method) {
	 case breadth_first:
	    EnQueueFloat(leftright * new.x, leftright * new.y);
	    return 1;
	 case depth_first:
	    PushFloat   (leftright * new.x, leftright * new.y);
	    return 1;
	 case random_run:
	 case random_walk:
	    break;
      }
   }

   /*
    * Read the pixel's color:
    * For MIIM, if color >= maxhits, discard the point
    *           else put the point's children onto the queue
    */
   color  = getcolor(newcol, newrow);
   switch (major_method) {
      case breadth_first:
	 if (color < maxhits)
	 {
	    putcolor(newcol, newrow, color+1);
/*	    new = ComplexSqrtFloat(new.x - Cx, new.y - Cy); */
	    EnQueueFloat( new.x,  new.y);
	    EnQueueFloat(-new.x, -new.y);
	 }
	 break;
      case depth_first:
	 if (color < maxhits)
	 {
	    putcolor(newcol, newrow, color+1);
/*	    new = ComplexSqrtFloat(new.x - Cx, new.y - Cy); */
	    if (minor_method == left_first)
	    {
	       if (QueueFullAlmost())
		  PushFloat(-new.x, -new.y);
	       else
	       {
		  PushFloat( new.x,  new.y);
		  PushFloat(-new.x, -new.y);
	       }
	    }
	    else
	    {
	       if (QueueFullAlmost())
		  PushFloat( new.x,  new.y);
	       else
	       {
		  PushFloat(-new.x, -new.y);
		  PushFloat( new.x,  new.y);
	       }
	    }
	 }
	 break;
      case random_run:
	 if (random_len-- == 0)
	 {
	    random_len = RANDOM(run_length);
	    random_dir = RANDOM(3);
	 }
	 switch (random_dir) {
	    case 0:	/* left */
	       break;
	    case 1:	/* right */
	       new.x = -new.x;
	       new.y = -new.y;
	       break;
	    case 2:	/* random direction */
	       new.x = leftright * new.x;
	       new.y = leftright * new.y;
	       break;
	 }
	 if (color < colors-1)
	    putcolor(newcol, newrow, color+1);
	 break;
      case random_walk:
	 if (color < colors-1)
	    putcolor(newcol, newrow, color+1);
	 new.x = leftright * new.x;
	 new.y = leftright * new.y;
	 break;
   }
   return 1;

}

int
Linverse_julia_orbit()
{
   static int   random_dir = 0, random_len = 0;
   int    newrow, newcol;
   int    color;

   /*
    * First, compute new point
    */
   switch (major_method) {
      case breadth_first:
	 if (QueueEmpty())
	    return -1;
	 lnew = DeQueueLong();
	 break;
      case depth_first:
	 if (QueueEmpty())
	    return -1;
	 lnew = PopLong();
	 break;
      case random_walk:
	 lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
	 if (RANDOM(2))
	 {
	    lnew.x = -lnew.x;
	    lnew.y = -lnew.y;
	 }
	 break;
      case random_run:
	 lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
	 if (random_len == 0)
	 {
	    random_len = RANDOM(run_length);
	    random_dir = RANDOM(3);
	 }
	 switch (random_dir) {
	    case 0:	/* left */
	       break;
	    case 1:	/* right */
	       lnew.x = -lnew.x;
	       lnew.y = -lnew.y;
	       break;
	    case 2:	/* random direction */
	       if (RANDOM(2))
	       {
		  lnew.x = -lnew.x;
		  lnew.y = -lnew.y;
	       }
	       break;
	 }
   }

   /*
    * Next, find its pixel position
    *
    * Note: had to use a bitshift of 21 for this operation because
    * otherwise the values of lcvt were truncated.  Used bitshift
    * of 24 otherwise, for increased precision.
    */
   newcol = (multiply(lcvt.a, lnew.x >> (bitshift - 21), 21) +
	     multiply(lcvt.b, lnew.y >> (bitshift - 21), 21) + lcvt.e) >> 21;
   newrow = (multiply(lcvt.c, lnew.x >> (bitshift - 21), 21) +
	     multiply(lcvt.d, lnew.y >> (bitshift - 21), 21) + lcvt.f) >> 21;

   if (newcol < 1 || newcol >= xdots || newrow < 1 || newrow >= ydots)
   {
      /*
       * MIIM must skip points that are off the screen boundary,
       * since it cannot read their color.
       */
      if (RANDOM(2))
	 color =  1;
      else
	 color = -1;
      switch (major_method) {
	 case breadth_first:
	    lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
	    EnQueueLong(color * lnew.x, color * lnew.y);
	    break;
	 case depth_first:
	    lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
	    PushLong(color * lnew.x, color * lnew.y);
	    break;
	 case random_run:
	    random_len--;
	 case random_walk:
	    break;
      }
      return 1;
   }

   /*
    * Read the pixel's color:
    * For MIIM, if color >= maxhits, discard the point
    *           else put the point's children onto the queue
    */
   color  = getcolor(newcol, newrow);
   switch (major_method) {
      case breadth_first:
	 if (color < maxhits)
	 {
	    putcolor(newcol, newrow, color+1);
	    lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
	    EnQueueLong( lnew.x,  lnew.y);
	    EnQueueLong(-lnew.x, -lnew.y);
	 }
	 break;
      case depth_first:
	 if (color < maxhits)
	 {
	    putcolor(newcol, newrow, color+1);
	    lnew = ComplexSqrtLong(lnew.x - CxLong, lnew.y - CyLong);
	    if (minor_method == left_first)
	    {
	       if (QueueFullAlmost())
		  PushLong(-lnew.x, -lnew.y);
	       else
	       {
		  PushLong( lnew.x,  lnew.y);
		  PushLong(-lnew.x, -lnew.y);
	       }
	    }
	    else
	    {
	       if (QueueFullAlmost())
		  PushLong( lnew.x,  lnew.y);
	       else
	       {
		  PushLong(-lnew.x, -lnew.y);
		  PushLong( lnew.x,  lnew.y);
	       }
	    }
	 }
	 break;
      case random_run:
	 random_len--;
	 /* fall through */
      case random_walk:
	 if (color < colors-1)
	    putcolor(newcol, newrow, color+1);
	 break;
   }
   return 1;
}

#if 0
int inverse_julia_orbit(double *x, double *y, double *z)
{
   double r, xo, yo;
   int waste;
   if(*z >= 1.0) /* this assumes intial value is 1.0!!! */
   {
      waste = 20; /* skip these points at first */
      *z = 0;
   }
   else
      waste = 1;
   while(waste--)
   {
      xo = *x - param[0]; 
      yo = *y - param[1];
      r  = sqrt(xo*xo + yo*yo);
      *x  = sqrt((r + xo)/2);
      if (yo < 0)
         *x = - *x;
      *y = sqrt((r - xo)/2);
      if(RANDOM(10) > 4)
      {
		  *x = -(*x);
		  *y = -(*y);
      }
   }
   return(0);
}
#endif

int lorenz3dlongorbit(long *l_x, long *l_y, long *l_z)
{
      l_xdt = multiply(*l_x,l_dt,bitshift);
      l_ydt = multiply(*l_y,l_dt,bitshift);
      l_dx  = -multiply(l_adt,*l_x,bitshift) + multiply(l_adt,*l_y,bitshift);
      l_dy  =  multiply(l_bdt,*l_x,bitshift) -l_ydt -multiply(*l_z,l_xdt,bitshift);
      l_dz  = -multiply(l_cdt,*l_z,bitshift) + multiply(*l_x,l_ydt,bitshift);

      *l_x += l_dx;
      *l_y += l_dy;
      *l_z += l_dz;
      return(0);
}

int lorenz3d1floatorbit(double *x, double *y, double *z)
{
      double norm;

      xdt = (*x)*dt;
      ydt = (*y)*dt;
      zdt = (*z)*dt;

      /* 1-lobe Lorenz */
      norm = sqrt((*x)*(*x)+(*y)*(*y));
      dx   = (-adt-dt)*(*x) + (adt-bdt)*(*y) + (dt-adt)*norm + ydt*(*z);
      dy   = (bdt-adt)*(*x) - (adt+dt)*(*y) + (bdt+adt)*norm - xdt*(*z) -
	     norm*zdt;
      dz   = (ydt/2) - cdt*(*z);

      *x += dx;
      *y += dy;
      *z += dz;
      return(0);
}

int lorenz3dfloatorbit(double *x, double *y, double *z)
{
      xdt = (*x)*dt;
      ydt = (*y)*dt;
      zdt = (*z)*dt;

      /* 2-lobe Lorenz (the original) */
      dx  = -adt*(*x) + adt*(*y);
      dy  =  bdt*(*x) - ydt - (*z)*xdt;
      dz  = -cdt*(*z) + (*x)*ydt;

      *x += dx;