ParametricNDSolve Value and Manipulate with 3 parameters

Question

I am trying to solve a set of coupled differential equations, using ParametricNDSolveValue, where there are 3 parameters. The parameters are Rast, \[Rho] and one of the initial condition values, a[0]. I want to plot the result of this integration as a Manipulate plot, where each parameter can be independently varied.

My code is as follows:

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])/(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, a[0] == a0, e[0] == 0.3}, {a[t], e[t]}, {t, 0, 9*Gyr}, 
    {Rast, \[Rho]}];

The parameters are Rast, \[Rho], and the initial condition value, a0. I am unsure how to create a Manipulate plot with all three parameters.

PARAMETER DOMAINS: Rast from 0.001 to 0.01 \[rho] from 1000 to 7000 a0 from 3*au to 20*au and the Plot domain: t=0 to t=9 Gyr

Any help would be greatly appreciated.

Answer 1

$Version

(* "12.2.0 for Mac OS X x86 (64-bit) (December 12, 2020)" *)

Clear["Global`*"]

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 = 6/10*Msun;
Cst = 1.27;
U = 10^17;

Functions

L[t_] := (326/100*Lsun*(Mwd/(6/10*Msun)))/
   (1/10 + t/Myr)^(118/100);
Roche[dens_] := (65/100*Cst*Rsun*(Mwd/(6/10*Msun))^(1/3))/
   (dens/3000)^3^(-1);
Papsis[t_] := a[t]*(1 - e[t]);

Radiative Drag

RDdadtRρa = -((3*L[t]*Qpr*(2 + 3*e[t]^2))/
     (c^2*(16*Pi*ρ*Rast*a[t]*(1 - e[t]^2)^(3/2))));
RDdedtRρa = -((15*L[t]*e[t])/
     (c^2*(32*Pi*Rast*ρ*a[t]^2*Sqrt[1 - e[t]^2])));

Include a0 as a parameter

RDsolRρa = 
  ParametricNDSolveValue[{a'[t] == RDdadtRρa, e'[t] == RDdedtRρa, 
    a[0] == a0, e[0] == 3/10}, {a, e}, 
   {t, 0, 9*Gyr}, {Rast, ρ, a0}];

EDIT: Plot of a rescaled to units of au

Manipulate[
 factor = If[func == 1, 1/au, 1];
 Plot[factor*RDsolRρa[Rast, ρ, a0*au][[func]][t],
  {t, 0, 9 Gyr}],
 {{func, 1}, {1 -> "a", 2 -> "e"}},
 {{Rast, 0.005}, 0.001, 0.01, 0.001, Appearance -> "Labeled"},
 {{ρ, 3000}, 1000, 7000, 50, Appearance -> "Labeled"},
 {{a0, 10, "a0 (au)"}, 3, 20, 0.2, Appearance -> "Labeled"}]