MathGroup Archive 2005

[Date Index] [Thread Index] [Author Index]

Search the Archive

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


  • Prev by Date: Re: Add terms surrounded by zero together in matrix
  • Next by Date: Re: Add terms surrounded by zero together in matrix
  • Previous by thread: BlankSequence