(*  The following are examples of curvature analysis  on  individual tubercles in seeds of Silene species*)






	(*******************************  AJ300  *)

Needs["Splines`"]

yunnanensis[1]={{5.8922,14.8708},{7.0146,15.1514},{8.1369,15.7126},{9.2592,16.5543},{10.3815,17.1155},{11.5039,17.6767},{12.6262,18.5184},{13.7485,19.0796},{14.8708,20.2019},{15.9932,20.7631},{17.1155,21.8854},{18.2378,22.4466},{19.3602,23.5689},{20.4825,24.6912},{21.6048,25.8135},{22.7271,26.9359},{23.8495,28.6194},{24.9718,30.5834},{26.0941,32.5475},{27.2165,34.7922},{28.3388,36.7562},{29.4611,39.2815},{30.5834,40.4038},{31.7058,40.9650},{32.8281,41.2456},{33.9504,41.2456},{35.0728,41.2456},{36.1951,40.9650},{37.3174,40.4038},{38.4397,39.5621},{39.5621,38.7203},{40.6844,37.5980},{41.8067,36.1951},{42.9291,34.7922},{44.0514,33.1087},{45.1737,31.4252},{46.2960,29.7417},{47.4184,27.7776},{48.5407,25.8135},{49.6630,24.6912},{50.7853,22.7271},{51.9077,21.3242},{53.0300,19.6407},{54.1523,18.5184},{55.2747,17.6767},{56.3970,16.8349},{57.5193,15.7126},{58.6416,15.1514},{59.7640,14.5903},{60.8863,14.3097}}
yunnanensis[2]={{7.0146,12.0650},{8.1369,12.6262},{9.2592,12.9068},{10.3815,13.4679},{11.5039,14.3097},{12.6262,15.1514},{13.7485,15.7126},{14.8708,16.8349},{15.9932,17.9572},{17.1155,19.0796},{18.2378,19.6407},{19.3602,20.7631},{20.4825,21.8854},{21.6048,23.0077},{22.7271,23.5689},{23.8495,24.6912},{24.9718,25.8135},{26.0941,27.2165},{27.2165,28.3388},{28.3388,29.1805},{29.4611,30.8640},{30.5834,31.9864},{31.7058,33.6698},{32.8281,35.0728},{33.9504,36.7562},{35.0728,38.1592},{36.1951,40.1232},{37.3174,42.3679},{38.4397,44.6125},{39.5621,45.4543},{40.6844,46.5766},{41.8067,47.6989},{42.9291,48.5407},{44.0514,49.1019},{45.1737,49.6630},{46.2960,49.6630},{47.4184,49.6630},{48.5407,49.6630},{49.6630,49.6630},{50.7853,49.1019},{51.9077,48.5407},{53.0300,47.4184},{54.1523,46.5766},{55.2747,44.8931},{56.3970,43.4902},{57.5193,40.9650},{58.6416,39.0009},{59.7640,37.3174},{60.8863,35.3533},{62.0086,32.8281},{63.1310,31.4252},{64.2533,30.3029},{65.3756,28.8999},{66.4979,28.0582},{67.6203,26.9359},{68.7426,25.8135},{69.8649,24.6912},{70.9873,23.2883},{72.1096,22.1660},{73.2319,21.3242},{74.3542,20.2019},{75.4766,19.6407},{76.5989,18.7990},{77.7212,18.5184},{78.8436,17.9572},{79.9659,17.6767}}
yunnanensis[3]={{3.6476,0.2806},{4.7699,1.1223},{5.8922,2.2447},{7.0146,3.0864},{8.1369,3.9281},{9.2592,5.0505},{10.3815,6.1728},{11.5039,7.2951},{12.6262,7.8563},{13.7485,9.2592},{14.8708,10.1010},{15.9932,11.2233},{17.1155,11.7844},{18.2378,12.6262},{19.3602,13.7485},{20.4825,14.3097},{21.6048,15.1514},{22.7271,16.2738},{23.8495,16.8349},{24.9718,17.6767},{26.0941,18.5184},{27.2165,19.6407},{28.3388,20.7631},{29.4611,21.6048},{30.5834,22.7271},{31.7058,24.1301},{32.8281,25.5330},{33.9504,27.7776},{35.0728,30.3029},{36.1951,33.1087},{37.3174,35.6339},{38.4397,37.5980},{39.5621,38.7203},{40.6844,39.8426},{41.8067,40.4038},{42.9291,40.6844},{44.0514,41.2456},{45.1737,41.2456},{46.2960,41.2456},{47.4184,41.2456},{48.5407,40.6844},{49.6630,39.8426},{50.7853,39.2815},{51.9077,37.8786},{53.0300,36.4757},{54.1523,34.5116},{55.2747,31.9864},{56.3970,29.1805},{57.5193,26.6553},{58.6416,24.4106},{59.7640,23.0077},{60.8863,21.3242},{62.0086,20.2019},{63.1310,19.0796},{64.2533,18.2378},{65.3756,17.3961},{66.4979,16.8349},{67.6203,16.5543}}
yunnanensis[4]={{0.2806,8.6980},{1.4029,9.2592},{2.5252,9.5398},{3.6476,10.1010},{4.7699,10.9427},{5.8922,11.7844},{7.0146,12.6262},{8.1369,13.1874},{9.2592,14.0291},{10.3815,15.1514},{11.5039,15.7126},{12.6262,16.8349},{13.7485,17.3961},{14.8708,18.5184},{15.9932,19.6407},{17.1155,20.7631},{18.2378,22.1660},{19.3602,23.5689},{20.4825,24.9718},{21.6048,26.6553},{22.7271,28.6194},{23.8495,30.5834},{24.9718,32.8281},{26.0941,35.0728},{27.2165,37.3174},{28.3388,39.5621},{29.4611,40.9650},{30.5834,41.5261},{31.7058,42.6485},{32.8281,43.2096},{33.9504,43.7708},{35.0728,44.0514},{36.1951,44.0514},{37.3174,44.0514},{38.4397,43.4902},{39.5621,42.9291},{40.6844,42.3679},{41.8067,41.2456},{42.9291,39.5621},{44.0514,37.8786},{45.1737,35.6339},{46.2960,33.1087},{47.4184,30.5834},{48.5407,28.6194},{49.6630,26.6553},{50.7853,24.6912},{51.9077,23.0077},{53.0300,21.3242},{54.1523,19.9213},{55.2747,18.2378},{56.3970,17.1155},{57.5193,16.2738},{58.6416,15.4320},{59.7640,15.1514},{60.8863,14.5903},{62.0086,14.3097}}
yunnanensis[5]={{7.0146,17.1155},{8.1369,17.3961},{9.2592,17.9572},{10.3815,19.0796},{11.5039,20.2019},{12.6262,21.6048},{13.7485,23.2883},{14.8708,24.9718},{15.9932,27.2165},{17.1155,29.7417},{18.2378,32.8281},{19.3602,36.4757},{20.4825,39.2815},{21.6048,40.9650},{22.7271,42.0873},{23.8495,42.3679},{24.9718,42.9291},{26.0941,42.9291},{27.2165,42.9291},{28.3388,42.3679},{29.4611,41.5261},{30.5834,39.8426},{31.7058,36.1951},{32.8281,32.5475},{33.9504,29.7417},{35.0728,26.6553},{36.1951,23.8495},{37.3174,21.3242},{38.4397,19.6407},{39.5621,18.5184},{40.6844,17.6767},{41.8067,17.1155},{42.9291,16.2738},{44.0514,15.7126},{45.1737,15.4320},{46.2960,14.3097},{47.4184,12.9068},{48.5407,12.3456}}
yunnanensis[6]={{7.0146,8.1369},{8.1369,8.6980},{9.2592,8.9786},{10.3815,9.8204},{11.5039,10.3815},{12.6262,11.2233},{13.7485,11.7844},{14.8708,12.9068},{15.9932,14.0291},{17.1155,15.4320},{18.2378,17.1155},{19.3602,18.7990},{20.4825,20.4825},{21.6048,22.4466},{22.7271,24.1301},{23.8495,26.9359},{24.9718,32.2669},{26.0941,39.0009},{27.2165,42.6485},{28.3388,44.8931},{29.4611,46.0155},{30.5834,46.5766},{31.7058,47.1378},{32.8281,47.4184},{33.9504,47.4184},{35.0728,47.4184},{36.1951,47.4184},{37.3174,46.8572},{38.4397,46.0155},{39.5621,44.8931},{40.6844,43.7708},{41.8067,40.1232},{42.9291,33.6698},{44.0514,29.1805},{45.1737,23.2883},{46.2960,21.3242},{47.4184,20.2019}}

