(* James Lawrence University of Massachusetts at Amherst GANG Amherst, MA 01003 *) BeginPackage["EnvelopePackage`"] DrawEnvelope::usage = "DrawEnvelope[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr] evaluates and draws the numerically parameterized envelope function corresponding to the one-parameter plane curve f and the one-parameter circle-radius scalar function r along the domain tmin..tmax at increments of tincr. lincr is a small number which determines the accuracy of each evaluation point. DrawEnvelope[f[u,v],r[u,v],{u,umin,umax,uincr},{v,vmin,vmax,vincr},sqincr:defaultsqincr] evaluates and draws the numerically parameterized envelope function corresponding to the two-parameter surface f and the two-parameter sphere-radius scalar function r along the domain {umin..umax,vmin..vmax} at increments of {uincr,vincr}. sqincr is a small number which determines the accuracy of each evaluation point."; Envelope::usage = "Envelope[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr], Envelope[f[u,v],r[u,v],{u,umin,umax,uincr},{v,vmin,vmax,vincr},sqincr:defaultsqincr], is the unrendered graphics object of DrawEnvelope."; NEnvelopeFunction::usage = "NEnvelopeFunction[f[u,v],r[u,v],{u,v},{u0,v0},sqincr:defaultsqincr] is the numerically parameterized envelope function. It gives the point on the envelope corresponding to f[u0,v0] on original function."; DrawFunction::usage = "DrawFunction[f[t],{t,tmin,tmax}], DrawFunction[f[u,v],{u,umin,umax},{v,vmin,vmax}] draws f."; Funct::usage = "Funct[f[t],{t,tmin,tmax}], Funct[f[u,v],{u,umin,umax,uincr},{v,vmin,vmax,vincr}] is the unrendered graphics object of DrawFunction."; DrawBoth::usage = "DrawBoth[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr] draws the original function(blue) and the envelope(red) DrawBoth[f[u,v],r[u,v],{u,umin,umax,uincr},{v,vmin,vmax,vincr},sqincr:defaultsqincr], See DrawEnvelope for more information."; Both::usage = "Both[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr], Both[f[u,v],r[u,v],{u,umin,umax,uincr},{v,vmin,vmax,vincr},sqincr:defaultsqincr], is the unrendered graphics object of DrawBoth."; DrawEnvelopeSphere::usage = "DrawEnvelopeSphere[f[t],r[t],{t,t0}] draws one circle tangent to the curve f at f[t0] with radius r[t0]. DrawEnvelopeSphere[f[u,v],r[u,v],{u,v},{u0,v0}] draws one sphere tangent to the surface f at f[u0,v0] with radius r[u0,v0]."; EnvelopeSphere::usage = "EnvelopeSphere[f[t],r[t],{t,t0}], EnvelopeSphere[f[u,v],r[u,v],{u,v},{u0,v0}] is the unrendered graphics object of DrawEnvelopeSphere."; DrawAll::usage = "DrawAll[f[t],r[t],{t,tmin,tmax,tincr},t0,lincr:defaultlincr] draws the original function, the envelope, and one circle located at f[t0]. DrawAll[f[u,v],r[u,v],{u,umin,umax,uincr},{v,vmin,vmax,vincr},{u0,v0},sqincr:defaultsqincr], draws the original function, the envelope, and one sphere located at f[u0,v0]. When the sphere is very large, only a portion is drawn. See DrawEnvelope for more information."; Allll::usage = "Allll[f[t],r[t],{t,tmin,tmax,tincr},t0,lincr:defaultlincr], Allll[f[u,v],r[u,v],{u,umin,umax,uincr},{v,vmin,vmax,vincr},{u0,v0},sqincr:defaultsqincr], is the unrendered graphics object of DrawAll."; DrawEnvelopeSpheres::usage = "DrawEnvelopeSpheres[f[t],r[t],{t,tmin,tmax,tincr}] draws circles tangent to the one-parameter plane curve f each with radius corresponding to the one-parameter circle-radius scalar function r. The circles range from tmin..tmax, drawn at increments of tincr. DrawEnvelopeSpheres[f[u,v],r[u,v],{u,umin,umax,uincr},{v,vmin,vmax,vincr}] Draws spheres tangent to the two-parameter surface f each with radius corresponding to the two-parameter sphere-radius function r. The spheres are placed along {umin..umax,vmin..vmax}, spaced along intervals of {uincr,vincr}."; EnvelopeSpheres::usage = "EnvelopeSpheres[f[t],r[t],{t,tmin,tmax,tincr}], EnvelopeSpheres[f[u,v],r[u,v]],{u,umin,umax,uincr},{v,vmin,vmax,vincr}] is the unrendered graphics object of DrawSpheres."; EnvelopePoints::usage = "EnvelopePoints[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr] evaluates the numerically parameterized envelope function corresponding to the one-parameter plane curve f and the one-parameter circle-radius scalar function r along the domain tmin..tmax at increments of tincr. lincr is a small number which determines the accuracy of each evaluation point. The output is a list of points on the envelope curve. EnvelopePoints[f[u,v],r[u,v],{u,umin,umax,uincr},{v,vmin,vmax,vincr},sqincr:defaultsqincr] evaluates the numerically parameterized envelope function corresponding to the two-parameter surface f and the two-parameter sphere-raduis scalar function r along the domain {umin..umax,vmin..vmax} at increments of {uincr,vincr}. sqincr is a small number which determines the accuracy of each evaluation point. The output is a matrix of points that correspond to adjacent points on the envelope surface."; (*----------------------------------------------------------------*) DrawChannelEnvelope::usage = "DrawChannelEnvelope[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr] evaluates and draws the two numerically parameterized channel envelope functions corresponding to the one-parameter circle-center plane curve f and the one-parameter radius function r along the domain tmin..tmax at increments of tincr. lincr is a small number which determines the accuracy of each evaluation point. DrawChannelEnvelope[f[t],r[t],{t,tmin,tmax,tincr},spokes:defaultspokes,lincr:defaultlincr] draws the channel surface corresponding to the one-parameter sphere-center space-curve f and the one-parameter scalar radius function r along the domain tmin..tmax at increments of tincr. If f[t] is a straight line then it must be oriented {0,0,z(t)}. spokes is the number of evaulation points spaced around the tube at a particular point on the curve. lincr is a small number that determines the accuracy of each evaluation point."; ChannelEnvelopeAbove::usage = "ChannelEnvelopeAbove[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr] (two-dimensional case) is one of the unrendered envelopes of DrawChannelEnvelope."; ChannelEnvelopeBelow::usage = "ChannelEnvelopeBelow[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr], (two-dimensional case) is one of the unrendered envelopes of DrawChannelEnvelope."; ChannelEnvelope::usage = "ChannelEnvelope[f[t],r[t],{t,tmin,tmax,tincr},spokes:defaultspokes,lincr:defaultlincr] (three-dimensional case) is the unrendered graphics object of DrawChannelEnvelope."; DrawChannelFunction::usage = "DrawChannelFunction[f[t],{t,tmin,tmax}] draws the sphere center function f."; ChannelFunct::usage = "ChannelFunct[f[t],{t,tmin,tmax}] (two-dimensional) ChannelFunct[f[t],{t,tmin,tmax}] (three-dimensional) is the unrendered graphics object of DrawChannelFunction."; DrawChannelBoth::usage = "DrawChannelBoth[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr], DrawChannelBoth[f[t],r[t],{t,tmin,tmax,tincr},spokes:defaultspokes,lincr:defaultlincr], draws the original function and the envelope. See DrawChannelEnvelope for more information."; ChannelBoth::usage = "ChannelBoth[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr], ChannelBoth[f[t],r[t],{t,tmin,tmax,tincr},spokes:defaultspokes,lincr:defaultlincr], is the unrendered graphics object of DrawChannelBoth."; DrawChannelAll::usage = "DrawChannelAll[f[t],r[t],{t,tmin,tmax,tincr},t0,lincr:defaultlincr], DrawChannelAll[f[t],r[t],{t,tmin,tmax,tincr},t0,spokes:defaultspokes,lincr:defaultlincr], draws the original function, the envelope, and one sphere/circle at t=t0. See DrawChannelEnvelope for more information."; ChannelAll::usage = "ChannelAll[f[t],r[t],{t,tmin,tmax,tincr},t0,lincr:defaultlincr], ChannelAll[f[t],r[t],{t,tmin,tmax,tincr},spokes:defaultspokes,t0,lincr:defaultlincr], is the unrendered graphics object of DrawChannelAll."; DrawChannelSphere::usage = "DrawChannelSphere[f[t],r[t],{t,t0}], (two-dimensional) DrawChannelSphere[f[t],r[t],{t,t0}] (three-dimensional) draws one sphere/circle with center f[t0] and radius r[t0]"; ChannelSphere::usage = "ChannelSphere[f[t],r[t],{t,t0}], (two-dimensional), DrawChannelSphere[f[t],r[t],{t,t0}] (three-dimensional) is the unrendered graphics object of DrawChannelSphere."; DrawChannelSpheres::usage = "DrawChannelSpheres[f[t],r[t],{t,tmin,tmax,tincr}], (two-imensional) DrawChannelSpheres[f[t],r[t],{t,tmin,tmax,tincr}] (three-dimensional) draws spheres/circles along the the sphere-center function f[t] with radius r[t] along the domain tmin..tmax at increments of tincr."; ChannelSpheres::usage = "ChannelSpheres[f[t],r[t],{t,tmin,tmax,tincr}] (two-imensional) ChannelSpheres[f[t],r[t],{t,tmin,tmax,tincr}] (three-imensional) is the unrendered graphics object of DrawChannelSpheres."; ChannelEnvelopePoints::usage = "ChannelEnvelopePoints[f[t],r[t],{t,tmin,tmax,tincr},spokes:defaultspokes,lincr:defaultlincr] (three-dimensional case) gives the points on the channel envelope corresponding to the one-parameter sphere-center function f and the one-parameter radius function r along the domain tmin..tmax at increments of tincr. spokes is the number of evaulation points around the tube at a particular point. lincr is a small number that determines the accuracy of each evaluation point. The output is a matrix of points that correspond to points on the channel envelope."; ChannelEnvelopePointsAbove::usage = "ChannelEnvelopePointsAbove[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr] (two-dimensional case) evaluates one of the two numerically parameterized channel envelope functions corresponding to the one-parameter circle-center function f and the one-parameter radius function r along the domain tmin..tmax at increments of tincr. lincr is a small number which determines the accuracy of each evaluation point. The output is a list of points on the envelope."; ChannelEnvelopePointsBelow::usage = "ChannelEnvelopePointsBelow[f[t],r[t],{t,tmin,tmax,tincr},lincr:defaultlincr] (two-dimensional case) evaluates one of the two numerically parameterized channel envelope functions corresponding to the one-parameter circle-center function f and the one-parameter radius function r along the domain tmin..tmax at increments of tincr. lincr is a small number which determines the accuracy of each evaluation point. The output is a list of points on the envelope."; (*---------------------------*) DrawInsertSphere::usage = "DrawInsertSphere[f[u,v],{u,v},{u0,v0},{px,py,pz}] draws the sphere tangent to the surface f[u,v] at {u0,v0} and through the point {px,py,pz}. For use when the radius function is unknown."; InsertSphere::usage = "InsertSphere[f[u,v],{u,v},{u0,v0},{px,py,pz}] is the unrendered graphics object of DrawInsertSphere."; InterpolateSphereRadius::usage = "InterpolateSphereRadius[f[u,v],{u,v},{u0,v0},{px,py,pz}] gives the numerical radius of the interpolated sphere. See InsertSphere for more information"; InterpolateSphereCenter::usage = "InterpolateSphereCenter[f[u,v],{u,v},{u0,v0},{px,py,pz}] gives the numerical sphere center coordinate of the interpolated sphere. See InsertSphere for more information."; (*--------------------------*) SphereCenterFunction::usage = "SphereCenterFunction[f[u,v],r[u,v],{u,v}]" Begin["`Private`"] (*----constants-----*) defaultlincr = 0.01; defaultsqincr = 0.005; defaultspokes = 10; rot90 = {{0,-1},{1,0}}; (*------------------*) (*-----------------------some standard functions-----------------------*) Magnitude[{x_,y_}] = Sqrt[x^2 + y^2]; Magnitude[{x_,y_,z_}] = Sqrt[x^2 + y^2 + z^2]; Cross[{x1_,x2_,x3_},{y1_,y2_,y3_}] = { x2 y3 - x3 y2, -(x1 y3 - x3 y1), x1 y2 - x2 y1 }; NormalFunction[{fx_,fy_},t_] := ( (rot90 . D[{fx,fy},t])/Magnitude[D[{fx,fy},t]] ); NormalFunction[{fx_,fy_,fz_},{u_,v_}] := ( Cross[D[{fx,fy,fz},u],D[{fx,fy,fz},v]] / Magnitude[Cross[D[{fx,fy,fz},u],D[{fx,fy,fz},v]]] ); Sphere[{x_,y_,z_},r_] := ParametricPlot3D[ {x + r Cos[u] Sin[v], y + r Sin[u] Sin[v], z + r Cos[v]}, {u,0,2 Pi}, {v,0,Pi}, DisplayFunction->Identity ]; (*-----------------------------Envelopes 2D and 3D-------------------------*) SphereCenterFunction[{fx_,fy_},r_,t_] := ( {fx,fy} + r NormalFunction[{fx,fy},t] ); SphereCenterFunction[{fx_,fy_,fz_},r_,{u_,v_}] := ( {fx,fy,fz} + r NormalFunction[{fx,fy,fz},{u,v}] ); SingleEnvelopePoint[a1_,n1_,r1_,a2_,n2_,r2_] := Module[{c1,c2,d,sinalpha,cosalpha,sinbeta,cosbeta,b1,b2}, c1 = a1 + r1 n1; c2 = a2 + r2 n2; d = Magnitude[c2 - c1]; cosalpha = (r1^2 + d^2 - r2^2)/(2 r1 d); sinalpha = Sqrt[Abs[1 - cosalpha^2]]; cosbeta = (c2 - c1) . (a1 - c1) / (d r1); sinbeta = Sqrt[Abs[1 - cosbeta^2]]; b1 = r1 (sinalpha/sinbeta); b2 = r1 (sinalpha cosbeta + cosalpha sinbeta) / sinbeta; c1 + b1 n1 + b2 (c2 - c1)/d ]; ThreeSpheresIntersection[c1_,r1_,c2_,r2_,c3_,r3_] := Module[{cosalpha,cosbeta,magnitudec2c1, magnitudec3c2,basenormal,h1,h2,eq1,eq2,eq3,k,h}, magnitudec2c1 = Magnitude[c2 - c1]; magnitudec3c2 = Magnitude[c3 - c2]; basenormal = Cross[c2 - c1,c3 - c1]; cosalpha = (r1^2 + magnitudec2c1^2 - r2^2)/(2 r1 magnitudec2c1); cosbeta = (r2^2 + magnitudec3c2^2 - r3^2)/(2 r2 magnitudec3c2); h1 = c1 + r1 cosalpha (c2 - c1)/magnitudec2c1; h2 = c2 + r2 cosbeta (c3 - c2)/magnitudec3c2; eq1 = Expand[({x,y,z} - h1) . (c2 - c1)]; eq2 = Expand[({x,y,z} - h2) . (c3 - c2)]; eq3 = Expand[({x,y,z} - c1) . basenormal]; k = ThreePlanesIntersection[ {Coefficient[eq1,x], Coefficient[eq1,y], Coefficient[eq1,z], eq1 /. {x->0,y->0,z->0}}, {Coefficient[eq2,x], Coefficient[eq2,y], Coefficient[eq2,z], eq2 /. {x->0,y->0,z->0}}, {Coefficient[eq3,x], Coefficient[eq3,y], Coefficient[eq3,z], eq3 /. {x->0,y->0,z->0}} ]; h = Sqrt[Abs[r1^2 - Magnitude[k - c1]^2]]; k + h basenormal/Magnitude[basenormal] ]; ThreePlanesIntersection[{a1_,b1_,c1_,d1_},{a2_,b2_,c2_,d2_},{a3_,b3_,c3_,d3_}] = { (-(b3 c2 d1) + b2 c3 d1 + b3 c1 d2 - b1 c3 d2 - b2 c1 d3 + b1 c2 d3)/ ( a3 b2 c1 - a2 b3 c1 - a3 b1 c2 + a1 b3 c2 + a2 b1 c3 - a1 b2 c3), (-(a3 c2 d1) + a2 c3 d1 + a3 c1 d2 - a1 c3 d2 - a2 c1 d3 + a1 c2 d3)/ (-(a3 b2 c1) + a2 b3 c1 + a3 b1 c2 - a1 b3 c2 - a2 b1 c3 + a1 b2 c3), ( a3 b2 d1 - a2 b3 d1 - a3 b1 d2 + a1 b3 d2 + a2 b1 d3 - a1 b2 d3 )/ (-(a3 b2 c1) + a2 b3 c1 + a3 b1 c2 - a1 b3 c2 - a2 b1 c3 + a1 b2 c3) }; NEnvelopeFunctionHelper[{fx_,fy_},r_,norm_,{t_,t0_},lincr_:defaultlincr] := ( ( SingleEnvelopePoint[N[{fx,fy}/.t->t0],N[norm/.t->t0],N[r/.t->t0], N[{fx,fy}/.t->t0+lincr],N[norm/.t->t0+lincr],N[r/.t->t0+lincr]] + SingleEnvelopePoint[N[{fx,fy}/.t->t0],N[norm/.t->t0],N[r/.t->t0], N[{fx,fy}/.t->t0-lincr],N[norm/.t->t0-lincr],N[r/.t->t0-lincr]] )/2 ); NEnvelopeFunctionHelper[{fx_,fy_,fz_},r_,scf_,{u_,v_},{u0_,v0_},sqincr_:defaultlincr] := Module[{c,c1,c2,c3,c4,r1,r2,r3,r4,rc}, c1 = N[scf /. {u->u0+sqincr,v->v0}]; c2 = N[scf /. {u->u0,v->v0+sqincr}]; c3 = N[scf /. {u->u0-sqincr,v->v0}]; c4 = N[scf /. {u->u0,v->v0-sqincr}]; c = N[scf /. {u->u0,v->v0}]; r1 = N[r /. {u->u0+sqincr,v->v0}]; r2 = N[r /. {u->u0,v->v0+sqincr}]; r3 = N[r /. {u->u0-sqincr,v->v0}]; r4 = N[r /. {u->u0,v->v0-sqincr}]; rc = N[r /. {u->u0,v->v0}]; ( ThreeSpheresIntersection[c,rc,c1,r1,c2,r2] + ThreeSpheresIntersection[c,rc,c2,r2,c3,r3] + ThreeSpheresIntersection[c,rc,c3,r3,c4,r4] + ThreeSpheresIntersection[c,rc,c4,r4,c1,r1] )/4 ]; NEnvelopeFunction[{fx_,fy_,fz_},r_,{u_,v_},{u0_,v0_},sqincr_:defaultlincr] := Module[{scf}, scf = SphereCenterFunction[{fx,fy,fz},r,{u,v}]; NEnvelopeFunctionHelper[{fx,fy,fz},r,scf,{u,v},{u0,v0},sqincr] ] (*----*) EnvelopePoints[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := Module[{norm}, norm = NormalFunction[{fx,fy},t]; Table[ NEnvelopeFunctionHelper[{fx,fy},r,norm,{t,x},lincr], {x,tmin,tmax,tincr} ] ]; EnvelopePoints[{fx_,fy_,fz_},r_,{u_,umin_,umax_,uincr_},{v_,vmin_,vmax_,vincr_},sqincr_:defaultsqincr] := Module[{scf}, scf = SphereCenterFunction[{fx,fy,fz},r,{u,v}]; Table[ NEnvelopeFunctionHelper[{fx,fy,fz},r,scf,{u,v},{pp,qq},sqincr], {pp,umin,umax,uincr}, {qq,vmin,vmax,vincr} ] ]; (*-------*) EnvelopeSphere[{fx_,fy_},r_,{t_,t0_}] := Show[ Graphics[Circle[SphereCenterFunction[{fx,fy},r,t]/.{t->t0},r/.t->t0]], DisplayFunction->Identity]; EnvelopeSphere[{fx_,fy_,fz_},r_,{u_,v_},{u0_,v0_}] := Show[ InsertSphere[{fx,fy,fz},{u,v},{u0,v0},NEnvelopeFunction[{fx,fy,fz},r,{u,v},{u0,v0}]], DisplayFunction->Identity]; DrawEnvelopeSphere[{fx_,fy_},r_,{t_,t0_}] := Show[ EnvelopeSphere[{fx,fy},r,{t,t0}], DisplayFunction->$DisplayFunction]; DrawEnvelopeSphere[{fx_,fy_,fz_},r_,{u_,v_},{u0_,v0_}] := Show[ EnvelopeSphere[{fx,fy,fz},r,{u,v},{u0,v0}], DisplayFunction->$DisplayFunction]; EnvelopeSpheres[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_}] := Module[{g}, g = SphereCenterFunction[{fx,fy},r,t]; Show[ Table[ Graphics[ Circle[g /.t->x, Abs[r/.t->x] ]], {x,tmin,tmax,tincr}], DisplayFunction->Identity] ]; EnvelopeSpheres[{fx_,fy_,fz_},r_,{u_,umin_,umax_,uincr_},{v_,vmin_,vmax_,vincr_}] := Show[ Table[ EnvelopeSphere[{fx,fy,fz},r,{u,v},{xx,yy}], {xx,umin,umax,uincr}, {yy,vmin,vmax,vincr}], DisplayFunction->Identity ]; DrawEnvelopeSpheres[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_}] := Show[ EnvelopeSpheres[{fx,fy},r,{t,tmin,tmax,tincr}], DisplayFunction->$DisplayFunction ]; DrawEnvelopeSpheres[{fx_,fy_,fz_},r_,{u_,umin_,umax_,uincr_},{v_,vmin_,vmax_,vincr_}] := Show[ EnvelopeSpheres[{fx,fy,fz},r,{u,umin,umax,uincr},{v,vmin,vmax,vincr}], DisplayFunction->$DisplayFunction ]; (*-------*) Allll[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},t0_,lincr_:defaultlincr] := Show[ Both[{fx,fy},r,{t,tmin,tmax,tincr},lincr], EnvelopeSphere[{fx,fy},r,{t,t0}], DisplayFunction->Identity ]; Allll[{fx_,fy_,fz_},r_,{u_,umin_,umax_,uincr_},{v_,vmin_,vmax_,vincr_},{u0_,v0_},sqincr_:defaultsqincr] := Show[ Both[{fx,fy,fz},r,{u,umin,umax,uincr},{v,vmin,vmax,vincr},sqincr], EnvelopeSphere[{fx,fy,fz},r,{u,v},{u0,v0}], DisplayFunction->Identity ]; DrawAll[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},t0_,lincr_:defaultlincr] := Show[ Allll[{fx,fy},r,{t,tmin,tmax,tincr},t0,lincr], DisplayFunction->$DisplayFunction]; DrawAll[{fx_,fy_,fz_},r_,{u_,umin_,umax_,uincr_},{v_,vmin_,vmax_,vincr_},{u0_,v0_},sqincr_:defaultsqincr] := Show[ Allll[{fx,fy,fz},r,{u,umin,umax,uincr},{v,vmin,vmax,vincr},{u0,v0},sqincr], DisplayFunction->$DisplayFunction]; (*-------*) Both[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := Show[ Funct[{fx,fy},{t,tmin,tmax}], Envelope[{fx,fy},r,{t,tmin,tmax,tincr},lincr]]; Both[{fx_,fy_,fz_},r_,{u_,umin_,umax_,uincr_},{v_,vmin_,vmax_,vincr_},sqincr_:defaultsqincr] := Show[ Funct[{fx,fy,fz},{u,umin,umax},{v,vmin,vmax}], Envelope[{fx,fy,fz},r,{u,umin,umax,uincr},{v,vmin,vmax,vincr},sqincr]]; DrawBoth[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := Show[ Both[{fx,fy},r,{t,tmin,tmax,tincr},lincr], DisplayFunction->$DisplayFunction ]; DrawBoth[{fx_,fy_,fz_},r_,{u_,umin_,umax_,uincr_},{v_,vmin_,vmax_,vincr_},sqincr_:defaultsqincr] := Show[ Both[{fx,fy,fz},r,{u,umin,umax,uincr},{v,vmin,vmax,vincr},sqincr], DisplayFunction->$DisplayFunction ]; (*-------*) Funct[{fx_,fy_},{t_,tmin_,tmax_}] := ParametricPlot[{fx,fy},{t,tmin,tmax}, DisplayFunction->Identity, PlotStyle->{RGBColor[0,0,1]}]; Funct[{fx_,fy_,fz_},{u_,umin_,umax_},{v_,vmin_,vmax_}] := ParametricPlot3D[{fx,fy,fz},{u,umin,umax},{v,vmin,vmax}, DisplayFunction->Identity]; DrawFunction[{fx_,fy_},{t_,tmin_,tmax_}] := Show[ Funct[{fx,fy},{t,tmin,tmax}], DisplayFunction->$DisplayFunction]; DrawFunction[{fx_,fy_,fz_},{u_,umin_,umax_},{v_,vmin_,vmax_}] := Show[ Funct[{fx,fy,fz},{u,umin,umax},{v,vmin,vmax}], DisplayFunction->$DisplayFunction]; (*-------*) Envelope[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := ListPlot[ EnvelopePoints[{fx,fy},r,{t,tmin,tmax,tincr},lincr], PlotJoined->True, PlotStyle->{RGBColor[1,0,0]}, DisplayFunction->Identity ]; Envelope[{fx_,fy_,fz_},r_,{u_,umin_,umax_,uincr_}, {v_,vmin_,vmax_,vincr_},sqincr_:defaultsqincr] := Module[{pts}, pts = EnvelopePoints[{fx,fy,fz},r,{u,umin,umax,uincr},{v,vmin,vmax,vincr},sqincr]; Print["Done calculating the points on the envelope. Now forming polygons..."]; Show[ Graphics3D[ Table[ Polygon[ {pts[[i,j]],pts[[i+1,j]],pts[[i+1,j+1]],pts[[i,j+1]]} ], {i,1,Length[pts] - 1}, {j,1,Length[Part[pts,1]] - 1} ] ], DisplayFunction->Identity] ]; DrawEnvelope[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := Show[ Envelope[{fx,fy},r,{t,tmin,tmax,tincr},lincr], DisplayFunction->$DisplayFunction ]; DrawEnvelope[{fx_,fy_,fz_},r_,{u_,umin_,umax_,uincr_}, {v_,vmin_,vmax_,vincr_},sqincr_:defaultsqincr] := Show[ Envelope[{fx,fy,fz},r,{u,umin,umax,uincr},{v,vmin,vmax,vincr},sqincr], DisplayFunction->$DisplayFunction]; (*---------------------------Channel Envelopes 2D--------------------------*) SingleChannelEnvelopePointAbove[c1_,r1_,c2_,r2_] := Module[{d,cosalpha,sinalpha}, d = Magnitude[c1 - c2]; cosalpha = (r1^2 + d^2 - r2^2)/(2 r1 d); sinalpha = Sqrt[Abs[1 - cosalpha^2]]; c1 + (r1 cosalpha / d) (c2 - c1) + (r1 sinalpha / d) ( rot90 . (c2 - c1) ) ]; SingleChannelEnvelopePointBelow[c1_,r1_,c2_,r2_] := Module[{d,cosalpha,sinalpha}, d = Magnitude[c1 - c2]; cosalpha = (r1^2 + d^2 - r2^2)/(2 r1 d); sinalpha = Sqrt[Abs[1 - cosalpha^2]]; c1 + (r1 cosalpha / d) (c2 - c1) - (r1 sinalpha / d) ( rot90 . (c2 - c1) ) ]; (*-----*) NChannelEnvelopeFunctionAbove[{fx_,fy_},r_,{t_,t0_},lincr_] := ( ( SingleChannelEnvelopePointAbove[N[{fx,fy}/.t->t0], N[r/.t->t0], N[{fx,fy}/.t->t0+lincr], N[r/.t->t0+lincr]] + SingleChannelEnvelopePointBelow[N[{fx,fy}/.t->t0], N[r/.t->t0], N[{fx,fy}/.t->t0-lincr], N[r/.t->t0-lincr]] )/2 ); NChannelEnvelopeFunctionBelow[{fx_,fy_},r_,{t_,t0_},lincr_] := ( ( SingleChannelEnvelopePointBelow[N[{fx,fy}/.t->t0], N[r/.t->t0], N[{fx,fy}/.t->t0+lincr], N[r/.t->t0+lincr]] + SingleChannelEnvelopePointAbove[N[{fx,fy}/.t->t0], N[r/.t->t0], N[{fx,fy}/.t->t0-lincr], N[r/.t->t0-lincr]] )/2 ); (*-----*) ChannelEnvelopePointsAbove[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := ( Table[ NChannelEnvelopeFunctionAbove[{fx,fy},r,{t,x},lincr], {x,tmin,tmax,tincr} ] ); ChannelEnvelopePointsBelow[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := ( Table[ NChannelEnvelopeFunctionBelow[{fx,fy},r,{t,x},lincr], {x,tmin,tmax,tincr} ] ); (*-----*) ChannelEnvelopeAbove[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := Show[ ListPlot[ ChannelEnvelopePointsAbove[{fx,fy},r,{t,tmin,tmax,tincr},lincr], DisplayFunction->Identity, PlotJoined->True, PlotStyle->{RGBColor[1,0,0]} ] ]; ChannelEnvelopeBelow[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := Show[ ListPlot[ ChannelEnvelopePointsBelow[{fx,fy},r,{t,tmin,tmax,tincr},lincr], DisplayFunction->Identity, PlotJoined->True, PlotStyle->{RGBColor[1,0,0]} ] ]; (*-----*) DrawChannelEnvelope[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := Show[ ChannelEnvelopeAbove[{fx,fy},r,{t,tmin,tmax,tincr},lincr], ChannelEnvelopeBelow[{fx,fy},r,{t,tmin,tmax,tincr},lincr], DisplayFunction->$DisplayFunction ]; (*-----*) ChannelFunct[{fx_,fy_},{t_,tmin_,tmax_}] := ParametricPlot[ {fx,fy}, {t,tmin,tmax}, PlotStyle->{RGBColor[0,0,1]}, DisplayFunction->Identity]; DrawChannelFunction[{fx_,fy_},{t_,tmin_,tmax_}] := Show[ ChannelFunct[{fx,fy},{t,tmin,tmax}], DisplayFunction->$DisplayFunction ]; (*-----*) ChannelBoth[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := Show[ ChannelFunct[{fx,fy},{t,tmin,tmax}], ChannelEnvelopeAbove[{fx,fy},r,{t,tmin,tmax,tincr},lincr], ChannelEnvelopeBelow[{fx,fy},r,{t,tmin,tmax,tincr},lincr], DisplayFunction->Identity ]; DrawChannelBoth[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},lincr_:defaultlincr] := Show[ ChannelBoth[{fx,fy},r,{t,tmin,tmax,tincr},lincr], DisplayFunction->$DisplayFunction]; (*-----*) ChannelAll[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},t0_,lincr_:defaultlincr] := Show[ ChannelBoth[{fx,fy},r,{t,tmin,tmax,tincr},lincr], ChannelSphere[{fx,fy},r,{t,t0}], DisplayFunction->Identity]; DrawChannelAll[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_},t0_,lincr_:defaultlincr] := Show[ ChannelAll[{fx,fy},r,{t,tmin,tmax,tincr},t0,lincr], DisplayFunction->$DisplayFunction ]; (*-----*) ChannelSphere[{fx_,fy_},r_,{t_,t0_}] := Graphics[Circle[ {fx,fy}/.t->t0, r/.t->t0 ]]; DrawChannelSphere[{fx_,fy_},r_,{t_,t0_}] := Show[ ChannelSphere[{fx,fy},r,{t,t0}], DisplayFunction->$DisplayFunction ]; ChannelSpheres[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_}] := Show[ Table[ Graphics[Circle[{fx,fy}/.t->x, Abs[r/.t->x]]], {x,tmin,tmax,tincr}], DisplayFunction->Identity]; DrawChannelSpheres[{fx_,fy_},r_,{t_,tmin_,tmax_,tincr_}] := Show[ ChannelSpheres[{fx,fy},r,{t,tmin,tmax,tincr}], DisplayFunction->$DisplayFunction]; (*----------------------------channel envelopes 3D-----------------------------*) LocalCoordFunctions[{fx_,fy_,fz_},t_] := Module[{v,a,bi}, v = D[{fx,fy,fz},t]; a = D[v,t]; bi = Cross[v,a]; {v/Magnitude[v], a/Magnitude[a], bi/Magnitude[bi]} ]; EnvelopePivot[c1_,r1_,c2_,r2_] := Module[{cosalpha}, d = Magnitude[c2 - c1]; cosalpha = (r1^2 + d^2 - r2^2)/(2 r1 d); c1 + (r1 cosalpha) (c2 - c1)/Magnitude[c2 - c1] ]; EnvelopeRadius[c1_,r1_,c2_,r2_] := Module[{cosalpha,sinalpha}, d = Magnitude[c2 - c1]; cosalpha = (r1^2 + d^2 - r2^2)/(2 r1 d); sinalpha = Sqrt[Abs[1 - cosalpha^2]]; r1 sinalpha ]; NEnvelopePivotFunction[{fx_,fy_,fz_},r_,{t_,t0_},lincr_] := ( ( EnvelopePivot[N[{fx,fy,fz}/.t->t0],N[r/.t->t0],N[{fx,fy,fz}/.t->t0+lincr],N[r/.t->t0+lincr]] + EnvelopePivot[N[{fx,fy,fz}/.t->t0],N[r/.t->t0],N[{fx,fy,fz}/.t->t0-lincr],N[r/.t->t0-lincr]] )/2 ); NEnvelopeRadiusFunction[{fx_,fy_,fz_},r_,{t_,t0_},lincr_] := ( ( EnvelopeRadius[N[{fx,fy,fz}/.t->t0],N[r/.t->t0],N[{fx,fy,fz}/.t->t0+lincr],N[r/.t->t0+lincr]] + EnvelopeRadius[N[{fx,fy,fz}/.t->t0],N[r/.t->t0],N[{fx,fy,fz}/.t->t0-lincr],N[r/.t->t0-lincr]] )/2 ); ChannelEnvelopePoints[{fx_,fy_,fz_},r_,{t_,tmin_,tmax_,tincr_},spokes_:defaultspokes,lincr_:defaultlincr] := Module[{coords}, coords = If[(fx == 0) && (fy == 0), (* straight lines don't have a canonical frame *) {{0,0,1},{0,1,0},{1,0,0}}, LocalCoordFunctions[{fx,fy,fz},t], LocalCoordFunctions[{fx,fy,fz},t] ]; Table[ NEnvelopePivotFunction[{fx,fy,fz},r,{t,x},lincr] + NEnvelopeRadiusFunction[{fx,fy,fz},r,{t,x},lincr] N[Cos[theta]] N[Part[coords,2]/.t->x] + NEnvelopeRadiusFunction[{fx,fy,fz},r,{t,x},lincr] N[Sin[theta]] N[Part[coords,3]/.t->x], {x,tmin,tmax,tincr}, {theta,0,2Pi,2Pi/spokes} ] ]; ChannelEnvelope[{fx_,fy_,fz_},r_,{t_,tmin_,tmax_,tincr_},spokes_:defaultspokes,lincr_:defaultlincr] := Module[{pts}, pts = ChannelEnvelopePoints[{fx,fy,fz},r,{t,tmin,tmax,tincr},spokes,lincr]; Print["Done calculating the points on the envelope. Now forming polygons..."]; Show[ Graphics3D[ Table[ Polygon[ {pts[[i,j]],pts[[i+1,j]],pts[[i+1,j+1]],pts[[i,j+1]]} ], {i,1,Length[pts] - 1}, {j,1,Length[Part[pts,1]] - 1} ] ], DisplayFunction->Identity] ]; DrawChannelEnvelope[{fx_,fy_,fz_},r_,{t_,tmin_,tmax_,tincr_},spokes_:defaultspokes,lincr_:defaultlincr] := Show[ ChannelEnvelope[{fx,fy,fz},r,{t,tmin,tmax,tincr},spokes,lincr], DisplayFunction->$DisplayFunction ]; ChannelFunct[{fx_,fy_,fz_},{t_,tmin_,tmax_}] := ParametricPlot3D[{fx,fy,fz},{t,tmin,tmax},DisplayFunction->Identity]; DrawChannelFunction[{fx_,fy_,fz_},{t_,tmin_,tmax_}] := Show[ ChannelFunct[{fx,fy,fz},{t,tmin,tmax}], DisplayFunction->$DisplayFunction]; ChannelBoth[{fx_,fy_,fz_},r_,{t_,tmin_,tmax_,tincr_},spokes_:defaultspokes,lincr_:defaultlincr] := Show[ ChannelFunct[{fx,fy,fz},{t,tmin,tmax}], ChannelEnvelope[{fx,fy,fz},r,{t,tmin,tmax,tincr},spokes,lincr], DisplayFunction->Identity]; DrawChannelBoth[{fx_,fy_,fz_},r_,{t_,tmin_,tmax_,tincr_},spokes_:defaultspokes,lincr_:defaultlincr] := Show[ ChannelBoth[{fx,fy,fz},r,{t,tmin,tmax,tincr},spokes,lincr], DisplayFunction->$DisplayFunction ]; ChannelSphere[{fx_,fy_,fz_},r_,{t_,t0_}] := Show[Sphere[{fx,fy,fz}/.t->t0, r/.t->t0], DisplayFunction->Identity]; DrawChannelSphere[{fx_,fy_,fz_},r_,{t_,t0_}] := Show[ ChannelSphere[{fx,fy,fz},r,{t,t0}], DisplayFunction->$DisplayFunction]; ChannelAll[{fx_,fy_,fz_},r_,{t_,tmin_,tmax_,tincr_},t0_,spokes_:defaultspokes,lincr_:defaultlincr] := Show[ ChannelSphere[{fx,fy,fz},r,{t,t0}], ChannelBoth[{fx,fy,fz},r,{t,tmin,tmax,tincr},spokes,lincr], DisplayFunction->Identity ]; DrawChannelAll[{fx_,fy_,fz_},r_,{t_,tmin_,tmax_,tincr_},t0_,spokes_:defaultspokes,lincr_:defaultlincr] := Show[ ChannelAll[{fx,fy,fz},r,{t,tmin,tmax,tincr},t0,spokes,lincr], DisplayFunction->$DisplayFunction ]; ChannelSpheres[{fx_,fy_,fz_},r_,{t_,tmin_,tmax_,tincr_}] := Show[ Table[ Sphere[{fx,fy,fz}/.t->x,r/.t->x], {x,tmin,tmax,tincr}], DisplayFunction->Identity]; DrawChannelSpheres[{fx_,fy_,fz_},r_,{t_,tmin_,tmax_,tincr_}] := Show[ ChannelSpheres[{fx,fy,fz},r,{t,tmin,tmax,tincr}], DisplayFunction->$DisplayFunction]; (*---------------------------------*) InsertSphere[{fx_,fy_,fz_},{u_,v_},{u0_,v0_},b_] := Module[{a,r0,c,e1,e2,e3}, a = {fx,fy,fz} /. {u->u0,v->v0}; normvect = NormalFunction[{fx,fy,fz},{u,v}] /. {u->u0,v->v0}; costheta = (normvect . (b - a))/(Magnitude[normvect] Magnitude[b - a]); r0 = Magnitude[(b - a)/2]/costheta; c = {fx,fy,fz} + r0 NormalFunction[{fx,fy,fz},{u,v}] /. {u->u0,v->v0}; If[costheta == 1, (*change coordinates when envelope point lies on normal*) e1 = D[{fx,fy,fz},u]/Magnitude[D[{fx,fy,fz},u]] /. {u->u0,v->v0}; e2 = D[{fx,fy,fz},v]/Magnitude[D[{fx,fy,fz},v]] /. {u->u0,v->v0}; e3 = normvect;, e3 = N[(c - (a + b)/2)/Magnitude[(c - (a + b)/2)]]; e2 = N[(b - a)/Magnitude[b - a]]; e1 = N[Cross[e2,e3]/Magnitude[Cross[e2,e3]]]; ]; If[N[costheta] >= N[Sqrt[2]/2], ParametricPlot3D[ c + (r0 Cos[s] Sin[t]) e1 + (r0 Sin[s] Sin[t]) e2 + (r0 Cos[t]) e3, {s,0,2Pi},{t,0,Pi},DisplayFunction->Identity], ParametricPlot3D[ c + (r0 Cos[s] Sin[t]) e1 + (r0 Sin[s] Sin[t]) e2 + (r0 Cos[t]) e3, {s,0,2Pi},{t,Pi/2 + ArcCos[Abs[costheta]],Pi},DisplayFunction->Identity] ] ]; DrawInsertSphere[{fx_,fy_,fz_},{u_,v_},{u0_,v0_},b_] := Show[InsertSphere[{fx,fy,fz},{u,v},{u0,v0},b],DisplayFunction->Identity]; InterpolateSphereRadius[{fx_,fy_,fz_},{u_,v_},{u0_,v0_},b_] := Module[{a,normvect,costheta}, a = N[{fx,fy,fz} /. {u->u0,v->v0}]; normvect = N[NormalFunction[{fx,fy,fz},{u,v}] /. {u->u0,v->v0}]; costheta = N[(normvect . (b - a))/(Magnitude[normvect] Magnitude[b - a])]; N[Magnitude[(b - a)/2]/costheta] ]; InterpolateSphereCenter[{fx_,fy_,fz_},{u_,v_},{u0_,v0_},b_] := Module[{r0}, r0 = InterpolateSphereRadius[{fx,fy,fz},{u,v},{u0,v0},b]; N[{fx,fy,fz} + r0 NormalFunction[{fx,fy,fz},{u,v}] /. {u->u0,v->v0}] ]; End [ ] EndPackage[ ]