(*This Mathematica code was written by John Zweck (zweck@utdallas.edu) and generates a mathematical curve/surface that can be printed using a 3D printer. This model can (probably) be printed by Shapeways.*)

(*If you evaluate this Mathematica notebook,the curve/surface will be rendered using computer graphics.*)

(*If you uncomment the export line at the end of the code and set the directory path appropriately for your system you will obtain an STL file that can be sent to a 3D printer.*)

(*The STL file can also be viewed using the open source software Meshlab.*)

(*You may modify this code in any way you like.*)

(*This project was supported by the National Science Foundation through a Supplement to the UTeach Dallas Robert Noyce Scholarship Program (DUE 1544159).*)

(* This was one of the first surfaces I made and as a consequence the code is probably more complicated than it needs to be. The resultiung model is however of very good quality. *)

upperSurface[u_,v_]:={u*Cosh[v],u*Sinh[v],u^2};
lowerSurface[x_,y_]:={x*Sinh[y],x*Cosh[y],-x^2};
negX[t_]:={t,-t,0};
posX[t_]:={t,t,0};
intersectionCylin[r_,t_]={r*(Cos[t]),r*(Sin[t]),r^2*((Cos[t])^2-(Sin[t])^2)};

scale = 16;
radius = 2;

RadiusCylinder=2;

numPoints=24;
gridstepsLevel=8;

umin=-2;
umax=2;
vmin=-2;
vmax=2;

levelSteps={-4^(1/2),-3^(1/2),-2^(1/2),-1,1,2^(1/2),3^(1/2),4^(1/2)};

(*For vertices where sin curve on cylinder intersects level curves we need times on intersectionCylin where these intersection points occur*) 

CylinderIntersectionTimes1=Table[0.5*ArcCos[(k/RadiusCylinder)^2], {k, levelSteps}];
CylinderIntersectionTimes2=Table[0.5*ArcCos[-(k/RadiusCylinder)^2],{k, levelSteps}];
CylinderIntersectionTimes3=Table[0.5*(Pi+ArcCos[(k/RadiusCylinder)^2]),{k, levelSteps}];
CylinderIntersectionTimes4=Table[0.5*(Pi+ArcCos[-(k/RadiusCylinder)^2]),{k, levelSteps}];

curvesUpperYplane=Table[scale*upperSurface[u,i],{i,0,0,1}]; (*One upward parabola at v=0*) 
curvesLowerYplane=Table[scale*lowerSurface[u,i],{i,0,0,1}];  (*One downward parabola at v=0*)

curvesUpperXplane=Table[scale*upperSurface[j,v],{j,levelSteps}]; (*Hyperbolae above xy-plane*)
curvesLowerXplane=Table[scale*lowerSurface[j,v],{j,levelSteps}];  (*Hyperbolae below xy-plane*)


tubesUpperYplane=ParametricPlot3D[curvesUpperYplane,{u,umin,umax},RegionFunction->Function[{x,y,z},(x^2+y^2)^0.5<scale*RadiusCylinder],
PlotStyle->Tube[radius,PlotPoints->numPoints],PlotRange->All];

tubesLowerYplane=ParametricPlot3D[curvesLowerYplane,{u,umin,umax},
RegionFunction->Function[{x,y,z},(x^2+y^2)^0.5<scale*RadiusCylinder],
PlotStyle->Tube[radius,PlotPoints->numPoints],PlotRange->All];


tubesUpperXplane=ParametricPlot3D[curvesUpperXplane,{v,vmin,vmax},RegionFunction->Function[{x,y,z},(x^2+y^2)^0.5<scale*RadiusCylinder],
PlotStyle->Tube[radius,PlotPoints->numPoints],PlotRange->All];

tubesLowerXplane=ParametricPlot3D[curvesLowerXplane,{v,vmin,vmax},
RegionFunction->Function[{x,y,z},(x^2+y^2)^0.5<scale*RadiusCylinder],
PlotStyle->Tube[radius,PlotPoints->numPoints],PlotRange->All];

