0
$\begingroup$

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

$\endgroup$
7
  • $\begingroup$ "parameter values such that the function crosses a certain point" - have you already looked at WhenEvent[]? $\endgroup$ – J. M.'s torpor Feb 16 at 17:28
  • $\begingroup$ @J.M. yes but it didn't stop the integration. I'm trying to get it to stop when Papsis[t] i.e. q in the plot, is less than or equal to Roche[\[rho]]. $\endgroup$ – testing09 Feb 16 at 17:33
  • $\begingroup$ Did you notice the "StopIntegration" property of WhenEvent[] mentioned in its docs? $\endgroup$ – J. M.'s torpor Feb 16 at 17:34
  • 2
    $\begingroup$ @testing09 Please show us the code you tried using WhenEvent. The combination of WhenEvent and StopIntegration would be the obvious answer here, so it should work, so maybe there's something simple we can fix if you show us the code you tried. $\endgroup$ – MarcoB Feb 16 at 19:09
  • 1
    $\begingroup$ It looks like the WhenEvent condition is triggered at $t$=0 for all values throughout the range of the Manipulate 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 the WhenEvent at $t$=0. $\endgroup$ – LouisB Feb 17 at 10:30

Your Answer

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

Browse other questions tagged or ask your own question.