Mathematica at Work
Introduzione ad un Software estremamente potente 
e ad alcune sue applicazioni nella ricerca e nel mondo del lavoro

INCA - ABACO Infrastutture di Calcolo a Basso Costo
Università degli Studi di Roma Tor Vergata

Corrado Falcolini
Dipartimento di Architettura
Università degli Studi di Roma Tre

Poligoni regolari ed approssimazione di π .

Rilievo, nuvole di punti e modelli in Architettura.

Equazioni differenziali. 

Coesistenza di attrattori.

Poliedri semi-regolari.

Poligoni regolari ed approssimazione di π .

Sia P(n) un poligono regolare di n lati.

Domande:

1) Qual è il minimo numero n tale che il rapporto tra il perimetro di P(n) e il diametro della circonferenza circoscritta  sia maggiore di 3.14?

2) Qual è il minimo valore di n tale che P(n)  non sia distinguibile da un cerchio?

Manipulate[Graphics[{{Red,Circle[{0,0},1]},Line[Table[{Cos[k*2Pi/n],Sin[k*2Pi/n]},{k,0,n,1}]]},ImageSize->500,PlotRange->{{-1.1,1.1},{-1.1,1.1}}],{n,3,100,1}]

Manipulate[Graphics[{Line[Table[{Cos[k*2Pi/n],Sin[k*2Pi/n]},{k,0,n,1}]]},ImageSize->500,PlotRange->{{-1.1,1.1},{-1.1,1.1}}],{n,3,100,1}]

(* (Subscript[ℓ, n])=2r*sin((π/n)) , (Subscript[π, n])=n*((Subscript[ℓ, n])/(2r))=n*sin((π/n)) *)

TableForm[Table[{n,N[n*Sin[Pi/n],7]},{n,3,60,1}],TableHeadings->{None,{"n","(Subscript[π, n])"}}]

appk[k_]:=Block[{n=3},While[n*Sin[Pi/n]<ToExpression[StringJoin[Drop[Characters[ToString[N[Pi,k+1]]],-1]]],n++];Return[n];];

TableForm[Table[{p,appk[p]},{p,1,10,1}]]

Poligoni regolari ed approssimazione di π .

Sia P(n) un poligono regolare di n lati.

Domande:

3) Per quali interi n il numero  cos(((2π)/n))  è rappresentabile per radicali?

4) Per quali interi n il poligono P(n) è costruibile con riga e compasso?

Cosa si può fare con Mathematica:

The function Cos[ ]

Table[Cos[2Pi/k],{k,1,20}]

The function ToRadicals[ ]

Table[ToRadicals[Cos[2Pi/k]],{k,1,11}]

note

c[n_]:=HoldForm[Cos[2Pi/n]]

Table[{c[k],Cos[2Pi/k]},{k,1,8}]

TableForm[Table[{k,c[k],Cos[2Pi/k]},{k,1,8}],TableHeadings->{None,{"k","HoldForm","Cos[((2π)/k)]"}}]

TableForm[Table[{n,c[n],"=",Cos[2Pi/n],ToRadicals[Cos[2Pi/n]]},{n,1,50}],TableHeadings->{None,{"k","Hold","","Cos[((2π)/k)]","ToRadicals"}}]

"constructible" numbers

supscr[a_,b_]:=If[b==1,a,Superscript[a,b]]

prod[n_]:=Block[{com=Apply[supscr,FactorInteger[n],1],res=n},
If[Length[com]>1,com=Apply[CenterDot,com],com=com[[1]]];
Return[com]]

costr[n_]:=Block[{com=Complement[FactorInteger[n],{{3,1},{5,1},{17,1},{257,1},{65537,1}}],res},
res=If[Length[com]==0||Length[com]==1&&com[[1,1]]==2,True,False,False];
Return[res];
];

a table with divisors

TableForm[Table[{Style[n,If[costr[n],FontColor->Red,FontColor->Black]],"=",prod[n],c[n],If[costr[n],"=",""],If[costr[n],ToRadicals[Cos[2Pi/n]],""]},{n,2,99}]]

constructible numbers: an error?

