0
$\begingroup$

I wrote the following true anomaly function based on a Fortran program that I found in a computational space flight dynamics book. The function uses an iterative method to calculate the true anomaly of and object based on its orbital period and eccentricity.

The real problem is that the function stops working after producing 13,923 data points -- meaning that the argument simLength cannot be bigger that 13923. It doesn't crash, it just displays the "Running..." message indefinitely. Is it a problem with memory? Or is it a problem with Mathematica's loop functions?

I am using version 10.3.0 and I first noticed the problem by when trying run [3*58944, 58944, 0.8]. I found the exact limit by having the data points exported to a text file and then importing it back and checking the length of the list.

The number of data points that the function should return is determined by the input "simLength" and the local variable "dT" at "kEnd = Round[simLength/dT]". So if dT is set to 10, and want the simLength to be 10,000 I should get 1000 data points in return.

trueAnomaly[simLength_, period_, e_] := 
  Module[
    {tp = 0., endPoint = 1.*10^-5, dT = 10., k = 1., kEnd, theta, τ, dE, e0, e1, t},
    kEnd = Round[simLength/dT];
    theta = Table[0, kEnd];
    t = Table[0, kEnd];
    For[k, k <= kEnd, k++,
      t[[k]] = k dT;
      τ = 2 π (t[[k]] - tp)/period;
      dE = .1;
      e0 = .1;
      While[dE > endPoint,
        e1 = e0 - (e0 - e Sin[e0] - τ)/(1. - e Cos[e0]);
        dE = Abs[e1 - e0];
        e0 = e1;];
      theta[[k]] = 2.*ArcTan[Sqrt[(1 + e)/(1 - e)]*Tan[e0/2.]];
      If[theta[[k]] < 0, theta[[k]] = 2.*π + theta[[k]]];];
    Return[{t, theta}];]
$\endgroup$

closed as off-topic by Simon Woods, MarcoB, Bob Hanlon, m_goldberg, user9660 Jan 5 '16 at 6:13

  • The question does not concern the technical computing software Mathematica by Wolfram Research. Please see the help center to find out about the topics that can be asked here.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Michael E2 Jan 3 '16 at 6:11
  • 1
    $\begingroup$ You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this this meta Q&A helpful. (You have formatted the code as a "quote.") $\endgroup$ – Michael E2 Jan 3 '16 at 6:12
  • 1
    $\begingroup$ I can not reproduce your problem. I ran data = trueAnomaly[150000, 365., .01]; without any problem. I am using V10.3 on OS X 10.10.2. $\endgroup$ – m_goldberg Jan 3 '16 at 7:33
  • 1
    $\begingroup$ This is not a problem with memory or with Mathematica's loop functions. The problem is that your code contains a While loop which terminates only when dE <= endPoint, a condition which is never true for the k=13924 point in trueAnomaly[3*58944, 58944, 0.8]. At this point instead of converging to a single fixed value, e1 stabilises on a cycle of 3 values. $\endgroup$ – Simon Woods Jan 4 '16 at 21:49
  • 4
    $\begingroup$ I'm voting to close this question as off-topic because the code implements an infinite loop and therefore running indefinitely is the correct behaviour. $\endgroup$ – Simon Woods Jan 4 '16 at 21:55
1
$\begingroup$

This is not an answer, but and extended comment.

I ran your code with some minor modifications without any problems. Here is the code I ran:

trueAnomaly[simLength_, period_, e_] := 
  Module[
     {tp = 0., endPoint = 1.*10^-5, dT = 10., k, kEnd, theta, τ, dE, e0, e1, t, a},
    a = Sqrt[(1. + e)/(1. - e)]; (* only compute this once *)
    kEnd = Round[simLength/dT];
    t = Table[0, kEnd];
    theta = Table[0, kEnd];
    For[k = 1, k <= kEnd, k++,
      t[[k]] = k dT;
      τ = 2 π (t[[k]] - tp)/period;
      dE = .1; e0 = .1;
      While[dE > endPoint,
        e1 = e0 - (e0 - e Sin[e0] - τ)/(1. - e Cos[e0]);
        dE = Abs[e1 - e0];
        e0 = e1];
      theta[[k]] = 2.*ArcTan[a*Tan[e0/2.]];
      If[theta[[k]] < 0, theta[[k]] = 2.*π + theta[[k]]]];
    {t, theta}] (* no need for Return *)

I was able to evaluate

data = trueAnomaly[150000, 365., .01];
$\endgroup$
  • $\begingroup$ Hey, I tried using this code and noticed that since the variable "k" isn't given a starting value, the For loop doesn't actually run, so all that is returned are 2 lists of zeros haha. I set k = 1 and ran into the same problem. But thanks for the tip on calculating "a" first. $\endgroup$ – astroDynamic Jan 4 '16 at 20:08
  • $\begingroup$ @astroDynamic. Thanks for pointing the typo in my code and also for the adding additional information to your post. Using your new info, I found the problem is not in the number of data points you used, but in the value for e. The algorithm appears to fail for e >= .5. Unfortunately, I don't have time to investigate further, but my finding should give you some new ideas on where to look for what's going wrong. $\endgroup$ – m_goldberg Jan 4 '16 at 22:08

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