Have you tried:

plotParametric(x1, x2)

The label of the x-axis is wrong in my example, it's written "time" but the values are ok.

I found this command in OpenModelicaUsersGuide.pdf on page 31.  You can download this document from openmodelica.org under Home / User Documentation / Users Guide.

Regards,

Rolf

bboymq wrote:

I have a nonlinear Differential Equation for an Orbit of a Satellite
Consider the problem of an orbit of satellite, whose position and velocity are obtained as the solution of the following state equation:

x'₁(t) = x₃(t)
x'₂(t) = x₄(t)
x'₃(t) = -GM_Ex₁(t) / (x²₁(t) + x²₂(t))^3/2
x'₄(t) = -GM_Ex₂(t) / (x²₁(t) + x²₂(t))^3/2

where G = 6.672 x 10^-11 Nm²/kg² is the gravitational constant, and M_E =  6.97 x 10²⁴ kg is the mass of the earth. Note (x₁, x₂) and (x₃, x₄) denote the position and velocity, respectively, of the satellite on the plane having the earth at its origin.

How can I plot them in openModelica? I have tried, but I can't plot them.

Thanks a lot.

Thanks Rolf, can you check my program

(x10,x20) = (4.223 × 107, 0)[m]and (x30,x40) = (v10,v20) =
(0, 3071)[m/s].

Here is the code in modelica

model A
  constant Real G = 6.672e-11;
  constant Real ME = 5.97e24;
  Real x1(start = 4.223e7);
  Real x2(start = 0);
  Real x3(start = 0);
  Real x4(start = 3071);
equation
  der(x1) = x3;
  der(x2) = x4;
  der(x3) = -G * ME * x1 / ((x1^2 + x2^2)^(3/2));
  der(x4) = -G * ME * x2 / ((x1^2 + x2^2)^(3/2));
end A;

Is this program wrong?

Thanks a lot.

M_E should be 6.97e24, not 5.97e24, the rest of code seems to be ok.

I used x1(start = 6.37E06) and x4(start = 10.4E03)
and 20000 time steps. With these settings I get a perfect ellipse but the simulation takes some time (30s).

I guess that your start values should work too.

oh, thanks Rolf once again for guiding me,

You have got a perfect ellipse and that's correct, but, my program can't show a perfect ellipse, only a part of ellipse, I will show you my result:



Once again, may you show to me your "simulate code", how to get a perfect ellipse?

Thanks Rolf

As you can see, you have simulated 1/4 from an ellipse. So you have to simulate 4 times longer or 5 times longer to be sure that you simulate more than 360 degrees.

rolffankhauser wrote:

As you can see, you have simulated 1/4 from an ellipse. So you have to simulate 4 times longer or 5 times longer to be sure that you simulate more than 360 degrees.

Thanks Rolf, you can guide me concretely because I'm a beginner, How to simulate, I dont know how to simulate 4 times longer or more.

You used a command simulate(... stopTime = 100 ...) it seems. Simply change that to 400.

Here is the code I used:

model orbit
 
  Real x1(start = 6.37E06);
  Real x2(start = 0.0);
  Real x3(start = 0.0);
  Real x4(start = 10.4E03);
  parameter Real m_E = 6.97E24;
  parameter Real G = 6.672E-11;

equation

  der(x1) = x3;
  der(x2) = x4;
  der(x3) = -G*m_E*x1/(x1^2 + x2^2)^1.5;
  der(x4) = -G*m_E*x2/(x1^2 + x2^2)^1.5;

end orbit;

simulate(orbit, startTime = 0.0, stopTime = 20000.0)

plotParametric(x1, x2)

Maybe you have chosen a height and velocity of the satellite which did not result in elliptical curve.

The total energy of the satellite is E = 0.5*m*v^2 - G*M*m/(R + h). It must be negativ for an elliptical curve.

So, velocity v < sqrt(G*M/(R + h)). But if you choose v too small the elliptical curve will lie inside the earth, that means the satellite will crash into the earth.

For example for h = 10 km you get v < 5330 m/s. If you choose v = 5000 m/s you will get an elliptical curve which is near a circle. You can see that if you plot x1 and x2 vs. time ( plot({x1, x2}) )

The theoretical value of the semimajor axis should be a = G*M (R + h) / (2*G*M - v^2*(R + h)) = 14.6 E06. Thus, the width of the ellipse is 29.2 E06 which is close to the value in the plot.

Oh, Rolf, I can't use any word to my show deep gratitude to you. You're very great !!!