units.c 15.5 KB
Newer Older
1
2
/*******************************************************************************
 * This file is part of SWIFT.
3
 * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
4
 *                    Matthieu Schaller (matthieu.schaller@durham.ac.uk).
5
 *
6
7
8
9
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
10
 *
11
12
13
14
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
15
 *
16
17
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
18
 *
19
20
21
22
23
 ******************************************************************************/

/* Config parameters. */
#include "../config.h"

24
/* Some standard headers. */
25
26
#include <math.h>
#include <stddef.h>
27
28
29
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
30
31

/* MPI headers. */
32
33
34
#ifdef WITH_MPI
#include <mpi.h>
#endif
35

36
37
38
39
/* This object's header. */
#include "units.h"

/* Includes. */
40
41
#include "const.h"
#include "error.h"
42

43
/**
44
 * @brief Initialises the UnitSystem structure with the constants given in
45
46
47
48
 * rhe parameter file.
 *
 * @param us The UnitSystem to initialize.
 * @param params The parsed parameter file.
49
 * @param category The section of the parameter file to read from.
50
 */
51
52
void units_init(struct UnitSystem* us, const struct swift_params* params,
                const char* category) {
53

54
55
56
57
58
59
60
  char buffer[200];
  sprintf(buffer, "%s:UnitMass_in_cgs", category);
  us->UnitMass_in_cgs = parser_get_param_double(params, buffer);
  sprintf(buffer, "%s:UnitLength_in_cgs", category);
  us->UnitLength_in_cgs = parser_get_param_double(params, buffer);
  sprintf(buffer, "%s:UnitVelocity_in_cgs", category);
  const double unitVelocity = parser_get_param_double(params, buffer);
61
  us->UnitTime_in_cgs = us->UnitLength_in_cgs / unitVelocity;
62
63
64
65
  sprintf(buffer, "%s:UnitCurrent_in_cgs", category);
  us->UnitCurrent_in_cgs = parser_get_param_double(params, buffer);
  sprintf(buffer, "%s:UnitTemp_in_cgs", category);
  us->UnitTemperature_in_cgs = parser_get_param_double(params, buffer);
66
67
}

68
69
70
/**
 * @brief Returns the base unit conversion factor for a given unit system
 * @param us The UnitSystem used
71
 * @param baseUnit The base unit
72
 */
73
74
double units_get_base_unit(const struct UnitSystem* us,
                           enum BaseUnits baseUnit) {
75
76
77
78
79
80
81
82
83
84
85
86
87
88
  switch (baseUnit) {
    case UNIT_MASS:
      return us->UnitMass_in_cgs;
    case UNIT_LENGTH:
      return us->UnitLength_in_cgs;
    case UNIT_TIME:
      return us->UnitTime_in_cgs;
    case UNIT_CURRENT:
      return us->UnitCurrent_in_cgs;
    case UNIT_TEMPERATURE:
      return us->UnitTemperature_in_cgs;
    default:
      error("Invalid base Unit");
  }
89
  return 0.0;
90
91
92
}

/**
93
 * @brief Returns the base unit symbol used internally
94
 * @param baseUnit The base unit
95
 */
96
const char* units_get_base_unit_internal_symbol(enum BaseUnits baseUnit) {
97
98
99
100
101
102
103
104
105
106
107
108
109
110
  switch (baseUnit) {
    case UNIT_MASS:
      return "U_M";
    case UNIT_LENGTH:
      return "U_L";
    case UNIT_TIME:
      return "U_t";
    case UNIT_CURRENT:
      return "U_I";
    case UNIT_TEMPERATURE:
      return "U_T";
    default:
      error("Invalid base Unit");
  }
111
  return "";
112
113
114
115
}

/**
 * @brief Returns the base unit symbol in the cgs system
116
 * @param baseUnit The base unit
117
 */
118
const char* units_get_base_unit_cgs_symbol(enum BaseUnits baseUnit) {
119
120
121
122
123
124
125
126
127
128
129
130
131
132
  switch (baseUnit) {
    case UNIT_MASS:
      return "g";
    case UNIT_LENGTH:
      return "cm";
    case UNIT_TIME:
      return "s";
    case UNIT_CURRENT:
      return "A";
    case UNIT_TEMPERATURE:
      return "K";
    default:
      error("Invalid base Unit");
  }
133
  return "";
134
135
}

