1
$\begingroup$

I am wanting to pull multiple pieces of data from a while loop and put it in a list. The example I have added below is a table with inputted data that I found. Im pretty sure the syntax is unrelated to Mathematica, but it gives you an idea of how I want to display it once I run all of my simulations. For mine I want it to display the following along the top similar to the example.

  • Planet
  • Mass
  • Perihelion (Distance/ Velocity)
  • Aphelion (Distance/Velocity)
  • Semimajor Axis
  • Sidereal Orbit period
  • Orbit eccentricity
  • Mean orbital velocity

From my simulation, I am obtaining the distance and velocity of the aphelion and the sidereal period from the While loop and I want to extract that data so that I can insert it into a table.

I didn't know how to give just example to where it made sense, but I have uploaded my mathematica sheet as a pdf here and the notebook here.

enter image description here

EDIT

Clear["Global`*"];
mSun = 1.98855*10^30; (* Mass of Sun, kg *) 
G = 6.67384*10^-11;     (*Universal Gravitational Constant,N m^2 kg^-2*)
(* Create list of planet names *)
planets = {"Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"};
(* Create list of planet masses *)
mass = {0.3301, 4.8676, 5.9726, 0.64174, 1898.3, 568.36, 86.816, 102.42};
mass = 10^24*mass;

Print["EARTH"];     (* EARTH *)

vini = {0., 30.29*10^3};     (* Initial Velocity of Earth{x, y} , m/s *)
Print["vini= " , vini];
rini = {147.09*10^9, 0.};     (*  Initial Position of Earth {x, y}, m *)
Print["rini= " , rini];     (* copy rini \[Rule] rc, vini \[Rule] vc, where 'c' means "current" *)

rc = rini; vc = vini;
(* time in s *)
\[CapitalDelta]t = 60.;
t = 0;
f = 0;
i = 0;
KeepGoing = True;
nsp = 10000;

While[KeepGoing,
 {
  t = t + \[CapitalDelta]t;
  i = i + 1;
  fold = f;
  If[Mod[i, nsp] == 0, 
   Print["\nstep = ", i, ", time = ", t, " s = ", t/60., " min = ", t/
    3600., " h = ", t/86400., " days"]];
  (* update rn, where 'n' means "new" *)
  rn = rc + vc*\[CapitalDelta]t;
  If[Mod[i, nsp] == 0, Print["\trn = ", rn, ", |rn| = ", Norm[rn]]];
  (* update v *)
  vn = vc - (G*mSun)/Norm[rc]^3*rc*\[CapitalDelta]t;
  (* print *)
  If[Mod[i, nsp] == 0, Print["\tvn = ", vn, " |vMn| = ", Norm[vn]]];
  (* copy vn \[Rule] vc and rn \[Rule] rc *)
  rc = rn; vc = vn;
  (* calculate true anomaly *)
  polar = ToPolarCoordinates[rc];
  f = polar[[2]]/Degree;
  If[Mod[i, nsp] == 0, Print["\tf (true anomaly) = ", f, "\[Degree]"]];
  If[fold < 0 && f >= 0, {KeepGoing = False; 
    Print["\nNOW WILL STOP - f changed from - to +"]}];
  If[Abs[Abs[f] - 180.] <= 0.005 || ! KeepGoing,
   {
    Print[
     "\n\tTrue anomaly is close to 180\[Degree] or completed full \
orbit, f = ", f, "\[Degree]"];
    Print["\tstep = ", i, ", time = ", t, " s = ", t/60., " min = ", 
     t/3600., " h = ", t/86400., " days"];
    Print["\trn = ", rn, ", |rMn| = ", Norm[rn]*10^-3*10^-6, 
     "\[Times]\!\(\*SuperscriptBox[\(10\), \(6\)]\) km"];
    Print["\tvn = ", vn, " |vMn| = ", Norm[vn]*10^-3, " km/s"]
    }];


  }
 ]
$\endgroup$
  • 2
    $\begingroup$ Take a look at Sow/Reap. Thanks for providing a background but could you create a toy example to mimic your procedure and the result you would like to get? $\endgroup$ – Kuba Apr 17 '18 at 7:55
  • $\begingroup$ Welcome! To make the most of Mma.SE start by taking the tour now. It will help us to help you if you write an excellent question. Edit if improvable, show due diligence, give brief context, include minimal working example of code and data in formatted form. As you receive give back, vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. $\endgroup$ – rhermans Jul 7 '18 at 10:29
1
$\begingroup$

Your main problem is that you have not migrated to Mathematica's programming paradigm. Your "C" programming style is not the best way to do things in Mathematica, you should read What are the most common pitfalls awaiting new users? to help you change paradigm.

Some examples of things that could have been done differently in your code:

NDSolve

The integration of the equations of motion you doesn't need to be done by hand, you could use NDSolve. If you need to stop the integration at specic events you can use "EventLocator" Method for NDSolve

Constructing Lists

At any point of your iteration you could AppendTo a list or use Sow and Reap. You could also use Table for your iterations.

Using Print into a notebook make little sense in Mathematica.

PlanetData

You can take advantage of Mathematica curated data

Planet names:

PlanetData[All, "Name"]
(* {"Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"} *)

Planet mass:

QuantityMagnitude@
 UnitConvert[PlanetData[PlanetData[], "Mass"], "Kilograms"]
(* {3.30104*10^23, 4.86732*10^24, 5.9721986*10^24, 6.41693*10^23,
  1.89813*10^27, 5.68319*10^26, 8.68103*10^25, 1.02410*10^26} *)
$\endgroup$
1
$\begingroup$

For an iterated function system like this, Nest and relatives (in particular NestListWhile here) are a fine toolkit. First, write your time step as a function of the state variables, yielding a new list of state variables. Here's a simple example, constant 1D acceleration. Here, t is time, h is height, v is velocity.

timeStep[{t_, h_, v_}] := {t + deltat, h + deltat (v + a deltat/2 ), v + a deltat}

To decide when to end, we need a test function.

test[{t_, h_, v_}] := h < hlimit

For injecting constant parameters, Block is useful for avoiding global definitions. TableForm is good for formatting output.

Block[{a = 30., deltat = 1., hlimit = 100000.},
  NestWhileList[timeStep, {0., 0., 0.}, test]] // TableForm

Yielding:

0.  0.      0.
1.  15.     30.
2.  60.     60.
3.  135.    90.
...
81. 98415.  2430.
82. 100860. 2460.
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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