Units consistency between model units and real units

AgustinVallejo commented 1 year ago

Original flash sim uses unitless values, nice and round. We decided MSS will use AU, km/s and years. This entails a bunch of potential problems.

A teacher said that their astronomy lab now is broken! Answer: The sim is 10 years old, perhaps that's okay.
Measuring G inside of the sim might yield dubious results (close but not exact due to significant units).
We need to figure out the exact conversion rate between Model units and View units, or change the model units to MKS altogether.

Here's what we wrote to the teacher:

There's multiple items regarding the simulation units change:

As Diana said, we re-adjusted the scales for them to better portray real physical systems, with units 10^28 kg, AU, and km/s. Because of that, the values have re-scaled.

The scaling factors are 0.01 for distance and 0.2109 for velocity. i.e. if an old distance was 100 old units, it's now 1.00AU. Again, the unfortunate scaling factor for velocities is to allow for real G calculations.

Only the 'Sun and Planet' preset had its values changed for design reasons, the rest of them are just re-scaled.

The number display currently shows 1 significant unit, so old values of 75 old distance units now show as 0.7AU. We should probably change that to 2. Maybe add a queryParameter for that preference?

We're deciding if the Center of Mass should be centered and followed by default. If that's the case, presets would show a slight position shift.

samreid commented 1 year ago

From discussion with @AgustinVallejo and @jonathanolson, about possibly using mks.

