Code for dynamical pictures

zpt[z_]:={Re[z],Im[z]}

triple[{x_,y_}]:=Drop[Flatten[{zpt[y],zpt[x]}],-1]

F[{z0_,z1_},b_,c_]:={z1,z1^2+c-b*z0}

Finv[{z0_,z1_},b_,c_]:={(z0^2+c-z1)/b,z0}

i2[z_]:={slice,z}

slice=0;

p2[{z_,w_}]:=w

gg[n_,z_?NumericQ,b_,c_]:=p2[Nest[F[#,b,c]&,i2[z],n]]

ff[n_,z_,b_,c_]:=Nest[F[#,b,c]&,i2[z],n]

rr[b_,c_]:=(Abs[b]+2+Sqrt[(Abs[b]+2)^2+4*Abs[c]])/2
r0[b_,c_]:=rr[b,c]+0.5

gl[lst_]:=Line[Map[zpt,lst,{-1}]]

gl2[prpth_]:=Line[Map[triple,prpth]]

frn[n_,z1_,z2_,b_,c_]:=z/.FindRoot[gg[n,z,b,c]==z2,{z,z1}]

sspath[{z_,w_}]:=Table[(1-s)z+s*w,{s,0,1,0.001}]

py:=Print["This is a y-plane picture with parameters ","b0=",b0," c0=",c0]

px:=Print["This is an x-plane picture with parameters ","b0=",b0," c0=",c0]

x[1,b_,c_]:=Table[r0[b,c]*Exp[2*Pi*I*th],{th,0,1,0.005}]

y[j_,b_,c_]:=y[j,b,c]=Map[gg[j,#,b,c]&,x[j,b,c]]

cpath[deg_,len_,b_,c_]:=Table[r0[b,c]*Exp[(2 Pi I  2^deg j)/len],{j,0,len}]

homy[j_,b_,c_]:=homy[j,b,c]=Map[sspath,Transpose[{y[j,b,c],cpath[j,Length[y[j,b,c]]-1,b,c]}]]

frn[n_,z1_,z2_,b_,c_]:=z/.FindRoot[gg[n,z,b,c]==z2,{z,z1}]

newlift[n_,{z_,path1_},b_,c_]:=Drop[FoldList[frn[n,#1,#2,b,c]&,z,path1],1]

liftpths[n_,initptslst_,pthlst_,b_,c_]:=Map[newlift[n,#,b,c]&,Transpose[{initptslst,pthlst}]]

homx[j_,b_,c_]:=homx[j,b,c]=liftpths[j,x[j,b,c],homy[j,b,c],b,c]

double[pth_]:=Table[If[OddQ[k],pth[[(k+1)/2]],(pth[[k/2]]+pth[[k/2+1]])/2],{k,1,2Length[pth]-1}]

x[j_?OddQ,b_,c_]:=x[j,b,c]=Map[Last,homx[j-1,b,c]]
x[j_?EvenQ,b_,c_]:=x[j,b,c]=double[Map[Last,homx[j-1,b,c]]]

homz[j_,k_,b_,c_]:=homz[j,k,b,c]=Map[ff[j,#,b,c]&,homx[k,b,c],{2}]

z[j_,k_,b_,c_]:=z[j,k,b,c]=Map[ff[j,#,b,c]&,x[k,b,c]]

xx[1,th_,b_,c_]:=r0[b,c]*Exp[2*Pi*I*th];

yy[j_,th_,b_,c_]:=y[j,th,b,c]=gg[j,xx[j,th,b,c],b,c]

ccpath[deg_,th_,b_,c_]:=r0[b,c]*Exp[(2 Pi I  2^deg)*th]

hhomy[j_,th_,b_,c_]:=hhomy[j,th,b,c]=sspath[{yy[j,th,b,c],ccpath[j,th,b,c]}]

lliftpth[n_,initpt_,pth_,b_,c_]:=newlift[n,{initpt,pth},b,c]

hhomx[j_,th_,b_,c_]:=hhomx[j,th,b,c]=lliftpth[j,xx[j,th,b,c],hhomy[j,th,b,c],b,c]

xx[j_,th_,b_,c_]:=xx[j,th,b,c]=Last[hhomx[j-1,th,b,c]]

raypath[th_,lev_,b_,c_]:=Join[Table[hhomx[j,th,b,c],{j,lev}]]

raypath[th_,lev_,b_,c_]:=Flatten[Join[Table[hhomx[j,th,b,c],{j,lev}]]]
raypath2[th_,lev1_,lev2_,b_,c_]:=Flatten[Join[Table[hhomx[j,th,b,c],{j,lev1,lev2}]]]


gpt[z_]:=Graphics[{PointSize[Tiny],Point[z]}];
parameterlist=Table[(0.14)*E^({j+4}*I*Pi/32) + (-1.9-0.15I),{j,0,64}] /.{x_}->x;

bseq={0,0};

lstfrm[n_,l_]:=IntegerDigits[n-1,2,l]

even[{n_,l_}]:=If[l>0,EvenQ[n],False]

odd[{n_,l_}]:=If[l>0,OddQ[n],False]

extend[{n_,l_}]:={2n,l+1}

div[{n_,l_}]:={Quotient[n,2],l-1}

inc[{n_,l_}]:={n+1,l}

dec[{n_,l_}]:={n-1,l}

len[{n_,l_}]:=l

res=5;

circle[b_,c_]:=Table[r0[b,c]*Exp[2 Pi I th/2^res],{th,0,2^res}];

pt[seq_,b_,c_]:=pt[seq,b,c]=If[len[seq]==0,r0[b,c],If[even[seq],Last[down[div[seq],b,c]],Last[around[dec[seq],b,c]]]]

around[seq_,b_,c_]:=around[seq,b,c]=newlift[len[seq],{pt[seq,b,c],circle[b,c]},b,c]

down[seq_,b_,c_]:=down[seq,b,c]=newlift[len[seq]+1,{pt[seq,b,c],sspath[{gg[len[seq+1],pt[seq,b,c],b,c],r0[b,c]}]},b,c]

pf:=Print["frame=",frame," c0=",c0, " scale=",scale]

level[k_]:=Table[{j,k},{j,0,2^k-1}]

shlev[l_,b_,c_]:=Map[gl[around[#,b,c]]&,level[l]]

shdownlev[l_,b_,c_]:=Map[gl[down[#,b,c]]&,level[l]]

aend[seq_,b_,c_]:=Last[around[seq,b,c]]

Permutation functions

near[z_,w_]:=Abs[z-w]<.0001

permutation given as a function

apermfcn[seq1_,b_,c_]:=Select[level[len[seq1]],near[aend[seq1,b,c],pt[#,b,c]]&][[1]]

permutation viewed as a list of numbers. given in the form of cycles.

aperm[l_,b_,c_]:=FindPermutation[Map[apermfcn[#,b,c]&,level[l]],level[l]]

aperminv[l_,b_,c_]:=FindPermutation[level[l],Map[apermfcn[#,b,c]&,level[l]]]

List of components

comps[l_,b_,c_]:=aperm[l,b,c][[1]]

Number of components

compno[l_,b_,c_]:=Length[comps[l,b,c]]

indfrm[n_,l_]:={n-1,l}

indcomps[l_,b_,c_]:=Map[indfrm[#,l]&,comps[l,b,c],{2}]

mp[{n_,l_}]:={Mod[n,2^(l-1)],l-1}

mplist[list_]:=Map[mp,list]

r3pic[b_,c_,level_,n_]:=Graphics3D[Map[Line,Map[triple[F[i2[#],b,c]]&,Map[around[#,b,c]&,indcomps[level,b,c],{2}],{3}]]]

nestedf[z_,b_,c_,n_]:=Nest[F[#1,b,c]&,z,n]

Below is the command which draws a picture in R^3, parameters b,c are obvious from the Henon map. the level corresponds to better and better approximations. n determines the number of times we map the slice approximation forward by the Henon map.

mappingforwardpicture[b_,c_,level_,n_]:=mappingforwardpicture[b,c,level,n]=Graphics3D[Map[Line,Map[triple[nestedf[#,b,c,n]]&,Map[i2[#]&,Map[around[#,b,c]&,indcomps[level+n-1,b,c],{2}],{3}],{3}]]]

mappingforwardpictureblue[b_,c_,level_,n_]:=Graphics3D[{Blue,Map[Line,Map[triple[nestedf[#,b,c,n]]&,Map[i2[#]&,Map[around[#,b,c]&,indcomps[level+n-1,b,c],{2}],{3}],{3}]]}]

mappedrayshowray[th_,b_,c_,n_]:=Graphics3D[gl2[Map[nestedf[#,b,c,n]&,Map[i2[#]&,raypath2[th,2+n-1,10,b,c]]]]]

Code for parameter space pictures

p2[{z_,w_}]:=w

i2[z_]:={0,z}

F2[{z0_,z1_},b_,c_]:={z1,z1^2+c-b*z0}

ggs2[n_,z_,b_,c_]:=p2[Nest[F2[#,b,c]&,i2[z],n]]

Critical points:

cp22[its_,b_,c_,r_]:=z/.FindRoot[D[ggs2[its,z,b,c],z]==0,{z,r}]

Test for NestWhile loop

test22[{p_,b_,c_,ct_}]:=(Abs[p[[2]]]<3) && (ct<30)

Iteration for NestWhile loop

f22[{p_,b_,c_,ct_}]:={F2[p,b,c],b,c,ct+1}

Loops with different critical points:

esc20[b_,c_]:=Last[NestWhile[f22,{i2[cp22[3,b,c,0]],b,c,0},test22]]

esc21[b_,c_]:=Last[NestWhile[f22,{i2[cp22[3,b,c,Sqrt[-c]]],b,c,0},test22]]

esc22[b_,c_]:=Last[NestWhile[f22,{i2[cp22[3,b,c,-Sqrt[-c]]],b,c,0},test22]]

Test with different stopping point:

test22[{p_,b_,c_,ct_}]:=(Abs[p[[2]]]<3) && (ct<20)

Matrix corresponding to a square in parameter space. Arguments correspond to: center, magnification and resolution.

sq[re_,im_,mag_,res_]:=Table[x+y*I,{y,re-1/mag,re+1/mag,1/(res*mag)},{x,im-1/mag,im+1/mag,1/(res*mag)}]

Used for data range function.

dr[re_,im_,mag_,res_]:={{im-1/mag,im+1/mag},{re-1/mag,re+1/mag}}

Display functions. Incorporates Jacobian parameter.

disp20[cen_,mag_,res_,b_]:=ArrayPlot[array20[cen,mag,res,b],DataRange->dr[-Im[cen],Re[cen],mag,res],FrameTicks->Automatic]

disp21[cen_,mag_,res_,b_]:=ArrayPlot[array21[cen,mag,res,b],DataRange->dr[-Im[cen],Re[cen],mag,res],FrameTicks->Automatic]

disp22[cen_,mag_,res_,b_]:=ArrayPlot[array22[cen,mag,res,b],DataRange->dr[-Im[cen],Re[cen],mag,res],FrameTicks->Automatic]

Array computation functions.

array20[cen_,mag_,res_,b_]:=array20[cen,mag,res,b]=Reverse[Map[esc20[b,#]&,sq[-Im[cen],Re[cen],mag,res],{2}]]

array21[cen_,mag_,res_,b_]:=array21[cen,mag,res,b]=Reverse[Map[esc21[b,#]&,sq[-Im[cen],Re[cen],mag,res],{2}]]

array22[cen_,mag_,res_,b_]:=array22[cen,mag,res,b]=Reverse[Map[esc22[b,#]&,sq[-Im[cen],Re[cen],mag,res],{2}]]

Saddle=Compile[{{b,_Complex},{c,_Complex},{init,_Complex},{n,_Integer}},
Module[{z0,z1,z2,j},
z0=0.0+0.0 I;
z1=init;
For[j=0,j<=n,j++,
If[Abs[z1]>=3.0,Return[j],
z2=z1^2+c-b*z0;
z0=z1;
z1=z2;]];
Return[0]],
RuntimeOptions->"Speed"
];

sq[re_,im_,mag_,res_]:=Table[x+y*I,{y,re-1/mag,re+1/mag,1/(res*mag)},{x,im-1/mag,im+1/mag,1/(res*mag)}]

dr[re_,im_,mag_,res_]:={{im-1/mag,im+1/mag},{re-1/mag,re+1/mag}}

cp[its_,b_,c_,r_]:=z/.FindRoot[D[ggs2[its,z,b,c],z]==0,{z,r}];

array1[cen_,mag_,res_,b_,maxit_]:=Reverse[ParallelMap[Saddle[b,#,cp[3,b,#,Sqrt[-#]],maxit]&,sq[Im[cen],Re[cen],mag,res],{2}]]

array2[cen_,mag_,res_,b_,maxit_]:=Reverse[ParallelMap[Saddle[b,#,cp[3,b,#,-Sqrt[-#]],maxit]&,sq[Im[cen],Re[cen],mag,res],{2}]]

ctab={RGBColor[66/256 ,30/256, 15/256],RGBColor[25/256, 7/256, 26/256],RGBColor[9/256 ,1/256, 47/256],RGBColor[4/256 ,4/256 ,73/256],RGBColor[0/256 ,7/256, 100/256],RGBColor[12/256, 44/256, 138/256],RGBColor[24/256, 82/256 ,177/256],RGBColor[57/256 ,125/256, 209/256],RGBColor[134/256, 181/256, 229/256],RGBColor[211 /256,236/256 ,248/256],RGBColor[241/256, 233/256, 191/256],RGBColor[248/256, 201/256, 95/256],RGBColor[255/256, 170/256, 0/256],RGBColor[204/256, 128/256, 0/256],RGBColor[153/256, 87/256, 0],RGBColor[106 /256,52/256 ,3/256]};

Wikicolor[j_]:=ctab[[Mod[j,16]+1]]

disp1[cen_,mag_,res_,b_,maxit_]:=ArrayPlot[array1[cen,mag,res,b,maxit],ColorFunctionScaling->False,ColorFunction->Wikicolor,DataRange->dr[Im[cen],Re[cen],mag,res],FrameTicks->Automatic]

disp2[cen_,mag_,res_,b_,maxit_]:=ArrayPlot[array2[cen,mag,res,b,maxit],ColorFunctionScaling->False,ColorFunction->Wikicolor,DataRange->dr[Im[cen],Re[cen],mag,res],FrameTicks->Automatic]

Critical value commands

gg1[n_,z_,b_,c_]:=p2[Nest[F[#,b,c]&,i2[z],n]]

criticalpoint[n_,b_,c_]:=z/.Solve[D[gg1[n,z,b,c],z]==0]

cv[n_,b_,c_]:=Map[gg1[n,#,b,c]&,z/.Solve[D[gg1[n,z,b,c],z]==0]]

critpoint[n_,b_,c_]:=z/.Solve[D[gg1[n,z,b,c],z]==0]

criticalvaluesofalist[n_,b_,list_,critptnum_]:=Show[Table[Graphics3D[Point[triple[F[i2[cv[n,b,list[[i]]][[critptnum]]],b,list[[i]]]]]],{i,1,Length[list]}]]

criticalpointsofalist[n_,b_,list_,critptnum_]:=Show[Table[Graphics3D[Point[triple[F[i2[critpoint[n,b,list[[i]]][[critptnum]]],b,list[[i]]]]]],{i,1,Length[list]}]]

classicalcritptandvalpic[level_ ,mapforward_,critptlevel_,b_,list_,critptnum_]:=Show[mappingforwardpicture[b,list[[1]],level,mapforward],Graphics3D[Point[triple[F[i2[critpoint[critptlevel,b,list[[1]]][[critptnum]]],b,list[[1]]]]]],Table[Graphics3D[Point[triple[F[i2[cv[critptlevel,b,list[[i]]][[critptnum]]],b,list[[i]]]]]],{i,1,Length[list]}]]

lipaslice=;

Critical value angles for spinning b herd, Lipa b herd. Lie in the 3/7,4/7 piece

BherdBcritangle={0.5708176,0.4466,0.448,0.45,0.452,0.4536,0.4556,0.459,0.46,0.4628,
0.46580, 0.4694,0.4724,0.4764,0.482,0.489,0.498,0.501,0.504,0.51,
0.516,0.52,0.525,0.5296,0.5305,0.534,0.5362,0.5398,0.5436,0.5468,
0.550,0.553,0.556,0.558,0.561,0.5638,0.5662,0.5676,0.569,0.5694,
0.5701,0.5706};

Spinning A herd angles

BherdAcritangle={0.5668,0.565,0.563,0.561,0.559,0.5576,0.5562,0.554,0.553,0.551,
0.549, 0.547,0.5456,0.5442,0.5424,0.541,0.54,0.5396,0.5412,0.5425,
0.544,0.5456,0.5474,0.5494,0.55,0.5526,0.554,0.5562,0.5588,0.561,
0.563,0.565,0.5667,0.5685,0.5704,0.5724,0.5736, 0.5736,0.5734,0.5724,
0.5705,0.5684};

Critical value angles for spinning A herd, Lipa A herd . Lie in the 13/31, 18/31 piece

AherdAcritangle={0.574,0.5755, 0.5758, 0.41925, 0.4195, 0.4208,0.4229,0.4258,0.4288,0.4314,
0.434,0.441,0.4442,0.4486,0.453,0.458,0.4615,0.4667,0.468,0.4725,
0.4778,0.483,0.487,0.489,0.492,0.495,0.498,0.503,0.517,0.525,
0.5316,0.5358,0.5408,0.5444,0.546,0.55,0.5536,0.556,0.561,0.5665,
0.57,0.5726};

Critical value angles for spinning B herd, Lipa A herd

AherdBcritangle={0.4245,0.4208,0.4184,0.4168,0.4167,0.4166,0.4184,0.421,0.4232,0.4244,
0.4269,0.4292,0.4318,0.4346,0.438,0.4406,0.4432,0.4466,0.4474,0.4498,
0.4518,0.454,0.456,0.458, 0.4604,0.4609,0.4604,0.4596,0.458,0.4562,
0.454,0.4516,0.449,0.4462,0.445,0.4426,0.441,0.4395, 0.437,0.434,
0.4304,0.427};

lipacoordpic=Show[RegionPlot[Rectangle[{-2.1,-0.4},{-1,0.4}],BoundaryStyle->None,PlotStyle->Texture[lipaslice]]];

parameterlist=Table[(0.14)*E^({j+4}*I*Pi/32) + (-1.9-0.15I),{j,0,64}] /.{x_}->x;

c[i_]:=parameterlist[[i]];

newparametersplitpair={{c[1]-0.01-0.01I,5},{c[2]-0.02-0.015I,5},{c[3]-0.03-0.02I,5},{c[4]-0.04-0.026I,5},{c[5]-0.049-0.034I,5},{c[6]-0.05-0.042I},{c[7]-0.052-0.052I,5},{c[8]-0.058-0.06I,5},{c[9]-0.052-0.062I,5},{c[10]-0.052-0.066I,5},{c[11]-0.052-0.068I,5},{c[12]-0.052-0.07I,5},{c[13]-0.052-0.072I,5},{c[14]-0.052-0.074I,5}, {c[15]-0.052-0.074I,5},{c[16]-0.052-0.074I,5},{c[17]-0.052-0.074I,5},{c[18]-0.05-0.074I,5},{c[19]-0.03-0.074I,5}, {c[20]-0.01-0.07I,5},{c[21]+0.01-0.065I,5},{c[22]+0.028-0.06I,5},{c[23]+0.048-0.052I,5},{c[24]+0.068-0.044I,5},{c[25]+0.076-0.034I,5},{c[26]+0.096-0.028I,5},{c[27]+0.110-0.021I,5},{c[28]+0.126-0.012I,5},{c[29]+0.144-0.006I,5},{c[30]+0.158-0.00I,5},{c[31]+0.172+0.008I,5},{c[32]+0.184+0.016I,5},{c[33]+0.192+0.022I,5}, {c[34]+0.2+0.032I,5},{c[35]+0.21+0.038I,5},{c[36]+0.22+0.044I,5},{c[37]+0.23+0.06I,5}, {c[38]+0.23+0.08I,5},{c[39]+0.226+0.1I,5},{c[40]+0.212+0.12I,5},{c[41]+0.194+0.144I,5},{c[42]+0.168+0.165I,5}};

newparameterlist=Map[First,newparametersplitpair];

aherdlist=Reverse[1.6E^(I*(-Pi/12))*Conjugate[newparameterlist]+1.02-0.846I];

img05i=;

M05i=Show[RegionPlot[Rectangle[{-2.1,-0.4},{-1,0.4}],BoundaryStyle->None,PlotStyle->Texture[img05i]]];

S[25]:=;

R1[i_]=Show[RegionPlot[Rectangle[{-0.3378086898,0.8019635263},{0.1,1.13}],BoundaryStyle->None,PlotStyle->Texture[S[i]]]];
R2[i_]=Show[RegionPlot[Rectangle[{-0.3378086898,0.8020639372},{0.1,1.13}],BoundaryStyle->None,PlotStyle->Texture[S[i]]]];

M[i_]:= If[1<i<54 ,R1[i],R2[i]]

kochparams=Table[(0.066)*E^({j+4}*I*Pi/16) + (-0.27+1.04I)/.{x_}->x,{j,0,32}];

b[25]:=0.0524857+0.049287I;

twostringherd=;

twostringherdb=0.2+0.3I;

bbherdlist=newparameterlist-0.46-0.51I;

offbasilicaherd=Table[(0.02)*E^({j+4}*I*Pi/16) + (-1.41-.23I)/.{x_}->x,{j,0,32}];

offbasilicaherdper5=Table[(0.06)*E^({j+4}*I*Pi/16) + (-1.26-.52I)/.{x_}->x,{j,0,32}];

mappedrayshowray[th_,b_,c_,n_]:=Graphics3D[gl2[Map[nestedf[#,b,c,n]&,Map[i2[#]&,raypath2[th,2+n-1,10,b,c]]]]]

familymappedray[family_,b_,c_,n_]:=Map[mappedrayshowray[#,b,c,n]&,family]

cmappedrayshowray[th_,b_,c_,n_,colour_]:=Graphics3D[{colour,gl2[Map[nestedf[#,b,c,n]&,Map[i2[#]&,raypath2[th,2+n-1,18,b,c]]]]}]

cfamilymappedray[family_,b_,c_,n_,colour_]:=Map[cmappedrayshowray[#,b,c,n,colour]&,family]

exportgif[name_,table_]:=Export[name,table,"AnimationRepetitions"-> ∞]

Show[mappingforwardpicture[0.125,-1.124,10,2]]