I'm plotting a function, which is the result of solving a set of coupled differential equations. I am using the Manipulate
function to allow the variation of 3 defined parameters. I am trying to find suitable parameter values such that the function crosses a certain point, i.e. when the plotted function, q
i.e. Papsis[t]
, is less than or equal to Roche[ρ]
(see code for definition).
After this point, the equations become nonsensical, hence the plots no longer display a valid input. Is there a way of stopping the integration/plotting at the point when the plotted, q
i.e. Papsis[t]
, has value smaller than Roche[ρ]
?
Constants
au = QuantityMagnitude[UnitConvert[Quantity[1, "AstronomicalUnit"], "Meters"]];
c = QuantityMagnitude[UnitConvert[Quantity[1, "SpeedOfLight"], "MetersPerSecond"]];
Qpr = 1;
Lsun = QuantityMagnitude[UnitConvert[Quantity[1, "SolarLuminosity"], "Watts"]];
Rsun = QuantityMagnitude[UnitConvert[Quantity[1, "SolarRadius"], "Meters"]];
Msun = QuantityMagnitude[UnitConvert[Quantity[1, "SolarMass"], "Kilograms"]];
G = QuantityMagnitude[UnitConvert[Quantity[1, "GravitationalConstant"], ("Meters"^2*"Newtons")/"Kilograms"^2]];
year = QuantityMagnitude[UnitConvert[Quantity[1, "Years"], "Seconds"]];
Myr = year*10^6;
Gyr = year*10^9;
Mwd = 0.6*Msun;
Cst = 1.27;
U = 1*10^17;
Functions
L[t_] := (3.26*Lsun*(Mwd/(0.6*Msun)))/(0.1 + t/Myr)^1.18;
Roche[dens_] := (0.65*Cst*Rsun*(Mwd/(0.6*Msun))^(1/3))/(dens/3000)^3^(-1);
Papsis[t_] := a[t]*(1 - e[t]);
Radiative Drag
RDdadtR\[Rho]a = -((3*L[t]*Qpr*(2 + 3*e[t]^2))/(c^2*(16*Pi*\[Rho]*Rast*a[t]*(1 - e[t]^2)^(3/2))));
RDdedtR\[Rho]a = -((15*L[t]*e[t]*Qpr)/(c^2*(32*Pi*Rast*\[Rho]*a[t]^2*Sqrt[1 - e[t]^2])));
RDsolR\[Rho]a = ParametricNDSolveValue[{Derivative[1][a][t] == RDdadtR\[Rho]a, Derivative[1][e][t] == RDdedtR\[Rho]a,
WhenEvent[Evaluate[Roche[\[Rho]] >= Papsis[t]],"StopIntegration"],
a[0] == a0, e[0] == 0.3},
{a, e}, {t, 0, 9*Gyr},
{Rast, \[Rho], a0}];
fRDticks = {{Automatic, Automatic}, {Charting`FindTicks[{0, 1}, {0, 1/Myr}], Automatic}};
Manipulate[Grid[{Style["Radiative Drag Working Plot", Bold], Plot[fun[func, t]/scale[func], {t, 0, 9*Gyr}, FrameTicks -> fRDticks,
Epilog -> {Red, Dashed, InfiniteLine[{{0, Roche[\[Rho]]}, {10, Roche[\[Rho]]}}]}, PlotStyle -> {Directive[Blue, Thickness[0.01]]}, PlotRange -> All],
Style["Compiled Plot", Bold], If[comp === {}, Plot[fun[func, t]/scale[func], {t, 0, 9*Gyr}, FrameTicks -> fRDticks,
Epilog -> {Red, Dashed, InfiniteLine[{{0, Roche[\[Rho]]}, {10, Roche[\[Rho]]}}]}, PlotStyle -> {Directive[Blue, Thickness[0.01]]}, PlotRange -> All],
Plot[comp/scale[func], {t, 0, 9*Gyr}, FrameTicks -> fRDticks, Epilog -> {Red, Dashed, InfiniteLine[{{0, Roche[\[Rho]]}, {10, Roche[\[Rho]]}}]},
PlotStyle -> {Directive[Blue, Thickness[0.01]], Directive[Red, Thickness[0.01]]}, PlotRange -> All]]}], {{func, 1}, {1 -> "a", 2 -> "e", 3 -> "q"}},
{{Rast, 0.05}, 0.0001, 0.1, 0.001, Appearance -> "Labeled"}, {{\[Rho], 2000}, 1000, 7000, 50, Appearance -> "Labeled"},
{{a0, 5, "a0 (au)"}, 0.5, 20, 0.2, Appearance -> "Labeled"}, Button["Append", AppendTo[comp, fun[func, t]]], Button["Reset", comp = {}],
TrackedSymbols -> {func, Rast, \[Rho], a0}, Initialization :> {comp = {}, fun[sel_, t_] := Switch[sel, 1, RDsolR\[Rho]a[Rast, \[Rho], a0*au][[1]][t], 2,
RDsolR\[Rho]a[Rast, \[Rho], a0*au][[2]][t], 3, RDsolR\[Rho]a[Rast, \[Rho], a0*au][[1]][t]*(1 - RDsolR\[Rho]a[Rast, \[Rho], a0*au][[2]][t])],
scale[sel_] := Switch[sel, 1 | 3, au, 2, 1]}]
Any help would be appreciated
WhenEvent[]
? $\endgroup$ – J. M.'s ennui♦ Feb 16 at 17:28Papsis[t]
i.e.q
in the plot, is less than or equal toRoche[\[rho]]
. $\endgroup$ – testing09 Feb 16 at 17:33"StopIntegration"
property ofWhenEvent[]
mentioned in its docs? $\endgroup$ – J. M.'s ennui♦ Feb 16 at 17:34WhenEvent
condition is triggered at $t$=0 for all values throughout the range of theManipulate
parameters. For example, with $\rho$=2000, $a$(0)=5 and $e$(0)=0.3 we get Roche[$\rho$]=6.6x10^8 and Papsis(0)=3.5.FunctionRange
may be useful in choosing values for $\rho$ and $a_0$ that do not trigger theWhenEvent
at $t$=0. $\endgroup$ – LouisB Feb 17 at 10:30