Question: Implementing Modulo 2Pi

Hello,

My code records the values I need, however, I need to implement a modulo of 2*Pi on my result for theta. But this leads to a graph with no plots and I'm not sure how to fix it. Any help is greatly aprreciated! Thank you in advance!

Kind regards,

Gam

with(plots):

a := 1.501*10^9:

Th := sqrt(4*Pi^2*a^3/(G*(Mh+Msat))):

HyperionOrbit := proc (`θIC`, `ωIC`, n) local a, Mh, Msat, G, e, beta, M, Eqns, ICs, soln; option remember; global `ωH`, Th; a := 1.501*10^9; Mh := 5.5855*10^18; Msat := 5.6832*10^26; G := 6.67259/10^11; e := .232; beta := .89; M := Mh+Msat; Eqns := diff(theta(t), t) = omega(t), diff(omega(t), t) = -G*Msat*beta^2*(xH(t)*sin(theta(t))-yH(t)*cos(theta(t)))*(xH(t)*cos(theta(t))+yH(t)*sin(theta(t)))/(xH(t)^2+yH(t)^2)^2.5, diff(xH(t), t) = vxH(t), diff(vxH(t), t) = -G*M*xH(t)/(xH(t)^2+yH(t)^2)^(3/2), diff(yH(t), t) = vyH(t), diff(vyH(t), t) = -G*M*yH(t)/(xH(t)^2+yH(t)^2)^(3/2); ICs := xH(0) = a*(1+e), yH(0) = 0, vxH(0) = 0, vyH(0) = sqrt(G*M*(1-e)/(a*(1+e))), theta(0) = `θIC`, omega(0) = `ωIC`; soln := dsolve({Eqns, ICs}, numeric, maxfun = 0, output = array([seq(i, i = 0 .. n*Th, Th)])); plots:-odeplot(soln, [modp(theta(t), 2*Pi), omega(t)/`ωH`], 0 .. n*Th, labels = ["θ(t)","ω(t)/ωH"], axes = boxed, style = plottools:-point, size = [.25, .75]) end proc:

plots:-display(HyperionOrbit(.5, 1.8*`ωH`, 10));

Download Poincare_section_Boyd_plot_fixing_theta.mw

bia Man

Please Wait...