There appears to be a bug with Maple 2023.2 which will remove units from physical constants.

AddConstant(Solar_equatorial_radius, symbol = r[e,Sol], value = 696342., uncertainty = 65., units = km) :

AddConstant(Solar_flattening, symbol = f[Sol], value = 0.000009) :

AddConstant(Solar_polar_radius, symbol = r[p,Sol], derive = -r[e,Sol]*(f[Sol] - 1)) :

AddConstant(Solar_nonradius, symbol = x[Sol], derive = f[Sol]*r[e,Sol]) :

In the example here, "GetUnit(Constant(r[p,Sol]))" will return "1" and "GetUnit(Constant(x[Sol]))" will return "m".

The only workaround is to not save dimensionless quantities as physical constants or otherwise at least some constants derived from a dimensionless constant will also be dimensionless.

Adding "units = 1" to "AddConstant" does not help.