```diff Subject: [PATCH] Remove unnecessary calls to setScreenSummaryIntroAndTitle, see https://github.com/phetsims/my-solar-system/issues/103 --- Index: main/solar-system-common/js/SolarSystemCommonConstants.ts IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/main/solar-system-common/js/SolarSystemCommonConstants.ts b/main/solar-system-common/js/SolarSystemCommonConstants.ts --- a/main/solar-system-common/js/SolarSystemCommonConstants.ts (revision 5e236f59ebcc772c8754bea259f60b3fecde8a66) +++ b/main/solar-system-common/js/SolarSystemCommonConstants.ts (date 1679005664907) @@ -62,9 +62,9 @@ DEFAULT_SOUND_OUTPUT_LEVEL: 0.1, // Multipliers that modify the numeric value shown in Number Displays - POSITION_MULTIPLIER: 0.01, - VELOCITY_MULTIPLIER: 0.2109, - TIME_MULTIPLIER: 0.218, + POSITION_MULTIPLIER: 1, + VELOCITY_MULTIPLIER: 1, + TIME_MULTIPLIER: 1, MAX_ORBITAL_DIVISIONS: 6, MIN_ORBITAL_DIVISIONS: 2 Index: main/my-solar-system/js/lab/model/LabModel.ts IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/main/my-solar-system/js/lab/model/LabModel.ts b/main/my-solar-system/js/lab/model/LabModel.ts --- a/main/my-solar-system/js/lab/model/LabModel.ts (revision 7d97b291e5e8ae0fa5f233b52929582cb0354c46) +++ b/main/my-solar-system/js/lab/model/LabModel.ts (date 1679007354571) @@ -88,71 +88,86 @@ } public setModesToMap(): void { + + const toMass = ( mass: number ) => mass * 1e28; + const toPosition = v => { + return new Vector2( v.x * 1.496e11 / 100, v.y * 1.496e+11 / 100 ); // 1AU = 1.496e+11m + }; + const toVelocity = v => { + return new Vector2( v.x / 1e9, v.y / 1e9 ); + }; + // F = Gm1m2/r^2 = 6.6743E-11 * t + + const distance = toPosition( new Vector2( 100, 0 ) ).magnitude; + const force = 6.6743E-11 * toMass( 200.0000 ) * toMass( 0.1 ) / Math.pow( distance, 2 ); + console.log( force ); this.modeMap.set( LabMode.SUN_PLANET, [ - { active: true, mass: 200, position: new Vector2( 0, 0 ), velocity: new Vector2( 0, -5 ) }, - { active: true, mass: 10, position: new Vector2( 200, 0 ), velocity: new Vector2( 0, 100 ) } + { active: true, mass: toMass( 200.0000 ), position: toPosition( new Vector2( 0, 0 ) ), velocity: toVelocity( new Vector2( 0, -5 ) ) }, + // { active: true, mass: toMass( 10.0000 ), position: toPosition( new Vector2( 200, 0 ) ), velocity: toVelocity( new Vector2( 0, 100 ) ) } + + { active: true, mass: toMass( 0.1 ), position: toPosition( new Vector2( 100, 0 ) ), velocity: toVelocity( new Vector2( 0, 141.42135 ) ) } ] ); this.modeMap.set( LabMode.SUN_PLANET_MOON, [ - { active: true, mass: 200, position: new Vector2( 0, 0 ), velocity: new Vector2( 0, 0 ) }, - { active: true, mass: 10, position: new Vector2( 160, 0 ), velocity: new Vector2( 0, 120 ) }, - { active: true, mass: 0.000001, position: new Vector2( 140, 0 ), velocity: new Vector2( 0, 53 ) } + { active: true, mass: toMass( 200.0000 ), position: toPosition( new Vector2( 0, 0 ) ), velocity: toVelocity( new Vector2( 0, 0 ) ) }, + { active: true, mass: toMass( 10.0000 ), position: toPosition( new Vector2( 160, 0 ) ), velocity: toVelocity( new Vector2( 0, 120 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( 140, 0 ) ), velocity: toVelocity( new Vector2( 0, 53 ) ) } ] ); this.modeMap.set( LabMode.SUN_PLANET_COMET, [ - { active: true, mass: 200, position: new Vector2( 0, 0 ), velocity: new Vector2( 0, 0 ) }, - { active: true, mass: 1, position: new Vector2( 150, 0 ), velocity: new Vector2( 0, 120 ) }, - { active: true, mass: 0.000001, position: new Vector2( -220, 130 ), velocity: new Vector2( -20, -35 ) } + { active: true, mass: toMass( 200.0000 ), position: toPosition( new Vector2( 0, 0 ) ), velocity: toVelocity( new Vector2( 0, 0 ) ) }, + { active: true, mass: toMass( 1.000 ), position: toPosition( new Vector2( 150, 0 ) ), velocity: toVelocity( new Vector2( 0, 120 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( -220, 130 ) ), velocity: toVelocity( new Vector2( -20, -35 ) ) } ] ); this.modeMap.set( LabMode.TROJAN_ASTEROIDS, [ - { active: true, mass: 200, position: new Vector2( 0, 0 ), velocity: new Vector2( 0, 0 ) }, - { active: true, mass: 5, position: new Vector2( 150, 0 ), velocity: new Vector2( 0, 119 ) }, - { active: true, mass: 0.000001, position: new Vector2( 75, -130 ), velocity: new Vector2( 103, 60 ) }, - { active: true, mass: 0.000001, position: new Vector2( 75, 130 ), velocity: new Vector2( -103, 60 ) } + { active: true, mass: toMass( 200.0000 ), position: toPosition( new Vector2( 0, 0 ) ), velocity: toVelocity( new Vector2( 0, 0 ) ) }, + { active: true, mass: toMass( 5.000 ), position: toPosition( new Vector2( 150, 0 ) ), velocity: toVelocity( new Vector2( 0, 119 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( 75, -130 ) ), velocity: toVelocity( new Vector2( 103, 60 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( 75, 130 ) ), velocity: toVelocity( new Vector2( -103, 60 ) ) } ] ); this.modeMap.set( LabMode.ELLIPSES, [ - { active: true, mass: 250, position: new Vector2( -200, 0 ), velocity: new Vector2( 0, 0 ) }, - { active: true, mass: 0.000001, position: new Vector2( -115, 0 ), velocity: new Vector2( 0, 151 ) }, - { active: true, mass: 0.000001, position: new Vector2( 50, 0 ), velocity: new Vector2( 0, 60 ) }, - { active: true, mass: 0.000001, position: new Vector2( 220, 0 ), velocity: new Vector2( 0, 37 ) } + { active: true, mass: toMass( 250.0000 ), position: toPosition( new Vector2( -200, 0 ) ), velocity: toVelocity( new Vector2( 0, 0 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( -115, 0 ) ), velocity: toVelocity( new Vector2( 0, 151 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( 50, 0 ) ), velocity: toVelocity( new Vector2( 0, 60 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( 220, 0 ) ), velocity: toVelocity( new Vector2( 0, 37 ) ) } ] ); this.modeMap.set( LabMode.HYPERBOLIC, [ - { active: true, mass: 250, position: new Vector2( 0, 25 ), velocity: new Vector2( 0, 0 ) }, - { active: true, mass: 0.000001, position: new Vector2( -250, -70 ), velocity: new Vector2( 120, 0 ) }, - { active: true, mass: 0.000001, position: new Vector2( -250, -140 ), velocity: new Vector2( 120, 0 ) }, - { active: true, mass: 0.000001, position: new Vector2( -250, -210 ), velocity: new Vector2( 120, 0 ) } + { active: true, mass: toMass( 250.0000 ), position: toPosition( new Vector2( 0, 25 ) ), velocity: toVelocity( new Vector2( 0, 0 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( -250, -70 ) ), velocity: toVelocity( new Vector2( 120, 0 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( -250, -140 ) ), velocity: toVelocity( new Vector2( 120, 0 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( -250, -210 ) ), velocity: toVelocity( new Vector2( 120, 0 ) ) } ] ); this.modeMap.set( LabMode.SLINGSHOT, [ - { active: true, mass: 200, position: new Vector2( 1, 0 ), velocity: new Vector2( 0, -1 ) }, - { active: true, mass: 10, position: new Vector2( 131, 55 ), velocity: new Vector2( -55, 115 ) }, - { active: true, mass: 0.000001, position: new Vector2( -6, -128 ), velocity: new Vector2( 83, 0 ) } + { active: true, mass: toMass( 200.0000 ), position: toPosition( new Vector2( 1, 0 ) ), velocity: toVelocity( new Vector2( 0, -1 ) ) }, + { active: true, mass: toMass( 10.0000 ), position: toPosition( new Vector2( 131, 55 ) ), velocity: toVelocity( new Vector2( -55, 115 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( -6, -128 ) ), velocity: toVelocity( new Vector2( 83, 0 ) ) } ] ); this.modeMap.set( LabMode.DOUBLE_SLINGSHOT, [ - { active: true, mass: 200, position: new Vector2( 0, 0 ), velocity: new Vector2( 0, -1 ) }, - { active: true, mass: 5, position: new Vector2( 0, -112 ), velocity: new Vector2( 134, 0 ) }, - { active: true, mass: 5, position: new Vector2( 186, -5 ), velocity: new Vector2( 1, 111 ) }, - { active: true, mass: 0.000001, position: new Vector2( 70, 72 ), velocity: new Vector2( -47, 63 ) } + { active: true, mass: toMass( 200.0000 ), position: toPosition( new Vector2( 0, 0 ) ), velocity: toVelocity( new Vector2( 0, -1 ) ) }, + { active: true, mass: toMass( 5.00000 ), position: toPosition( new Vector2( 0, -112 ) ), velocity: toVelocity( new Vector2( 134, 0 ) ) }, + { active: true, mass: toMass( 5.00000 ), position: toPosition( new Vector2( 186, -5 ) ), velocity: toVelocity( new Vector2( 1, 111 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( 70, 72 ) ), velocity: toVelocity( new Vector2( -47, 63 ) ) } ] ); this.modeMap.set( LabMode.BINARY_STAR_PLANET, [ - { active: true, mass: 150, position: new Vector2( -100, 0 ), velocity: new Vector2( 0, -60 ) }, - { active: true, mass: 120, position: new Vector2( 100, 0 ), velocity: new Vector2( 0, 50 ) }, - { active: true, mass: 0.000001, position: new Vector2( -50, 0 ), velocity: new Vector2( 0, 120 ) } + { active: true, mass: toMass( 150.0000 ), position: toPosition( new Vector2( -100, 0 ) ), velocity: toVelocity( new Vector2( 0, -60 ) ) }, + { active: true, mass: toMass( 120.0000 ), position: toPosition( new Vector2( 100, 0 ) ), velocity: toVelocity( new Vector2( 0, 50 ) ) }, + { active: true, mass: toMass( 0.000001 ), position: toPosition( new Vector2( -50, 0 ) ), velocity: toVelocity( new Vector2( 0, 120 ) ) } ] ); this.modeMap.set( LabMode.FOUR_STAR_BALLET, [ - { active: true, mass: 120, position: new Vector2( -100, 100 ), velocity: new Vector2( -50, -50 ) }, - { active: true, mass: 120, position: new Vector2( 100, 100 ), velocity: new Vector2( -50, 50 ) }, - { active: true, mass: 120, position: new Vector2( 100, -100 ), velocity: new Vector2( 50, 50 ) }, - { active: true, mass: 120, position: new Vector2( -100, -100 ), velocity: new Vector2( 50, -50 ) } + { active: true, mass: toMass( 120.0000 ), position: toPosition( new Vector2( -100, 100 ) ), velocity: toVelocity( new Vector2( -50, -50 ) ) }, + { active: true, mass: toMass( 120.0000 ), position: toPosition( new Vector2( 100, 100 ) ), velocity: toVelocity( new Vector2( -50, 50 ) ) }, + { active: true, mass: toMass( 120.0000 ), position: toPosition( new Vector2( 100, -100 ) ), velocity: toVelocity( new Vector2( 50, 50 ) ) }, + { active: true, mass: toMass( 120.0000 ), position: toPosition( new Vector2( -100, -100 ) ), velocity: toVelocity( new Vector2( 50, -50 ) ) } ] ); this.modeMap.set( LabMode.DOUBLE_DOUBLE, [ - { active: true, mass: 60, position: new Vector2( -115, -3 ), velocity: new Vector2( 0, -154 ) }, - { active: true, mass: 70, position: new Vector2( 102, 0 ), velocity: new Vector2( 1, 150 ) }, - { active: true, mass: 55, position: new Vector2( -77, -2 ), velocity: new Vector2( -1, 42 ) }, - { active: true, mass: 62, position: new Vector2( 135, 0 ), velocity: new Vector2( -1, -52 ) } + { active: true, mass: toMass( 60.0000 ), position: toPosition( new Vector2( -115, -3 ) ), velocity: toVelocity( new Vector2( 0, -154 ) ) }, + { active: true, mass: toMass( 70.0000 ), position: toPosition( new Vector2( 102, 0 ) ), velocity: toVelocity( new Vector2( 1, 150 ) ) }, + { active: true, mass: toMass( 55.0000 ), position: toPosition( new Vector2( -77, -2 ) ), velocity: toVelocity( new Vector2( -1, 42 ) ) }, + { active: true, mass: toMass( 62.0000 ), position: toPosition( new Vector2( 135, 0 ) ), velocity: toVelocity( new Vector2( -1, -52 ) ) } ] ); this.modeMap.set( LabMode.CUSTOM, [ - { active: true, mass: 120, position: new Vector2( -100, 100 ), velocity: new Vector2( -50, -50 ) }, - { active: true, mass: 120, position: new Vector2( 100, 100 ), velocity: new Vector2( -50, 50 ) }, - { active: true, mass: 120, position: new Vector2( 100, -100 ), velocity: new Vector2( 50, 50 ) }, - { active: true, mass: 120, position: new Vector2( -100, -100 ), velocity: new Vector2( 50, -50 ) } + { active: true, mass: toMass( 120.0000 ), position: toPosition( new Vector2( -100, 100 ) ), velocity: toVelocity( new Vector2( -50, -50 ) ) }, + { active: true, mass: toMass( 120.0000 ), position: toPosition( new Vector2( 100, 100 ) ), velocity: toVelocity( new Vector2( -50, 50 ) ) }, + { active: true, mass: toMass( 120.0000 ), position: toPosition( new Vector2( 100, -100 ) ), velocity: toVelocity( new Vector2( 50, 50 ) ) }, + { active: true, mass: toMass( 120.0000 ), position: toPosition( new Vector2( -100, -100 ) ), velocity: toVelocity( new Vector2( 50, -50 ) ) } ] ); } } Index: main/my-solar-system/js/common/model/NumericalEngine.ts IDEA additional info: Subsystem: com.intellij.openapi.diff.impl.patch.CharsetEP <+>UTF-8 =================================================================== diff --git a/main/my-solar-system/js/common/model/NumericalEngine.ts b/main/my-solar-system/js/common/model/NumericalEngine.ts --- a/main/my-solar-system/js/common/model/NumericalEngine.ts (revision 7d97b291e5e8ae0fa5f233b52929582cb0354c46) +++ b/main/my-solar-system/js/common/model/NumericalEngine.ts (date 1679006541997) @@ -24,7 +24,7 @@ const CHI = -0.06626458266981849; // Gravitational constant -const G = 10000; +const G = 6.6743E-11; export default class NumericalEngine extends Engine { @@ -195,6 +195,7 @@ const mass2 = body2.massProperty.value; assert && assert( mass2 > 0, 'mass2 should not be 0' ); const force: Vector2 = this.getForce( body1, body2 ); + console.log( force ); // REVIEW: Can this scratchVector be removed? body1.forceProperty.value = body1.forceProperty.value.plus( scratchVector.set( force ) ); ```