136
137
void units_get_base_unit_exponants_array(float baseUnitsExp[5],
                                         enum UnitConversionFactor unit) {
138
  switch (unit) {
139
140
141
    case UNIT_CONV_NO_UNITS:
      break;

142
143
144
    case UNIT_CONV_MASS:
      baseUnitsExp[UNIT_MASS] = 1.f;
      break;
145

146
147
148
    case UNIT_CONV_LENGTH:
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      break;
149

150
151
152
    case UNIT_CONV_TIME:
      baseUnitsExp[UNIT_TIME] = 1.f;
      break;
153

154
155
156
    case UNIT_CONV_FREQUENCY:
      baseUnitsExp[UNIT_TIME] = -1.f;
      break;
157

158
159
160
161
    case UNIT_CONV_DENSITY:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = -3.f;
      break;
162

163
164
165
166
    case UNIT_CONV_SPEED:
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      baseUnitsExp[UNIT_TIME] = -1.f;
      break;
167

168
169
170
171
    case UNIT_CONV_ACCELERATION:
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
172

173
174
175
176
177
    case UNIT_CONV_FORCE:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
178

179
180
181
182
183
    case UNIT_CONV_ENERGY:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
184

185
186
187
188
    case UNIT_CONV_ENERGY_PER_UNIT_MASS:
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
189

190
191
192
193
194
    case UNIT_CONV_ENTROPY:
      baseUnitsExp[UNIT_MASS] = 1.f - const_hydro_gamma;
      baseUnitsExp[UNIT_LENGTH] = 3.f * const_hydro_gamma - 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
195

196
197
198
199
200
    case UNIT_CONV_ENTROPY_PER_UNIT_MASS:
      baseUnitsExp[UNIT_MASS] = -const_hydro_gamma;
      baseUnitsExp[UNIT_LENGTH] = 3.f * const_hydro_gamma - 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
201

202
203
204
205
206
    case UNIT_CONV_POWER:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -3.f;
      break;
207

208
209
210
211
212
    case UNIT_CONV_PRESSURE:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = -1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
213

214
    case UNIT_CONV_ELECTRIC_CHARGE:
215
216
217
      baseUnitsExp[UNIT_TIME] = 1.f;
      baseUnitsExp[UNIT_CURRENT] = 1.f;
      break;
218
219

    case UNIT_CONV_ELECTRIC_VOLTAGE:
220
221
222
223
224
225
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -3.f;
      baseUnitsExp[UNIT_CURRENT] = -1.f;
      break;

226
    case UNIT_CONV_ELECTRIC_CAPACITANCE:
227
228
229
230
231
      baseUnitsExp[UNIT_MASS] = -1.f;
      baseUnitsExp[UNIT_LENGTH] = -2.f;
      baseUnitsExp[UNIT_TIME] = 4;
      baseUnitsExp[UNIT_CURRENT] = 2.f;
      break;
232
233

    case UNIT_CONV_ELECTRIC_RESISTANCE:
234
235
236
237
238
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -3.f;
      baseUnitsExp[UNIT_CURRENT] = -2.f;
      break;
239
240

    case UNIT_CONV_ELECTRIC_CONDUCTANCE:
241
242
243
244
245
246
      baseUnitsExp[UNIT_MASS] = -1.f;
      baseUnitsExp[UNIT_LENGTH] = -2.f;
      baseUnitsExp[UNIT_TIME] = 3.f;
      baseUnitsExp[UNIT_CURRENT] = 2.f;
      break;

247
    case UNIT_CONV_MAGNETIC_FLUX:
248
249
250
251
252
253
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      baseUnitsExp[UNIT_CURRENT] = -1.f;
      break;

254
    case UNIT_CONV_MAGNETIC_FIELD:
255
256
257
258
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      baseUnitsExp[UNIT_CURRENT] = -1.f;
      break;
259
260

    case UNIT_CONV_MAGNETIC_INDUCTANCE:
261
262
263
264
265
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      baseUnitsExp[UNIT_CURRENT] = -2.f;
      break;
266
267

    case UNIT_CONV_TEMPERATURE:
268
      baseUnitsExp[UNIT_TEMPERATURE] = 1.f;
269
  }
270
271
272
}

/**
273
274
 * @brief Returns the conversion factor for a given unit in the chosen unit
 * system
275
276
277
 * @param us The system of units in use
 * @param unit The unit to convert
 */
278
279
double units_cgs_conversion_factor(const struct UnitSystem* us,
                                   enum UnitConversionFactor unit) {
280
  float baseUnitsExp[5] = {0.f};
281

282
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
283

284
  return units_general_cgs_conversion_factor(us, baseUnitsExp);
285
286
287
288
289
290
291
}

/**
 * @brief Returns the h factor exponentiation for a given unit
 * @param us The system of units in use
 * @param unit The unit to convert
 */
