(* scegliere il file con la lista delle coordinate dei punti della "nuvola di punti" denominata nuvola1tr utilizzando il Menù Insert/FilePath *)

<<"/Users/corradofalcolini/Desktop/curve e superfici 2018/nuvola1tr.txt";

(* numero di punti della lista nuvola1tr *)

Length[nuvola1tr]

(* seleziona i primi 5 punti della lista nuvola1tr *)

Take[nuvola1tr,5]

(* Immagine con i punti della lista nuvola1tr. Notare l'orientamento del rocchio di colonna rappresentato *)

Graphics3D[{PointSize[.0005],Point[nuvola1tr]},Axes->True]

(* DETERMINAZIONE DELL'ASSE DELLA COLONNA: METODO 1 *)

(* seleziona i punti della lista nuvola1tr che distano dal piano y=-0.4 meno di 1/30. 
#[[2]]rappresenta la seconda coordinata del generico elemento della lista nuvola1tr *)

sel1=Select[nuvola1tr,Abs[#[[2]]+0.4]<1/30&];

(* seleziona i punti della lista nuvola1tr che distano dal piano 
y=0.5 meno di 1/30.  *)

sel2=Select[nuvola1tr,Abs[#[[2]]-0.5]<1/30&];

(* Immagine con i punti della lista nuvola1tr ed i punti, in blù, delle liste sel1 e sel2 selezionati. *)

Graphics3D[{PointSize[.0005],Point[nuvola1tr],Blue,PointSize[.003],Point[sel1],Point[sel2]},Axes->True]

(* rappresenta i punti della selezione sel2 sul piano y=0.5 . 
La proiezione dei punti sul piano parallelo al piano (x,z) si ottiene eliminando la coordinata y .
 *)

Graphics[{Blue,PointSize[.003],Point[Map[Delete[#,{2}]&,sel2]]},Axes->True,PlotRange->{{-0.6,0.7},{-0.9,0.4}}]

(* Cerco la circonferenza che meglio approssima dall'esterno i punti selezionati in sel2.
Il comando Manipulate permette la manipolazione dei valori (xc,yc) ed r 
che determinano l'equazione parametrica della circonferenza di centro (xc,yc) e raggio r.
I valori iniziali sono xc=0.1, yc=-0.158, r=0.588 .
 *)

Manipulate[Graphics[{Blue,PointSize[.003],Point[Map[Delete[#,{2}]&,sel2]],Black,Circle[{xc,yc},r]},Axes->True,PlotRange->{{-0.6,0.7},{-0.9,0.4}}],{{xc,0.1},-0.7,0.5},{{yc,-0.158},-0.5,0.3},{{r,0.588},0.01,1},ControlPlacement->Right]

(* Cerco la circonferenza che meglio approssima dall'esterno i punti selezionati in sel1.
I valori iniziali sono xc=-0.13, yc=-0.264, r=0.6 .
 *)

Manipulate[Graphics[{Blue,PointSize[.003],Point[Map[Delete[#,{2}]&,sel1]],Black,Circle[{xc,yc},r]},Axes->True,PlotRange->{{-0.8,0.5},{-0.9,0.4}}],{{xc,-0.13},-0.7,0.5},{{yc,-0.264},-0.5,0.3},{{r,0.6},0.01,1},ControlPlacement->Right]

(* I centri delle due sezioni individuate sono c1 e c2 *)

c1={-0.13,-0.4,-0.264};

c2={0.1,0.5,-0.158};

(* Lo pseudo asse della colonna passante per c1 e c2 *)

Graphics3D[{PointSize[.0005],Point[nuvola1tr],Blue,PointSize[.005],Point[sel1],Point[sel2],PointSize[.03],Point[{c1,c2}],Thickness[.01],Line[{c1,c2}]},Axes->True]

(* il vettore per c1 e c2 *)

c2-c1

(* la matrice di rotazione m che allinea il vettore c2-c1 con il versore dell'asse z *)

m=RotationMatrix[{c2-c1,{0,0,1}}]

(* la matrice di rotazione m in una forma più leggibile *)

MatrixForm[m]

(* nuvola rappresenta la nuvola raddrizzata applicando la matrice a tutti i punti di nuvola1tr *)

nuvola=Map[m.#&,nuvola1tr];

(* le due nuvole ed i loro pseudo assi *)

Graphics3D[{PointSize[.0005],Point[nuvola1tr],Point[nuvola],Thickness[.01],Line[{c1,c2}],Line[{m.c1,m.c2}]},Axes->True]

(* la proiezione dell'intera nuvola sul piano (x,y) *)

Graphics[{PointSize[.0005],Point[Map[Delete[#,{3}]&,nuvola]]},Axes->True]

(* ripetiamo la procedura precedente con sezioni orizzontali z=-0.4 e z=0.5 *)

sel3=Select[nuvola,Abs[#[[3]]+0.4]<1/30&];

sel4=Select[nuvola,Abs[#[[3]]-0.5]<1/30&];

(* liste di punti nel piano parallelo al piano (x,y) ottenute eliminando la terza coordinata z *)

sel3xy=Map[Delete[#,{3}]&,sel3];

sel4xy=Map[Delete[#,{3}]&,sel4];

Manipulate[Graphics[{Blue,PointSize[.003],Point[Map[Delete[#,{3}]&,sel3]],Black,Circle[{xc,yc},r]},Axes->True,PlotRange->{{-0.6,0.7},{-0.5,0.8}}],{{xc,0.008},-0.6,0.7},{{yc,0.183},-0.5,0.8},{{r,0.56},0.01,1},ControlPlacement->Right]

Manipulate[Graphics[{Blue,PointSize[.003],Point[Map[Delete[#,{3}]&,sel4]],Black,Circle[{xc,yc},r]},Axes->True,PlotRange->{{-0.6,0.7},{-0.5,0.8}}],{{xc,0.023},-0.6,0.7},{{yc,0.207},-0.5,0.8},{{r,0.56},0.01,1},ControlPlacement->Right]

c3={0.008,0.183,-0.4};

c4={.023,0.207,0.5};

m2=RotationMatrix[{c4-c3,{0,0,1}}]

(* nuvola2 rappresenta la nuvola raddrizzata nuovamente applicando la matrice a tutti i punti di nuvola *)

nuvola2=Map[m2.#&,nuvola];

(* la proiezione dell'intera nuova nuvola2 sul piano (x,y). Notare la maggiore precisione *)

Graphics[{PointSize[.0005],Point[Map[Delete[#,{3}]&,nuvola2]]},Axes->True]

(* DETERMINAZIONE DELL'ASSE DELLA COLONNA: METODO 2 *)

(* funzione distanza tra i punti P e Q nel piano *)

Expand[((dist[sel3xy[[2]],{xc,yc}])^2-r^2)^2]

(* funzione che esprime la distanza complessiva dei punti della selezione sel3xy
dal centro {xc,yc} *)

fun=Expand[Total[Table[((dist[sel3xy[[k]],{xc,yc}])^2-r^2)^2,{k,1,Length[sel3xy]}]]]

(* Cerco il minimo della funzione fun(xc,yc,r) *)

FindMinimum[fun,{{xc,0},{yc,0},{r,0.5}}]

(* chiamo val3 il valore minimo di fun(xc,yc,r) e sost3 i valori di xc, yc ed r da sostituire nell'equazione parametrica della circonferenza per ottenere quella ottimale *)

{val3,sost3}=FindMinimum[fun,{{xc,0},{yc,0},{r,0.5}}]

(* grafico della circonferenza ottimale e di quella esterna *)

Graphics[{PointSize[.005],Point[sel3xy],Red,Circle[{xc,yc},r]/.sost3,Blue,Circle[{xc,yc},r+.033]/.sost3,PointSize[.02],Point[{xc,yc}]/.sost3,Text[Style[(Subscript[C, 3]),Large],{xc+.1,yc}]/.sost3},Axes->True]

c3={xc,yc,-0.4}/.sost3

(* funzione che esprime la distanza complessiva dei punti della selezione sel4xy
dal centro {xc,yc} *)

fun=Expand[Total[Table[((dist[sel4xy[[k]],{xc,yc}])^2-r^2)^2,{k,1,Length[sel4xy]}]]]

(* chiamo val4 il valore minimo di fun(xc,yc,r) e sost4 i valori di xc, yc ed r da sostituire nell'equazione parametrica della circonferenza per ottenere quella ottimale *)

{val4,sost4}=FindMinimum[fun,{{xc,0},{yc,0},{r,0.5}}]

(* grafico della circonferenza ottimale e di quella esterna *)

Graphics[{PointSize[.005],Point[sel4xy],Red,Circle[{xc,yc},r]/.sost4,Blue,Circle[{xc,yc},r+.033]/.sost4,PointSize[.02],Point[{xc,yc}]/.sost4,Text[Style[(Subscript[C, 4]),Large],{xc+.1,yc}]/.sost4},Axes->True]

c4={xc,yc,0.5}/.sost4

(* selxy è la proiezione sul piano (x,y) dei punti di nuvola2 con 
la coordinata -0.4<z<0.5 *)

selxy=Map[Delete[#,{3}]&,Select[nuvola2,-.4<#[[3]]<.5&]];

(* equazione parametrica della circonferenza esterna ottimale gxy[t] *)

Manipulate[Graphics[{PointSize[.0005],Point[selxy],Blue,Circle[{xc,yc},r]},Axes->True,PlotRange->{{-0.6,0.7},{-0.5,0.8}}],{{xc,0.022},-0.6,0.7,.001},{{yc,0.205},-0.5,0.8,.001},{{r,0.57},0.01,1,.001},ControlPlacement->Right]

gxy[t_]:={0.022+0.57Cos[t],0.21+0.57Sin[t]}

(* raggio delle 24 scanalature: 0.057 *)

Show[{Graphics[{PointSize[.0005],Point[selxy],Red,Table[Circle[gxy[2k*Pi/24+Pi/23],.057],{k,1,24}],Table[Point[gxy[2k*Pi/24+Pi/23]],{k,1,24}]},PlotRange->{{-0.7,0.7},{-0.5,0.9}}],ParametricPlot[gxy[t],{t,0,2Pi},PlotStyle->Blue]}]

(* modello della sezione di una colonna scanalata *)

Manipulate[Module[{r0=.57,r,arc,tc,n=48,p=9},
r=r0/p;
tc=ArcCos[1-1/(2*p^2)];
arc=1/2ArcTan[1-1/(2*p^2),1/(2*p^2)(Sqrt[4*p^2-1])];
Graphics[{PointSize[.0005],Point[selxy],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}],t,{x0,y0}],{k,2,n,2}],Table[Rotate[Circle[{x0,y0},r0,{tc+k*2Pi/n,-tc+(k+2)2Pi/n}],t,{x0,y0}],{k,2,n,2}]},Axes->True,PlotRange->{{-.6,.6},{-.4,.8}}]],{{x0,0.02},-.5,.5,.001},{{y0,0.206},-.5,.5,.001},{{t,0.134},0,2Pi,.001},ControlPlacement->Right]