(*******************************************************************************)
(****** Procedure qui calcule les voisinages et prépare le tracé des liens *****)
 (*******************************************************************************)

voisins[ppn_,d_]:=Module[{i,j,ns,dd1,dd2,vois,a,dist},
dd1=1.05*d;
dd2=d*0.8;
ns=Length[ppn];
a=Reap[
Do[
Do[
{dist=Norm[ppn[[i]]-ppn[[j]]];
If[dist<dd1∧dist>dd2,Sow[{i,j}]]
},
{j,i+1,ns}],{i,1,ns-1}];
];(*fin reap*)
vois=a[[2,1]];
vois];

preparegraph[pp_,vois_]:=Module[{i,j,nv,liens,a},
nv=Length[vois];
a=Reap[
Do[
Sow[{Line[{pp[[vois[[i,1]]]],pp[[vois[[i,2]]]]}]}],{i,1,nv}]
]; (** fin reap **)
liens=a[[2,1]];
;liens];

 (*************************************************************)
(********* Entrée dans programme principal *******************)
(*************************************************************)
(*************************************************************)

(** serie d'insructions décidant ce qui est imprimé ou dessiné **)
printinfo=False;
plotgraph=False;

(**************************** nombre d'inflations ********************************************)
(**** c'est ce nombre qu'il faut varier pour construire les approximants de taille variable***)
(*********************************************************************************************)
nbinfl=2;
(********************************************************************************************)
(********************************************************************************************)

(**************************quelques paramètre d'initialisation ****************************)
startcpu=TimeUsed[];
startreel=SessionTime[];
g=Sqrt[2]/2;
eperp={{1,g, 0,-g},{0,g,1,g}}; (*projecteur sur l'espace E_perp *)
epar={{0,-g, 1,-g},{1,-g, 0, g}};(*projecteur sur l'espace E_par *)
Print["e_par : ",MatrixForm[epar]];
gg=1/2;
epariso={{0,-gg, 1,-gg},{1,-gg, 0,gg}};(*projecteur sur l'espace E_par version "isométrique *)

e=(({{1, 0, 0, 0, -1, 0, 0, 0}, {0, 1, 0, 0, 0, -1, 0, 0}, {0, 0, 1, 0, 0, 0, -1, 0}, {0, 0, 0, 1, 0, 0, 0, -1}})); (*e : les 8 vecteurs pointant sur les 1ers voisins dans Z4 *)
Print["e : ",MatrixForm[e]];
M={{1,-1,0,1},{-1, 1,-1, 0},{0,-1, 1,-1},{1,0,-1,1}};(* La matrice d'inflation - version Pavel*)
MatrixForm[M]

(* 4 points initiaux:*)
f={{0,0,1,1},{-1,-1,-1,-1},{0,1,1,0},{0,0,0,0}};
Print["f : ",MatrixForm[f]];

(*couleur fonction du niveau d'inflation *)
typecoul={Black,Red,Orange,Green,Blue,Cyan,Yellow}

(****** projection des  4 points sur eperp et epar ******)
Print["projection des  4 points sur eperp et epar"];
Print["sur eperp :",MatrixForm[N[eperp.f]]];
Print["sur epar :",MatrixForm[N[epar.f]]];
ftriangle=Transpose[N[epar.f]];
coul=Table[typecoul[[1]],{i,1,4}];

(*************** paramètres pour le tracé ***************)
npts=4;
rperp=0.1/nbinfl;(*rayon des disques sur eperp *)
rpar=0.5/nbinfl;(*rayon des disques sur epar *)
If[nbinfl>4,{rperp=0.05/nbinfl;rpar=1.5/nbinfl}];
If[nbinfl<2,{rperp=0.02/nbinfl;rpar=0.25/nbinfl}];
(************************)
(* debut de l'itération *)
(************************)

Do[{   (* do iter **)
Print["iteration : ",iter];
(** inflation sur les npts points déjà créés **)
f=M.f;
nf=Length[nf];
If[printinfo,Print["f : ",MatrixForm[f]]];
If[printinfo,Print[" transpose de f : ",MatrixForm[Transpose[f]]]];
Print["Apres entree dans  la boucle , projection des  npts points sur eperp et epar"];

fperpnpts=Transpose[N[eperp.f]];(*projection sur eperp des points anciennement créés*)
fparnpts=Transpose[N[epar.f]];(*projection sur epar des points anciennement créés*)
fcarre=Transpose[N[epar.f[[;;,1;;4]]]]; (* projection des quatre sommets du carré *)

If[printinfo,Print[" fcarré : ",ftriangle]];
If[printinfo,Print["sur eperp :",MatrixForm[fperpnpts]]];
If[printinfo,Print["sur epar :",MatrixForm[fparnpts]]];

p=Table[0,{i,1,4},{j,1,8*npts}];
Uns=Transpose[ConstantArray[1,{1,npts}]];

(*contruction des nouveaux points*)
Do[{
p[[;;,(m-1)*npts+1;;m*npts]]=ArrayReshape[Flatten[f+Transpose[{e[[;;,m]]}].{Uns}],{4,npts}];
If[printinfo,Print["m :",m, "  p  :",MatrixForm[p]]]
},{m,1,8}];

(*****************************************)
(* Find unique points among the new ones*)
(*****************************************)
nbtemp=8*npts;
CC={p[[;;,1]]};
Do[
If[¬MemberQ[CC,p[[;;,j]]],
CC=AppendTo[CC,p[[;;,j]]]]
,{j,2,nbtemp}];
pu=Transpose[CC];
If[printinfo,Print[" pu : ",MatrixForm[pu]]];

(***********************************************)
(* Préparation pour la selection dans E_perp **)
(* ********************************************)

v=N[Transpose[eperp.pu]]; (*les nouveaux points projetés sur E_perp *)
If[printinfo,Print[" v (eperp) : ",MatrixForm[v]]];
nbv=Length[v];
idxdom={};

(*listes (à venir) correspondantes aux site acceptés sur eprp*)
pua={};zza={};ppn={}; ppniso={};
v=Abs[v];(* pour simplifier la sélection *)


Print["*************************"];
Print["entrée dans la boucle de sélection e_perp"];
Print["*************************"];

(*************************************************************)
(* selection des points dans la zone d'acceptance octogonale *)
(*************************************************************)

eps=0.000001;(* Tolerance*)
Do[If[v[[i,1]]<(0.5+g)∧v[[i,2]]<(0.5+g)∧(v[[i,1]]+v[[i,2]])<(1+g),
{ccc=N[epar.pu[[;;,i]]];
ccci=N[epariso.pu[[;;,i]]];
(****** selection des points dans le carré ******)
If[((ccc[[2]]+eps)>fcarre[[1,2]])∧ccc[[1]]<(fcarre[[2,1]]+eps)∧ ccc[[2]]<(fcarre[[4,2]]+eps)∧(ccc[[1]]+eps)>fcarre[[1,1]],
{
ppn=AppendTo[ppn,ccc];
ppniso=AppendTo[ppniso,ccci];
idxdom=AppendTo[idxdom,i];(*numéro des sites sélectionnés*)
pua=AppendTo[pua,pu[[;;,i]]];
}]
}],{i,1,nbv}];
If[printinfo,Print[" index des points sélectionnés : ",idxdom]];

(*pua contient les points sélectionnés à 4d *)
npua=Length[pua];
Print[" nb de points sélectionnés : ",npua];
pua=Transpose[pua];

If[printinfo,Print["pua : ",MatrixForm[pua]]];
If[printinfo,Print["pua : ",MatrixForm[Transpose[pua]]]];
If[printinfo,Print["ppn : ",MatrixForm[ppn]]];

If[printinfo,Print[" transpose de f : ",Transpose[f]]];
f=Transpose[f];
pua=Transpose[pua];
ninit=Length[f];
f=Join[f,pua];
npts=Length[f];
Do[coul=AppendTo[coul,typecoul[[iter+1]]],{j,ninit+1,npts}];
f=Transpose[f];
If[printinfo,Print["f en fin de boucle  : ",MatrixForm[f]]];
Print["npts en fin de boucle : ",npts]
},{iter,1,nbinfl}];

Print["************************"];
Print[" Fin du calcul itératif, les sites sont créés " ];
Print["***********************"];
Print["  durée CPU depuis le début : ",TimeUsed[]-startcpu, "  durée réelle : ",SessionTime[] -startreel]  ;
Print[];

(*nf=Length[f];*)
v=N[Transpose[eperp.f]];
za=Table[{coul[[i]],Disk[v[[i]],rperp]},{i,1,npts}];

(***********************************************)
(*********** calcul des voisinages *************)
(***********************************************)

ppn=N[Transpose[epar.f]];
ppniso=N[Transpose[epariso.f]];
ns=Length[ppn];

(** translation de la maille pour l'accrocher à l'origine **)
Do[{
ppn[[i,1]]=ppn[[i,1]]-ppn[[1,1]],ppn[[i,2]]=ppn[[i,2]]-ppn[[1,2]],
ppniso[[i,1]]=ppniso[[i,1]]-ppniso[[1,1]],ppniso[[i,2]]=ppniso[[i,2]]-ppniso[[1,2]]}
,{i,2,ns}
];
ppn[[1]]={0,0};
ppniso[[1]]={0,0};
(******************************************************)

(*construction des aretes *)

Print["entrée dans la procédure de calcul des liens "];
vois=voisins[ppn,1];
Print["sortie de la procédure de calcul des liens "];
Print["  durée CPU depuis le début : ",TimeUsed[]-startcpu, "  durée réelle : ",SessionTime[] -startreel]  ;

If[plotgraph,

{
(****************************************)
(***** trace des points sur E_parellele ********)
(****************************************)
liens=preparegraph[ppn,vois];
liensiso=preparegraph[ppniso,vois];

Print["*********************************************************************"];
Print[" plot de la maille carrée pour un nombre d'iteration : ",nbinfl ];
Print["**********************************************************************"];

zza=Table[{coul[[i]],Disk[ppn[[i]],rpar]},{i,1,npts}];
If[printinfo,Print[MatrixForm[ppn]]];
If[printinfo,Print[MatrixForm[f]]];

zzaiso=Table[{coul[[i]],Disk[ppniso[[i]],rpar]},{i,1,npts}];
If[printinfo,Print[MatrixForm[ppn]]];
If[printinfo,Print[MatrixForm[f]]];

maille={Dashed, Thin, Red,Line[{{ppn[[1,1]],ppn[[1,2]]},{ppn[[2,1]],ppn[[2,2]]},{ppn[[3,1]],ppn[[3,2]]},{ppn[[4,1]],ppn[[4,2]]},{ppn[[1,1]],ppn[[1,2]]}}]};

mailleiso={Dashed, Thin, Red,Line[{{ppniso[[1,1]],ppniso[[1,2]]},{ppniso[[2,1]],ppniso[[2,2]]},{ppniso[[3,1]],ppniso[[3,2]]},{ppniso[[4,1]],ppniso[[4,2]]},{ppniso[[1,1]],ppniso[[1,2]]}}]};

Graphics[{maille,zza,liens}]//Print;
Graphics[{mailleiso,zzaiso,liensiso}]//Print
}];
Print["  durée CPU depuis le début : ",TimeUsed[]-startcpu, "  durée réelle : ",SessionTime[] -startreel]  ;

l=Length[ppniso];
lv=Length[vois];
ppniso=Sqrt[2]*ppniso;

Length[ppn]

vois[[46]]

Norm[ppniso[[9]]-ppniso[[250]]]

ppnisonorm3=Sqrt[2]*ppniso;

Export["ppnisonorm3.dat",ppnisonorm3]

vois

Norm[ppnisonor[[1]]-ppnisonor[[13]]]

ppniso[[13]]