lorenz.c

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

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

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

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

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

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

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

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

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

int henonfloatorbit(double *x, double *y, double *z)
{
      double newx,newy;
      newx  = 1 + *y - a*(*x)*(*x);
      newy  = b*(*x);
      *x = newx;
      *y = newy;
      return(0);
}

int henonlongorbit(long *l_x, long *l_y, long *l_z)
{
      long newx,newy;
      newx = multiply(*l_x,*l_x,bitshift);
      newx = multiply(newx,l_a,bitshift);
      newx  = fudge + *l_y - newx;
      newy  = multiply(l_b,*l_x,bitshift);
      *l_x = newx;
      *l_y = newy;
      return(0);
}

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

      dx = -ydt - (*z)*dt;
      dy = xdt + (*y)*adt;
      dz = bdt + (*z)*xdt - (*z)*cdt;

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

int pickoverfloatorbit(double *x, double *y, double *z)
{
      double newx,newy,newz;
      newx = sin(a*(*y)) - (*z)*cos(b*(*x));
      newy = (*z)*sin(c*(*x)) - cos(d*(*y));
      newz = sin(*x);
      *x = newx;
      *y = newy;
      *z = newz;
      return(0);
}
/* page 149 "Science of Fractal Images" */
int gingerbreadfloatorbit(double *x, double *y, double *z)
{
      double newx;
      newx = 1 - (*y) + fabs(*x);
      *y = *x;
      *x = newx;
      return(0);
}

int rosslerlongorbit(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  = -l_ydt - multiply(*l_z,l_dt,bitshift);
      l_dy  =  l_xdt + multiply(*l_y,l_adt,bitshift);
      l_dz  =  l_bdt + multiply(*l_z,l_xdt,bitshift)
		     - multiply(*l_z,l_cdt,bitshift);

      *l_x += l_dx;
      *l_y += l_dy;
      *l_z += l_dz;

      return(0);
}

/* OSTEP  = Orbit Step (and inner orbit value) */
/* NTURNS = Outside Orbit */
/* TURN2  = Points per orbit */
/* a	  = Angle */


int kamtorusfloatorbit(double *r, double *s, double *z)
{
   double srr;
   if(t++ >= l_d)
   {
      orbit += b;
      (*r) = (*s) = orbit/3;
      t = 0;
      *z = orbit;
      if(orbit > c)
	 return(1);
   }
   srr = (*s)-(*r)*(*r);
   (*s)=(*r)*sinx+srr*cosx;
   (*r)=(*r)*cosx-srr*sinx;
   return(0);
}

int kamtoruslongorbit(long *r, long *s, long *z)
{
   long srr;
   if(t++ >= l_d)
   {
      l_orbit += l_b;
      (*r) = (*s) = l_orbit/3;
      t = 0;
      *z = l_orbit;
      if(l_orbit > l_c)
	 return(1);
   }
   srr = (*s)-multiply((*r),(*r),bitshift);
   (*s)=multiply((*r),l_sinx,bitshift)+multiply(srr,l_cosx,bitshift);
   (*r)=multiply((*r),l_cosx,bitshift)-multiply(srr,l_sinx,bitshift);
   return(0);
}
/*#define sign(x) ((x)>0?1:((x)<0?(-1):0))*/
#define sign(x) ((x)>=0?1:-1)
int hopalong2dfloatorbit(double *x, double *y, double *z)
{
   double tmp;
   tmp = *y - sign(*x)*sqrt(fabs(b*(*x)-c));
   *y = a - *x;
   *x = tmp;
   return(0);
}

int martin2dfloatorbit(double *x, double *y, double *z)
{
   double tmp;
   tmp = *y - sin(*x);
   *y = a - *x;
   *x = tmp;
   return(0);
}

int mandelcloudfloat(double *x, double *y, double *z)
{
    double newx,newy,x2,y2;
    x2 = (*x)*(*x);
    y2 = (*y)*(*y);
    if (x2+y2>2) return 1;
    newx = x2-y2+a;
    newy = 2*(*x)*(*y)+b;
    *x = newx;
    *y = newy;
    return(0);
}