292
293
float units_h_factor(const struct UnitSystem* us,
                     enum UnitConversionFactor unit) {
294
  float baseUnitsExp[5] = {0.f};
295

296
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
297

298
  return units_general_h_factor(us, baseUnitsExp);
299
300
301
302
303
304
305
}

/**
 * @brief Returns the scaling factor exponentiation for a given unit
 * @param us The system of units in use
 * @param unit The unit to convert
 */
306
307
float units_a_factor(const struct UnitSystem* us,
                     enum UnitConversionFactor unit) {
308
  float baseUnitsExp[5] = {0.f};
309

310
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
311

312
  return units_general_a_factor(us, baseUnitsExp);
313
314
}

315
/**
316
 * @brief Returns a string containing the exponents of the base units making up
317
 * the conversion factors
318
 */
319
320
void units_cgs_conversion_string(char* buffer, const struct UnitSystem* us,
                                 enum UnitConversionFactor unit) {
321
  float baseUnitsExp[5] = {0.f};
322

323
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
324

325
  units_general_cgs_conversion_string(buffer, us, baseUnitsExp);
326
327
}

328
/**
329
330
 * @brief Returns the conversion factor for a given unit (expressed in terms of
 * the 5 fundamental units) in the chosen unit system
331
 * @param us The unit system used
332
 * @param baseUnitsExponants The exponent of each base units required to form
333
 * the desired quantity. See conversionFactor() for a working example
334
 */
335
336
double units_general_cgs_conversion_factor(const struct UnitSystem* us,
                                           const float baseUnitsExponants[5]) {
337
338
339
  double factor = 1.;
  int i;

340
341
  for (i = 0; i < 5; ++i)
    if (baseUnitsExponants[i] != 0)
342
      factor *= pow(units_get_base_unit(us, i), baseUnitsExponants[i]);
343
  return factor;
344
345
346
}

/**
347
348
 * @brief Returns the h factor exponentiation for a given unit (expressed in
 * terms of the 5 fundamental units)
349
 * @param us The unit system used
350
 * @param baseUnitsExponants The exponent of each base units required to form
351
 * the desired quantity. See conversionFactor() for a working example
352
 */
353
float units_general_h_factor(const struct UnitSystem* us,
354
                             const float baseUnitsExponants[5]) {
355
  float factor_exp = 0.f;
356

357
358
359
  factor_exp += -baseUnitsExponants[UNIT_MASS];
  factor_exp += -baseUnitsExponants[UNIT_LENGTH];
  factor_exp += -baseUnitsExponants[UNIT_TIME];
360

361
  return factor_exp;
362
363
364
}

/**
365
366
 * @brief Returns the scaling factor exponentiation for a given unit (expressed
 * in terms of the 5 fundamental units)
367
 * @param us The unit system used
368
 * @param baseUnitsExponants The exponent of each base units required to form
369
 * the desired quantity. See conversionFactor() for a working example
370
 */
371
float units_general_a_factor(const struct UnitSystem* us,
372
                             const float baseUnitsExponants[5]) {
373
  float factor_exp = 0.f;
374

375
  factor_exp += baseUnitsExponants[UNIT_LENGTH];
376
377

  return factor_exp;
378
}
379
380

/**
381
 * @brief Returns a string containing the exponents of the base units making up
382
383
384
 * the conversion factors (expressed in terms of the 5 fundamental units)
 * @param buffer The buffer in which to write (The buffer must be long enough,
 * 140 chars at most)
385
386
 * @param us The UnitsSystem in use.
 * @param baseUnitsExponants The exponent of each base units required to form
387
 * the desired quantity. See conversionFactor() for a working example
388
 */
