## Define Function(s) to Retrieve System of Equations

```
Clear(P,V,n,R,T);
Rval=QuantityMagnitude@UnitConvert@Quantity(1, "MolarGasConstant");
idealGasEqn := Module({R=Rval,eqns}, eqns = {P*V == n*R*T})
```

## Known Variables

### Case 1: P, V, and n are knowns (Solve for T)

```
Pval1 = Quantity(1.5, "Atmospheres");
Vval1 = Quantity(3, "Liters");
nval1 = Quantity(1, "Moles");
```

### Case 2: V, T, and n are knowns (Solve for P)

```
Vval2 = Quantity(3, "Liters");
nval2 = Quantity(1, "Moles");
Tval2 = Quantity(55,"Kelvins");
```

## Procedure

### Setup

Equations, solve variables, and inputs

- get system of equations based on input argument (e.g. type = “IdealGas”) using a Switch statement.
- define list of solve variables (Symbols that are left unassigned)
- define list of input variables (mixture of unassigned and assigned)

Units

- Get output unit and SI unit Quantities both with magnitude 1
- find positions in solve variable and input lists based on variable type (Symbol or Quantity) using Position
- replace quantities with magnitude of SI-converted quantities

### Solve

- solve for unknowns using SI magnitudes, no units in output (i.e. unitless) using Solve
- attach SI magnitudes to unitless solution and convert to output units

### Output

- output a unitless or unit-containing solution

## Module

```
idealGasSolver(P1_,V1_,n1_,T1_,type_:"IdealGas",unitlessQ_:False) :=
Module(
{eqns,vars},
(*get system of equations*)
eqns = Switch(type,"IdealGas",idealGasEqn);
vars = {P,V,n,T}; (*Symbols for solve, keep unassigned throughout*)
valsTmp = {P1,V1,n1,T1}; (*input values, some are Symbols, some are Quantities*)
(*units with magnitude 1*)
outUnits = Quantity(1,#)&/@{"Atmospheres","Liters","Moles","DegreesCelsius"};
SIunits = Quantity(1,#)&/@QuantityUnit@UnitConvert@outUnits;
(*find positions based on variable type*)
getIDs(head_) := Position(Head@#===head&/@valsTmp,True)//Flatten;
quantityIDs = getIDs(Quantity);
symbolIDs = getIDs(Symbol);
(*replace quantities with magnitude of SI - converted quantities*)
rules1 = MapThread(#1->#2&,
{quantityIDs,QuantityMagnitude@UnitConvert@valsTmp((quantityIDs))});
vals = ReplacePart(valsTmp,rules1);
(*solve for unknowns using SI magnitudes, no units in output*)
rules2 = MapThread(#1->#2&,{vars((quantityIDs)),vals((quantityIDs))});
unitlessSoln = Solve(eqns/.rules2,vars((symbolIDs)))((1));
(*convert solution to output units and include units*)
rules3 = MapThread(#1 -> #2 &, {vars((symbolIDs)),
vals((symbolIDs))*SIunits((symbolIDs))});
outVals = MapThread(UnitConvert(#1,#2)&,
{vars((symbolIDs))/.rules3/.unitlessSoln,outUnits((symbolIDs))});
unitSoln = MapThread(#1->#2&,{vals((symbolIDs)),outVals});
(*output a solution based on unitlessQ argument*)
outsoln = If(unitlessQ,unitlessSoln,unitSoln)
)
```

## Case 1

Make sure that temperature input is unassigned and specifically use the variable “T”.

```
Clear(T);
idealGasSolver(Pval1, Vval1, nval1, T) (*output in units based on outUnits (deg C)*)
idealGasSolver(Pval1, Vval1, nval1, T, "IdealGas", True) (*output temperature SI unit (K) magnitude*)
```

`{T -> Quantity(-218.31031631383098, "DegreesCelsius")}`

`{T -> 54.83968368616898}`

We get units with the first output, and an SI magnitude with the second.

## Case 2

Make sure that pressure input is unassigned and specifically use the variable “P”.

```
Clear(P);
idealGasSolver(P, Vval2, nval2, Tval2) (*output in units based on outUnits (atm)*)
idealGasSolver(P, Vval2, nval2, Tval2, "IdealGas", True) (*output pressure SI unit (Pa) magnitude*)
```

`{P -> Quantity(2286477219992141/1519875000000000, "Atmospheres")}`

`{P -> 2286477219992141/15000000000}`

Exact arithmetic is preserved in this case.

## Case 3 (additional case, underdetermined system of equations)

Make sure that all inputs except temperature are unassigned and specifically use the variables “P”, “V”, and “n”.

```
Clear(P, V, n)
idealGasSolver(P, V, n, Tval2) // N
idealGasSolver(P, V, n, Tval2, "IdealGas", True) // N (*output SI magnitude*)
```

`{P -> UnitConvert(P*Quantity(1., "Kilograms"/("Meters"*"Seconds"^2)), Quantity(1., "Atmospheres")), V -> UnitConvert(V*Quantity(1., "Meters"^3), Quantity(1., "Liters")), n -> UnitConvert(P*V*Quantity(0.002186770091685928, "Moles"), Quantity(1., "Moles"))}`

`{n -> 0.002186770091685928*P*V}`

The second output (SI magnitude) is more parsable and less subject to issues if you were to apply this process successively (i.e. use the outputs as inputs to the next system of equations).