int dynamfloat(double *x, double *y, double *z)
{
      _CMPLX cp,tmp;
      double newx,newy;
      cp.x = b* *x;
      cp.y = 0;
      CMPLXtrig0(cp, tmp);
      newy = *y + dt*sin(*x + a*tmp.x);
      if (euler) {
	  *y = newy;
      }

      cp.x = b* *y;
      cp.y = 0;
      CMPLXtrig0(cp, tmp);
      newx = *x - dt*sin(*y + a*tmp.x);
      *x = newx;
      *y = newy;
      return(0);
}

/* dmf */
#undef  LAMBDA
#define LAMBDA  param[0]
#define ALPHA   param[1]
#define BETA    param[2]
#define GAMMA   param[3]
#define OMEGA   param[4]
#define DEGREE  param[5]

int iconfloatorbit(double *x, double *y, double *z)
{

    double oldx, oldy, zzbar, zreal, zimag, za, zb, zn, p;
    unsigned char i;

    oldx = *x;
    oldy = *y;

    zzbar = oldx * oldx + oldy * oldy;
    zreal = oldx;
    zimag = oldy;

    for(i=1; i <= DEGREE-2; i++) {
        za = zreal * oldx - zimag * oldy;
        zb = zimag * oldx + zreal * oldy;
        zreal = za;
        zimag = zb;
    }
    zn = oldx * zreal - oldy * zimag;
    p = LAMBDA + ALPHA * zzbar + BETA * zn;
    *x = p * oldx + GAMMA * zreal - OMEGA * oldy;
    *y = p * oldy - GAMMA * zimag + OMEGA * oldx;

/*    if (fractype==ICON3D) */
    *z = zzbar; 
    return(0);
}
#ifdef LAMBDA  /* Tim says this will make me a "good citizen" */
#undef LAMBDA  /* Yeah, but you were bad, Dan - LAMBDA was already */
#undef ALPHA   /* defined! <grin!> TW */
#undef BETA
#undef GAMMA
#endif
/**********************************************************************/
/*   Main fractal engines - put in fractalspecific[fractype].calctype */
/**********************************************************************/

int inverse_julia_per_image()
{
   int color = 0;

   if (resuming)            /* can't resume */
      return -1;

   while (color >= 0)       /* generate points */
   {
      if (check_key())
      {
	 Free_Queue();
	 return -1;
      }
      color = curfractalspecific->orbitcalc();
      old = new;
   }
   Free_Queue();
   return 0;
}

