Using units in NBodySimulation function

I’m trying to use NBodySimulation to run a simulation of objects such as stellar masses.
I’ve tried using the InverseSquare and Newtonian settings on NBodySimulation for this, but InverseSquare doesn’t have units and I can’t find any examples for how to use my own/randomized data for the Newtonian setting.

Here’s an example of a unitless simulation I ran with the InverseSquare setting:

data = NBodySimulation[
  "InverseSquare", {<|"Mass" -> 100, "Position" -> {0, 0}, 
    "Velocity" -> {0, -2/25}|>,
   <|"Mass" -> 2, "Position" -> {8, 0}, "Velocity" -> {0, 4}|>}, 20]
ParametricPlot[Evaluate[data[All, "Position", t]], {t, 0, 20}]

And here’s an example of a simulation I ran with the Newtonian setting, of the sun-earth system:

sun = Entity["Star", "Sun"];
earth = Entity["Planet", "Earth"];
date = DateObject[{2018, 1, 1}, "Day", "Gregorian", -4.`];
sunPosition = 
 Table[Quantity[0, "AstronomicalUnits"], 3]; earthPosition = 
 EntityValue[earth, 
  EntityProperty["Planet", "HelioCoordinates", {"Date" -> date}]];
day = Quantity[1, "Days"];
sunVelocity = sunPosition/day;
earthVelocity = (EntityValue[earth, 
      EntityProperty["Planet", 
       "HelioCoordinates", {"Date" -> date + day}]] - earthPosition)/
   day;
dataplanets = 
 NBodySimulation[
  "Newtonian", <|
   sun -> <|"Position" -> sunPosition, "Velocity" -> sunVelocity|>, 
   earth -> <|"Position" -> earthPosition, 
     "Velocity" -> earthVelocity|>|>, Quantity[6, "Months"]]
tend = dataplanets["SimulationTime"]
ParametricPlot3D[
 Evaluate[dataplanets[sun, "Position", t]], {t, 0, tend}]

To input my own data into the Newtonian setting, I’ve tried setups such as the following, but they’ve only produced completely flat graphs so far (just an x-axis with values such as 1.0 * 10^11), even if I try to imitate the sun-earth system:

data = NBodySimulation[
  "Newtonian", {<|"Mass" -> Quantity[2000000, "Kilograms"], 
    "Position" -> Quantity[{0, 0}, "AstronomicalUnit"], 
    "Velocity" -> Quantity[{2, 0}, "AstronomicalUnit"/"Hours"]|>,
   <|"Mass" -> Quantity[2000000, "Kilograms"], 
    "Position" -> Quantity[{10, 0}, "AstronomicalUnit"], 
    "Velocity" -> Quantity[{-2, 0}, "AstronomicalUnit"/"Hours"]|>}, 20]
ParametricPlot[Evaluate[data[All, "Position", t]], {t, 0, 20}]

data  = NBodySimulation[
  "Newtonian", {<|"Mass" -> Quantity[1.989*10^30, "Kilograms"], 
    "Position" -> Quantity[{0, 0}, "Miles"], 
    "Velocity" -> Quantity[{0, -0.008}, "Miles"/"Days"]|>,
   <|"Mass" -> Quantity[5.973*10^24, "Kilograms"], 
    "Position" -> Quantity[{93*10^6, 0}, "Miles"], 
    "Velocity" -> Quantity[{0, 2792}, "Miles"/"Days"]|>}, 100]
ParametricPlot[Evaluate[data[All, "Position", t]], {t, 0, 100}]

So, does anyone know how I could fix this so I can use my own data with units?