Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I have an interplanetary trajectory solver that calculates a spacecraft route from Earth to Mars, and I am trying to zoom into the spacecraft when it reaches its Mars intercept orbit but am having some trouble. I originally asked a question of a similar nature on these forums (Zooming into animated plot), but used a simple parametric function as an example, and although I was given a good answer, I was unable to translate that solution to the trajectory solver as it is a bit more involved. Here is a stripped down version of the code which shows the orbits of Earth (inner circle), Mars (outer circle) and the spacecraft (line joining the inner and outer circles). Although what I tried hasn't worked, if you look at the PlotRange code I used when making the trajectory animation, I used an if function to change the PlotRange and essentially "zoom-in" to the spacecraft when it got close to Mars; unfortunately this code did not work, but I have left it in as it is essentially what I want to do in pseudo-code if you will. Any help would be appreciated, thanks guys.

Remove["Global`*"]
(*Gravitational Constant*)
G = 6.672*10^-11;
(*Simulation running time*)
tmax = (800) (86400);
(*Spacecraft TOF between Earth and Mars*)
dt = (254) (86400);
(*Time at which Mid-course correction is applied*)
dtmcc = dt/2;
(*Mass of Sun, Earth, Mars and spacecraft*)
m[0] = 1.988544*10^30 ;(*Mass of Sun*)
m[1] = 5.97219*10^24 ;(*Mass of Earth*)
m[2] = 6.4185*10^23 ;(*Mass of Mars*)
m[3] = 1000;
(*Planetary radii of Sun, Earth and Mars*)
r[0] = 6.963*10^8 ;(*Mean radius of Sun*) 
r[1] = 6.37101 *10^6;(*Mean radius of Earth*) 
r[2] = 3.3899*10^6 ;(*Mean radius of Mars*)
(*Heliocentric positions of Earth and Mars on 26 November 2011*)
p[1] = 149597870700 {4.416639858432274*10^-1, 
   8.798967504313304*10^-1} ;
p[2] = 149597870700 {-8.159382724017646*10^-1, 1.414986880765974*10^0};
p[3] = {6.6078363676586365`*^10, 1.3162870491119447`*^11}
(*Heliocentric velocities of Earth and Mars on 26 November 2011*)
v[1] = 149597870700/
   86400 {-1.563974110293042*10^-2, 7.690252775639107*10^-3} ;
v[2] = 149597870700/
   86400 {-1.160326991502370*10^-2, -5.778933879736245*10^-3} ;

vLambert[1] = {-8473.224022968232`, -23353.77562466311`};
vLambert[2] = {17642.608238929563`, -11842.147739330356`};
mindr = 4.660029187323423`*10^6 + 1000000;
theta = 2.1223566349593104`;
v2x = 20874.193923950195`;
v2y = -11635.115234375`;

