From 593200ab1fa4d14e22147d945730fb21f49ce8ce Mon Sep 17 00:00:00 2001 From: RoryBarnes Date: Fri, 23 Jan 2026 12:12:19 -0800 Subject: [PATCH] Add unit conversion functions from v3.0 branch Port 16 unit function pairs (fd* and fs*) from the v3.0 branch to enable proper unit handling for quantities like velocity, pressure, mass rate, acceleration, and thermal properties. Also add Python unit aliases for mass rate variants (Earth Masses/Myr, EarthMasses/Myr, Mearth/Myr). This resolves vplot test failures caused by undefined unit strings. Co-Authored-By: Claude Opus 4.5 --- src/control.c | 289 ++++++++++++++++++++++++++++++++++++++++ src/control.h | 33 ++++- vplanet/custom_units.py | 3 + 3 files changed, 324 insertions(+), 1 deletion(-) diff --git a/src/control.c b/src/control.c index e7a7e01e9..78af9869a 100644 --- a/src/control.c +++ b/src/control.c @@ -1265,6 +1265,295 @@ void fsUnitsTempRate(int iType, char **cUnit) { } } +/* Return proper velocity conversion (numeric version) */ +double fdUnitsVel(UNITS *units) { + double dConversion = + fdUnitsLength(units->iLength) / fdUnitsTime(units->iTime); + return dConversion; +} + +/* Particle flux units */ +void fsUnitsParticleFlux(UNITS *units, char **cUnit) { + char *cUnitLength = NULL, *cUnitTime = NULL; + + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsTime(units->iTime, &cUnitTime); + fvFormattedString(cUnit, "%s^-2 %s^-1", cUnitLength, cUnitTime); + free(cUnitLength); + free(cUnitTime); +} + +double fdUnitsParticleFlux(UNITS *units) { + double dTmp; + + dTmp = 1. / (fdUnitsLength(units->iLength) * fdUnitsLength(units->iLength) * + fdUnitsTime(units->iTime)); + return dTmp; +} + +/* Volume rate (meter cubed per second) units */ +void fsUnitsMeterCubedPerSecond(UNITS *units, char **cUnit) { + char *cUnitLength = NULL, *cUnitTime = NULL; + + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsTime(units->iTime, &cUnitTime); + fvFormattedString(cUnit, "%s^3/%s", cUnitLength, cUnitTime); + free(cUnitLength); + free(cUnitTime); +} + +double fdUnitsMeterCubedPerSecond(UNITS *units) { + double dConversion = + pow(fdUnitsLength(units->iLength), 3) / fdUnitsTime(units->iTime); + return dConversion; +} + +/* Pressure units */ +double fdUnitsPressure(UNITS *units) { + double dConversion = fdUnitsMass(units->iMass) / + fdUnitsLength(units->iLength) / + (fdUnitsTime(units->iTime) * fdUnitsTime(units->iTime)); + return dConversion; +} + +void fsUnitsPressure(UNITS *units, char **cUnit) { + char *cUnitLength = NULL, *cUnitMass = NULL, *cUnitTime = NULL; + + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsMass(units->iMass, &cUnitMass); + fsUnitsTime(units->iTime, &cUnitTime); + fvFormattedString(cUnit, "%s/%s/%s^2", cUnitMass, cUnitLength, cUnitTime); + free(cUnitLength); + free(cUnitMass); + free(cUnitTime); +} + +/* Charge units */ +void fsUnitsCharge(char **cUnit) { + fvFormattedString(cUnit, "C"); +} + +double fdUnitsCharge(void) { + return 1; +} + +/* Mass rate units */ +void fsUnitsMassRate(UNITS *units, char **cUnit) { + char *cUnitMass = NULL, *cUnitTime = NULL; + + fsUnitsMass(units->iMass, &cUnitMass); + fsUnitsTime(units->iTime, &cUnitTime); + fvFormattedString(cUnit, "%s/%s", cUnitMass, cUnitTime); + free(cUnitTime); + free(cUnitMass); +} + +double fdUnitsMassRate(UNITS *units) { + double dConversion = fdUnitsMass(units->iMass) / fdUnitsTime(units->iTime); + return dConversion; +} + +/* Length per temperature units */ +void fsUnitsLengthPerTemperature(UNITS *units, char **cUnit) { + char *cUnitLength = NULL, *cUnitTemp = NULL; + + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsTemp(units->iTemp, &cUnitTemp); + fvFormattedString(cUnit, "%s/%s", cUnitLength, cUnitTemp); + free(cUnitTemp); + free(cUnitLength); +} + +double fdUnitsLengthPerTemperature(UNITS *units, double dTemp) { + double dTempFactor; + + if (units->iTemp == U_CELSIUS || units->iTemp == U_KELVIN) { + dTempFactor = 1; + } else if (units->iTemp == U_FARENHEIT) { + dTempFactor = 5. / 9; + } else { + fprintf(stderr, "Unknown temperature unit.\n"); + exit(EXIT_UNITS); + } + double dConversion = fdUnitsLength(units->iLength) / dTempFactor; + return dConversion; +} + +/* Force units */ +void fsUnitsForce(UNITS *units, char **cUnit) { + char *cUnitLength = NULL, *cUnitMass = NULL, *cUnitTime = NULL; + + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsMass(units->iMass, &cUnitMass); + fsUnitsTime(units->iTime, &cUnitTime); + fvFormattedString(cUnit, "%s*%s/%s^2", cUnitMass, cUnitLength, cUnitTime); + free(cUnitMass); + free(cUnitTime); + free(cUnitLength); +} + +double fdUnitsForce(UNITS *units) { + double dConversion = fdUnitsLength(units->iLength) * + fdUnitsMass(units->iMass) / + (fdUnitsTime(units->iTime) * fdUnitsTime(units->iTime)); + return dConversion; +} + +/* Power per length per temperature (thermal conductivity) units */ +void fsUnitsPowerPerLengthPerTemperature(UNITS *units, char **cUnit) { + char *cUnitLength = NULL, *cUnitPower = NULL, *cUnitTemp = NULL; + + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsPower(units, &cUnitPower); + fsUnitsTemp(units->iTemp, &cUnitTemp); + fvFormattedString(cUnit, "%s/%s/%s", cUnitPower, cUnitLength, cUnitTemp); + free(cUnitPower); + free(cUnitTemp); + free(cUnitLength); +} + +double fdUnitsPowerPerLengthPerTemperature(UNITS *units, double dTemp) { + double dConversion = + fdUnitsPower(units->iTime, units->iMass, units->iLength) / + fdUnitsLength(units->iLength) / + fdUnitsTemp(dTemp, U_KELVIN, units->iTemp); + return dConversion; +} + +/* Length squared per time cubed units */ +void fsUnitsLengthSquaredPerTimeCubed(UNITS *units, char **cUnit) { + char *cUnitLength = NULL, *cUnitTime = NULL; + + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsTime(units->iTime, &cUnitTime); + fvFormattedString(cUnit, "%s^2/%s^3", cUnitLength, cUnitTime); + free(cUnitTime); + free(cUnitLength); +} + +double fdUnitsLengthSquaredPerTimeCubed(UNITS *units) { + double dConversion = fdUnitsLength(units->iLength) * + fdUnitsLength(units->iLength) / + pow(fdUnitsTime(units->iTime), 3); + return dConversion; +} + +/* Acceleration units */ +void fsUnitsAcceleration(UNITS *units, char **cUnit) { + char *cUnitLength = NULL, *cUnitTime = NULL; + + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsTime(units->iTime, &cUnitTime); + fvFormattedString(cUnit, "%s/%s^2", cUnitLength, cUnitTime); + free(cUnitTime); + free(cUnitLength); +} + +double fdUnitsAcceleration(UNITS *units) { + double dConversion = fdUnitsLength(units->iLength) / + (fdUnitsTime(units->iTime) * fdUnitsTime(units->iTime)); + return dConversion; +} + +/* Charge per time per length squared (current density) units */ +void fsUnitsChargePerTimePerLengthSquared(UNITS *units, char **cUnit) { + char *cUnitCharge = NULL, *cUnitLength = NULL, *cUnitTime = NULL; + + fsUnitsCharge(&cUnitCharge); + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsTime(units->iTime, &cUnitTime); + fvFormattedString(cUnit, "%s/%s/%s^2", cUnitCharge, cUnitLength, cUnitTime); + free(cUnitTime); + free(cUnitLength); + free(cUnitCharge); +} + +double fdUnitsChargePerTimePerLengthSquared(UNITS *units) { + double dConversion = fdUnitsCharge() / fdUnitsLength(units->iLength) / + (fdUnitsTime(units->iTime) * fdUnitsTime(units->iTime)); + return dConversion; +} + +/* Energy per temperature per mass (specific heat) units */ +void fsUnitsEnergyPerTemperaturePerMass(UNITS *units, char **cUnit) { + char *cUnitMass = NULL, *cUnitEnergy = NULL, *cUnitTemp = NULL; + + fsUnitsMass(units->iMass, &cUnitMass); + fsUnitsEnergy(units, &cUnitEnergy); + fsUnitsTemp(units->iTemp, &cUnitTemp); + fvFormattedString(cUnit, "%s/%s/%s", cUnitEnergy, cUnitTemp, cUnitMass); + free(cUnitEnergy); + free(cUnitTemp); + free(cUnitMass); +} + +double fdUnitsEnergyPerTemperaturePerMass(UNITS *units, double dTemp) { + double dConversion = + fdUnitsEnergy(units->iTime, units->iMass, units->iLength) / + fdUnitsMass(units->iMass) / fdUnitsTemp(dTemp, U_KELVIN, units->iTemp); + return dConversion; +} + +/* Angular momentum units (expanded version) */ +void fsUnitsAngularMomentum(UNITS *units, char **cUnit) { + char *cUnitMass = NULL, *cUnitLength = NULL, *cUnitTime = NULL; + + fsUnitsMass(units->iMass, &cUnitMass); + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsTime(units->iTime, &cUnitTime); + fvFormattedString(cUnit, "%s*%s^2/%s", cUnitMass, cUnitLength, cUnitTime); + free(cUnitMass); + free(cUnitLength); + free(cUnitTime); +} + +double fdUnitsAngularMomentum(UNITS *units) { + double dConversion = + fdUnitsMass(units->iMass) * fdUnitsLength(units->iLength) * + fdUnitsLength(units->iLength) / fdUnitsTime(units->iTime); + return dConversion; +} + +/* Mass per area (surface density) units */ +void fsUnitsMassPerArea(UNITS *units, char **cUnit) { + char *cUnitLength = NULL, *cUnitMass = NULL; + + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsMass(units->iMass, &cUnitMass); + fvFormattedString(cUnit, "%s/%s^2", cUnitMass, cUnitLength); + free(cUnitMass); + free(cUnitLength); +} + +double fdUnitsMassPerArea(UNITS *units) { + double dConversion = + fdUnitsMass(units->iMass) / + (fdUnitsLength(units->iLength) * fdUnitsLength(units->iLength)); + return dConversion; +} + +/* Energy flux per temperature units */ +void fsUnitsEnergyFluxPerTemperature(UNITS *units, char **cUnit) { + char *cUnitPower = NULL, *cUnitLength = NULL, *cUnitTemp = NULL; + + fsUnitsPower(units, &cUnitPower); + fsUnitsLength(units->iLength, &cUnitLength); + fsUnitsTemp(units->iTemp, &cUnitTemp); + fvFormattedString(cUnit, "%s/%s^2/%s", cUnitPower, cUnitLength, cUnitTemp); + free(cUnitPower); + free(cUnitLength); + free(cUnitTemp); +} + +double fdUnitsEnergyFluxPerTemperature(UNITS *units, int iOldTemp, int iNewTemp, + double dTemp) { + double dConversion = + fdUnitsPower(units->iTime, units->iMass, units->iLength) / + (fdUnitsLength(units->iLength) * fdUnitsLength(units->iLength)) / + fdUnitsTemp(dTemp, iOldTemp, iNewTemp); + return dConversion; +} + /* * FILES struct functions */ diff --git a/src/control.h b/src/control.h index eb5fc2713..a45cd4d0a 100644 --- a/src/control.h +++ b/src/control.h @@ -17,7 +17,7 @@ void InitializeControl(CONTROL *, MODULE *); void InitializeControlEvolve(BODY *, CONTROL *, MODULE *, UPDATE *); void InitializeControlVerifyProperty(CONTROL *); -void InitializeFiles(FILES *, OPTIONS *, char *,char **, int); +void InitializeFiles(FILES *, OPTIONS *, char *, char **, int); void WriteHelpOption(OPTIONS *, int); void WriteHelpOutput(OUTPUT *, int); @@ -38,6 +38,22 @@ double fdUnitsAngle(int); double fdUnitsPower(int, int, int); double fdUnitsEnergy(int, int, int); double fdUnitsEnergyFlux(int, int, int); +double fdUnitsVel(UNITS *); +double fdUnitsPressure(UNITS *); +double fdUnitsMeterCubedPerSecond(UNITS *); +double fdUnitsMassRate(UNITS *); +double fdUnitsLengthPerTemperature(UNITS *, double); +double fdUnitsForce(UNITS *); +double fdUnitsPowerPerLengthPerTemperature(UNITS *, double); +double fdUnitsLengthSquaredPerTimeCubed(UNITS *); +double fdUnitsAcceleration(UNITS *); +double fdUnitsChargePerTimePerLengthSquared(UNITS *); +double fdUnitsEnergyPerTemperaturePerMass(UNITS *, double); +double fdUnitsAngularMomentum(UNITS *); +double fdUnitsMassPerArea(UNITS *); +double fdUnitsEnergyFluxPerTemperature(UNITS *, int, int, double); +double fdUnitsParticleFlux(UNITS *); +double fdUnitsCharge(void); void fsUnitsRateSquared(int, char **); // double fdUnitsRate(int); @@ -56,6 +72,21 @@ void fsUnitsPower(UNITS *, char **); void fsUnitsEnergy(UNITS *, char **); void fsUnitsEnergyFlux(UNITS *, char **); void fsUnitsViscosity(UNITS *, char **); +void fsUnitsPressure(UNITS *, char **); +void fsUnitsMeterCubedPerSecond(UNITS *, char **); +void fsUnitsMassRate(UNITS *, char **); +void fsUnitsLengthPerTemperature(UNITS *, char **); +void fsUnitsForce(UNITS *, char **); +void fsUnitsPowerPerLengthPerTemperature(UNITS *, char **); +void fsUnitsLengthSquaredPerTimeCubed(UNITS *, char **); +void fsUnitsAcceleration(UNITS *, char **); +void fsUnitsChargePerTimePerLengthSquared(UNITS *, char **); +void fsUnitsEnergyPerTemperaturePerMass(UNITS *, char **); +void fsUnitsAngularMomentum(UNITS *, char **); +void fsUnitsMassPerArea(UNITS *, char **); +void fsUnitsEnergyFluxPerTemperature(UNITS *, char **); +void fsUnitsParticleFlux(UNITS *, char **); +void fsUnitsCharge(char **); void InfileCopy(INFILE *, INFILE *); diff --git a/vplanet/custom_units.py b/vplanet/custom_units.py index 4b59c22b0..dcf47ed34 100644 --- a/vplanet/custom_units.py +++ b/vplanet/custom_units.py @@ -21,6 +21,9 @@ u.def_unit("Earthradii", u.Rearth), u.def_unit("Earth Masses", u.Mearth), u.def_unit("Degrees", u.deg), + u.def_unit("Earth Masses/Myr", u.Mearth / u.Myr), + u.def_unit("EarthMasses/Myr", u.Mearth / u.Myr), + u.def_unit("Mearth/Myr", u.Mearth / u.Myr), # # Non-standard quantities we'll interpret as unitless: #