fit a tube around a trajectory
- To: mathgroup at smc.vnet.net
- Subject: [mg59173] fit a tube around a trajectory
- From: Georgios Oikonomou <georgios.oikonomou at gmail.com>
- Date: Sun, 31 Jul 2005 01:30:45 -0400 (EDT)
- Reply-to: Georgios Oikonomou <georgios.oikonomou at gmail.com>
- Sender: owner-wri-mathgroup at wolfram.com
Dear forum
below is one of the problems I am trying to solve lately. Actually is
a simplified version of my problem but the structure is the same. So I
am having three trajectories plotted from the program below. Now I
want to fit a tube around the resulted trajectories. I was trying to
use the a package that its code is illustrated below. However it didnt
work and I was trying to change the m.File in order to fit my plots.
Unfortunately without any luck.
If anyone has a clue how I can short it out it, I would be greatful. I
suspect that the problem has to do with the derivative in the m.file.
(* .nb file for trajectories*)
(* parameters *)
aa = 1.0;
bb = 2.0;
(* starting points *)
x0 = aa;
y0 = bb;
z0 = 0.0;
(* angular velocity in orbit *)
w = 2.0*Pi;
cc = 0.5;
c1 = 0.0;
k1 = 0.0;
k2 = k1;
k3 = k1;
(* velocity *)
vx[xx_, yy_, zz_, tt_, c1_,
c2_, c3_, d1_, d2_, d3_] := -aa*Sin[w*tt] + Sin[xx];
vy[xx_, yy_, zz_, tt_, c1_, c2_, c3_, d1_, d2_, d3_] := bb*Cos[w*tt] + \
Sin[zz] + Sin[yy];
vz[xx_, yy_, zz_, tt_, c1_, c2_, c3_, d1_, d2_, d3_] := cc +
Sin[yy] + Exp[-zz];
(* solve for trajectory : outputs interpolating function *)
sol0 = NDSolve[
{x'[t] == N[vx[x[t], y[t], z[t], t, k1, k2, k3, 1, 2, 3]],
y'[t] == N[vy[x[t], y[t], z[t], t, k1, k2, k3, 1, 2, 3]],
z'[t] == N[vz[x[t], y[t], z[t], t, k1, k2, k3, 1, 2, 3]],
x[0] == x0, y[0] == y0, z[0] == z0},
{x, y, z}, {t, 0, 10}, Method -> ExplicitRungeKutta];
(*Plot trajectory in real space *)
pic0 = ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol0], {t, 0, 5}, \
PlotRange -> All];
(* do this for three orbits with different starting points *)
(* number of trajectories *)
Ntrajs = 3;
(* initialise a set of initial xx, yy, zz values *)
xxstart = Table[0.0, {jj, 1, Ntrajs}];
yystart = Table[0.0, {jj, 1, Ntrajs}];
zzstart = Table[0.0, {jj, 1, Ntrajs}];
(*Choose initial starting points to lie in x = aa, z =
0 plane uniformily in range - bb < yy < bb *)
Do[
xxstart[[jj]] = aa;
yystart[[jj]] = -bb + 2*bb*(jj - 1)/(Ntrajs - 1);
zzstart[[jj]] = 0.0,
{jj, 1, Ntrajs}];
(*initialise array of trajectories *)
sol = Table[sol0, {jj, 1, Ntrajs}];
(*initialise graphics array of trajectories *)
pic = Table[pic0, {jj, 1, Ntrajs}];
(* Do loop to calculate Njtrajs trajectories *)
Do[
sol[[jj]] = NDSolve[
{x'[t] == N[vx[x[t], y[t], z[t], t, k1, k2, k3, 1, 2, 3]],
y'[t] == N[vy[x[t], y[t], z[t], t, k1, k2, k3,
1, 2, 3]], z'[t] == N[vz[x[t], y[t], z[t], t, k1, k2, k3, 1, 2, 3]],
x[0] == xxstart[[jj]], y[0] == yystart[[jj]], z[0] == zzstart[[
jj]]},
{x, y, z}, {t, 0, 10}, Method -> ExplicitRungeKutta];
pic[[jj]] = ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. sol[[jj]]], {
t, 0, 1}, PlotRange -> All, AspectRatio -> 1], {jj, 1, Ntrajs}];
Show[pic[[1]], pic[[2]], pic[[3]]]
(*--------------------------------------------*)
(*Tubes.m file*)
(* code is taken from the web,
http://www.unca.edu/~mcmcclur/mathematicaGraphics/index.html *)
TubePlot[curve_List, {var_, min_, max_}, radius_,
crossVector_List:{1, 1, 1}, opts___] :=
Module[{tangent, unitTangent, normal,
unitNormal, biNormal},
tangent = D[curve, t];
unitTangent = tangent/Sqrt[tangent.tangent];
normal = Cross[tangent, crossVector];
unitNormal = normal/Sqrt[normal.normal];
biNormal = Cross[unitTangent, unitNormal];
ParametricPlot3D[
curve + radius Cos[s] unitNormal + radius Sin[s] biNormal //
Evaluate,
{var, min, max}, {s, 0, 2Pi}, opts]]
(*-------Example from the author of Tubes.m-----The Klein Bottle---*)
d = .00000001;
r = If[-Pi/2 <= t && t <= Pi/2,
(1 - Sqrt[Cos[t]]/2), (1 + Sqrt[-Cos[t]]/2)];
g = TubePlot[{ 2Cos[t] (1 + Sin[t]), 6 Sin[t], 0},
{t, -Pi/2 + d, 3Pi/2 - d}, r, {0, 0, 1},
PlotPoints -> {48, 12}, Axes -> False, Boxed -> False,
LightSources -> {{{1., 0., 1.}, RGBColor[1, 0, 0]},
{{1., 1., 1.}, RGBColor[0, 1, 0]}, {{0., 1., 1.},
RGBColor[0, 0, 1]}, {{-1., 0., -1.}, RGBColor[1, 0, 0]},
{{-1., -1., -1.}, RGBColor[0, 1, 0]}, {{0., -1., -1.},
RGBColor[0, 0, 1]}}]
Thanx for your time.
Ragards
Georgios