Line data Source code
1 : #ifndef UNIT_UNIT_H
2 : #define UNIT_UNIT_H
3 :
4 :
5 : #include <array>
6 : #include <charconv>
7 : #include <cstdlib>
8 : #include <iostream>
9 : #include <system_error>
10 : #include <unordered_map>
11 : #include <string>
12 : #include <stdexcept>
13 : #include <cmath>
14 : #include <regex>
15 :
16 : //#include "Set/Base.H"
17 : //#include "String.H"
18 :
19 : struct Unit : std::pair<double,std::array<int,7>>
20 : {
21 :
22 : //
23 : // STATIC CONST DEFINITIONS, TRUE EVERYWHERE
24 : //
25 : enum class Type {
26 : Length=0, Time=1, Mass=2, Temperature=3, Current=4, Amount=5, LuminousIntensity=6
27 : };
28 : // Fundamental units
29 : inline static const std::map<std::string, std::pair<double,std::array<int,7>>> base_units =
30 : {
31 : {"1", {1.0, {0,0,0,0,0,0,0}}},
32 : {"m", {1.0, {1,0,0,0,0,0,0}}},
33 : {"s", {1.0, {0,1,0,0,0,0,0}}},
34 : {"kg", {1.0, {0,0,1,0,0,0,0}}},
35 : {"K", {1.0, {0,0,0,1,0,0,0}}},
36 : {"A", {1.0, {0,0,0,0,1,0,0}}},
37 : {"mol",{1.0, {0,0,0,0,0,1,0}}},
38 : {"cd", {1.0, {0,0,0,0,0,0,1}}}
39 : };
40 : //
41 : // Compound unit relationships and nomenclature
42 : // Any desired unit relationship can be added to this list.
43 : // The unit conversion can be specified in any desired format, as long as
44 : // all of the units on the right hand side are resolvable.
45 : //
46 : inline static const std::map<std::string, std::pair<double,std::string>> compound = {
47 : // Lengths
48 : {"km", {1e+3, "m" }},
49 : {"hm", {1e+2, "m" }},
50 : {"dam", {1e+1, "m" }},
51 : {"dm", {1e-1, "m" }},
52 : {"cm", {1e-2, "m" }},
53 : {"mm", {1e-3, "m" }},
54 : {"um", {1e-6, "m" }},
55 : {"nm", {1e-9, "m" }},
56 : {"pm", {1e-12, "m" }},
57 : {"in", {0.0254, "m" }},
58 : {"ft", {0.3048, "m" }},
59 : {"yd", {0.9144, "m" }},
60 : {"mi", {5280.0, "ft" }},
61 : {"furlong", {660.0, "ft" }},
62 : {"au", {1.495978707e11, "m" }}, // astronomical unit
63 : {"ly", {9.4607304725808e15, "m"}}, // lightyear
64 : {"pc", {3.08567758149137e16, "m"}}, // parsec
65 :
66 : // Times
67 : {"ms", {1e-3, "s" }},
68 : {"us", {1e-6, "s" }},
69 : {"ns", {1e-9, "s" }},
70 : {"ps", {1e-12, "s" }},
71 : {"min", {60.0, "s" }},
72 : {"hr", {60.0, "min" }},
73 : {"day", {24.0, "hr" }},
74 : {"week", {7.0, "day" }},
75 : {"fortnight", {14.0, "day" }},
76 : {"year", {365.25, "day" }},
77 : {"decade", {10.0, "year" }},
78 :
79 : // Mass
80 : {"g", {1e-3, "kg" }},
81 : {"mg", {1e-6, "kg" }},
82 : {"ug", {1e-9, "kg" }},
83 : {"tonne", {1000.0, "kg" }},
84 : {"lb", {0.45359237, "kg" }},
85 : {"oz", {0.0283495, "kg" }},
86 : {"slug", {14.5939, "kg" }},
87 :
88 : // Speed
89 : {"mph", {1.0, "mi/hr" }},
90 : {"kph", {1.0, "km/hr" }},
91 : {"mps", {1.0, "m/s" }},
92 : {"fps", {1.0, "ft/s" }},
93 : {"knot", {1852.0 / 3600.0, "m/s"}},
94 :
95 : // Area
96 : {"ha", {1e4, "m^2" }},
97 : {"acre", {4046.85642, "m^2" }},
98 : {"sqft", {1.0, "ft^2" }},
99 : {"sqin", {1.0, "in^2" }},
100 : {"sqkm", {1.0, "km^2" }},
101 : {"sqmi", {1.0, "mi^2" }},
102 :
103 : // Volume
104 : {"L", {1e-3, "m^3" }},
105 : {"mL", {1e-6, "m^3" }},
106 : {"cc", {1e-6, "m^3" }},
107 : {"gal", {3.78541e-3, "m^3" }},
108 : {"qt", {0.946353e-3, "m^3" }},
109 : {"pt", {0.473176e-3, "m^3" }},
110 : {"cup", {0.24e-3, "m^3" }},
111 : {"floz", {29.5735e-6, "m^3" }},
112 : {"tbsp", {14.7868e-6, "m^3" }},
113 : {"tsp", {4.92892e-6, "m^3" }},
114 :
115 : // Force
116 : {"N", {1.0, "kg*m/s^2" }},
117 : {"kN", {1e3, "N" }},
118 : {"lbf", {4.44822162, "N" }},
119 : {"dyn", {1e-5, "N" }},
120 :
121 : // Energy
122 : {"J", {1.0, "N*m" }},
123 : {"kJ", {1e3, "J" }},
124 : {"MJ", {1e6, "J" }},
125 : {"eV", {1.602176634e-19, "J"}},
126 : {"cal", {4.184, "J" }},
127 : {"kcal", {4184.0, "J" }},
128 : {"BTU", {1055.06, "J" }},
129 : {"ft*lbf", {1.35582, "J" }},
130 :
131 : // Power
132 : {"W", {1.0, "J/s" }},
133 : {"kW", {1e3, "W" }},
134 : {"MW", {1e6, "W" }},
135 : {"hp", {745.7, "W" }},
136 :
137 : // Pressure
138 : {"Pa", {1.0, "N/m^2" }},
139 : {"kPa", {1e3, "Pa" }},
140 : {"MPa", {1e6, "Pa" }},
141 : {"GPa", {1e9, "Pa" }},
142 : {"bar", {1e5, "Pa" }},
143 : {"atm", {101325.0,"Pa" }},
144 : {"psi", {6894.76, "Pa" }},
145 : {"mmHg", {133.322, "Pa" }},
146 : {"torr", {133.322, "Pa" }},
147 :
148 : // Temperature differences (not absolute temperatures)
149 : {"K", {1.0, "K" }},
150 : {"degR", {5.0/9.0, "K" }},
151 :
152 : // Charge
153 : {"C", {1.0, "A*s" }},
154 : {"mAh", {3.6, "C" }},
155 :
156 : // Voltage
157 : {"V", {1.0, "J/C" }},
158 :
159 : // Capacitance
160 : {"F", {1.0, "C/V" }},
161 :
162 : // Resistance
163 : {"ohm", {1.0, "V/A" }},
164 : {"kOhm", {1e3, "ohm" }},
165 : {"MOhm", {1e6, "ohm" }},
166 :
167 : // Conductance
168 : {"S", {1.0, "1/ohm" }},
169 :
170 : // Magnetic
171 : {"T", {1.0, "kg/(A*s^2)"}},
172 : {"G", {1e-4, "T" }}, // gauss
173 :
174 : // Frequency
175 : {"Hz", {1.0, "1/s" }},
176 : {"kHz", {1e3, "Hz" }},
177 : {"MHz", {1e6, "Hz" }},
178 : {"GHz", {1e9, "Hz" }},
179 : };
180 :
181 :
182 : //
183 : // These static factory functions return unit types.
184 : // They are used to check unit compliance, mainly at the parsing stage.
185 : //
186 : // Primary
187 3419 : static Unit Less() { return Unit{1.0, {0,0,0,0,0,0,0}}; } // 1
188 282 : static Unit Length() { return Unit{1.0, {1,0,0,0,0,0,0}}; } // m
189 409 : static Unit Time() { return Unit{1.0, {0,1,0,0,0,0,0}}; } // s
190 49 : static Unit Mass() { return Unit{1.0, {0,0,1,0,0,0,0}}; } // kg
191 25 : static Unit Temperature() { return Unit{1.0, {0,0,0,1,0,0,0}}; } // K
192 : static Unit Current() { return Unit{1.0, {0,0,0,0,1,0,0}}; } // A
193 : static Unit Amount() { return Unit{1.0, {0,0,0,0,0,1,0}}; } // mol
194 : static Unit LuminousIntensity() { return Unit{1.0, {0,0,0,0,0,0,1}}; } // cd
195 :
196 : // Geometry and kinematics
197 24 : static Unit Area() { return Length()*Length(); }
198 7 : static Unit Volume() { return Length()*Length()*Length(); }
199 40 : static Unit Velocity() { return Length() / Time();}
200 38 : static Unit Acceleration(){ return Velocity() / Time(); }
201 : // Mechanics
202 38 : static Unit Force() { return Mass() * Acceleration(); } // kg·m/s² = N
203 : static Unit Momentum() { return Mass() * Velocity(); }
204 : static Unit Impulse() { return Force() * Time(); }
205 13 : static Unit Pressure() { return Force() / Area(); } // N/m² = Pa
206 25 : static Unit Energy() { return Force() * Length(); } // N·m = J
207 6 : static Unit Power() { return Energy() / Time(); } // J/s = W
208 7 : static Unit Density() { return Mass() / Volume(); }
209 : static Unit SpecificWeight() { return Force() / Volume(); }
210 : static Unit Work() { return Energy(); }
211 : // Thermodynamics
212 4 : static Unit SpecificHeatCapacity() {return Energy() / Mass() / Temperature() ;}
213 4 : static Unit ThermalConductivity() {return Power()/Length()/Temperature();}
214 2 : static Unit ThermalDiffusivity() {return ThermalConductivity() / Density() / SpecificHeatCapacity();}
215 : static Unit HeatFlux() { return Power() / Area(); }
216 : static Unit HeatTransferCoefficient() { return HeatFlux() / Temperature(); }
217 : // Electromagnetism
218 :
219 : static Unit Charge() { return Current() * Time(); } // A·s = C
220 : static Unit Voltage() { return Power() / Current(); } // W/A = V
221 : static Unit Resistance() { return Voltage() / Current(); } // V/A = Ω
222 : static Unit Capacitance() { return Charge() / Voltage(); } // C/V = F
223 : static Unit Conductance() { return Less() / Resistance(); } // S = 1/Ω
224 : static Unit MagneticFlux(){ return Voltage() * Time(); } // V·s = Wb
225 : static Unit MagneticField(){ return MagneticFlux() / Area(); } // Wb/m² = T
226 : static Unit Inductance() { return MagneticFlux() / Current(); } // Wb/A = H
227 : // Chemistry
228 : static Unit MolarMass() { return Mass() / Amount(); } // kg/mol
229 : static Unit Concentration(){ return Amount() / Volume(); } // mol/m³
230 : static Unit MolarEnergy() { return Energy() / Amount(); } // J/mol
231 : // Optics and radiation
232 : static Unit LuminousFlux() { return LuminousIntensity(); } // for scalar usage (lm)
233 : static Unit Illuminance() { return LuminousFlux() / Area(); } // lux = lm/m²
234 :
235 :
236 : //
237 : // This is where SYSTEM NORMALIZATION is stored.
238 : // Should usually be set once (at the beginning of the simulation) and then
239 : // left constant throughout the simulation.
240 : // All normalizations are in terms of standard SI units: m,s,kg,K,A,mol,cd.
241 : //
242 : inline static std::array<std::pair<std::string,double>, 7> normalization = {
243 : std::pair{"m", 1.0},
244 : std::pair{"s", 1.0},
245 : std::pair{"kg", 1.0},
246 : std::pair{"K", 1.0},
247 : std::pair{"A", 1.0},
248 : std::pair{"mol",1.0},
249 : std::pair{"cd", 1.0}
250 : };
251 :
252 : using std::pair<double, std::array<int,7>>::pair;
253 3759 : Unit(const std::pair<double, std::array<int, 7>>& p) : std::pair<double, std::array<int, 7>>(p) {}
254 : static std::map<std::string, int> ParseUnitString(const std::string& input);
255 :
256 :
257 679 : static Unit Parse(double val, std::string unitstring)
258 : {
259 679 : Unit ret = StringParse(unitstring);
260 679 : ret.first *= val;
261 679 : return ret;
262 : }
263 1283 : static Unit Parse(std::string unit)
264 : {
265 1283 : if (unit.size() == 0) return Unit::Less();
266 1260 : std::vector<std::string> split;
267 1260 : std::stringstream ss(unit);
268 1260 : std::string item;
269 3199 : while (std::getline(ss, item, '_')) {
270 1939 : split.push_back(item);
271 : }
272 1260 : if (split.size() == 2)
273 : {
274 679 : return Parse(std::stod(split[0]),split[1]);
275 : }
276 581 : else if (split.size() == 1)
277 : {
278 : // Attempt to parse the string as if it's a unitless number
279 581 : double value = NAN;
280 :
281 : // auto first = unit.data();
282 : // auto last = unit.data() + unit.size();
283 : // auto result = std::from_chars(first, last, value); // apparently this fails on mac
284 :
285 581 : char *end = nullptr;
286 581 : value = std::strtod(unit.c_str(), &end);
287 :
288 : // If it is a untless number, then return a unitless result
289 : //if (result.ec == std::errc() && result.ptr == last) // fails on mac
290 581 : if (end != unit.data() && *end == '\0')
291 : {
292 581 : Unit ret = Unit::Less();
293 581 : ret.first = value;
294 581 : return ret;
295 : }
296 : // Otherwise, it must be a pure unit, so return that.
297 : else
298 : {
299 0 : return StringParse(unit);
300 : }
301 : }
302 : else
303 : {
304 0 : throw std::runtime_error("Invalid unit string "+unit);
305 : }
306 1260 : }
307 :
308 :
309 1983 : static Unit StringParse(std::string unitstring)
310 : {
311 1983 : Unit type = Unit::Less();
312 1983 : std::map<std::string, int> map = ParseUnitString(unitstring);
313 6968 : for (auto it = map.begin(); it != map.end(); it++)
314 : {
315 4985 : auto [str, exp] = *it;
316 :
317 4985 : if (Unit::base_units.find(str) != base_units.end())
318 : {
319 3759 : Unit t = Unit::base_units.at(str);
320 3759 : type = type * (t^exp);
321 : }
322 1226 : else if (compound.find(it->first) != compound.end())
323 : {
324 1226 : auto [subfac, substr] = compound.at(str);
325 1226 : Unit t = StringParse(substr);
326 1226 : t.first *= subfac;
327 1226 : type = type * t^exp;
328 1226 : }
329 : else
330 : {
331 0 : throw std::runtime_error("Unknown unit: " + unitstring);
332 : }
333 4985 : }
334 3966 : return type;
335 1983 : }
336 :
337 :
338 39 : static void setLengthUnit(std::string unit)
339 : {
340 39 : Unit parsed = StringParse(unit);
341 39 : if (parsed.second != Length().second) throw std::runtime_error("Not a unit of length: " + unit);
342 39 : normalization[int(Type::Length)].first = unit;
343 39 : normalization[int(Type::Length)].second = 1.0 / parsed.first;
344 39 : }
345 39 : static void setTimeUnit(std::string unit)
346 : {
347 39 : Unit parsed = StringParse(unit);
348 39 : if (parsed.second != Time().second) throw std::runtime_error("Not a unit of time: " + unit);
349 39 : normalization[int(Type::Time)].first = unit;
350 39 : normalization[int(Type::Time)].second = 1.0 / parsed.first;
351 39 : }
352 :
353 1180 : bool isType(const Unit &test) const
354 : {
355 1180 : return this->second == test.second;
356 : }
357 :
358 5086 : Unit operator * (const Unit &rhs)
359 : {
360 5086 : Unit ret;
361 5086 : ret.first = this->first * rhs.first;
362 45774 : for (unsigned int i = 0; i < this->second.size(); i++) ret.second[i] = (this->second)[i] + rhs.second[i];
363 5086 : return ret;
364 : }
365 : friend Unit operator*(double lhs, const Unit& rhs) {
366 : Unit ret;
367 : ret.first = lhs * rhs.first;
368 : ret.second = rhs.second;
369 : return ret;
370 : }
371 3 : friend Unit operator / (double lhs, const Unit& rhs) {
372 3 : Unit ret = Unit::Less() / rhs;
373 3 : ret.first = lhs * rhs.first;
374 3 : return ret;
375 : }
376 150 : Unit operator / (const Unit &rhs)
377 : {
378 150 : Unit ret;
379 150 : ret.first = this->first / rhs.first;
380 1350 : for (unsigned int i = 0; i < this->second.size(); i++) ret.second[i] = (this->second)[i] - rhs.second[i];
381 150 : return ret;
382 : }
383 4985 : Unit operator ^ (const int &rhs)
384 : {
385 4985 : Unit ret;
386 4985 : ret.first = std::pow(this->first,rhs);
387 79760 : for (unsigned int i = 0; i < this->second.size(); i++)
388 : {
389 34895 : ret.second[i] = (this->second)[i]*rhs;
390 : }
391 4985 : return ret;
392 : }
393 :
394 : friend std::ostream& operator<<(std::ostream &out, const Unit &a)
395 : {
396 : out << a.normalized_value() << "_";
397 : out << a.normalized_unitstring();
398 : return out;
399 : }
400 :
401 700 : double normalized_value() const
402 : {
403 700 : double fac = 1.0;
404 11200 : for (unsigned int i = 0 ; i < this->second.size(); i++)
405 4900 : fac *= std::pow(normalization[i].second, this->second[i]);
406 :
407 700 : return this->first * fac;
408 : }
409 :
410 0 : std::string normalized_unitstring() const
411 : {
412 : // little utility to join strings with delimiter
413 0 : auto join = [] (std::vector<std::string>vec,std::string sep){
414 0 : std::ostringstream oss;
415 0 : for (size_t i = 0; i < vec.size(); ++i) {
416 0 : oss << vec[i];
417 0 : if (i < vec.size() - 1) {
418 0 : oss << sep;
419 : }
420 : }
421 0 : return oss.str();
422 0 : };
423 :
424 0 : std::vector<std::string> num, den;
425 0 : for (unsigned int i = 0 ; i < this->second.size(); i++)
426 : {
427 0 : if (this->second[i] > 1)
428 0 : num.push_back(normalization[i].first + "^" + std::to_string(this->second[i]));
429 0 : else if (this->second[i] == 1)
430 0 : num.push_back(normalization[i].first);
431 0 : else if (this->second[i] == -1)
432 0 : den.push_back(normalization[i].first);
433 0 : else if (this->second[i] < -1)
434 0 : den.push_back(normalization[i].first + "^" + std::to_string(std::abs(this->second[i])));
435 : }
436 :
437 0 : std::string ret;
438 0 : if (num.size()) ret = join(num,"*");
439 0 : else ret = "1";
440 0 : if (den.size()) ret += "/" + join(den,"/");
441 :
442 0 : return ret;
443 0 : }
444 :
445 :
446 :
447 :
448 :
449 :
450 : };
451 :
452 :
453 : inline std::map<std::string, int>
454 1983 : Unit::ParseUnitString(const std::string& input) {
455 1983 : std::map<std::string, int> result;
456 : enum class Op { Mul = 1, Div = -1 };
457 1983 : Op currentOp = Op::Mul;
458 1983 : std::string token;
459 1983 : size_t i = 0;
460 4985 : auto parse_unit = [&](const std::string& t, Op op) {
461 4985 : if (t.empty()) return;
462 :
463 4985 : size_t caret = t.find('^');
464 4985 : std::string base = t;
465 4985 : int exp = 1;
466 :
467 4985 : if (caret != std::string::npos) {
468 616 : base = t.substr(0, caret);
469 616 : std::string exp_str = t.substr(caret + 1);
470 : try {
471 616 : exp = std::stoi(exp_str);
472 0 : } catch (...) {
473 0 : throw std::runtime_error("Invalid exponent: " + exp_str);
474 0 : }
475 616 : }
476 4985 : int signed_exp = (op == Op::Mul ? +1 : -1) * exp;
477 4985 : result[base] += signed_exp;
478 4985 : };
479 :
480 13008 : while (i < input.size()) {
481 11025 : char c = input[i];
482 :
483 11025 : if (c == '*')
484 : {
485 1195 : parse_unit(token, currentOp);
486 1195 : token.clear();
487 1195 : currentOp = Op::Mul;
488 1195 : ++i;
489 : }
490 9830 : else if (c == '/')
491 : {
492 1807 : parse_unit(token, currentOp);
493 1807 : token.clear();
494 1807 : currentOp = Op::Div;
495 1807 : ++i;
496 : }
497 : else
498 : {
499 8023 : token += c;
500 8023 : ++i;
501 : }
502 : }
503 :
504 : // Parse final token
505 1983 : parse_unit(token, currentOp);
506 :
507 3966 : return result;
508 1983 : }
509 :
510 :
511 :
512 :
513 :
514 :
515 :
516 :
517 :
518 :
519 : #endif
|
