# Selecting initial conditions from a list

Let's define a couple of equations:

Clear["Global*"];
d1 = Sqrt[(x - x1)^2 + y^2];
d2 = Sqrt[(x - x2)^2 + y^2];
R = Sqrt[x^2 + y^2]; U = Sqrt[xp^2 + yp^2]; A = -xp y + x yp;
m1 = 1 - μ;
m2 = μ;
x1 = -μ;
x2 = 1 - μ;
G = 1;
c = 10^4;
μ = 1/2;
ϵ = 1;
ω1 = μ/2 ((1 - μ) - 3);

JC = m1/d1 + m2/d2 + 1/2 (R^2 - U^2) + ϵ/
c^2 (-(A^2/2) - m1^2/(2 d1^2) - (1/d1 + 1/d2) m1 m2 - (m1 m2)/(
d1 d2) - m2^2/(2 d2^2) + R^4/8 - A U^2 - (R^2 U^2)/4 - (3 U^4)/
8 - 3/2 (m1/d1 + m2/d2) (-R^2 + U^2) -
7/2 x ((m1 x1)/d1 + (m2 x2)/d2) +
3/2 ((m1 x1^2)/d1 + (m2 x2^2)/d2) -
1/2 ((m1 x1^2)/d1^3 + (m2 x2^2)/d2^3) y^2 + R^2 ω1);


Then we create a two-dimensional grid of initial conditions as follows

ICs = Flatten[Table[{i, j}, {i, -6, 6, 0.2}, {j, -6, 6, 0.2}], 1];
Ntot = Length[ICs]


and then we use a loop procedure in order to obtain some additional information:

data = (x0 = #[[1]];
y0 = 0;
xp0 = 0;
C0 = #[[2]];
JC0 = JC /. {x -> x0, y -> y0, xp -> xp0};
sol = Solve[JC0 == C0, yp, Reals];
yp0 = Abs[yp /. sol[[1]]];
{x0, yp0, C0}) & /@ ICs;


As we can see our aim is to determine the value of yp0. However for some values of the energy C0 there is no real solution of yp0 (energetically forbidden initial conditions). So I want the following: whenever there is no real solution of yp0 the code should jump this energy level, do not print anything in the data list and proceed to the next energetically allowed set of initial conditions.

• Use Nothing; something like yp0 = If[Length@sol > 0, Abs[yp /. sol[[1]]], Nothing] (not tested). – corey979 Nov 16 '16 at 10:35
• @Unfortunately Nothing is not recognized in v9.0 which I use. – Vaggelis_Z Nov 16 '16 at 10:39
• Unevaluated@Sequence[] should work like Nothing. – corey979 Nov 16 '16 at 10:54
• @corey979 Indeed it works! However when there is no real solution the code prints in the data list only two elements x0 and C0, instead of three. But I want nothing to be printed in this case. So the IF should be somehow expanded also in the printing format. – Vaggelis_Z Nov 16 '16 at 10:58
• It won't be costly to just do Select[data, Length@# == 3 &] after data is completed. – corey979 Nov 16 '16 at 11:14

The output of the Solve determines what should be done:

yp0 = If[Length@sol > 0, Abs[yp /. sol[[1]]], Nothing]


In case of pre-Nothing versions of Mathematica, Unevaluated@Sequence[] works like Nothing:

yp0 = If[Length@sol > 0, Abs[yp /. sol[[1]]], Unevaluated@Sequence[]]


This produces elements of data that can contain three (incuding yp0) and two (only x0 and C0) elements. I guess some kind of If procedure might be constructed to keep only the three-element parts, but it won't be costly in any manner to filter out the unwanted two-element results after the whole data is produced:

data = Select[data, Length@# == 3 &]
`