int orbit2dfloat()
{
   FILE *fp;
   double *soundvar;
   double x,y,z;
   int color,col,row;
   int count;
   int oldrow, oldcol;
   double *p0,*p1,*p2;
   struct affine cvt;
   int ret,start;
   start = 1;
   
   fp = open_orbitsave();
   /* setup affine screen coord conversion */
   setup_convert_to_screen(&cvt);

   /* set up projection scheme */
   if(projection==0)
   {
      p0 = &z;
      p1 = &x;
      p2 = &y;
   }
   else if(projection==1)
   {
      p0 = &x;
      p1 = &z;
      p2 = &y;
   }
   else if(projection==2)
   {
      p0 = &x;
      p1 = &y;
      p2 = &z;
   }
   if(soundflag==1)
      soundvar = &x;
   else if(soundflag==2)
      soundvar = &y;
   else if(soundflag==3)
      soundvar = &z;

   count = 0;
   if(inside > 0)
      color = inside;
   if(color >= colors)
      color = 1;
   oldcol = oldrow = -1;
   x = initorbit[0];
   y = initorbit[1];
   z = initorbit[2];

   if (resuming)
   {
      start_resume();
      get_resume(sizeof(count),&count,sizeof(color),&color,
          sizeof(oldrow),&oldrow,sizeof(oldcol),&oldcol,
          sizeof(x),&x,sizeof(y),&y,sizeof(z),&z,
          sizeof(t),&t,sizeof(orbit),&orbit,
          0);
      end_resume();
   }

   ret = 0;
   while(1)
   {
      if(check_key())
      {
         nosnd();
         alloc_resume(100,1);
         put_resume(sizeof(count),&count,sizeof(color),&color,
             sizeof(oldrow),&oldrow,sizeof(oldcol),&oldcol,
             sizeof(x),&x,sizeof(y),&y,sizeof(z),&z,
             sizeof(t),&t,sizeof(orbit),&orbit,
             0);
         ret = -1;
         break;
      }

      if (++count > 1000)
      {        /* time to switch colors? */
         count = 0;
         if (++color >= colors)   /* another color to switch to? */
	    color = 1;  /* (don't use the background color) */
      }

      col = cvt.a*x + cvt.b*y + cvt.e;
      row = cvt.c*x + cvt.d*y + cvt.f;
      if ( col >= 0 && col < xdots && row >= 0 && row < ydots )
      {
         if (soundflag > 0)
            snd((int)(*soundvar*100+basehertz));
/* dmf */
         if(fractype!=ICON)
         {
         if(oldcol != -1 && connect)
            draw_line(col,row,oldcol,oldrow,color&(colors-1));
            else
            (*plot)(col,row,color&(colors-1));
         } else {
            /* should this be using plothist()? */
            color = getcolor(col,row)+1;
            if( color < colors ) /* color sticks on last value */
               (*plot)(col,row,color);

         }

         oldcol = col;
         oldrow = row;
         start = 0;
      }
      else if((long)abs(row) + (long)abs(col) > BAD_PIXEL) /* sanity check */
         return(ret);
      else   
         oldrow = oldcol = -1;

      if(curfractalspecific->orbitcalc(p0, p1, p2))
         break;
      if(fp)
         fprintf(fp,orbitsave_format,*p0,*p1,0.0);
   }
   if(fp)
      fclose(fp);
   return(ret);
}

int orbit2dlong()
{
   FILE *fp;
   long *soundvar;
   long x,y,z;
   int color,col,row;
   int count;
   int oldrow, oldcol;
   long *p0,*p1,*p2;
   struct l_affine cvt;
   int ret,start;
   start = 1;
   fp = open_orbitsave();

   /* setup affine screen coord conversion */
   l_setup_convert_to_screen(&cvt);

   /* set up projection scheme */
   if(projection==0)
   {
      p0 = &z;
      p1 = &x;
      p2 = &y;
   }
   else if(projection==1)
   {
      p0 = &x;
      p1 = &z;
      p2 = &y;
   }
   else if(projection==2)
   {
      p0 = &x;
      p1 = &y;
      p2 = &z;
   }
   if(soundflag==1)
      soundvar = &x;
   else if(soundflag==2)
      soundvar = &y;
   else if(soundflag==3)
      soundvar = &z;
   count = 0;
   if(inside > 0)
      color = inside;
   if(color >= colors)
      color = 1;
   oldcol = oldrow = -1;
   x = initorbitlong[0];
   y = initorbitlong[1];
   z = initorbitlong[2];

   if (resuming)
   {
      start_resume();
      get_resume(sizeof(count),&count,sizeof(color),&color,
          sizeof(oldrow),&oldrow,sizeof(oldcol),&oldcol,
          sizeof(x),&x,sizeof(y),&y,sizeof(z),&z,
          sizeof(t),&t,sizeof(l_orbit),&l_orbit,
          0);
      end_resume();
   }

   ret = 0;
   while(1)
   {
      if(check_key())
      {
         nosnd();
         alloc_resume(100,1);
         put_resume(sizeof(count),&count,sizeof(color),&color,
	     sizeof(oldrow),&oldrow,sizeof(oldcol),&oldcol,
             sizeof(x),&x,sizeof(y),&y,sizeof(z),&z,
             sizeof(t),&t,sizeof(l_orbit),&l_orbit,
             0);
         ret = -1;
         break;
      }
      if (++count > 1000)
      {        /* time to switch colors? */
         count = 0;