Skip to content

Commit 7ce742a

Browse files
pdpollypdpolly
authored
Add files via upload
1 parent 401d276 commit 7ce742a

2 files changed

Lines changed: 200 additions & 0 deletions

File tree

‎Guide to Mathematica Snails 1.1.pdf‎

388 KB
Binary file not shown.

‎Snails for Mathematica 1.1.m‎

Lines changed: 200 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,200 @@
1+
(* ::Package:: *)
2+
3+
(* ::Title:: *)
4+
(*Snail functions for Mathematica*)
5+
6+
7+
(* ::Text:: *)
8+
(*Version 1.1*)
9+
(*15 October 2025*)
10+
(**)
11+
(*This notebook contains functions for simulating and processing virtual gastropod (and other mollusc) shells. *)
12+
(**)
13+
(*P. David Polly*)
14+
(*pdpolly@pollylab.org*)
15+
(*https://pollylab.org/*)
16+
17+
18+
(* This function prints the version number for the installed verison of this package. *)
19+
VersionNumber="1.1";
20+
VersionDate="15 October 2025";
21+
SnailsVersion[]:=Print["Snails for Mathematica "<>VersionNumber<>"\n(c) P. David Polly, "<>VersionDate<>"\n"];
22+
SnailsVersion[]
23+
24+
25+
(* This function generates a three-dimensional rendered shell using the functions
26+
developed by David Raup (1966).
27+
28+
Usage: RaupCoil3D[shape, W, T, D, TurnNum]
29+
where shape is a matrix of 3D semilandmarks that circumscribe the aperture
30+
W is Raup's whorl expansion rate parameter
31+
T is Raup's vertical translation rate parameter
32+
D is Raup's parameter for distance of aperture from columella
33+
TurnNum is the desired number of whorls
34+
35+
Created 16 March 2016.
36+
Updated 5 February 2023 with adjustment to aperture position relative to coiling.
37+
*)
38+
39+
RaupCoil3D[S_,W_,T_,Di_, turns_,CoordsOnly_:False]:=Quiet[Block[{r,z,x,CoiledCoords,\[Theta],Sprime,S3D,yaspect,rc},
40+
41+
Sprime=#-Mean[S]&/@S;
42+
Sprime=(Sprime/(Max[Sprime[[1;;,1]]]-Min[Sprime[[1;;,1]]]));
43+
Sprime[[1;;,1]]=Sprime[[1;;,1]]-Min[Sprime[[1;;,1]]]+Di;
44+
yaspect=Max[Sprime[[1;;,3]]]-Min[Sprime[[1;;,3]]];
45+
rc=Mean[Sprime[[1;;,1]]];
46+
S3D=Sprime;
47+
CoiledCoords=Table[Table[
48+
r=S3D[[x,1]]*(W^(\[Theta]/(2Pi)));
49+
z=(S3D[[x,3]]*(W^(\[Theta]/(2Pi))))+(T*(W^(\[Theta]/(2Pi))-0.5));
50+
{r*Cos[\[Theta]],r*Sin[\[Theta]]+S3D[[x,2]]*(W^(\[Theta]/(2Pi))),z}
51+
,{x,Length[S3D]}]
52+
,{\[Theta],0,turns*2Pi,10Degree}];
53+
If[CoordsOnly==True,Return[CoiledCoords],
54+
Return[
55+
Graphics3D[
56+
Table[{
57+
Table[{RGBColor[0.997681, 0.980789, 0.842313],EdgeForm[None],Opacity[1],Polygon[{CoiledCoords[[x,y]],CoiledCoords[[x-1,y]],CoiledCoords[[x-1,y+1]],CoiledCoords[[x,y+1]]}]},{y,Length[CoiledCoords[[x]]]-1}],
58+
{GrayLevel[0.5],Opacity[0],AbsolutePointSize[1.75],Point[CoiledCoords[[x]]]}},
59+
{x,2,Length[CoiledCoords]}],Boxed->False]]];
60+
]];
61+
62+
63+
64+
(* This function generates a three-dimensional rendered shell using the functions
65+
developed by David Raup (1966) as a three-dimensional mesh suitable for printing
66+
with a 3D preinter. It is nearly idential ato RaupCoil3D execpt it creates a shell
67+
with two layers that is capped at apex and aperture.
68+
69+
Usage: RaupCoil3DForPrint[shape, W, T, D, TurnNum]
70+
where shape is a matrix of 3D semilandmarks that circumscribe the aperture
71+
W is Raup's whorl expansion rate parameter
72+
T is Raup's vertical translation rate parameter
73+
D is Raup's parameter for distance of aperture from columella
74+
TurnNum is the desired number of whorls
75+
76+
Created 16 March 2016.
77+
*)
78+
79+
RaupCoil3DForPrint[S_,W_,T_,Di_, turns_]:=Quiet[Block[{r,z,x,CoiledCoords,\[Theta],Sprime,S3D,yaspect,rc},
80+
81+
Sprime=#-Mean[S]&/@S;
82+
Sprime=(Sprime/(Max[Sprime[[1;;,1]]]-Min[Sprime[[1;;,1]]]));
83+
Sprime[[1;;,1]]=Sprime[[1;;,1]]-Min[Sprime[[1;;,1]]]+Di;
84+
yaspect=Max[Sprime[[1;;,3]]]-Min[Sprime[[1;;,3]]];
85+
rc=Mean[Sprime[[1;;,1]]];
86+
S3D=Sprime;
87+
CoiledCoords=Table[Table[
88+
r=S3D[[x,1]]*(W^(\[Theta]/(2Pi)));
89+
z=(S3D[[x,3]]*(W^(\[Theta]/(2Pi))))+rc*T*yaspect*((W^(\[Theta]/(2Pi)))-1);
90+
{r*Cos[\[Theta]],r*Sin[\[Theta]]+S3D[[x,2]]*(W^(\[Theta]/(2Pi))),z}
91+
,{x,Length[S3D]}]
92+
,{\[Theta],0,turns*2Pi,10Degree}];
93+
Return[CoiledCoords];
94+
]];
95+
96+
97+
(* This function tilts the aperture shape at a specified number of degrees.
98+
99+
Usage: TiltAperture[shape, angle]
100+
where shape is a matrix of 3D semilandmarks that circumscribe the aperture
101+
and angle is an angle from vertical given in degrees
102+
103+
Created 16 March 2016.
104+
*)
105+
106+
TiltAperture[aperture_,angle_]:=Module[{tiltedaperture},
107+
tiltedaperture=aperture . RotationMatrix[angle Degree,{1,0,0}];
108+
tiltedaperture[[1;;,2]]=tiltedaperture[[1;;,2]]+Mean[tiltedaperture[[1;;,2]]];
109+
Return[tiltedaperture];
110+
];
111+
112+
113+
(* This function generates three-dimensional elliptical Fourier coefficiences from
114+
a closed curve (like an aperture shape). It returns Fourier cofficients for a
115+
single shape.
116+
117+
Usage: EllipticalFourierCoefficients[points]
118+
where points is a matrix of 3D semilandmarks that circumscribe the aperture
119+
120+
Created 16 March 2016.
121+
*)
122+
123+
EllipticalFourierCoefficients[standardizedoutline_]:=Module[{deltaXY,deltaT,TotalT,FourierCoeffs},
124+
deltaXY=Table[standardizedoutline[[k+1]]-standardizedoutline[[k]],{k,Length[standardizedoutline]-1}];
125+
deltaT=Sqrt[Plus@@(#^2)]&/@deltaXY;
126+
TotalT=Plus@@deltaT;
127+
FourierCoeffs=Table[{TotalT/(2*(n^2)*(Pi^2))*(Plus@@Table[deltaXY[[k,1]]/deltaT[[k]]*(Cos[(2Pi*n*Plus@@deltaT[[1;;k+1]])/TotalT]-Cos[(2Pi*n*Plus@@deltaT[[1;;k]])/TotalT]),{k,Length[deltaXY]-1}]),TotalT/(2*(n^2)*(Pi^2))*(Plus@@Table[deltaXY[[k,1]]/deltaT[[k]]*(Sin[(2Pi*n*Plus@@deltaT[[1;;k+1]])/TotalT]-Sin[(2Pi*n*Plus@@deltaT[[1;;k]])/TotalT]),{k,Length[deltaXY]-1}]),TotalT/(2*(n^2)*(Pi^2))*(Plus@@Table[deltaXY[[k,2]]/deltaT[[k]]*(Cos[(2Pi*n*Plus@@deltaT[[1;;k+1]])/TotalT]-Cos[(2Pi*n*Plus@@deltaT[[1;;k]])/TotalT]),{k,Length[deltaXY]-1}]),TotalT/(2*(n^2)*(Pi^2))*(Plus@@Table[deltaXY[[k,2]]/deltaT[[k]]*(Sin[(2Pi*n*Plus@@deltaT[[1;;k+1]])/TotalT]-Sin[(2Pi*n*Plus@@deltaT[[1;;k]])/TotalT]),{k,Length[deltaXY]-1}])},{n,Length[deltaXY]/2}];
128+
Return[{FourierCoeffs,deltaT,TotalT}];
129+
];
130+
131+
132+
(* This function generates a shape in X, Y coordinates from elliptical Fourier coefficients.
133+
134+
Usage: HarmonicReconstruction[coefficients, deltaT, T, n]
135+
where coefficients is a matrix of Fourier coefficients
136+
deltaT is the offset
137+
T is the number
138+
n is
139+
140+
Created 16 March 2016.
141+
*)
142+
HarmonicReconstruction[coeffs_,deltaT_,T_,n_]:=Module[{xy},
143+
xy=Table[Flatten[{Plus@@Table[coeffs[[j,1]]*Cos[(2Pi*j*Plus@@deltaT[[1;;i]])/T]+coeffs[[j,2]]*Sin[(2Pi*j*Plus@@deltaT[[1;;i]])/T],{j,n}],Plus@@Table[coeffs[[j,3]]*Cos[(2Pi*j*Plus@@deltaT[[1;;i]])/T]+coeffs[[j,4]]*Sin[(2Pi*j*Plus@@deltaT[[1;;i]])/T],{j,n}]}],{i,Length[deltaT]}];
144+
Return[xy]
145+
];
146+
147+
148+
(* This function saves a 3D shell suitable for printing as an OBJ mesh file.
149+
150+
Usage: SaveShell[shellobject, parameters, title, path]
151+
where shellobject is a shell generated with the RaupCoil3DForPrint function
152+
parameters are the shell coiling parameters used to generate the shell(which are included as metadata in the OBJ file)
153+
title is a title for the shell object (included as metadata in the OBJ file)
154+
path is the path to the location where the file will be saved.
155+
156+
Created 16 March 2016.
157+
*)
158+
159+
SaveShell[myShell_,parameters_,title_,path_]:=Module[{s,vertexnumbers,verts,outer,inner,outerfaces,innerfaces,endfaces,innerprotoconchcenter,innerprotofaces, outerprotofaces,outerprotoconchcenter,W,T,shp,turn,Di},
160+
{W,T,Di,turn}=parameters;
161+
162+
(* calculating the faces *)
163+
164+
vertexnumbers=Partition[Table[x,{x,Dimensions[myShell][[1]]*Dimensions[myShell][[2]]}],Dimensions[myShell][[2]]];
165+
verts=StringRiffle[Prepend[#,"v"]]&/@Flatten[myShell,1];
166+
167+
outer={1,Dimensions[vertexnumbers][[2]]/2};
168+
inner={(Dimensions[vertexnumbers][[2]]/2)+1,Dimensions[vertexnumbers][[2]]};
169+
170+
outerfaces=Flatten[{Table[Table[{StringRiffle[Prepend[
171+
{vertexnumbers[[i+1,j]],vertexnumbers[[i,j+1]],vertexnumbers[[i,j]]},"f"]],StringRiffle[Prepend[{vertexnumbers[[i+1,j+1]],vertexnumbers[[i,j+1]],vertexnumbers[[i+1,j]]},"f"]]},{j,outer[[1]],outer[[2]]-1}],{i,Dimensions[vertexnumbers][[1]]-1}],Table[StringRiffle[Prepend[{vertexnumbers[[i,outer[[1]]]],vertexnumbers[[i,outer[[2]]]],vertexnumbers[[i+1,outer[[2]]]]},"f"]],{i,Dimensions[vertexnumbers][[1]]-1}],
172+
Table[StringRiffle[Prepend[{vertexnumbers[[i,outer[[1]]]],vertexnumbers[[i+1,outer[[2]]]],vertexnumbers[[i+1,outer[[1]]]]},"f"]],{i,Dimensions[vertexnumbers][[1]]-1}]}];
173+
174+
innerfaces=Flatten[{Table[Table[{StringRiffle[Prepend[
175+
{vertexnumbers[[i,j]],vertexnumbers[[i,j+1]],vertexnumbers[[i+1,j]]},"f"]],StringRiffle[Prepend[{vertexnumbers[[i+1,j]],vertexnumbers[[i,j+1]],vertexnumbers[[i+1,j+1]]},"f"]]},{j,inner[[1]],inner[[2]]-1}],{i,Dimensions[vertexnumbers][[1]]-1}],Table[StringRiffle[Prepend[{vertexnumbers[[i+1,inner[[2]]]],vertexnumbers[[i,inner[[2]]]],vertexnumbers[[i,inner[[1]]]]},"f"]],{i,Dimensions[vertexnumbers][[1]]-1}],
176+
Table[StringRiffle[Prepend[{vertexnumbers[[i+1,inner[[1]]]],vertexnumbers[[i+1,inner[[2]]]],vertexnumbers[[i,inner[[1]]]]},"f"]],{i,Dimensions[vertexnumbers][[1]]-1}]}];
177+
178+
endfaces=Flatten[{Table[{StringRiffle[{"f",vertexnumbers[[-1,i+1]],vertexnumbers[[-1,i]],vertexnumbers[[-1,i+inner[[1]]-1]]}],StringRiffle[{"f",vertexnumbers[[-1,i+inner[[1]]]],vertexnumbers[[-1,i+1]],vertexnumbers[[-1,i+inner[[1]]-1]]}]},{i,outer[[2]]-1}],
179+
StringRiffle[{"f",vertexnumbers[[-1,outer[[1]]]],vertexnumbers[[-1,outer[[2]]]],vertexnumbers[[-1,inner[[2]]]]}],
180+
StringRiffle[{"f",vertexnumbers[[-1,outer[[1]]]],vertexnumbers[[-1,inner[[2]]]],vertexnumbers[[-1,inner[[1]]]]}]}];
181+
innerprotoconchcenter=Mean[myShell[[1,inner[[1]];;inner[[2]]]]];
182+
verts=Append[verts,StringRiffle[Prepend[innerprotoconchcenter,"v"]]];
183+
innerprotofaces=Append[Table[StringRiffle[{"f",Length[verts],i+1,i}],{i,inner[[1]],inner[[2]]-1,1}],StringRiffle[{"f",Length[verts],inner[[1]],inner[[2]]}]];
184+
outerprotoconchcenter=innerprotoconchcenter+{0,-1*Norm[innerprotoconchcenter-myShell[[1,1]]]/2,0};
185+
verts=Append[verts,StringRiffle[Prepend[outerprotoconchcenter,"v"]]];
186+
outerprotofaces=Append[Table[StringRiffle[{"f",Length[verts],i,i+1}],{i,outer[[1]],outer[[2]]-1,1}],StringRiffle[{"f",Length[verts],outer[[2]],outer[[1]]}]];
187+
(* Writing the data *)
188+
s=OpenWrite[path<>title<>".obj"];
189+
WriteLine[s,"# Shell Mesh created by SaveShell[], P. David Polly "<>DateString[]];
190+
WriteLine[s, "# Snails for Mathematica, Version "<>VersionNumber];
191+
WriteLine[s,"# Parameters: W="<>ToString[W]<>" T="<>ToString[T]<>" D="<>ToString[Di]<>" turns="<>ToString[turn]];
192+
WriteLine[s,#]&/@verts;
193+
WriteLine[s,#]&/@outerfaces;
194+
WriteLine[s,#]&/@innerfaces;
195+
WriteLine[s,#]&/@endfaces;
196+
WriteLine[s,#]&/@innerprotofaces;
197+
WriteLine[s,#]&/@outerprotofaces;
198+
WriteLine[s,"# End of File \n"];
199+
Close[s]
200+
];

0 commit comments

Comments
 (0)