(*Poner aqui la familia*)
var=yunnanensis;
(* el número de curvas estudiadas/semillas de la familia*)
m=6;
For[j=1,j≤m,j++,
puntos=var[j];
(* pru es el spline cubico
d es el vector que contiene la información de pru 
n es el número de trozos del spline, es decir, el número de puntos menos 1
*)
pru=SplineFit[puntos,Cubic];
d=InputForm[pru];
n=IntegerPart[d[[1,2,2]]];
spli=ParametricPlot[pru[t],{t,0,n},PlotRange->All, DisplayFunction->Identity];
(* Con los coeficientes se calculan los polinomios cúbicos que dan el spline. Habrá n polinomios para la primera coordenada (X) y otros tantos para la segunda
*)
For[i=0,i<n,i++,
SX[i,t_]=d[[1,4,i+1,1,1]]+d[[1,4,i+1,1,2]]*(t-i)+d[[1,4,i+1,1,3]]*(t-i)^2+d[[1,4,i+1,1,4]]*(t-i)^3];
For[i=0,i<n,i++,SY[i,t_]=d[[1,4,i+1,2,1]]+d[[1,4,i+1,2,2]]*(t-i)+d[[1,4,i+1,2,3]]*(t-i)^2+d[[1,4,i+1,2,4]]*(t-i)^3];
(* Derivada primera y segunda de cada trozo del spline
*)
For[i=0,i<n,i++,DX[i,t_]=d[[1,4,i+1,1,2]]+2*d[[1,4,i+1,1,3]]*(t-i)+3*d[[1,4,i+1,1,4]]*(t-i)^2];
For[i=0,i<n,i++,DY[i,t_]=d[[1,4,i+1,2,2]]+2*d[[1,4,i+1,2,3]]*(t-i)+3*d[[1,4,i+1,2,4]]*(t-i)^2];
For[i=0,i<n,i++,DDX[i,t_]=2*d[[1,4,i+1,1,3]]+6*d[[1,4,i+1,1,4]]*(t-i)];
For[i=0,i<n,i++,DDY[i,t_]=2*d[[1,4,i+1,2,3]]+6*d[[1,4,i+1,2,4]]*(t-i)];
(*Cálculo de la curvatura a trozos
*)
For[i=0,i<n,i++,cur[i,t_]=(DX[i,t]DDY[i,t]-DDX[i,t]DY[i,t])/((DX[i,t])^2+(DY[i,t])^2)^(3/2)];
(* Dibujos de los trozos de curvatura 
*)
For[i=0,i<n,i++,dib[i]=Plot[cur[i,t],{t,i,i+1},DisplayFunction->Identity]];
pp=dib[0];
For[i=1,i<n,i++,pp=Show[pp,dib[i],PlotRange->All]];
(* bez es el spline de Bezier
b es el vector que contiene la información de bez
nb es el número de trozos del spline, es decir, el número de puntos menos 1 (Como nb=n, no se define nuevo)
*)
bez=SplineFit[puntos,Bezier];
b=InputForm[bez];
(*Calculo de las funciones de Bezier para la coordenada X y la Y; la curva es la imagen del intervalo[0,1]*)
For[i=1,i≤n+1,i++,BX[t_]=bez[[4,1,1]]/.#1->t];
For[i=1,i≤n+1,i++,BY[t_]=bez[[4,1,2]]/.#1->t];
dibez[j]=ParametricPlot[{BX[t],BY[t]},{t,0,1},PlotRange->All,DisplayFunction->Identity];
curbe[j,t_]=(BX'[t]BY''[t]-BX''[t]BY'[t])/((BX'[t])^2+(BY'[t])^2)^(3/2);
cur[j]=Abs[NMinValue[{curbe[j,t],0.30<t<0.7},t]];
vmed[j]=Abs[NIntegrate[curbe[j,t],{t,0.3,0.7}]/0.4];
dibcurbe[j]=Plot[{curbe[j,t],-vmed[j]},{t,0,1},PlotRange->{-0.4,0.1},DisplayFunction->Identity,PlotStyle->{RGBColor[0,1,0],{Dotted,Dashed}}];
total=GraphicsGrid[{
(*{spli,pp},*)
{dibez[j],dibcurbe[j]}},Spacings->0.3];
Print[var,j, ":", "Máximo valor de la curvatura (en valor absoluto)",cur[j]];
Print[total]
]
Print["Curvaturas de cada ", var ,":  ", Table[cur[j],{j,1,6}]];
Print["Valores medios de curvatura de cada ", var ," en [0.3,0.7]:  ", Table[vmed[j],{j,1,6}]];

Show[Table[dibcurbe[j],{j,1,6}]]

Table[NMinValue[{curbe[j,t],0.25<t<0.75},t],{j,1,6}]

dibcurbe[j],{j,1,6}

Table[dibcurbe[j],{j,1,6}]

curbe[j,t_]
dibcurbe[j]