Soln = NDSolve[{
   x[1]''[t] == -((G m[0] x[1][t])/(x[1][t]^2 + y[1][t]^2)^(3/2)),
   y[1]''[t] == -((G m[0] y[1][t])/(x[1][t]^2 + y[1][t]^2)^(3/2)),
   x[2]''[t] == -((G m[0] x[2][t])/(x[2][t]^2 + y[2][t]^2)^(3/2)),
   y[2]''[t] == -((G m[0] y[2][t])/(x[2][t]^2 + y[2][t]^2)^(3/2)),
   x[3]''[t] == -((G m[0] x[3][t])/(x[3][t]^2 + y[3][t]^2)^(3/2))
   - (G m[1] (x[3][t]- x[1][t]))/((x[3][t] - x[1][t])^2 + (y[3][t] - y[1][t])^2)^(3/2)
   - (G m[2] (x[3][t] - x[2][t]))/((x[3][t] - x[2][t])^2 + (y[3][t] - y[2][t])^2)^(3/2),
   y[3]''[t] == -((G m[0] y[3][t])/(x[3][t]^2 + y[3][t]^2)^(3/2))
   - (G m[1] (y[3][t] - y[1][t]))/((x[3][t] - x[1][t])^2 + (y[3][t] - y[1][t])^2)^(3/2) 
   - (G m[2] (y[3][t] - y[2][t]))/((x[3][t] - x[2][t])^2 + (y[3][t] - y[2][t])^2)^(3/2),

   x[1][0] == p[1][[1]], y[1][0] == p[1][[2]], x[2][0] == p[2][[1]], 
   y[2][0] == p[2][[2]], x[3][0] == p[3][[1]], y[3][0] == p[3][[2]], 
   x[1]'[0] == v[1][[1]], y[1]'[0] == v[1][[2]], 
   x[2]'[0] == v[2][[1]], y[2]'[0] == v[2][[2]], 
   x[3]'[0] == v[1][[1]] + 3373.0439329044348`, 
   y[3]'[0] == v[1][[2]] + 10880.11479034743`, 
   WhenEvent[
    t == dtmcc, {x[3]'[t] -> vLambert[1][[1]], 
     y[3]'[t] -> vLambert[1][[2]]}], 
   WhenEvent[
    t == ((253*86400) + 81460), {x[3]'[t] -> 
      v2x - Sqrt[(G m[2])/mindr] Sin[theta], 
     y[3]'[t] -> v2y + Sqrt[(G m[2])/mindr] Cos[theta]}]}, {x[1][
    t], y[1][t], x[2][t], y[2][t], x[3][t], y[3][t]}, {t, 0, tmax}, 
  StartingStepSize -> 0.001, AccuracyGoal -> 15, PrecisionGoal -> 15, 
  Method -> "StiffnessSwitching", MaxSteps -> 10000000]

Orbits = ParametricPlot[{{x[1][t], y[1][t]}, {x[2][t], 
     y[2][t]}, {x[3][t], y[3][t]}} /. Soln, {t, 0, tmax}, 
  AxesLabel -> {x, y}, ImageSize -> Large, 
  PlotRange -> {{-0.25*10^12, 0.25*10^12}, {-0.25*10^12, 0.25*10^12}}]

OrbitAnimation = 
 Animate[ParametricPlot[{{x[1][t], y[1][t]}, {x[2][t], y[2][t]}, 
 {x[3][t], y[3][t]}}/. Soln /. t -> a, {t, 0, a}, 
   AxesLabel -> {x, y}, ImageSize -> Large, PlotRange ->If[t>=253*86400, {{x[3][t]-1000000,x[
   3][t]+1000000},{y[3][t]-1000000,y[3][t]+1000000}},{{-0.25*10^12,
   0.25*10^12},{-0.25*10^12,0.25*10^12}}]], {a, 0, tmax}, 
  AnimationRate -> 1000000]

enter image description here

share|improve this question
2  
Unless I overlooked something, everything is fine, there is just a scoping issue. Symbols t, x, and y are not available where you set the PlotRange option value. a and Soln are, use them like you have done in the first argument of your ParametricPlot. –  Theo Tiger May 20 at 19:26

1 Answer 1

up vote 4 down vote accepted
+50

You need to apply the /. Soln /. t->a rule that you apply inside ParametricPlot to your parameters in PlotRange as well. The error occurred because the expression given to PlotRange didn't evaluate. It would just say x[3][somevalue] because x had not been replaced by the corresponding interpolating function. After that I realized that the list had been wrapped with an extra List in the process of doing those replacements, so I wrapped the expression in First.

OrbitAnimation = 
 Animate[ParametricPlot[{{x[1][t], y[1][t]}, {x[2][t], 
       y[2][t]}, {x[3][t], y[3][t]}} /. Soln /. t -> a, {t, 0, a}, 
   AxesLabel -> {x, y}, ImageSize -> Large, PlotRange -> If[
     a >= 253*86400,
     First[{
         {x[3][t] - 1000000, x[3][t] + 1000000},
         {y[3][t] - 1000000, y[3][t] + 1000000}
         } /. Soln /. t -> a],
     {
      {-0.25*10^12, 0.25*10^12},
      {-0.25*10^12, 0.25*10^12}
      }
     ]
   ], {a, 0, tmax}, AnimationRate -> 1000000]

In response to the comment about transitioning into the new Mars centric view:

switchViewLimit = 253*86400;
transitionPeriod = 2 1000000;
transitionRatio[a_] := Min[(a - switchViewLimit)/transitionPeriod, 1]
OrbitAnimation = Animate[
  ParametricPlot[
   {{x[1][t], y[1][t]}, {x[2][t], y[2][t]}, {x[3][t], y[3][t]}} /. 
     Soln /. t -> a, {t, 0, a},
   AxesLabel -> {x, y},
   ImageSize -> Large,
   PlotRange -> If[
     a >= switchViewLimit,
     First[
      {
         {transitionRatio[a] x[3][t] - 0.25*10^12, 
          transitionRatio[a] x[3][t] + 0.25*10^12},
         {transitionRatio[a] y[3][t] - 0.25*10^12, 
          transitionRatio[a] y[3][t] + 0.25*10^12}
         } /. Soln /. t -> a
      ], {
      {-0.25*10^12, 0.25*10^12},
      {-0.25*10^12, 0.25*10^12}
      }]], {a, 0, tmax}, AnimationRate -> 1000000]
share|improve this answer
2  
+1 for the effort :D I didn't even try to read the OP's code.. :( –  Öskå May 21 at 11:18
    
Thanks very much Pickett. I apologise for the amount of code, but I tried to cut it down as much as possible and that's about 20% of it. –  user7388 May 22 at 21:31
    
@user7388 No worries, it actually wasn't a problem for me because I realized from the beginning where the error was. Theo Tiger did too, which proves your question was manageable after all. –  Pickett May 23 at 0:11
    
@Pickett is there a way to make the zoom transition from the heliocentric view to the areocentric view (mars-centered) instead of it occuring instantaneously? –  user7388 May 23 at 17:17
    
@user7388 I've added a piece of code that addresses this. transitionPeriod determines how fast it transitions to the new view. –  Pickett May 23 at 17:46

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.