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 : // Angles
67 : {"rad", {1.0, "1"}},
68 : {"deg", {0.01745329251, "rad" }},
69 :
70 : // Times
71 : {"ms", {1e-3, "s" }},
72 : {"us", {1e-6, "s" }},
73 : {"ns", {1e-9, "s" }},
74 : {"ps", {1e-12, "s" }},
75 : {"min", {60.0, "s" }},
76 : {"hr", {60.0, "min" }},
77 : {"day", {24.0, "hr" }},
78 : {"week", {7.0, "day" }},
79 : {"fortnight", {14.0, "day" }},
80 : {"year", {365.25, "day" }},
81 : {"decade", {10.0, "year" }},
82 :
83 : // Mass
84 : {"g", {1e-3, "kg" }},
85 : {"mg", {1e-6, "kg" }},
86 : {"ug", {1e-9, "kg" }},
87 : {"tonne", {1000.0, "kg" }},
88 : {"lb", {0.45359237, "kg" }},
89 : {"oz", {0.0283495, "kg" }},
90 : {"slug", {14.5939, "kg" }},
91 :
92 : // Speed
93 : {"mph", {1.0, "mi/hr" }},
94 : {"kph", {1.0, "km/hr" }},
95 : {"mps", {1.0, "m/s" }},
96 : {"fps", {1.0, "ft/s" }},
97 : {"knot", {1852.0 / 3600.0, "m/s"}},
98 :
99 : // Area
100 : {"ha", {1e4, "m^2" }},
101 : {"acre", {4046.85642, "m^2" }},
102 : {"sqft", {1.0, "ft^2" }},
103 : {"sqin", {1.0, "in^2" }},
104 : {"sqkm", {1.0, "km^2" }},
105 : {"sqmi", {1.0, "mi^2" }},
106 :
107 : // Volume
108 : {"L", {1e-3, "m^3" }},
109 : {"mL", {1e-6, "m^3" }},
110 : {"cc", {1e-6, "m^3" }},
111 : {"gal", {3.78541e-3, "m^3" }},
112 : {"qt", {0.946353e-3, "m^3" }},
113 : {"pt", {0.473176e-3, "m^3" }},
114 : {"cup", {0.24e-3, "m^3" }},
115 : {"floz", {29.5735e-6, "m^3" }},
116 : {"tbsp", {14.7868e-6, "m^3" }},
117 : {"tsp", {4.92892e-6, "m^3" }},
118 :
119 : // Force
120 : {"N", {1.0, "kg*m/s^2" }},
121 : {"kN", {1e3, "N" }},
122 : {"lbf", {4.44822162, "N" }},
123 : {"dyn", {1e-5, "N" }},
124 :
125 : // Energy
126 : {"J", {1.0, "N*m" }},
127 : {"mJ", {1e-3, "J" }},
128 : {"kJ", {1e3, "J" }},
129 : {"MJ", {1e6, "J" }},
130 : {"eV", {1.602176634e-19, "J"}},
131 : {"cal", {4.184, "J" }},
132 : {"kcal", {4184.0, "J" }},
133 : {"BTU", {1055.06, "J" }},
134 : {"ft*lbf", {1.35582, "J" }},
135 :
136 : // Power
137 : {"W", {1.0, "J/s" }},
138 : {"mW", {1e-3, "W" }},
139 : {"kW", {1e3, "W" }},
140 : {"MW", {1e6, "W" }},
141 : {"hp", {745.7, "W" }},
142 :
143 : // Pressure
144 : {"Pa", {1.0, "N/m^2" }},
145 : {"kPa", {1e3, "Pa" }},
146 : {"MPa", {1e6, "Pa" }},
147 : {"GPa", {1e9, "Pa" }},
148 : {"bar", {1e5, "Pa" }},
149 : {"atm", {101325.0,"Pa" }},
150 : {"psi", {6894.76, "Pa" }},
151 : {"mmHg", {133.322, "Pa" }},
152 : {"torr", {133.322, "Pa" }},
153 :
154 : // Temperature differences (not absolute temperatures)
155 : {"K", {1.0, "K" }},
156 : {"degR", {5.0/9.0, "K" }},
157 :
158 : // Charge
159 : {"C", {1.0, "A*s" }},
160 : {"mAh", {3.6, "C" }},
161 :
162 : // Voltage
163 : {"V", {1.0, "J/C" }},
164 :
165 : // Capacitance
166 : {"F", {1.0, "C/V" }},
167 :
168 : // Resistance
169 : {"ohm", {1.0, "V/A" }},
170 : {"kOhm", {1e3, "ohm" }},
171 : {"MOhm", {1e6, "ohm" }},
172 :
173 : // Conductance
174 : {"S", {1.0, "1/ohm" }},
175 :
176 : // Magnetic
177 : {"T", {1.0, "kg/(A*s^2)"}},
178 : {"G", {1e-4, "T" }}, // gauss
179 :
180 : // Frequency
181 : {"Hz", {1.0, "1/s" }},
182 : {"kHz", {1e3, "Hz" }},
183 : {"MHz", {1e6, "Hz" }},
184 : {"GHz", {1e9, "Hz" }},
185 : };
186 :
187 :
188 : //
189 : // These static factory functions return unit types.
190 : // They are used to check unit compliance, mainly at the parsing stage.
191 : //
192 : // Primary
193 4225 : static Unit Less() { return Unit{1.0, {0,0,0,0,0,0,0}}; } // 1
194 282 : static Unit Length() { return Unit{1.0, {1,0,0,0,0,0,0}}; } // m
195 409 : static Unit Time() { return Unit{1.0, {0,1,0,0,0,0,0}}; } // s
196 49 : static Unit Mass() { return Unit{1.0, {0,0,1,0,0,0,0}}; } // kg
197 25 : static Unit Temperature() { return Unit{1.0, {0,0,0,1,0,0,0}}; } // K
198 : static Unit Current() { return Unit{1.0, {0,0,0,0,1,0,0}}; } // A
199 : static Unit Amount() { return Unit{1.0, {0,0,0,0,0,1,0}}; } // mol
200 : static Unit LuminousIntensity() { return Unit{1.0, {0,0,0,0,0,0,1}}; } // cd
201 :
202 : static Unit Angle() {return Unit::Less();}
203 :
204 : // Geometry and kinematics
205 24 : static Unit Area() { return Length()*Length(); }
206 7 : static Unit Volume() { return Length()*Length()*Length(); }
207 40 : static Unit Velocity() { return Length() / Time();}
208 38 : static Unit Acceleration(){ return Velocity() / Time(); }
209 : // Mechanics
210 38 : static Unit Force() { return Mass() * Acceleration(); } // kg·m/s² = N
211 : static Unit Momentum() { return Mass() * Velocity(); }
212 : static Unit Impulse() { return Force() * Time(); }
213 13 : static Unit Pressure() { return Force() / Area(); } // N/m² = Pa
214 25 : static Unit Energy() { return Force() * Length(); } // N·m = J
215 6 : static Unit Power() { return Energy() / Time(); } // J/s = W
216 7 : static Unit Density() { return Mass() / Volume(); }
217 : static Unit SpecificWeight() { return Force() / Volume(); }
218 : static Unit Work() { return Energy(); }
219 : // Thermodynamics
220 4 : static Unit SpecificHeatCapacity() {return Energy() / Mass() / Temperature() ;}
221 4 : static Unit ThermalConductivity() {return Power()/Length()/Temperature();}
222 2 : static Unit ThermalDiffusivity() {return ThermalConductivity() / Density() / SpecificHeatCapacity();}
223 : static Unit HeatFlux() { return Power() / Area(); }
224 : static Unit HeatTransferCoefficient() { return HeatFlux() / Temperature(); }
225 : // Electromagnetism
226 :
227 : static Unit Charge() { return Current() * Time(); } // A·s = C
228 : static Unit Voltage() { return Power() / Current(); } // W/A = V
229 : static Unit Resistance() { return Voltage() / Current(); } // V/A = Ω
230 : static Unit Capacitance() { return Charge() / Voltage(); } // C/V = F
231 : static Unit Conductance() { return Less() / Resistance(); } // S = 1/Ω
232 : static Unit MagneticFlux(){ return Voltage() * Time(); } // V·s = Wb
233 : static Unit MagneticField(){ return MagneticFlux() / Area(); } // Wb/m² = T
234 : static Unit Inductance() { return MagneticFlux() / Current(); } // Wb/A = H
235 : // Chemistry
236 : static Unit MolarMass() { return Mass() / Amount(); } // kg/mol
237 : static Unit Concentration(){ return Amount() / Volume(); } // mol/m³
238 : static Unit MolarEnergy() { return Energy() / Amount(); } // J/mol
239 : // Optics and radiation
240 : static Unit LuminousFlux() { return LuminousIntensity(); } // for scalar usage (lm)
241 : static Unit Illuminance() { return LuminousFlux() / Area(); } // lux = lm/m²
242 :
243 :
244 : //
245 : // This is where SYSTEM NORMALIZATION is stored.
246 : // Should usually be set once (at the beginning of the simulation) and then
247 : // left constant throughout the simulation.
248 : // All normalizations are in terms of standard SI units: m,s,kg,K,A,mol,cd.
249 : //
250 : inline static std::array<std::pair<std::string,double>, 7> normalization = {
251 : std::pair{"m", 1.0},
252 : std::pair{"s", 1.0},
253 : std::pair{"kg", 1.0},
254 : std::pair{"K", 1.0},
255 : std::pair{"A", 1.0},
256 : std::pair{"mol",1.0},
257 : std::pair{"cd", 1.0}
258 : };
259 :
260 : using std::pair<double, std::array<int,7>>::pair;
261 4405 : Unit(const std::pair<double, std::array<int, 7>>& p) : std::pair<double, std::array<int, 7>>(p) {}
262 : static std::map<std::string, int> ParseUnitString(const std::string& input);
263 :
264 :
265 913 : static Unit Parse(double val, std::string unitstring, bool verbose = false)
266 : {
267 913 : Unit ret = StringParse(unitstring, verbose);
268 913 : ret.first *= val;
269 913 : return ret;
270 : }
271 1519 : static Unit Parse(std::string unit, bool verbose = false)
272 : {
273 1519 : if (unit.size() == 0) return Unit::Less();
274 1496 : std::vector<std::string> split;
275 1496 : std::stringstream ss(unit);
276 1496 : std::string item;
277 3905 : while (std::getline(ss, item, '_')) {
278 2409 : split.push_back(item);
279 : }
280 1496 : if (split.size() == 2)
281 : {
282 913 : return Parse(std::stod(split[0]),split[1], verbose);
283 : }
284 583 : else if (split.size() == 1)
285 : {
286 : // Attempt to parse the string as if it's a unitless number
287 583 : double value = NAN;
288 :
289 : // auto first = unit.data();
290 : // auto last = unit.data() + unit.size();
291 : // auto result = std::from_chars(first, last, value); // apparently this fails on mac
292 :
293 583 : char *end = nullptr;
294 583 : value = std::strtod(unit.c_str(), &end);
295 :
296 : // If it is a untless number, then return a unitless result
297 : //if (result.ec == std::errc() && result.ptr == last) // fails on mac
298 583 : if (end != unit.data() && *end == '\0')
299 : {
300 581 : Unit ret = Unit::Less();
301 581 : ret.first = value;
302 581 : return ret;
303 : }
304 : // Otherwise, it must be a pure unit, so return that.
305 : else
306 : {
307 2 : return StringParse(unit);
308 : }
309 : }
310 : else
311 : {
312 0 : throw std::runtime_error("Invalid unit string "+unit);
313 : }
314 1496 : }
315 :
316 :
317 2553 : static Unit StringParse(std::string unitstring, bool verbose = false)
318 : {
319 2553 : Unit type = Unit::Less();
320 :
321 : // Get the map of base tokens with signed exponents from ParseUnitString
322 2553 : std::map<std::string,int> parsed_units = ParseUnitString(unitstring);
323 2553 : if (verbose)
324 : {
325 0 : for (auto parsed_unit : parsed_units)
326 : {
327 0 : std::cout << parsed_unit.first << " " << parsed_unit.second << std::endl;
328 0 : }
329 : }
330 :
331 8518 : for (auto &p : parsed_units)
332 : {
333 5965 : const std::string &token = p.first;
334 5965 : int exp = p.second;
335 :
336 5965 : if (Unit::base_units.find(token) != base_units.end())
337 : {
338 : // Base unit: apply exponent
339 4405 : Unit t = Unit(Unit::base_units.at(token)) ^ exp;
340 4405 : type = type * t;
341 : }
342 1560 : else if (compound.find(token) != compound.end())
343 : {
344 : // Compound unit: expand recursively
345 1560 : auto pair = compound.at(token);
346 1560 : double factor = pair.first;
347 1560 : const std::string &subunit_str = pair.second;
348 :
349 1560 : Unit t = StringParse(subunit_str); // recursive expansion
350 1560 : t = t ^ exp; // apply exponent to both factor and dimension
351 1560 : t.first *= std::pow(factor, exp); // scale numeric factor correctly
352 :
353 1560 : type = type * t;
354 1560 : }
355 : else
356 : {
357 0 : throw std::runtime_error("Unknown unit token: " + token);
358 : }
359 : }
360 :
361 5106 : return type;
362 2553 : }
363 :
364 :
365 39 : static void setLengthUnit(std::string unit)
366 : {
367 39 : Unit parsed = StringParse(unit);
368 39 : if (parsed.second != Length().second) throw std::runtime_error("Not a unit of length: " + unit);
369 39 : normalization[int(Type::Length)].first = unit;
370 39 : normalization[int(Type::Length)].second = 1.0 / parsed.first;
371 39 : }
372 39 : static void setTimeUnit(std::string unit)
373 : {
374 39 : Unit parsed = StringParse(unit);
375 39 : if (parsed.second != Time().second) throw std::runtime_error("Not a unit of time: " + unit);
376 39 : normalization[int(Type::Time)].first = unit;
377 39 : normalization[int(Type::Time)].second = 1.0 / parsed.first;
378 39 : }
379 :
380 1298 : bool isType(const Unit &test) const
381 : {
382 1298 : return this->second == test.second;
383 : }
384 :
385 6066 : Unit operator * (const Unit &rhs)
386 : {
387 6066 : Unit ret;
388 6066 : ret.first = this->first * rhs.first;
389 54594 : for (unsigned int i = 0; i < this->second.size(); i++) ret.second[i] = (this->second)[i] + rhs.second[i];
390 6066 : return ret;
391 : }
392 : friend Unit operator*(double lhs, const Unit& rhs) {
393 : Unit ret;
394 : ret.first = lhs * rhs.first;
395 : ret.second = rhs.second;
396 : return ret;
397 : }
398 3 : friend Unit operator / (double lhs, const Unit& rhs) {
399 3 : Unit ret = Unit::Less() / rhs;
400 3 : ret.first = lhs * rhs.first;
401 3 : return ret;
402 : }
403 150 : Unit operator / (const Unit &rhs)
404 : {
405 150 : Unit ret;
406 150 : ret.first = this->first / rhs.first;
407 1350 : for (unsigned int i = 0; i < this->second.size(); i++) ret.second[i] = (this->second)[i] - rhs.second[i];
408 150 : return ret;
409 : }
410 5965 : Unit operator ^ (const int &rhs)
411 : {
412 5965 : Unit ret;
413 5965 : ret.first = std::pow(this->first,rhs);
414 95440 : for (unsigned int i = 0; i < this->second.size(); i++)
415 : {
416 41755 : ret.second[i] = (this->second)[i]*rhs;
417 : }
418 5965 : return ret;
419 : }
420 :
421 0 : friend std::ostream& operator<<(std::ostream &out, const Unit &a)
422 : {
423 0 : out << a.normalized_value() << "_";
424 0 : out << a.normalized_unitstring();
425 0 : return out;
426 : }
427 :
428 936 : double normalized_value() const
429 : {
430 936 : double fac = 1.0;
431 14976 : for (unsigned int i = 0 ; i < this->second.size(); i++)
432 6552 : fac *= std::pow(normalization[i].second, this->second[i]);
433 :
434 936 : return this->first * fac;
435 : }
436 :
437 0 : std::string normalized_unitstring() const
438 : {
439 : // little utility to join strings with delimiter
440 0 : auto join = [] (std::vector<std::string>vec,std::string sep){
441 0 : std::ostringstream oss;
442 0 : for (size_t i = 0; i < vec.size(); ++i) {
443 0 : oss << vec[i];
444 0 : if (i < vec.size() - 1) {
445 0 : oss << sep;
446 : }
447 : }
448 0 : return oss.str();
449 0 : };
450 :
451 0 : std::vector<std::string> num, den;
452 0 : for (unsigned int i = 0 ; i < this->second.size(); i++)
453 : {
454 0 : if (this->second[i] > 1)
455 0 : num.push_back(normalization[i].first + "^" + std::to_string(this->second[i]));
456 0 : else if (this->second[i] == 1)
457 0 : num.push_back(normalization[i].first);
458 0 : else if (this->second[i] == -1)
459 0 : den.push_back(normalization[i].first);
460 0 : else if (this->second[i] < -1)
461 0 : den.push_back(normalization[i].first + "^" + std::to_string(std::abs(this->second[i])));
462 : }
463 :
464 0 : std::string ret;
465 0 : if (num.size()) ret = join(num,"*");
466 0 : else ret = "1";
467 0 : if (den.size()) ret += "/" + join(den,"/");
468 :
469 0 : return ret;
470 0 : }
471 : };
472 :
473 :
474 : inline std::map<std::string, int>
475 2553 : Unit::ParseUnitString(const std::string& input) {
476 2553 : std::map<std::string, int> result;
477 : enum class Op { Mul = 1, Div = -1 };
478 2553 : Op currentOp = Op::Mul;
479 2553 : std::string token;
480 2553 : size_t i = 0;
481 5965 : auto parse_unit = [&](const std::string& t, Op op) {
482 5965 : if (t.empty()) return;
483 :
484 5965 : size_t caret = t.find('^');
485 5965 : std::string base = t;
486 5965 : int exp = 1;
487 :
488 5965 : if (caret != std::string::npos) {
489 820 : base = t.substr(0, caret);
490 820 : std::string exp_str = t.substr(caret + 1);
491 : try {
492 820 : exp = std::stoi(exp_str);
493 0 : } catch (...) {
494 0 : throw std::runtime_error("Invalid exponent: " + exp_str);
495 0 : }
496 820 : }
497 5965 : int signed_exp = (op == Op::Mul ? +1 : -1) * exp;
498 5965 : result[base] += signed_exp;
499 5965 : };
500 :
501 15702 : while (i < input.size()) {
502 13149 : char c = input[i];
503 :
504 13149 : if (c == '*')
505 : {
506 1355 : parse_unit(token, currentOp);
507 1355 : token.clear();
508 1355 : currentOp = Op::Mul;
509 1355 : ++i;
510 : }
511 11794 : else if (c == '/')
512 : {
513 2057 : parse_unit(token, currentOp);
514 2057 : token.clear();
515 2057 : currentOp = Op::Div;
516 2057 : ++i;
517 : }
518 : else
519 : {
520 9737 : token += c;
521 9737 : ++i;
522 : }
523 : }
524 :
525 : // Parse final token
526 2553 : parse_unit(token, currentOp);
527 :
528 5106 : return result;
529 2553 : }
530 :
531 :
532 :
533 :
534 :
535 :
536 :
537 :
538 :
539 :
540 : #endif
|
