(** isoglass.m -- function for an iso glass rendering.
 * @author Eric Laroche
 * @version @(#)$Id: isoglass.m,v 1.4 1998/01/17 14:42:38 laroche Exp $
 **)

(* package isoglass` *)
BeginPackage["isoglass`"]

(* exported variables and functions *)
isoglasslength::usage =
	"isoglasslength[]: iso glass length"
isoglassradius::usage =
	"isoglassradius[]: iso glass maximal radius"
isoglassprofile::usage =
	"isoglassprofile[z]: iso glass radius profile function"
isoglassplot::usage =
	"isoglassplot[]: plot of the iso glass"

(* internal part *)
Begin["`Private`"]

(* the parameters for the glass rendering functions are:
 * footradius -- maximum foot radius
 * stemradius -- minimum stem radius
 * stemlength -- stem length
 * cupradius -- maximum cup radius
 * lowercuplength -- length below cupradius
 * openingradius -- opening radius
 * uppercuplength -- length above cupradius
 *)
(* glass profile function *)
glassprofile[z_,
	{footradius_, stemradius_, stemlength_,
	cupradius_, lowercuplength_, openingradius_, uppercuplength_}] :=

	Module[{stem, lowercup, uppercup},

	(* stem: negative exponential for now, approximate stemradius *)
	stem[zz_] = (footradius-stemradius) Exp[
		-25 zz/stemlength]+stemradius;

	(* lower part: quarter of an ellipsoid: horizontal at lowest,
		vertical at lowercuplength *)
	lowercup[zz_] = r /.
		(* solve to r[zz], get the positive solution *)
		Solve[(zzz/lowercuplength)^2+(r/cupradius)^2 == 1,
			r][[-1, 1]] /.
		(* inverse zz *)
		zzz->lowercuplength-zz;
	(* assertions *)
	(* (D[lowercup[zz], zz] /. zz->0) == ComplexInfinity; *)

	(* upper part: paraboloid, vertical at 0 *)
	uppercup[zz_] = p zz^2 + cupradius /.
		(* solve the parable parameter *)
		Solve[(p zz^2 + cupradius == r) /.
			{zz->uppercuplength, r->openingradius},
			p][[1, 1]];
	(* assertions *)
	(uppercup[zz] /. zz->0) == (lowercup[zz] /. zz->lowercuplength);
	(D[uppercup[zz], zz] /. zz->0) ==
		(D[lowercup[zz], zz] /. zz->lowercuplength);

	(* putting it all together *)
	Which[
		z < 0,
			0,
		z < stemlength,
			stem[z],
		z < stemlength+lowercuplength,
			Max[stem[z], lowercup[z-stemlength]],
		z <= stemlength+lowercuplength+uppercuplength,
			uppercup[z-(stemlength+lowercuplength)],
		True,
			0]];

(* total glass length *)
glasslength[{footradius_, stemradius_, stemlength_,
	cupradius_, lowercuplength_, openingradius_, uppercuplength_}] :=
	stemlength+lowercuplength+uppercuplength;

(* maximal glass radius *)
glassradius[{footradius_, stemradius_, stemlength_,
	cupradius_, lowercuplength_, openingradius_, uppercuplength_}] :=
	Max[footradius, stemradius, cupradius, openingradius];

(* rotated profile plot *)
rotateplot[profile_, {zbegin_, zend_}] :=
	ParametricPlot3D[{
		Cos[phi] profile[z],
		Sin[phi] profile[z],
		z},
		{z, zbegin, zend}, {phi, 0, 2 Pi},
		DisplayFunction->Identity];

(* iso glass parameters, millimeters *)
isoglassparameters = {
	32.5, (* foot radius *)
	4.5, (* stem radius *)
	55, (* stem length *)
	32.5, (* maximal cup radius (==foot radius) *)
	32.5, (* lower cup length (==maximal cup radius) *)
	23, (* opening radius *)
	67.5}; (* upper cup length (==100-lower cup length) *)

(* function instances for the iso glass *)
isoglassprofile[z_] = glassprofile[z, isoglassparameters];
isoglasslength[] = glasslength[isoglassparameters];
isoglassradius[] = glassradius[isoglassparameters];

(* iso glass plot function *)
isoglassplot[] := rotateplot[isoglassprofile, {0, isoglasslength[]}];

(* package end *)
End[]
EndPackage[]

(* draw the profile *)
Plot[isoglassprofile[z], {z, 0, isoglasslength[]},
	PlotRange->{{0, isoglasslength[]}, {0, isoglassradius[]}},
	AxesOrigin->{0, 0},
	AspectRatio->isoglassradius[]/isoglasslength[]];

(* draw the glass *)
Show[isoglassplot[],
	Boxed->False,
	AxesEdge->None,
	DisplayFunction->$DisplayFunction];