(*Straight lines in xy-plane*)

tubesNegX=ParametricPlot3D[scale*negX[i],{i,-umax,umax},
RegionFunction->Function[{x,y,z},(x^2+y^2)^0.5<scale*RadiusCylinder],
PlotStyle->Tube[radius,PlotPoints->numPoints],PlotRange->All];     

tubesPosX=ParametricPlot3D[scale*posX[i],{i,-umax,umax},
RegionFunction->Function[{x,y,z},(x^2+y^2)^0.5<scale*RadiusCylinder],
PlotStyle->Tube[radius,PlotPoints->numPoints],PlotRange->All];

tubesIntersect=ParametricPlot3D[scale*intersectionCylin[RadiusCylinder,t],{t,0,2*Pi},PlotStyle->Tube[radius,PlotPoints->numPoints],PlotRange->All];

verticesUpperSurface = Graphics3D[Table[Sphere[scale*upperSurface[i,j],radius],{i, levelSteps},{j, 0,0,1}],PlotPoints->numPoints];

verticesLowerSurface = Graphics3D[Table[Sphere[scale*lowerSurface[i,j],radius],{i, levelSteps},{j, 0,0,1}],PlotPoints->numPoints];

verticesPosX = Graphics3D[Table[Sphere[scale*posX[i],radius],{i, -2^(1/2), 2^(1/2), 2*2^(1/2)}],PlotPoints->numPoints];

verticesNegX = Graphics3D[Table[Sphere[scale*negX[i],radius],{i, -2^(1/2), 2^(1/2), 2*2^(1/2)}],PlotPoints->numPoints];

verticesIntersect1 = Graphics3D[Table[Sphere[scale*intersectionCylin[RadiusCylinder,t],radius],{t,CylinderIntersectionTimes1}],PlotPoints->numPoints];

verticesIntersect2 = Graphics3D[Table[Sphere[scale*intersectionCylin[RadiusCylinder,t],radius],{t,CylinderIntersectionTimes2}],PlotPoints->numPoints];


verticesIntersect3 = Graphics3D[Table[Sphere[scale*intersectionCylin[RadiusCylinder,t],radius],{t,CylinderIntersectionTimes3}],PlotPoints->numPoints];

verticesIntersect4 = Graphics3D[Table[Sphere[scale*intersectionCylin[RadiusCylinder,t],radius],{t,CylinderIntersectionTimes4}],PlotPoints->numPoints];

verticesIntersect5 = Graphics3D[Table[Sphere[scale*intersectionCylin[RadiusCylinder,t],radius],{t,-CylinderIntersectionTimes1}],PlotPoints->numPoints];

verticesIntersect6 = Graphics3D[Table[Sphere[scale*intersectionCylin[RadiusCylinder,t],radius],{t,-CylinderIntersectionTimes2}],PlotPoints->numPoints];


verticesIntersect7 = Graphics3D[Table[Sphere[scale*intersectionCylin[RadiusCylinder,t],radius],{t,-CylinderIntersectionTimes3}],PlotPoints->numPoints];

verticesIntersect8 = Graphics3D[Table[Sphere[scale*intersectionCylin[RadiusCylinder,t],radius],{t,-CylinderIntersectionTimes4}],PlotPoints->numPoints];

origin = Graphics3D[Table[Sphere[scale*posX[t],radius],{t, 0,0}],PlotPoints->numPoints];

output=Show[tubesUpperXplane,tubesUpperYplane,tubesLowerXplane,tubesLowerYplane,verticesUpperSurface,verticesLowerSurface, tubesNegX,tubesPosX,  tubesIntersect, verticesPosX, verticesNegX, verticesIntersect1,verticesIntersect2,verticesIntersect3,verticesIntersect4,verticesIntersect5,verticesIntersect6,verticesIntersect7,verticesIntersect8,origin]

(*Export["SaddleSurfaceLevelCurves.stl",output]*)