Re: Mandelbrot
- To: mathgroup at smc.vnet.net
- Subject: [mg13616] Re: Mandelbrot
- From: daiyanh at mindspring.com (Daitaro Hagihara)
- Date: Fri, 7 Aug 1998 03:08:18 -0400
- Organization: MindSpring Enterprises
- References: <6q3r31$f29@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
In article <6q3r31$f29 at smc.vnet.net>, Jon Prudhomme
<prudhomj at elwha.evergreen.edu> wrote:
> Hello
>
> I was just curious if anyone had found a decent way to plot the
> Mandelbrot or Julia sets with Mathematica yet. I have been able to do
> it with DensityPlot and ListDensityPlot, but I can't help but wonder if
> there is an easier way than either of these:
>
>
> iterations=200;
> pointColor[c_]:=Module[{i,p},
> For[i=1;p=0,i<=iterations&&Sqrt[Re[p]^2+Im[p]^2]<=2,i++,p=p^2+c];i]
>
> DensityPlot[pointColor[Complex[x,y]],{x,-2.5,1.5},{y,-1.5,1.5},Mesh->False,
>
> ColorFunction->Hue,AspectRatio->Automatic,PlotPoints->1000]
>
> (* or this for a ListDensityPlot... *)
>
> manSet=Table[pointColor[Complex[x,y]],{y,-1.5,1.5,.01},{x,-2.5,1.5,.01}];
> ListDensityPlot[manSet,Mesh->False,ColorFunction->Hue,AspectRatio->Automatic]
>
> Anyone got any other ideas?
>
> Jon Prudhomme
> The Evergreen State College
> prudhomj at elwha.evergreen.edu
>
>
> PS - The algorithm for the Mandelbrot set is z[[n]]=z[[n-1]]^2+c where
> z[[0]]=0 and c is the point on the complex plain being tested as a
> member of the set. If after an arbitrary number of iterations the
> point is not 2 units away from the origin on the complex plain, the
> point c is a member of the set. Colors of non-member points are based
> on the iteration number that them as excluded from the set.
The following example from MathSource might be what you're looking for.
First compile rast.tm. The makefile for Mac's MPW-PPC is enclosed for
your convenience. Then eval the following in Mathematica. You'll get
a nice Mandelbrot movie. As you can see, the time it took to plot 34
complex frames was under 5 min on my 200MHz PPC. If you want to keep
movies under different params, then it's a good idea to compress them
to save disk space. I find 8-bit Animation codec with 100% quality
(the default for Convert to QT command in Mathematica) to be quite
adequate for archival purposes. (Don't even think of Cinepak, Indeo,
or even MPEG. Their performance is very poor when it comes to CG
animations.)
Install["Macintosh HD:rast"];
f[n_] := 0.8^n;
g[n_] := {{-1.40835915-f[n],-1.40835915+f[n]},
{0.13627737-f[n],0.13627737+f[n]}};
plot[n_]:=Show[Graphics[RasterArray[
FractalML[0.0, 0.0, g[n], 216]],
AspectRatio -> 1]];
Timing[Do[plot[n],{n,-3, 30}]] -->
- series of plots -
{262.733 Second, Null}
*************************** MAKEFILE ***************************
#MPFLAGS=-Application
CFLAGS=-sym on
LFLAGS=-sym on -d
CC=MrC
AS=PPCAsm
LN=PPCLink
LINKLIBS= "Macintosh HD:Mathematica 2.2.2:"MathLinkLibrary \
"{{SharedLibraries}}"InterfaceLib \
"{{SharedLibraries}}"MathLib \
"{{SharedLibraries}}"StdCLib \
"{{PPCLibraries}}"StdCRuntime.o \
"{{PPCLibraries}}"PPCCRuntime.o
.tm.c:
mprep $(MPFLAGS) -o $@ $*.tm
.c.o:
$(CC) $(CFLAGS) -o $@ $*.c
.asm.o:
$(AS) -o $@ $*.asm
rast: rast.o
$(LN) $(LFLAGS) rast.o $(LINKLIBS) -o rast
Rez mathlink.r -o rast -append
SetVersion -t 'vers' -i 1 -version ""1.0"" rast
SetFile rast -c '????' -t APPL
MakeSYM -sym 3.2 -w rast.xcoff
****************************************************************
**************************** rast.tm ****************************
:Begin:
:Function: fractalml
:Pattern: FractalML[a_Real,b_Real,{{c_Real,d_Real},{e_Real,f_Real}},g_Integer]
:Arguments: { a,b,c,d,e,f,g }
:ArgumentTypes: { Real,Real,Real,Real,Real,Real,Integer}
:ReturnType: Manual
:End:
#include "mathlink.h"
void fractalml(a,b,c,d,e,f,g)
double a,b,c,d,e,f;
int g;
{
double x0,y0,xmin,ymin,xmax,ymax;
double x,y,dist, xadd, yadd, temp;
int counts;
int xco,yco;
int divs;
x0 = (double) a;
y0 = (double) b;
xmin = (double) c;
xmax = (double) d;
ymin = (double) e;
ymax = (double) f;
divs = g;
MLPutFunction(stdlink,"List",divs);
for(yco = 0; yco <divs; yco++){
yadd = ymin + (ymax-ymin)*yco/(divs-1.0);
MLPutFunction(stdlink,"List",divs);
for(xco = 0; xco <divs ; xco++){
x = x0;
y = y0;
xadd = xmin + (xmax-xmin)*xco/(divs-1.0);
dist = x0*x0+y0*y0;
for(counts = 0; counts<100 && dist < 10000; counts++){
temp = x*x-y*y+xadd;
y = 2*x*y+yadd;
x = temp;
dist = x*x+y*y;
}
if(counts == 100){
MLPutFunction(stdlink,"RGBColor",3);
MLPutReal(stdlink,0.0);
MLPutReal(stdlink,0.0);
MLPutReal(stdlink,0.0);}
else{
MLPutFunction(stdlink,"Hue",1);
MLPutReal(stdlink,counts/120.0);
}
}
}
}
int main(argc,argv)
int argc; char* argv[];
{
return MLMain(argc, argv);
}
****************************************************************