samreid commented 1 year ago

MKS + adapters may be good enough, but brainstorming one more idea--If we get MKS working, maybe it would be straightforward to migrate that consistently/confidently to our own system of units? Not sure that's a great idea but just wanted to write it down.

AgustinVallejo commented 1 year ago

Value	Equation	Model	Conversion	View
M	Given	200 uM	x1e28 kg/uM	1.99e30 kg
D	Given	100 uD	x 0.01 AU/uD	1AU
G	Given	10000
V	sqrt(GM/R)	... uD/uT	-> AU/yr -> km/s	29.8km/s
T	2piR/V	...uT	-> yr / uT	1yr

This table briefly summarizes the reasoning behind the value for each conversion multiplier. G value is given by the model, and suppose a single circular orbit around M=200uM and with radius D=100uD for uM and uD arbitrary model units. Kepler's third law can be used to infer T in uT and circular velocity equation for V in uD/uT. From there, it's a matter of unit conversion:

uD x 0.01 = 1 AU
uT x 0.22... = 1 yr
200 uM x 0.994e28 = 1.99e30kg (1Msun) (This one should be further discussed because we probably want to display 200e28kg)
circular velocity in uD/uT x 0.21... = 29.8 km/s (Earth's circular velocity) (This one was tricky because the conversion is uD/uT -> AU/yr -> km/s )

With all those values, G can be calculated with an error on the 5th significant unit, which incidentally is about the same measurement of uncertainty we have for it (Source).

I also asked for external advice from a professor who does this kinds of simulations and he advised us to keep using model units internally and change them at the end, because these numbers in MKS range several orders of magnitude and it's easy for the software to make mistakes:

1-1e-9**2 = 1 (with double precision)

AgustinVallejo commented 1 year ago

17/03 meeting decision:

[x] Keep implementation as it is.
[x] Dynamicaly calculate the multipliers.
[x] Create unit units tests.
[ ] (Future) Show in PhetIO only Displayed Units (AU,kg,yr).

AgustinVallejo commented 1 year ago

@samreid Do you think the unit tests could/should be done for this version of the sim or that's more of a task to be done on par with the PhetIO implementation?

AgustinVallejo commented 1 year ago

I believe the unit tests don't have to be done for this release.

jonathanolson commented 1 year ago

I'm fine with unit tests being done for the next release.

samreid commented 1 year ago

How about a quick force calculation sanity check and establishing a pattern of how the units are set up and converted?

AgustinVallejo commented 1 year ago

@samreid I don't know if this is what you mean but now the unit conversions are always calculated at the beggining of the sim; no more seemingly arbitrary numeric constants. It's at the top of SolarSystemCommonConstants.ts

samreid commented 1 year ago

@AgustinVallejo and I wrote a harness for the unit tests, but did not enable the first test yet because we got unexpected values.

AgustinVallejo commented 1 year ago

This issue is now done for this release. Units are being calculated and we make sure the sim engine and unit conversions are working the same as SI. Leaving on hold because we might want to look at this in the future PhET-IO release.

samreid commented 1 year ago

I converted the unit tests to TypeScript. I also had a question about this part:

  const mass1SimUnits = mass1 / MASS_MULTIPLIER;
  const mass2SimUnits = mass2 / MASS_MULTIPLIER;
  const distanceSimUnits = distance / METERS_IN_AU / POSITION_MULTIPLIER;

There seems to be an asymmetry in that the MASS_MULTIPLIER converts between kg and uM (where the view unit is 10²⁸kg). But the POSITION_MULTIPLIER converts between m and AU (where the view unit is AU)? I'm not saying there is a problem, but maybe good to document it somewhere? Or maybe it is written somewhere I could not find it?

Maybe add a link from https://github.com/phetsims/solar-system-common/blob/670cc9ea98513bf05573646623904c86e524f908/js/SolarSystemCommonConstants.ts#L13-L20 to the appropriate documentation, or mention the units there?

arouinfar commented 1 year ago

@jonathanolson will take a look at https://github.com/phetsims/my-solar-system/issues/116#issuecomment-1477990090 and determine if anything is needed for this release.

We can defer PhET-iO considerations for now.

jonathanolson commented 1 year ago

Things look good here, the unit test referenced seems to do unit conversion correctly, and the model.md documents the conversions to my satisfaction.

phetsims / my-solar-system

Units consistency between model units and real units #116