(* Mathematica Package *)
(* :Title: NaturalUnits *)
(* :Author: Eemeli Tomberg *)
(* :Summary: This package can be used to convert physical quantities into natural units, defined by setting a group of natural constants to one. By default, it uses conventions common in early universe cosmology, but it also supports custom units and unit systems. *)
(* :Package Version: 2.0 *)
(* :Copyright: Copyright 2022, Eemeli Tomberg. Distribution and modification allowed and encouraged. *)
(* :Context: NaturalUnits` *)
(* :Keywords: unit, conversion, physics *)

BeginPackage["NaturalUnits`"];

Print["NaturalUnits 2.0"];

(************ Interface **************)

(* Functions, options, and other central variables *)
NaturalConvert::usage="NaturalConvert[quantity] converts `quantity' into a default unit under the current natural units.\nNaturalConvert[quantity, targetunit] attempts to convert the specified `quantity' into the specified `targetunit', or a power thereof, under the current system of natural units.\nNaturalConvert[quantity, targetunits] attempts to convert the specified `quantity' into powers of the units listed in `targetunits'.";
ToUnitDefault::usage="ToUnitDefault is an option that gives the unit to which `NaturalConvert' converts quantities if no target unit is specified.";
KeepUnit::usage="KeepUnit is an option to `NaturalConvert' telling whether the target unit should be returned with the result of a unit conversion (`True' or `False'). The special value `Defer' wraps the unit inside Defer[] to keep its product with the numerical coefficient unevaluated. This is useful if the unit itself contains numbers, e.g. 0.05 Mpc^-1."
NoUnit::usage="NoUnit represents no chosen unit, interpreted as the default natural unit or units when needed.";

NaturalUnitsSetup::usage="NaturalUnitsSetup[] sets up the natural unit system based on the options given (see Options[NaturalUnitsSetup]). If no options are given, the previous ones are reused.";
NewUnitSystem::usage="NewUnitSystem is an option to `NaturalUnitsSetup', giving the unit system to use as a list of equations containing units, typically setting a number of physical constants to one. Default value 'Null' keeps the current unit system.";
NewNatUnits::usage="NewNatUnits is an option to `NaturalUnitsSetup', a list of units to be regarded as fundamental. Doesn't need to be given in full: the rest of the natural units are inferred automatically from the basic units. Deafult value 'Null' keeps the current natural units, unless 'NewUnitSystem' is set, in which case natural units revert to basic units by default.";
UnitMaxIterations::usage="UnitMaxIterations is an option to `NaturalUnitsSetup', the maximum number of iterations to use when checking whether a set of units converges to a smaller set of basic units.";
NaturalUnitsSetStyle::usage="NaturalUnitsSetStyle[style] sets the style (color etc.) in which units are displayed in the output.\nNaturalUnitsSetStyle[] resets the styles of all units.";
NaturalUnitsAdd::usage="NaturalUnitsAdd[newunit->unitdefinition] adds 'newunit' the unit list with the given definition in terms of other units.";
Prefixable::usage="Prefixable is an option to `NaturalUnitsAdd'. If `True', the new unit is marked to take prefixes.";
UnitDefinition::usage="UnitDefinition is an option to `NaturalUnitsAdd', giving the usage message of the new unit.";
NaturalUnitsCheck::usage="NaturalUnitsCheck[unit] checks if `unit' is defined and, if so, gives information on its definition and relations to other units.";
NaturalUnitsList::usage="NaturalUnitsList[] prints out all the unit symbols and physical constants and their explanations (excluding those formed with prefixes).";

(* Unit lists *)
$FromBasicToNatUnits::usage="$FromBasicToNatUnits contains replacement rules for basic (by default: SI) units into natural units.";
$UnitRules::usage="$UnitRules contains replacement rules for units derived from other units. It also includes physical constants. This is the basis of all unit conversion data in NaturalUnits`.";
$PrefixedUnitRules::usage="$PrefixedUnitRules is a version of $UnitRules with prefixed units added.";
$PrefixableUnits::usage="$PrefixableUnits is a list of all units that take metric prefixes.";
$PrefixedUnitsDroplist::usage="$PrefixedUnitsDroplist is a list of units that should be dropped from final unit rules to not cause confusion (in practice, kg).";
$PrefixTable::usage="$PrefixTable is a list of metric prefixes.";
$ToNatUnits::usage="$ToNatUnits contains the collected and combined replacement rules from all basic and derived units into $NatUnits.";
$BasicUnits::usage="$BasicUnits is a list of basic unit symbols to which all other units converge.";
$AllUnits::usage="$AllUnits is a list of all unit symbols.";

(* Units and constants *)
c::usage="speed of light";
hbar::usage="reduced Planck constant";
kB::usage="Boltzmann constant";
G::usage="Newton's constant";
eV::usage="electronvolt";

m::usage="meter";
s::usage="second";
kg::usage="kilogram";
K::usage="kelvin";
min::usage="minute";
year::usage="year";
Hz::usage="hertz";
mP::usage="Planck mass";
MP::usage="reduced Planck mass";
tP::usage="Planck time";
TP::usage="Planck temperature";
lP::usage="Planck length";
EP::usage="Planck energy";
pc::usage="parsec";
ly::usage="lightyear";
AU::usage="astronomical unit";
g::usage="gram";
MSolar::usage="solar mass";
MEarth::usage="Earth mass";
MMoon::usage="Moon mass";
mProton::usage="proton mass";
mElectron::usage="elecron mass";
J::usage="joule";
Newt::usage="newton";
W::usage="watt";
h::usage="hour";
month::usage="month";
y::usage="year";
Coul::usage="coulomb";
e::usage="elementary charge";
A::usage="ampere";
V::usage="volt";
epsilon0::usage="vacuum permittivity";
mu0::usage="vacuum permeability";
deg::usage="degree";
rad::usage="radian";
ke::usage="Coulomb's constant";
b::usage="barn";


$NatUnits::usage="$NatUnits contains the current natural units.";

(* Systems of equations for different unit systems *)
$NaturalUnitSystem::usage="$NaturalUnitSystem is the unit system currently in use.";
$NaturalUnitSystemList::usage="$NaturalUnitSystemList is a list of all built-in natural unit systems.";

CosmologyUnits::usage="Equation system for units used in early universe cosmology: c = hbar = kB = 1.";
PlanckUnits::usage="Equation system with c = hbar = G = kB = 1, setting Planck mass, length, time, and temperature to one.";
GeometrizedUnits::usage="Equation system with c = G = 1.";
AtomicUnits::usage="Equation system with e = mElectron = hbar = ke = 1.";
ParticlePhysicsUnits::usage="Equation system used in particle physics and QFT: c = hbar = epsilon0 = 1.";
SuperNaturalUnits::usage="Fully dimensionless units, with  c = hbar = kB = MP = epsilon0 = rad = 1.";

(************ Conversion data **************)

$UnitRules = 
<|c -> 299792458 m/s,
hbar -> 6.62607015*10^-34 m^2*kg/s/(2*Pi),
kB -> 1.380649*10^-23 kg*m^2/s^2/K,
G -> 6.674*10^-11 m^3/(kg*s^2),
eV -> 1.602176634*10^-19 kg*m^2/s^2,
min-> 60 s,
year -> 31556926 s,
Hz -> 1/s,
mP -> Sqrt[hbar*c/G],
MP -> Sqrt[hbar*c/(8*Pi*G)],
tP -> Sqrt[hbar*G/c^5],
TP -> Sqrt[hbar*c^5/(G*kB^2)],
lP -> Sqrt[hbar*G/c^3],
EP -> Sqrt[hbar*c^5/G],
pc -> 3.08567758*10^16 m,
ly -> 9.4605284*10^15 m,
AU -> 149597871*10^3 m,
MSolar -> 1.989*10^30 kg,
MEarth -> 5.972*10^24 kg,
MMoon -> 7.348*10^22 kg,
mProton -> 1.6726219*10^-27 kg,
mElectron -> 9.10938356*10^-31 kg,
J -> kg*m^2/s^2,
Newt -> kg*m/s^2,
W -> J/s,
h -> 60 min,
month -> year/12,
y -> year,
e -> 1.602176634*10^-19 Coul,
A -> Coul/s, 
V -> J/Coul,
epsilon0 -> 8.8541878128*10^-12 s^4*A^2*kg^-1*m^-3,
mu0 -> 1.25663706212*10^-6 kg*m*s^-2*A^-2,
deg -> (2*Pi)/360 rad,
ke -> 1/(4*Pi*epsilon0),
b -> 10^-28 m^2,
g -> 10^-3 kg|>;

$PrefixedUnitRules=<||>;

(* Units that take matric prefixes (milli-, centi-, kilo-, ...) *)
$PrefixableUnits={s,Hz,m,pc,ly,eV,J,Newt,W,y,A,V,Coul,g};

(* Units that should be dropped after adding prefixes (this is a trick to make g and kg work right) *)
$PrefixedUnitsDroplist={kg};

(* Metric prefixes and the corresponding multipliers. *)
$PrefixTable=<|"f"->10^-15,"p"->10^-12,"n"->10^-9,"mu"->10^-6,"m"->10^-3,"c"->10^-2,"k"->10^3,"M"->10^6,"G"->10^9,"T"->10^12|>;

(* To be set later. *)
Unprotect[$ToNatUnits];
Clear[$ToNatUnits];
Protect[$ToNatUnits];


(************ Units systems **************)

CosmologyUnits={c == hbar == kB == 1};
PlanckUnits={c == hbar == kB == G == 1};
GeometrizedUnits={c == G == 1};
AtomicUnits={e == mElectron == hbar == ke == 1};
ParticlePhysicsUnits={c == hbar == epsilon0 == 1};
SuperNaturalUnits={c == hbar == kB == MP == epsilon0 == rad == 1};

$NaturalUnitSystemList=Hold[CosmologyUnits, PlanckUnits, GeometrizedUnits, AtomicUnits, ParticlePhysicsUnits, SuperNaturalUnits];


(************ Messages (used in functions below) **************)
General::wrongopts="Bad options to ``: ``";
NaturalUnitsSetup::badeq="Error: Expression `` in $NaturalUnitSystem is not an equation. Setup failed.";
NaturalUnitsSetup::badunit="Error: Expression `` in $NaturalUnitSystem is not a symbol. Setup failed.";
NaturalUnitsSetup::nobasics="Error: $UnitRules not converging to a set of basic units in `` iterations. Setup failed.";
NaturalUnitsSetup::nobasicsol="Error: can't solve basic units in terms of natural units from the unit system equations. Setup failed.";
NaturalUnitsSetup::nonatunits="Error: $ToNatUnits not converging to a set of basic units in `` iterations. Setup failed.";
NaturalUnitsSetup::wrongarg="NaturalUnitsSetup expects no arguments, only options.";
NatAssociateTo::unitdpl="Some units duplicated in `` and ``. Duplicated units and tallies: ``. Later instance takes precedence.";
NaturalUnitsAdd::badinput="Error: bad input. Expected input is NaturalUnitsAdd[newunit->unitdefinition, options], where 'newunit' is a symbol or a string.";
NaturalUnitsCheck::wrongarg="NaturalUnitsCheck expects one argument (symbol or corresponding string).";
NaturalConvert::nounits="Error: Natural units not successfully defined. Try calling NaturalUnitsSetup[].";
NaturalConvert::nounits2="Error: Natural unit list is empty. Can't fix the powers of any unit list. Try the non-list version NaturalConvert[quantity, targetunit] with a single target unit instead.";
NaturalConvert::abort="Conversion aborted without altering the input quantity.";
NaturalConvert::zerounit="Error: target unit seems to be zero.";
NaturalConvert::npowers="Error: quantity `` is not monomial in natural units.";
NaturalConvert::unitlistcomp="Could not match the units of `` and those from the list ``.";
NaturalConvert::unitcomp="Units of `` and `` not compatible! (Or unit not recognized.)";
NaturalConvert::corrpow="Correcting the power of target units.";


(************ Function implementations **************)

Begin["`Private`"];


(* Custom function for appending associations which warns for repeated keys. *)
Unprotect[NatAssociateTo];
SetAttributes[NatAssociateTo, HoldFirst];
NatAssociateTo[asTo_, asFrom_]:=Module[{duplicates=Select[Tally[Join[Keys[asTo],Keys[asFrom]]], (#[[2]]>1)&]},
    If[Length[duplicates]>0, Message[NatAssociateTo::unitdpl, asTo, asFrom, duplicates]];
    AssociateTo[asTo,asFrom];
];
Protect[NatAssociateTo];

(*** Parse unit lists ($FromBasicToNatUnits, $UnitRules, $PrefixableUnits, $PrefixableUnits, $PrefixTable) into usable form ($ToNatUnits). Called once when the package is loaded. Call again, if units edited or new unit system introduced. ***)
Unprotect[NaturalUnitsSetup];
Options[NaturalUnitsSetup]={NewUnitSystem->Null, NewNatUnits->Null, UnitMaxIterations->20};
NaturalUnitsSetup[opts:OptionsPattern[]]:=Quiet[#,{OptionValue::nodef}]&@Unevaluated@Module[{unprotectedVars, newUnits, newSymbol, basicUnitEqs, natUnitTemp, unitSystemPos},
    optionCheck[NaturalUnitsSetup,opts];
    unprotectedVars:=Hold[unitAssumptions,$NatUnits,$ToNatUnits,$AllUnits,Evaluate[$AllUnits]];
    Unprotect@@unprotectedVars;
    Clear[$ToNatUnits];

    (* Parse input parameters: for Null, maintain current settings. Check validity. *)
    If[OptionValue[NewUnitSystem]=!=Null,
        (* Set the new unit system. If possible, recognize it, and keep it in held form so that the user can see the name stored in $NaturalUnitSystem. *)
        $NaturalUnitSystem=OptionValue[NewUnitSystem];
        unitSystemPos=Position[{ReleaseHold[$NaturalUnitSystemList]}, $NaturalUnitSystem];
        If[Length[unitSystemPos]>0, $NaturalUnitSystem=Extract[$NaturalUnitSystemList,First[unitSystemPos],HoldForm]];
    ];
    Scan[
        If[Head[#]=!=Equal,
            Message[NaturalUnitsSetup::badeq,#];
            Protect@@unprotectedVars;
            Return[$Failed,Module];
        ]&,
        ReleaseHold[$NaturalUnitSystem]
    ];

    If[OptionValue[NewNatUnits]=!=Null,
        $NatUnits=OptionValue[NewNatUnits],
        If[OptionValue[NewUnitSystem]=!=Null,
            (* New unit system is being set without new natural units; avoid errors by reverting to defaults. *)
            $NatUnits={};
        ];
    ];
    Scan[
        If[Head[#]=!=Symbol,
            Message[NaturalUnitsSetup::badunit,#];
            Protect@@unprotectedVars;
            Return[$Failed,Module];
        ]&,
        $NatUnits
    ];

    (* Add prefixes for those units that need them. *)
    $PrefixedUnitRules=$UnitRules;
    Scan[
        Function[{unitNow},
        newUnits=Map[Function[{prefixNow}, newSymbol=Symbol["NaturalUnits`"<>prefixNow<>ToString[unitNow]];(Evaluate[newSymbol]::usage=ToString[StringForm["unit formed from `` by prefix ``",unitNow::usage,prefixNow]]);newSymbol->$PrefixTable[prefixNow]*unitNow], Keys[$PrefixTable]];
        NatAssociateTo[$PrefixedUnitRules,newUnits]],
        $PrefixableUnits
    ];

    (* Drop some units (in practice, kg, so that it's not replaced by grams). *)
    KeyDropFrom[$PrefixedUnitRules,$PrefixedUnitsDroplist];

    (* Check that all derived units reduce to some basic units within a few iterations. *)
    Check[$BasicUnits = Variables[ReplaceRepeated[Keys[$PrefixedUnitRules], $PrefixedUnitRules, MaxIterations->OptionValue[UnitMaxIterations]]],
        Message[NaturalUnitsSetup::nobasics,OptionValue[UnitMaxIterations]];
        Protect@@unprotectedVars;
        Return[$Failed];
    ];

    (* Solve basic units from the system of equations, in terms of the natural units. *)
    basicUnitEqs = Join[ReleaseHold[$NaturalUnitSystem], MapIndexed[natUnitTemp@@#2 == #1&, $NatUnits]]//.$PrefixedUnitRules;
    $FromBasicToNatUnits = Association[Quiet[Solve[basicUnitEqs, $BasicUnits]]]/.MapIndexed[natUnitTemp@@#2 -> #1&, $NatUnits];
    If[$FromBasicToNatUnits==<||> && basicUnitEqs!={},
        Message[NaturalUnitsSetup::nobasicsol];
        Protect@@unprotectedVars;
        Return[$Failed];
    ];

    (* Join basic and derived units into one association in terms of $NatUnits. Drop $NatUnits from the lists. *)
    $ToNatUnits=KeyDrop[$FromBasicToNatUnits,$NatUnits];
    NatAssociateTo[$ToNatUnits,KeyDrop[$PrefixedUnitRules,$NatUnits]];
    
    (* "Flatten" $ToNatUnits so that everything is transformed to $NatUnits in one iteration. *)
    $ToNatUnits=$ToNatUnits//.$ToNatUnits;

    (* Update natural units list as the true, final list of all independent final units. *)
    Check[$NatUnits = Variables[ReplaceRepeated[Keys[$ToNatUnits], $ToNatUnits, MaxIterations->OptionValue[UnitMaxIterations]]],
        Message[NaturalUnitsSetup::nonatunits,OptionValue[UnitMaxIterations]];
        Clear[$ToNatUnits];
        Protect@@unprotectedVars;
        Return[$Failed];
    ];

    (* Finalize variables. *)
    unitAssumptions=And@@Thread[Join[$NatUnits,Keys[$ToNatUnits]] > 0];
    $ToNatUnits=AssociationMap[N[Simplify[#,unitAssumptions]]&, $ToNatUnits];
    $AllUnits=Union[Keys[$ToNatUnits],$NatUnits];

    Protect@@unprotectedVars; (* Protect units to prevent accidental overwrites. *)
];
NaturalUnitsSetup[___]:=Message[NaturalUnitsSetup::wrongarg];
NaturalUnitsSetup[NewUnitSystem->{}, NewNatUnits->{}];
Protect[NaturalUnitsSetup,NewUnitSystem,NewNatUnits,UnitMaxIterations];


(*** Set formatting for all units. ***)
Unprotect[NaturalUnitsSetStyle];
NaturalUnitsSetStyle[newStyle_]:=Module[{},
    Unprotect[Evaluate[$AllUnits]];
    Scan[(#/:MakeBoxes[#, StandardForm]:=With[{box = ToBoxes[Style[ToString[#], newStyle]]}, 
        InterpretationBox[RowBox[{box}], #]])&, $AllUnits];
    Protect[Evaluate[$AllUnits]];
];
NaturalUnitsSetStyle[]:=Module[{},
    Unprotect[Evaluate[$AllUnits]];
    Scan[(#/:MakeBoxes[#, StandardForm]=.)&, $AllUnits];
    Protect[Evaluate[$AllUnits]];
];
NaturalUnitsSetStyle[pars___]:=Message[NaturalUnitsSetStyle::argt,"NaturalUnitsSetStyle",Length[{pars}],0,1];
Protect[NaturalUnitsSetStyle];


(*** Add an unit. Strings are processed into symbols under NaturalUnits`. ***)
Unprotect[NaturalUnitsAdd]
Options[NaturalUnitsAdd]={Prefixable->False, UnitDefinition->"custom unit"};
NaturalUnitsAdd[newUnit_?AtomQ->newResult_,opt:OptionsPattern[]]:=Quiet[#,{OptionValue::nodef}]&@Unevaluated@Module[{modNewUnit=newUnit/.str_String:>Symbol["NaturalUnits`"<>str], modNewResult=newResult/.str_String:>Symbol["NaturalUnits`"<>str]},
    optionCheck[NaturalUnitsAdd,opt];

    NatAssociateTo[$UnitRules,Association[modNewUnit->modNewResult]];
    Evaluate[modNewUnit]::usage=OptionValue[UnitDefinition];
    If[OptionValue[Prefixable],
        $PrefixableUnits=DeleteDuplicates[Append[$PrefixableUnits, modNewUnit]],
        $PrefixableUnits=DeleteCases[$PrefixableUnits, modNewUnit];
    ];
    
    NaturalUnitsSetup[];
];
NaturalUnitsAdd[___]:=Message[NaturalUnitsAdd::badinput];
Protect[NaturalUnitsAdd,Prefixable,UnitDefinition];


(*** Check whether a given unit exists. ***)
Unprotect[NaturalUnitsCheck];
NaturalUnitsCheck[checkedUnit_Symbol]:=Quiet[Module[{},
    Which[KeyExistsQ[$ToNatUnits, checkedUnit], Print[StringForm["Unit '``' exists as ``=``", checkedUnit, checkedUnit, $ToNatUnits[checkedUnit]]],
        MemberQ[$NatUnits,checkedUnit], Print[StringForm["Unit '``' exists as a natural unit.", checkedUnit]],
        True, Print[StringForm["Unit '``' does not exist in the final replacement rules.", checkedUnit]]];
    If[MemberQ[$BasicUnits, checkedUnit],
        Print["It is a basic unit."],
        If[KeyExistsQ[$PrefixedUnitRules, checkedUnit],
            Print[StringForm["Original definition: ``=``", checkedUnit, $PrefixedUnitRules[checkedUnit]]]
        ];
    ];
    If[MemberQ[$PrefixableUnits, checkedUnit], Print[StringForm["Unit can take metric prefixes."]]];
]];
NaturalUnitsCheck[checkedUnit_String]:=NaturalUnitsCheck[Symbol[checkedUnit]];
NaturalUnitsCheck[___]:=Message[NaturalUnitsCheck::wrongarg];
Protect[NaturalUnitsCheck];


(*** List all units (without prefixes). ***)
Unprotect[NaturalUnitsList];
NaturalUnitsList[]:=Scan[Print[#, " = ", Evaluate[#]::usage] &, Variables[$UnitRules/.Association[x___]->{x}/.Rule[x_, y_]->Sequence[x, y]]];
NaturalUnitsList[pars__]:=Message[NaturalUnitsList::argrx,"NaturalUnitsList",Length[{pars}],"no"];
Protect[NaturalUnitsList];

(*** Convert one unit to another. ***)
Unprotect[NaturalConvert];
SetAttributes[NaturalConvert, HoldRest];
Options[NaturalConvert]={ToUnitDefault->NoUnit, KeepUnit->True};

NaturalConvert[qty_, opts:OptionsPattern[]]:=Quiet[#,{OptionValue::nodef}]&@Unevaluated@Module[{qtyNat},
    optionCheck[NaturalConvert,opts];
    
    If[OptionValue[ToUnitDefault]=!=NoUnit, Return[NaturalConvert[qty,OptionValue[ToUnitDefault],opts]]];
    If[!ValueQ[$ToNatUnits],
        Message[NaturalConvert::nounits];
        Message[NaturalConvert::abort];
        Return[qty];
    ];

    qtyNat=antiRound[qty/.$ToNatUnits];
    If[OptionValue[KeepUnit]===False,
        (qtyNat/.Thread[$NatUnits->1]),
        antiRound[PowerExpand[FullSimplify[qtyNat,unitAssumptions]]]
    ]
];

NaturalConvert[qty_, toUnitIn_, opts:OptionsPattern[]]:=Quiet[#,{OptionValue::nodef}]&@Unevaluated@Module[{toUnit=toUnitIn, toUnitDef=Defer[toUnitIn], qtyNat, toUnitNat, qtyGood, qtyPows, toUnitGood, toUnitPows, powerSol, extraPower, qtyFinal},
    optionCheck[NaturalConvert,opts];

    If[!ValueQ[$ToNatUnits],
        Message[NaturalConvert::nounits];
        Message[NaturalConvert::abort];
        Return[qty];
    ];

    If[PossibleZeroQ[qty], Return[qty]];

    If[PossibleZeroQ[toUnit],
        Message[NaturalConvert::zerounit];
        Message[NaturalConvert::abort];
        Return[qty];
    ];

    (* First, convert quantity and target unit to natural units. *)
    qtyNat=qty/.$ToNatUnits;
    toUnitNat=toUnit/.$ToNatUnits;

    If[$NatUnits!={},
        (* Find the powers of natural units in quantity and target unit. *)
        {qtyGood,qtyPows}=ExtractNaturalPowers[qtyNat];
        If[!qtyGood,
            Message[NaturalConvert::npowers,qty];
            Message[NaturalConvert::abort];
            Return[qty];
        ];

        {toUnitGood,toUnitPows}=ExtractNaturalPowers[toUnitNat];
        If[!toUnitGood,
            Message[NaturalConvert::npowers,toUnit];
            Message[NaturalConvert::abort];
            Return[qty];
        ];

        (* Check if they are equal; if not, try to correct this. *)
        If[qtyPows!=toUnitPows,
            Message[NaturalConvert::unitcomp, qty, toUnit];
            powerSol=Solve[qtyPows==extraPower*toUnitPows,extraPower];
            If[powerSol=={},
                Message[NaturalConvert::abort];
                Return[qty],
                Message[NaturalConvert::corrpow];
            ];
            extraPower=extraPower/.First[powerSol];
            toUnitDef=Defer[toUnitIn^"extraPower"]/."extraPower"->extraPower;
            toUnit=antiRound[PowerExpand[Simplify[toUnit^extraPower,unitAssumptions]]];
            toUnitNat=toUnitNat^extraPower;
        ];
    ];

    (* Calculate the final number. *)
    qtyFinal=antiRound[qtyNat/toUnitNat/.Thread[$NatUnits->1]]; (* Set $NatUnits->1 to get rid of any errors from numerical inaccuracy. *)

    (* Add the unit, if wanted. *)
    Which[OptionValue[KeepUnit]===False,qtyFinal,
        OptionValue[KeepUnit]===Defer,qtyFinal*toUnitDef,
        True, antiRound[qtyFinal*toUnit]]
];

NaturalConvert[qty_, toUnits_?(MatchQ[#,_List]&), opts:OptionsPattern[]]:=Quiet[#,{OptionValue::nodef}]&@Unevaluated@Module[{toUnit, toUnitDef, qtyNat, toUnitsNat, qtyGood, qtyPows, toUnitsGood, toUnitsPows, powerSol, extraPowers},
    optionCheck[NaturalConvert,opts];

    If[!ValueQ[$ToNatUnits],
        Message[NaturalConvert::nounits];
        Message[NaturalConvert::abort];
        Return[qty];
    ];

    If[$NatUnits=={},
        Message[NaturalConvert::nounits2];
        Message[NaturalConvert::abort];
        Return[qty];
    ];

    If[PossibleZeroQ[qty], Return[qty]];

    (* Convert quantity and target units to natural units. *)
    qtyNat=qty/.$ToNatUnits;
    toUnitsNat=toUnits/.$ToNatUnits;

    (* Find their powers of natural units. *)
    {qtyGood,qtyPows}=ExtractNaturalPowers[qtyNat];
    If[!qtyGood,
        Message[NaturalConvert::npowers,qty];
        Message[NaturalConvert::abort];
        Return[qty];
    ];

    {toUnitsGood,toUnitsPows}=Transpose[ExtractNaturalPowers[toUnitsNat]];
    MapIndexed[
        If[!#1,
            Message[NaturalConvert::npowers,Extract[toUnits,#2]];
            Message[NaturalConvert::abort];
            Return[qty,Module]]&,
        toUnitsGood
    ];

    (* Try to solve the correct powers for the target units to get a match. *)
    extraPowers=Table[pw[i],{i,Length[toUnits]}];
    If[Check[powerSol=Solve[qtyPows==extraPowers.toUnitsPows,extraPowers],True]||powerSol=={},
        Message[NaturalConvert::unitlistcomp,qty,toUnits];
        Message[NaturalConvert::abort];
        Return[qty];
    ];
    extraPowers=extraPowers/.First[powerSol];
    toUnitDef=Times@@Thread[Power[Thread[Defer[toUnits]], extraPowers]];
    toUnit=antiRound[PowerExpand[Simplify[Times@@(toUnits^extraPowers),unitAssumptions]]];
    toUnitNat=Times@@(toUnitsNat^extraPowers);

    (* Calculate the final number. *)
    qtyFinal=antiRound[qtyNat/toUnitNat/.Thread[$NatUnits->1]];

    (* Add the unit, if wanted. *)
    Which[OptionValue[KeepUnit]===False,qtyFinal,
        OptionValue[KeepUnit]===Defer,qtyFinal*toUnitDef,
        True, antiRound[qtyFinal*toUnit]]
];

NaturalConvert[pars___]:=Message[NaturalConvert::argt,"NaturalConvert",Length[{pars}],1,2];

(* Make NaturalConvert "Listable" w.r.t. its first argument: *)
PrependTo[DownValues[NaturalConvert],HoldPattern[NaturalConvert[qties_List,rest:___]]:>(NaturalConvert[#,rest]&/@qties)];


Protect[NaturalConvert,ToUnitDefault,KeepUnit,NoUnit];

(* Helper function: Check if expression is proportional to a single power of the given variable (to the best of our ability). Give truth value and the power out. *)
SetAttributes[OnePowerQ,Listable];
OnePowerQ[expr_,x_]:=Module[{y,testExpr},
    testExpr=PowerExpand[Log[expr]/.x->Exp[y]];
    {PolynomialQ[testExpr,y]&&(testExpr==0||Exponent[testExpr,y]==0||Exponent[testExpr,y]==1),Exponent[expr,x]}
];

(* Helper function: extract a row of powers, one for each natural unit. This is the second part of the output; the first part is True/False depending on whether the expression is proportional to single powers of all the natural units. *)
SetAttributes[ExtractNaturalPowers,Listable];
ExtractNaturalPowers[expr_]:={And@@#1,#2}&@@Transpose[OnePowerQ[expr,$NatUnits]];

(* Hepler function: check for bad options. *)
optionCheck[func_,opts_]:=Module[{badrules=FilterRules[opts,Except[Options[func]]]},
    If[Length[badrules]>0, Message[func::wrongopts,func,badrules];False,True]
];

(* Helper function: transform small real values to integers, if appropriate. Makes output cleaner in some special cases. *)
antiRound[x_]:=x/.y_?(NumericQ[#]&&#==Round[#]&&Abs[#]<=10^3&)->Round[y];

End[];

EndPackage[];
