Hi,
I've got lisp error , not maxima error
Adam
Adam Majewski pisze:
> Hi,
>
> Here is the program which IMHO needs more memory
> Try change jMax to 50 or 100 to see effect.
> ( result image : http://commons.wikimedia.org/wiki/File:Jung50e.png
> theory : http://fraktal.republika.pl/mset_jungreis.html
> or
> http://en.wikipedia.org/wiki/External_ray
>
> Regards
>
> Adam
>
> -====================== code ===================
> /* batch file for maxima
> uses :
> - symmetry around horizontal ( 0X ) axis
> - Psi_M function to map conjugate plane to parameter plane
> - jungreis algorithm to
>
>
> time for jMax:30; iMax:200; iMaxBig:400; is 1191 sec
> */
>
>
> start:elapsed_run_time ();
>
>
> jMax:30; /* precision = proportional to details and time of
> computations */
> iMax:200; /* number of points to draw */
> iMaxBig:400;
>
> /* computes b coefficient of Jungreis function*/
> betaF[n,m]:=block
> (
> [nnn:2^(n+1)-1],
> if m=0
> then 1.0
> else if ((n>0) and (m < nnn))
> then 0.0
> else (betaF[n+1,m]-
> sum(betaF[n,k]*betaF[n,m-k],k,nnn,m-nnn)-betaF[0,m-nnn])/2.0
> )$
>
> b[m]:=betaF[0,m+1]$
>
>
> /* -------------------------------*/
>
>
> /* Power of w to j */
> wn[w,j]:= if j=0 then 1 else w*wn[w,j-1]$
>
>
> /* ---------Jungreis function ; c = Psi_M(w)
> ----------------------------- */
> Psi_M(w):=w + sum(b[j]/wn[w,j],j,0,jMax)$
>
>
>
>
>
>
>
> /* exponential for of complex number with angle in turns */
> GiveCirclePoint(t):=R*%e^(%i*t*2*%pi)$ /* gives point of unit circle for
> angle t in turns */
> GiveWRayPoint(R):=R*%e^(%i*tRay*2*%pi)$ /* gives point of external ray
> for radius R and angle tRay in turns */
>
> NmbrOfCurves:9;
> /* coordinate of w-plane not c-plane */
> R_max:1.5;
> R_min:1;
> dR:(R_max-R_min)/NmbrOfCurves; /* for circles */
> dRR:(R_max-R_min)/iMax; /* for rays */
>
>
>
> /* --------------------------------------f_0 plane = w-plane
> -----------------------------------------*/
> /*-------------- unit circle ------------*/
> R:1;
> circle_angles:makelist(i/iMax,i,0,iMax/2)$
> CirclePoints:map(GiveCirclePoint,circle_angles)$
>
> /* ---------------external circles ---------*/
> circle_radii: makelist(R_min+i*dR,i,1,NmbrOfCurves)$
> circle_angles:makelist(i/iMaxBig,i,0,iMaxBig/2)$
> WCirclesPoints:[]$
> for R in circle_radii do
> WCirclesPoints:append(WCirclesPoints,map(GiveCirclePoint,circle_angles))$
>
> /* -------------- external w rays -------------*/
> ray_radii:makelist(R_min+dRR*i,i,1,iMax);
> ray_angles:[0,1/3,1/7 , 2/7 ,3/7 ]; /* list of angles < 1/2 of root
> points */
> WRaysPoints:[];
> for tRay in ray_angles do
> WRaysPoints:append(WRaysPoints,map(GiveWRayPoint,ray_radii));
>
>
>
>
>
> /* -------------------------parameter plane = c plane
> -----------------------------------*/
>
>
> MPoints:map(Psi_M,CirclePoints); /* Mandelbrot set points */
> CRaysPoints:map(Psi_M,WRaysPoints); /* external z rays */
> Equipotentials:map(Psi_M,WCirclesPoints);
>
> /* add points below horizontal axis */
> for w in CirclePoints do CirclePoints:cons(conjugate(w),CirclePoints);
> for w in WRaysPoints do WRaysPoints:cons(conjugate(w),WRaysPoints);
> for w in WCirclesPoints do WCirclesPoints:cons(conjugate(w),WCirclesPoints);
>
> for c in MPoints do MPoints:cons(conjugate(c),MPoints);
> for c in CRaysPoints do CRaysPoints:cons(conjugate(c),CRaysPoints);
> for c in Equipotentials do Equipotentials:cons(conjugate(c),Equipotentials);
>
>
> /* time */
>
> stop:elapsed_run_time ();
> time:fix(stop-start);
>
> /* ---------------- draw
> *--------------------------------------------------------------------------*/
>
> load(draw); /* Mario Rodr?guez Riotorto
> http://www.telefonica.net/web2/biomates/maxima/gpdraw/index.html */
> draw(file_name = "jung30es_2",
> pic_width=1000,
> pic_height= 500,
> terminal = 'png,
> columns = 2,
> gr2d(title = " unit circle with external rays & circles ",
> point_type = filled_circle,
> points_joined =true,
> point_size = 0.34,
> color = red,
>
> points(map(realpart, CirclePoints),map(imagpart, CirclePoints)),
> points_joined =false,
> color = green,
> points(map(realpart, WCirclesPoints),map(imagpart, WCirclesPoints)),
> color = black,
> points(map(realpart, WRaysPoints),map(imagpart, WRaysPoints))
> ),
> gr2d(title = "Parameter plane : Image under Psi_M(w) ",
> points_joined =true,
> point_type = filled_circle,
> point_size =0.34,
> color = blue,
>
> points(map(realpart, MPoints),map(imagpart, MPoints)),
> points_joined =false,
>
> color = green,
> points(map(realpart, Equipotentials),map(imagpart, Equipotentials)),
> color = black,
> points(map(realpart, CRaysPoints),map(imagpart, CRaysPoints))
> )
>
> );