TableForm[Select[Table[{Style[n,If[costr[n],FontColor->Red,FontColor->Black]],prod[n],c[n],If[costr[n],"=",""],If[costr[n],ToRadicals[Cos[2Pi/n]],""]},{n,2,300}],!FreeQ[#,"="]&]]

the function costr[ ]

costr[n_]:=Block[{com=Complement[FactorInteger[n],{{3,1},{5,1},{17,1},{257,1},{65537,1}}],res},
res=If[Length[com]==0||Length[com]==1&&com[[1,1]]==2,True,False,False];
Return[res];
];

A theorem of Gauss

A famous theorem of Gauss states that for a constructible regular polygon n should be the product of any power of 2 and distinct prime Fermat’s numbers ( of type (2^(2^k))  ): 

so that, apart from simpler geometrical constructions, 

polygons with 3, 5, 17, 257, 65537 sides 
(the only prime Fermat’s numbers known up to know) 

are all constructible.

lsradic=Select[Table[k,{k,65540}],costr]

particular cases

TableForm[Table[{c[2^k],Cos[2Pi/2^k],ToRadicals[Cos[2Pi/2^k]]},{k,3,8}]]

Let’s go on ...

ToRadicals[Cos[2Pi/257]]

We want to represent  cos(((2π)/257)) by radicals, making use of some modular algebra as suggested in a very nice article of Christian Gottlieb (The Simple and Straightforward Construction of the Regular 257-gon, THE MATHEMATICAL INTELLIGENCER, VOLUME21, NUMBER1, 1999).

n = 5 :  Cos[2π/5]

definitions

listas[n_]:=Block[{ls,lscnt,i,j,k}, ls=Flatten[Table[Mod[3^(2i*2^(n-1))+3^((2k+1)*2^(n-1)),5],{i,0,2^(2-n)-1},{k,0,2^(2-n)-1}]];
lscnt=Table[Count[ls,j],{j,1,16}]; Return[Table[lscnt[[Mod[3^s,5]]],{s,0,2^(n-1)-1}]];];
vec[n_]:=Block[{res=Table[0,{2^n}],temp,k,h}, Subsuperscript[a,1,0]=-1; For[k=1,k≤2^(n-1),k++,
temp=Solve[x^2-Subsuperscript[a,k,n-1]x+Sum[listaprod[[n,h]]Subsuperscript[a,Mod[k+h-2,2^(n-1)]+1,n-1],{h,1,2^(n-1)}]==0,x];
res[[k]]=temp[[2-listarev[[n,k]],1,2]]; res[[k+2^(n-1)]]=temp[[1+listarev[[n,k]],1,2]];]; res];
sub[n_]:=Table[Subsuperscript[a,k,n]->vec[n][[k]],{k,1,2^n}]

lists

listaprod=Table[listas[k],{k,1,1}]

listarev=Table[Table[Boole[TrueQ[Chop[N[Sum[
Exp[2Pi I/5*3^(2^n*j+k)]-Exp[2Pi I/5*3^(2^n*j+k+2^(n-1))],
{j,0,2^(4-n)-1}]]]<0]],{k,0,2^(n-1)-1}],{n,1,1}]

result

Simplify[vec[1][[1]] ]/2

Cos[2Pi/5]

n = 17 :  Cos[2π/17]

definitions

listas[n_]:=Block[{ls,lscnt,i,j,k}, ls=Flatten[Table[Mod[3^(2i*2^(n-1))+3^((2k+1)*2^(n-1)),17],{i,0,2^(4-n)-1},{k,0,2^(4-n)-1}]];
lscnt=Table[Count[ls,j],{j,1,16}]; Return[Table[lscnt[[Mod[3^s,17]]],{s,0,2^(n-1)-1}]];];
vec[n_]:=Block[{res=Table[0,{2^n}],temp,k,h}, Subsuperscript[a,1,0]=-1; For[k=1,k≤2^(n-1),k++,
temp=Solve[x^2-Subsuperscript[a,k,n-1]x+Sum[listaprod[[n,h]]Subsuperscript[a,Mod[k+h-2,2^(n-1)]+1,n-1],{h,1,2^(n-1)}]==0,x];
res[[k]]=temp[[2-listarev[[n,k]],1,2]]; res[[k+2^(n-1)]]=temp[[1+listarev[[n,k]],1,2]];]; res];
sub[n_]:=Table[Subsuperscript[a,k,n]->vec[n][[k]],{k,1,2^n}]

lists

listaprod=Table[listas[k],{k,1,3}]

listarev=Table[Table[Boole[TrueQ[Chop[N[Sum[
Exp[2Pi I/17*3^(2^n*j+k)]-Exp[2Pi I/17*3^(2^n*j+k+2^(n-1))],
{j,0,2^(4-n)-1}]]]<0]],{k,0,2^(n-1)-1}],{n,1,3}]

result

duecos17=Simplify[Simplify[vec[3][[1]] /. sub[2] ] /. sub[1] ]

ToRadicals[Cos[2Pi/17]]

N[Cos[2Pi/17],300]

N[1/2  duecos17,300]

Costruzione di poligoni regolari.

Si parte da due punti dati, che possono essere scelti uno come origine e l'altro per definire l'unità di misura. 
Se con operazioni lecite riusciamo a costruire il poligono regolare di n lati saremo anche riusciti a trovare la rappresentazione del Cos(2Pi/n) per radicali, visto che le eventuali intersezioni saranno soluzioni di sistemi di primo, di secondo grado o riducibili a questi.

retta e cerchio

retta[{{x1_,y1_},{x2_,y2_}}]:=Simplify[(y-y1)(x2-x1)-(y2-y1)(x-x1)==0];
circ[{{x1_,y1_},{x2_,y2_}}]:=Expand[(x-x1)^2+(y-y1)^2-(x2-x1)^2-(y2-y1)^2==0];

graphical functions

lista[{{x1_,y1_},{x2_,y2_}}]:=Block[{res=Map[Map[#[[2]]&,Sort[#]]&,Solve[{retta[{{x1,y1},{x2,y2}}],circ[{{x1,y1},{x2,y2}}]}]]},
res=Union[res,Map[Map[#[[2]]&,Sort[#]]&,Solve[{retta[{{x1,y1},{x2,y2}}],circ[{{x2,y2},{x1,y1}}]}]]];
res=Union[res,Map[Map[#[[2]]&,Sort[#]]&,Solve[{circ[{{x1,y1},{x2,y2}}],circ[{{x2,y2},{x1,y1}}]}]]];
Return[res];
];
graf[{{x1_,y1_},{x2_,y2_}}]:=Show[Graphics[{
{PointSize[0.01],Blue,Point[Complement[lista[{{x1,y1},{x2,y2}}],{{x1,y1},{x2,y2}}]]},
{PointSize[0.01],Green,Point[{{x1,y1},{x2,y2}}]},
{Thickness[.002],Black,Line[{5{x1,y1}-4{x2,y2},5{x2,y2}-4{x1,y1}}]},
Circle[{x1,y1},Sqrt[(x2-x1)^2+(y2-y1)^2]],
Circle[{x2,y2},Sqrt[(x2-x1)^2+(y2-y1)^2]]
}]]

other graphical functions

passograf[ls_]:=Block[{i,j,res={Show[Graphics[{PointSize[0.02],Red,Point[ls]}]]}},
For[i=1,i<Length[ls],i++,
For[j=i+1,j<=Length[ls],j++,
AppendTo[res,graf[{ls[[i]],ls[[j]]}]];
];
];
Return[res];
];
listaeq[{{x1_,y1_},{x2_,y2_}}]:={retta[{{x1,y1},{x2,y2}}],circ[{{x1,y1},{x2,y2}}],circ[{{x2,y2},{x1,y1}}]};
passoeq[ls_]:=Block[{i,j,res={}},
For[i=1,i<Length[ls],i++,
For[j=i+1,j<=Length[ls],j++,
res=Union[res,listaeq[{ls[[i]],ls[[j]]}]];
];
];
Return[res];
];
listagraf[ls_]:=Manipulate[Show[{
Graphics[{PointSize[0.03],Red,Point[ls]}],
Take[passograf[ls],k]},PlotRange->{{-5,6},{-5,5}}],
{k,1,Length[passograf[ls]],1}]
listagraf[lsin_,lstep_,lsfin_]:=Manipulate[Show[{graf[lsin],Graphics[{Thickness[.004],Red,Circle[{0,0},1]}],
Graphics[{PointSize[0.03],Red,Point[lstep]}],
Graphics[{PointSize[0.03],Blue,Point[lsin]}],
Graphics[{PointSize[0.009],Brown,Point[lsfin]}],
Take[passograf[lstep],k+1]},PlotRange->{{-5,6},{-5,5}},ImageSize->600],
{k,0,Length[passograf[lstep]]-1,1}]

The set of constructible points

ls0={{0,0},{1,0}};
passoeq[ls0]
solve[{eq1_,eq2_}]:=Block[{sl=Solve[{eq1,eq2}],res={}},If[FreeQ[sl,Complex],res=Map[Map[#[[2]]&,Sort[#]]&,sl]];Return[res];] ;
passocomp[ls_]:=Block[{i,j,res={}},
For[i=1,i<Length[ls],i++,
For[j=i+1,j<=Length[ls],j++,
res=Union[res,solve[{ls[[i]],ls[[j]]}]];
];
];
Return[res];
];

data

graf[{{0,0},{1,0}}]

ls1=passocomp[passoeq[ls0]]

2 steps

listagraf[ls1]

all constructible points in 2 steps

passoeq[ls1]

ls2=passocomp[passoeq[ls1]]

listagraf[ls0,ls1,ls2]

Costruzione del pentagono regolare.

definitions

grafc[{{x1_,y1_},{x2_,y2_}}]:=Show[Graphics[{
Circle[{x1,y1},Sqrt[(x2-x1)^2+(y2-y1)^2]]
}]]
grafr[{{x1_,y1_},{x2_,y2_}}]:=Show[Graphics[{
{Thickness[.002],Black,Line[{3.2{x1,y1}-2.2{x2,y2},3.2{x2,y2}-2.2{x1,y1}}]}
}]]
listavertici[n_]:=Table[{Cos[2k*Pi/n],Sin[2k*Pi/n]},{k,0,n}]
retta[{{x1_,y1_},{x2_,y2_}}]:=Simplify[(y-y1)(x2-x1)-(y2-y1)(x-x1)==0];
circ[{{x1_,y1_},{x2_,y2_}}]:=Expand[(x-x1)^2+(y-y1)^2-(x2-x1)^2-(y2-y1)^2==0];
int[eq1_,eq2_,n_]:={x,y} /. Solve[{eq1,eq2}][[n]]
int[eq1_,eq2_]:={x,y} /. Solve[{eq1,eq2}]
lsgpent={
Graphics[{PointSize[0.02],Red,Point[{{0,0},{1,0}}]}],
Graphics[{
{Text[StyleForm["O",FontSize->14],{0,0},{1.5,2}]},
{Text[StyleForm["(Subscript[P, 0])≡(1,0)",FontSize->14],{1,0},{0,2}]}
}],
grafr[{{0,0},{1,0}}],
grafc[{{0,0},{1,0}}],
Graphics[{{PointSize[0.02],Blue,Point[P1]},
{Text[StyleForm["(Subscript[P, 1])",FontSize->14],P1,{1.5,2}]}}],
grafc[{{-1,0},{1,0}}],
grafc[{{1,0},{-1,0}}],
Graphics[{{PointSize[0.02],Blue,Point[P2]},
{Text[StyleForm["(Subscript[P, 2])",FontSize->14],P2,{1.5,-2}]}}],
grafr[{{0,0},P2}],
Graphics[{{PointSize[0.02],Blue,Point[P3]},
{Text[StyleForm["(Subscript[P, 3])",FontSize->14],P3,{1.5,-1}]}}],
grafc[{P3,{0,0}}],
Graphics[{{PointSize[0.02],Blue,Point[P41]},
{Text[StyleForm["(Subscript[P, (Subscript[4, 1])])",FontSize->14],P41,{1.5,-1}]}}],
Graphics[{{PointSize[0.02],Blue,Point[P42]},
{Text[StyleForm["(Subscript[P, (Subscript[4, 2])])",FontSize->14],P42,{-1.5,-1}]}}],
grafr[{P41,P42}],
Graphics[{{PointSize[0.02],Blue,Point[P5]},
{Text[StyleForm["(Subscript[P, 5])",FontSize->14],P5,{1.5,-1.5}]}}],
grafc[{P5,{1,0}}],
Graphics[{{PointSize[0.02],Blue,Point[P6]},
{Text[StyleForm["(Subscript[P, 6])",FontSize->14],P6,{1.5,-1.5}]}}],
grafc[{{1,0},P6}],
Graphics[{{PointSize[0.02],Blue,Point[P7]},
{Text[StyleForm["(Subscript[P, 7])",FontSize->14],P7,{1,-1.5}]}}],
Graphics[{{Thickness[.004],Red,Circle[{0,0},1]}}],
Graphics[{Thickness[.004],Red,Line[Take[listavertici[5],2]]}],
Graphics[{Thickness[.004],Green,Line[{{0,0},listavertici[5][[2]]}]}],
Graphics[{Thickness[.004],Green,Line[{{1,0},{0,0}}]}],
Graphics[{Thickness[.004],Purple,Circle[{0,0},.15,{0,2Pi/5}]}],
Graphics[{Text[StyleForm[2Pi/5,FontSize->11,FontColor->Purple],{0,0},{-3,-1.4}]}],
Graphics[{Dashed,Thickness[.004],Red,Line[Take[listavertici[5],{2,3}]]}],
Graphics[{Dashed,Thickness[.004],Red,Line[Take[listavertici[5],{3,4}]]}],
Graphics[{Dashed,Thickness[.004],Red,Line[Take[listavertici[5],{4,5}]]}],
Graphics[{Thickness[.004],Red,Line[Take[listavertici[5],{5,6}]]}],
Graphics[{Thickness[.004],Red,Line[Take[listavertici[5],{5,6}]]}],
Graphics[{Thickness[.004],Red,Line[Take[listavertici[5],{5,6}]]}],
Graphics[{Thickness[.004],Red,Line[Take[listavertici[5],{5,6}]]}],
Graphics[{Thickness[.004],Red,Line[Take[listavertici[5],{5,6}]]}],
Graphics[{Thickness[.004],Red,Line[Take[listavertici[5],{5,6}]]}],
Graphics[{Thickness[.004],Red,Line[Take[listavertici[5],{5,6}]]}],
Graphics[{Thickness[.004],Red,Line[Take[listavertici[5],{5,6}]]}],
Graphics[{Thickness[.004],Red,Line[Take[listavertici[5],{5,6}]]}],
Graphics[{Thickness[.004],Red,Line[Take[listavertici[5],{5,6}]]}]
};

One (simple) example via constructible points

Manipulate[Block[{P1,P2,P3,P41,P42,P5,P6,P7},
P1=int[retta[{{0,0},{1,0}}],circ[{{0,0},{1,0}}],1];
P2=int[circ[{{-1,0},{1,0}}],circ[{{1,0},{-1,0}}],2];
P3=int[retta[{{0,0},P2}],circ[{{0,0},{1,0}}],2];
{P41,P42}=int[circ[{{0,0},{1,0}}],circ[{P3,{0,0}}]];
P5=int[retta[{{0,0},P2}],retta[{P41,P42}],1];
P6=int[retta[{{0,0},P2}],circ[{P5,{1,0}}],1];
P7=int[circ[{{0,0},{1,0}}],circ[{{1,0},P6}],2];
Show[Take[lsgpent,k],PlotRange->{{-2.2,2.2},{-2.2,2.2}},Frame->True,ImageSize->600]],{k,1,Length[lsgpent],1}]

moviepent=Table[Block[{P1,P2,P3,P41,P42,P5,P6,P7},
P1=int[retta[{{0,0},{1,0}}],circ[{{0,0},{1,0}}],1];
P2=int[circ[{{-1,0},{1,0}}],circ[{{1,0},{-1,0}}],2];
P3=int[retta[{{0,0},P2}],circ[{{0,0},{1,0}}],2];
{P41,P42}=int[circ[{{0,0},{1,0}}],circ[{P3,{0,0}}]];
P5=int[retta[{{0,0},P2}],retta[{P41,P42}],1];
P6=int[retta[{{0,0},P2}],circ[{P5,{1,0}}],1];
P7=int[circ[{{0,0},{1,0}}],circ[{{1,0},P6}],2];
Show[Take[lsgpent,k],PlotRange->{{-2.2,2.2},{-2.2,2.2}},Frame->True,ImageSize->600]],{k,1,Length[lsgpent],1}];
Export["moviepent.gif",moviepent,"DisplayDurations"->2]

Rilievo,nuvole di punti e modelli in Architettura.

Circo Massimo

Rilievo: laser scanner

input: la nuvola di punti

nuvola=ReadList["/Users/corradofalcolini/Desktop/Nepal 2015/lezioni/f3.xyz",Number,RecordLists->True];

Length[nuvola]

TableForm[Take[nuvola,3]]

nuvola1=Map[Drop[#,-3]&,nuvola];

mean=Map[Mean,Transpose[nuvola1]]

nuvola1tr=Map[#-mean&,nuvola1];

nuvola1tr[[1]]

Object: portion of a column

Graphics3D[{PointSize[.0005],Point[nuvola1tr]},ImageSize->1000]

problema: sezione verticale della colonna

Graphics3D[{PointSize[.0005],Point[Select[nuvola1tr,0<#[[2]]<0.1&]]}]

orientamento: cilindro ottimale

vec={1,a,c};offset={0,b,d};

perp[vec1_,vec_,offset_]:=vec1-offset-Projection[vec1-offset,vec,Dot]

Together[(perp[{x,y,z},vec,offset].perp[{x,y,z},vec,offset]-r^2)^2]

perpvar2[{x_,y_,z_}]:={x,y,z}-{0,b,d}-((x+a y+c z-a b-c d)/(1+(a^2)+(c^2))){1,a,c}

distanza di un punto di coordinate (x,y,z)  da un generico cilindro di parametri a,b,c,d

Together[(perpvar2[{x,y,z}].perpvar2[{x,y,z}]-r^2)^2]

numer=Expand[Numerator[Together[(perpvar2[{x,y,z}].perpvar2[{x,y,z}]-r^2)^2]]]

expo=Map[First,CoefficientRules[numer,{x,y,z}]]

par=Map[Last,CoefficientRules[numer,{x,y,z}]]

funzione da minimizzare: distanza di tutti i punti della nuvola dal cilindro

Timing[coef=Map[Sum[Apply[Times,nuvola1tr[[k]]^#],{k,1,Length[nuvola1tr]}]&,expo]]

(Apply[Plus,Thread[coef*par]])/(1+a^2+c^2)^2

func=(Expand[Apply[Plus,Thread[coef*par]]])/(1+a^2+c^2)^2

trovare il minimo: l’asse del cilindro sarà anche l’asse della colonna

Timing[{m2,rt2}=FindMinimum[func,{a,1},{b,0},{c,0.4},{d,0},{r,0.5},MaxIterations->500]]

Show[{Graphics3D[{PointSize[.0005],Point[nuvola1tr]}],Graphics3D[{Opacity[.5],Cylinder[{{x,a*x+b,c*x+d}
/.x->-1,{x,a*x+b,c*x+d}/.x->1},r]}/.rt2]},PlotRange->{{-1,1},{-1,1},{-1,1}},ImageSize->800]

la colonna può ora essere orientata verticalmente

Show[{Graphics3D[{PointSize[.0005],Rotate[Point[nuvola1tr],{{1,a*1+b,c*1+d}/.rt2,{0,0,1}}]}],Graphics3D[{Opacity[.5],Rotate[Cylinder[{{x,a*x+b,c*x+d}
/.x->-1,{x,a*x+b,c*x+d}/.x->1},r],{{1,a*1+b,c*1+d}/.rt2,{0,0,1}}]}/.rt2]},PlotRange->{{-1,1},{-1,1},{-1,1}},ImageSize->800,Axes->True]

le sezioni

Show[{Graphics3D[{PointSize[.0005],Rotate[Point[nuvola1tr],{{1,a*1+b,c*1+d}/.rt2,{0,0,1}}]}],Graphics3D[{Opacity[.5],Rotate[Cylinder[{{x,a*x+b,c*x+d}
/.x->-1,{x,a*x+b,c*x+d}/.x->1},r],{{1,a*1+b,c*1+d}/.rt2,{0,0,1}}]}/.rt2]},PlotRange->{{-1,1},{-1,1},{-1,1}},ImageSize->800,Axes->True,ViewPoint->Top]

il modello matematico: proiezione

mrt=RotationMatrix[{{1,a*1+b,c*1+d}/.rt2,{0,0,1}}]

Show[{Graphics[{PointSize[.0005],Point[Map[Drop[mrt.#,{3}]&,nuvola1tr]]}],Graphics[{Red,Thick,Circle[Drop[mrt.{x,a*x+b,c*x+d},{3}]
/.x->0,r]}/.rt2]},ImageSize->600]

sezioni

nuvola1trrt=Map[mrt.#&,nuvola1tr];

Graphics[{PointSize[.0005],Point[Map[Drop[#,-1]&,Select[nuvola1trrt,0<#[[3]]<.03&]]]}]

Manipulate[Graphics[{PointSize[.0005],Point[Map[Drop[#,-1]&,Select[nuvola1trrt,z<#[[3]]<z+.03&]]]}],{z,-.2,.8,.01}]

nuvola1sez=Select[nuvola1trrt,-.2<#[[3]]<.8&];

nuvola1sezrt=Map[mrt.#&,nuvola1sez];

Length[nuvola1sez]

modello matematico

centro=Drop[mrtsz.{x,a*x+b,c*x+d},{3}]/.x->0/.rt2sz

Manipulate[
Show[{Graphics[{PointSize[.0005],Point[Map[Drop[mrtsz.#,{3}]&,nuvola1sez2]]}],Graphics[{Red,Thick,Circle[Drop[mrtsz.{x,a*x+b,c*x+d},{3}]
/.x->0,r]}/.rt2sz],Graphics[{Blue,Point[{x0,y0}],Circle[{x0,y0},r0]}]},ImageSize->600],{{x0,centro[[1]]},-1,1},{{y0,centro[[2]]},-1,1},{{r0,r/.rt2sz},0.0001,1}]

il modello matematico della colonna

Show[{Graphics[{PointSize[.0005],Point[Map[Drop[mrtsz.#,{3}]&,nuvola1sez2]]}],Graphics[{Red,Thick,Circle[Drop[mrtsz.{x,a*x+b,c*x+d},{3}]
/.x->0,r]}/.rt2sz],Graphics[{Blue,Circle[{.026,.19},.5057]}],Graphics[{Blue,Circle[{.026,.19},.57]}],Graphics[{Red,Point[{-.497,-.032}],Blue,Circle[{-.497,-.032},.057]}]},ImageSize->600]

modello matematico

rt[a_]:={{Cos[a],-Sin[a]},{Sin[a],Cos[a]}}

sez=Import["/Users/falco/Desktop/strasburgo 2013/colonna_sezione_tipo.jpg"]

Manipulate[Module[{n=48,r,arc,tc,sp=.004},
r=(c/n)r0;
tc=ArcCos[(r0-r^2/(2r0))/r0];
arc=1/2ArcTan[r0-r^2/(2r0),(1/((2r0))) (Sqrt[4 r^2 r0^2-r^4])];GraphicsRow[{sez,
Graphics[{Red,
Table[Rotate[Circle[{x0,y0}+r0{Cos[k*2Pi/n],Sin[k*2Pi/n]},r,{2k*Pi/n+Pi/2+arc,2k*Pi/n+3Pi/2-arc}],alfa,{x0,y0}],{k,2,n,2}],Table[Rotate[Circle[{x0,y0},r0,{tc+k*2Pi/n,-tc+(k+2)2Pi/n}],alfa,{x0,y0}],{k,2,n,2}],Dashed,Circle[{x0,y0},r0],Circle[{x0,y0},r0-r]},Axes->True,PlotRange->{{-1,1},{-435/562,435/562}}]},Spacings->-360]],{{x0,.524},-1,1},{{y0,.084},-1,1},{{r0,.379},.1,2},{{c,5.09117},1,10},{{alfa,.333},-1,1},ControlPlacement->Right]

il modello matematico della colonna

With[{r0=6/10,r=1/10},ContourPlot[{x^2+y^2==r0^2,(x-r0)^2+y^2==r^2},{x,-.7,.7},{y,-.7,.7},Epilog->{Line[{{0,0},r0{Cos[Pi/4],Sin[Pi/4]}}],Text[Style["(Subscript[r, 0])",15],{0.2,0.25}],Line[{{r0,0},{r0,0}+r{Cos[Pi/4],Sin[Pi/4]}}],Text[Style["r",15],{0.63,0.05}],Dashed,Circle[{0,0},r0-r],Red,PointSize[.01],Point[{x,y}/. Solve[{x^2+y^2==r0^2,(x-r0)^2+y^2==r^2},{x,y}][[2]]],Text[Style["P",15],{0.61,0.15}]},PlotPoints->30,Axes->True,ImageSize->400]]

(* coordinate del punto P (in figura) in funzione di r0 ed r *)

Solve[{x^2+y^2==r0^2,(x-r0)^2+y^2==r^2},{x,y}]

Manipulate[Module[{x0=0,y0=0,r0=1.188/2,c=5.4,r,arc,tc},
r=(c/n)r0;
tc=ArcCos[(r0-r^2/(2r0))/r0];
arc=1/2ArcTan[r0-r^2/(2r0),(1/((2r0))) (Sqrt[4 r^2 r0^2-r^4])];
Graphics[{
Table[Circle[{x0,y0}+r0{Cos[k*2Pi/n],Sin[k*2Pi/n]},r,{2k*Pi/n+Pi/2+arc,2k*Pi/n+3Pi/2-arc}],{k,2,n,2}],Table[Circle[{x0,y0},r0,{tc+k*2Pi/n,-tc+(k+2)2Pi/n}],{k,2,n,2}]},ImageSize->500,Axes->True,PlotRange->{{-.7,.7},{-.7,.7}}]],{{n,48},6,50,2}]

sez=Import["/Users/falco/Desktop/strasburgo 2013/colonna_sezione_tipo.jpg"]

ImageDimensions[sez]

Manipulate[GraphicsRow[{sez,Module[{n=48,r,arc,tc},
r=(c/n)r0;
tc=ArcCos[(r0-r^2/(2r0))/r0];
arc=1/2ArcTan[r0-r^2/(2r0),(1/((2r0))) (Sqrt[4 r^2 r0^2-r^4])];
Graphics[{
Table[Circle[{x0,y0}+r0{Cos[k*2Pi/n],Sin[k*2Pi/n]},r,{2k*Pi/n+Pi/2+arc,2k*Pi/n+3Pi/2-arc}],{k,2,n,2}],Table[Circle[{x0,y0},r0,{tc+k*2Pi/n,-tc+(k+2)2Pi/n}],{k,2,n,2}],Red,Dashed,Circle[{x0,y0},r0],Circle[{x0,y0},r0-r]},Axes->True,PlotRange->{{-1,1},{-1485/1432,1485/1432}}]]},Spacings->-360],{{x0,-.018},-1,1},{{y0,-.055},-1,1},{{r0,.982},.1,2},{{c,5.4},1,10}]

(*  per c=5.4 ed r0=1.188 si ha r = (c/n)r0 = (5.4/48)1.188 = 0.13365 come in figura *)

col[x0_,y0_,r0_,n_]:=Module[{c=5.4,tc,arc,r0ent,style=Opacity[1]},
tc=ArcCos[1-((c^2)/(2 (n^2)))];
arc=1/2ArcTan[1-((c^2)/(2 (n^2))),(1/(2 (n^2))) (Sqrt[4(c^2) (n^2)-(c^4)])];
r0ent=(650r0-9v)/600;
{Table[ParametricPlot3D[{x0+r0(Cos[k*2Pi/n]+(c/n)Cos[t]),y0+r0(Sin[k*2Pi/n]+(c/n)Sin[t]),v},{t,2k*Pi/n+Pi/2+arc,2k*Pi/n+3Pi/2-arc},{v,0,(8+1/3)2r0/3},Mesh->None,Lighting->"Neutral",PlotStyle->style],{k,2,n,2}],Table[ParametricPlot3D[{x0+r0ent(Cos[k*2Pi/n]+(c/n)Cos[t]),y0+r0ent(Sin[k*2Pi/n]+(c/n)Sin[t]),v},{t,2k*Pi/n+Pi/2+arc,2k*Pi/n+3Pi/2-arc},{v,(8+1/3)2r0/3,(8+1/3)2r0},Mesh->None,Lighting->"Neutral",PlotStyle->style],{k,2,n,2}],Table[ParametricPlot3D[{x0+r0*Cos[t],y0+r0*Sin[t],v},{t,tc+k*2Pi/n,-tc+(k+2)2Pi/n},{v,0,(8+1/3)2r0/3},Mesh->None,Lighting->"Neutral",PlotStyle->style],{k,2,n,2}],Table[ParametricPlot3D[{x0+r0ent*Cos[t],y0+r0ent*Sin[t],v},{t,tc+k*2Pi/n,-tc+(k+2)2Pi/n},{v,(8+1/3)2r0/3,(8+1/3)2r0},Mesh->None,Lighting->"Neutral",PlotStyle->style],{k,2,n,2}]}]

Show[{col[1,1,1.188/2,48],col[3,1,1.188/2,48],col[.5,3,1.188/2,48],col[2.5,3,1.188/2,48]},PlotRange->{{0,4},{0,4},{0,10}},ImageSize->500]

Show[col[0,0,1.188/2,48],PlotRange->{{-1,1},{-1,1},{0,10}},ImageSize->300]

Module[{n=48,c=1.6},RegionPlot3D[Apply[And,Prepend[Table[Norm[{x,y}-rt[k*2Pi/n].{Cos[Pi/n],0}]^2>(c*Sin[Pi/n])^2,{k,2,n,2}],x^2+y^2<1]],{x,-1,1},{y,-1,1},{z,0,10},Mesh->None,PlotPoints->60,BoxRatios->Automatic]]

Manipulate[Module[{n=48,r,arc,tc,sp=.004},
r=(c/n)r0;
tc=ArcCos[(r0-r^2/(2r0))/r0];
arc=1/2ArcTan[r0-r^2/(2r0),(1/((2r0))) (Sqrt[4 r^2 r0^2-r^4])];GraphicsRow[{Graphics[{PointSize[.0005],Point[Map[Drop[mrtsz.#,{3}]&,nuvola1sez2]]}],
Graphics[{Red,
Table[Rotate[Circle[{x0,y0}+r0{Cos[k*2Pi/n],Sin[k*2Pi/n]},r,{2k*Pi/n+Pi/2+arc,2k*Pi/n+3Pi/2-arc}],alfa,{x0,y0}],{k,2,n,2}],Table[Rotate[Circle[{x0,y0},r0,{tc+k*2Pi/n,-tc+(k+2)2Pi/n}],alfa,{x0,y0}],{k,2,n,2}],Dashed,Circle[{x0,y0},r0],Circle[{x0,y0},r0-r]},Axes->True,PlotRange->{{-1,1},{-600/562,600/562}}]},Spacings->-360]],{{x0,-.036},-1,1},{{y0,.158},-1,1},{{r0,.994},.1,2},{{c,5.24},1,10},{{alfa,.418},-1,1},ControlPlacement->Right]

Bibliografia

Grazie !