389
390
391
void units_general_cgs_conversion_string(char* buffer,
                                         const struct UnitSystem* us,
                                         const float baseUnitsExponants[5]) {
392
  char temp[14];
393
394
  double a_exp = units_general_a_factor(us, baseUnitsExponants);
  double h_exp = units_general_h_factor(us, baseUnitsExponants);
395
396
  int i;

397
398
  /* Check whether we are unitless or not */
  char isAllNonZero = 1;
399
400
  for (i = 0; i < 5; ++i)
    if (baseUnitsExponants[i] != 0.) isAllNonZero = 0;
401

402
403
404
405
  if (isAllNonZero) {
    sprintf(buffer, "[ - ] ");
    return;
  }
406

407
  /* Add a-factor */
408
  if (a_exp == 0)
409
    sprintf(buffer, " ");
410
  else if (a_exp == 1)
411
    sprintf(buffer, "a ");
412
413
  else if (remainder(a_exp, 1.) == 0)
    sprintf(buffer, "a^%d ", (int)a_exp);
414
415
416
417
  else
    sprintf(buffer, "a^%7.4f ", a_exp);

  /* Add h-factor */
418
  if (h_exp == 0)
419
    sprintf(temp, " ");
420
  else if (h_exp == 1)
421
    sprintf(temp, "h ");
422
423
  else if (remainder(h_exp, 1.) == 0)
    sprintf(temp, "h^%d ", (int)h_exp);
424
425
426
427
428
  else
    sprintf(temp, "h^%7.4f ", h_exp);
  strncat(buffer, temp, 12);

  /* Add conversion units */
429
430
431
432
433
  for (i = 0; i < 5; ++i)
    if (baseUnitsExponants[i] != 0) {
      if (baseUnitsExponants[i] == 0.)
        sprintf(temp, " ");
      else if (baseUnitsExponants[i] == 1.)
434
        sprintf(temp, "%s ", units_get_base_unit_internal_symbol(i));
435
      else if (remainder(baseUnitsExponants[i], 1.) == 0)
436
        sprintf(temp, "%s^%d ", units_get_base_unit_internal_symbol(i),
437
438
                (int)baseUnitsExponants[i]);
      else
439
        sprintf(temp, "%s^%7.4f ", units_get_base_unit_internal_symbol(i),
440
                baseUnitsExponants[i]);
441
442
      strncat(buffer, temp, 12);
    }
443
444

  /* Add CGS units */
445
  strncat(buffer, " [ ", 3);
446
447
448
449
450
451

  for (i = 0; i < 5; ++i) {
    if (baseUnitsExponants[i] != 0) {
      if (baseUnitsExponants[i] == 0.)
        continue;
      else if (baseUnitsExponants[i] == 1.)
452
        sprintf(temp, "%s ", units_get_base_unit_cgs_symbol(i));
453
      else if (remainder(baseUnitsExponants[i], 1.) == 0)
454
        sprintf(temp, "%s^%d ", units_get_base_unit_cgs_symbol(i),
455
456
                (int)baseUnitsExponants[i]);
      else
457
        sprintf(temp, "%s^%7.4f ", units_get_base_unit_cgs_symbol(i),
458
459
                baseUnitsExponants[i]);
      strncat(buffer, temp, 12);
460
    }
461
462
  }

463
464
  strncat(buffer, "]", 2);
}
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508

/**
 * @brief Are the two unit systems equal ?
 *
 * @param a The First #UnitSystem
 * @param b The second #UnitSystem
 * @return 1 if the systems are the same, 0 otherwise
 */
int units_are_equal(const struct UnitSystem* a, const struct UnitSystem* b) {

  if (a->UnitMass_in_cgs != b->UnitMass_in_cgs) return 0;
  if (a->UnitLength_in_cgs != b->UnitLength_in_cgs) return 0;
  if (a->UnitTime_in_cgs != b->UnitTime_in_cgs) return 0;
  if (a->UnitCurrent_in_cgs != b->UnitCurrent_in_cgs) return 0;
  if (a->UnitTemperature_in_cgs != b->UnitTemperature_in_cgs) return 0;

  return 1;
}

/**
 * @brief Return the unit conversion factor between two systems
 *
 * @param from The #UnitSystem we are converting from
 * @param to The #UnitSystem we are converting to
 * @param baseUnitsExponants The exponent of each base units required to form
 * the desired quantity. See conversionFactor() for a working example
 */
double units_general_conversion_factor(const struct UnitSystem* from,
                                       const struct UnitSystem* to,
                                       const float baseUnitsExponants[5]) {

  const double from_cgs =
      units_general_cgs_conversion_factor(from, baseUnitsExponants);
  const double to_cgs =
      units_general_cgs_conversion_factor(to, baseUnitsExponants);

  return from_cgs / to_cgs;
}

/**
 * @brief Return the unit conversion factor between two systems
 *
 * @param from The #UnitSystem we are converting from
 * @param to The #UnitSystem we are converting to
Matthieu Schaller's avatar
Matthieu Schaller committed
509
510
511
 * @param unit The unit we are converting
 * 
 * @return The conversion factor
512
513
514
515
516
517
518
519
520
521
522
 */
double units_conversion_factor(const struct UnitSystem* from,
                               const struct UnitSystem* to,
                               enum UnitConversionFactor unit) {

  float baseUnitsExp[5] = {0.f};

  units_get_base_unit_exponants_array(baseUnitsExp, unit);

  return units_general_conversion_factor(from, to, baseUnitsExp);
}