units.c 13.3 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
 */
50
void units_init(struct UnitSystem* us, const struct swift_params* params) {
51

52
53
54
55
56
57
  us->UnitMass_in_cgs =
      parser_get_param_double(params, "UnitSystem:UnitMass_in_cgs");
  us->UnitLength_in_cgs =
      parser_get_param_double(params, "UnitSystem:UnitLength_in_cgs");
  const double unitVelocity =
      parser_get_param_double(params, "UnitSystem:UnitVelocity_in_cgs");
58
  us->UnitTime_in_cgs = us->UnitLength_in_cgs / unitVelocity;
59
60
61
62
  us->UnitCurrent_in_cgs =
      parser_get_param_double(params, "UnitSystem:UnitCurrent_in_cgs");
  us->UnitTemperature_in_cgs =
      parser_get_param_double(params, "UnitSystem:UnitTemp_in_cgs");
63
64
}

65
66
67
/**
 * @brief Returns the base unit conversion factor for a given unit system
 * @param us The UnitSystem used
68
 * @param baseUnit The base unit
69
 */
70
71
double units_get_base_unit(const struct UnitSystem* us,
                           enum BaseUnits baseUnit) {
72
73
74
75
76
77
78
79
80
81
82
83
84
85
  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");
  }
86
  return 0.0;
87
88
89
90
}

/**
 * @brief Returns the base unit symbol
91
 * @param baseUnit The base unit
92
 */
93
const char* units_get_base_unit_symbol(enum BaseUnits baseUnit) {
94
95
96
97
98
99
100
101
102
103
104
105
106
107
  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");
  }
108
  return "";
109
110
111
112
}

/**
 * @brief Returns the base unit symbol in the cgs system
113
 * @param baseUnit The base unit
114
 */
115
const char* units_get_base_unit_CGS_symbol(enum BaseUnits baseUnit) {
116
117
118
119
120
121
122
123
124
125
126
127
128
129
  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");
  }
130
  return "";
131
132
}

133
134
void units_get_base_unit_exponants_array(float baseUnitsExp[5],
                                         enum UnitConversionFactor unit) {
135
  switch (unit) {
136
137
138
    case UNIT_CONV_NO_UNITS:
      break;

139
140
141
    case UNIT_CONV_MASS:
      baseUnitsExp[UNIT_MASS] = 1.f;
      break;
142

143
144
145
    case UNIT_CONV_LENGTH:
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      break;
146

147
148
149
    case UNIT_CONV_TIME:
      baseUnitsExp[UNIT_TIME] = 1.f;
      break;
150

151
152
153
    case UNIT_CONV_FREQUENCY:
      baseUnitsExp[UNIT_TIME] = -1.f;
      break;
154

155
156
157
158
    case UNIT_CONV_DENSITY:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = -3.f;
      break;
159

160
161
162
163
    case UNIT_CONV_SPEED:
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      baseUnitsExp[UNIT_TIME] = -1.f;
      break;
164

165
166
167
168
    case UNIT_CONV_ACCELERATION:
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
169

170
171
172
173
174
    case UNIT_CONV_FORCE:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
175

176
177
178
179
180
    case UNIT_CONV_ENERGY:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
181

182
183
184
185
    case UNIT_CONV_ENERGY_PER_UNIT_MASS:
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
186

187
188
189
190
191
    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;
192

193
194
195
196
197
    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;
198

199
200
201
202
203
    case UNIT_CONV_POWER:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -3.f;
      break;
204

205
206
207
208
209
    case UNIT_CONV_PRESSURE:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = -1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
210

211
    case UNIT_CONV_ELECTRIC_CHARGE:
212
213
214
      baseUnitsExp[UNIT_TIME] = 1.f;
      baseUnitsExp[UNIT_CURRENT] = 1.f;
      break;
215
216

    case UNIT_CONV_ELECTRIC_VOLTAGE:
217
218
219
220
221
222
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -3.f;
      baseUnitsExp[UNIT_CURRENT] = -1.f;
      break;

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

    case UNIT_CONV_ELECTRIC_RESISTANCE:
231
232
233
234
235
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -3.f;
      baseUnitsExp[UNIT_CURRENT] = -2.f;
      break;
236
237

    case UNIT_CONV_ELECTRIC_CONDUCTANCE:
238
239
240
241
242
243
      baseUnitsExp[UNIT_MASS] = -1.f;
      baseUnitsExp[UNIT_LENGTH] = -2.f;
      baseUnitsExp[UNIT_TIME] = 3.f;
      baseUnitsExp[UNIT_CURRENT] = 2.f;
      break;

244
    case UNIT_CONV_MAGNETIC_FLUX:
245
246
247
248
249
250
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      baseUnitsExp[UNIT_CURRENT] = -1.f;
      break;

251
    case UNIT_CONV_MAGNETIC_FIELD:
252
253
254
255
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      baseUnitsExp[UNIT_CURRENT] = -1.f;
      break;
256
257

    case UNIT_CONV_MAGNETIC_INDUCTANCE:
258
259
260
261
262
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      baseUnitsExp[UNIT_CURRENT] = -2.f;
      break;
263
264

    case UNIT_CONV_TEMPERATURE:
265
      baseUnitsExp[UNIT_TEMPERATURE] = 1.f;
266
  }
267
268
269
}

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

279
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
280

281
  return units_general_conversion_factor(us, baseUnitsExp);
282
283
284
285
286
287
288
}

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

293
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
294

295
  return units_general_h_factor(us, baseUnitsExp);
296
297
298
299
300
301
302
}

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

307
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
308

309
  return units_general_a_factor(us, baseUnitsExp);
310
311
}

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

320
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
321

322
  units_general_conversion_string(buffer, us, baseUnitsExp);
323
324
}

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

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

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

354
355
356
  factor_exp += -baseUnitsExponants[UNIT_MASS];
  factor_exp += -baseUnitsExponants[UNIT_LENGTH];
  factor_exp += -baseUnitsExponants[UNIT_TIME];
357

358
  return factor_exp;
359
360
361
}

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

372
  factor_exp += baseUnitsExponants[UNIT_LENGTH];
373
374

  return factor_exp;
375
}
376
377

/**
378
 * @brief Returns a string containing the exponents of the base units making up
379
380
381
 * 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)
382
383
 * @param us The UnitsSystem in use.
 * @param baseUnitsExponants The exponent of each base units required to form
384
 * the desired quantity. See conversionFactor() for a working example
385
 */
386
void units_general_conversion_string(char* buffer, const struct UnitSystem* us,
387
                                     const float baseUnitsExponants[5]) {
388
  char temp[14];
389
390
  double a_exp = units_general_a_factor(us, baseUnitsExponants);
  double h_exp = units_general_h_factor(us, baseUnitsExponants);
391
392
  int i;

393
394
  /* Check whether we are unitless or not */
  char isAllNonZero = 1;
395
396
  for (i = 0; i < 5; ++i)
    if (baseUnitsExponants[i] != 0.) isAllNonZero = 0;
397

398
399
400
401
  if (isAllNonZero) {
    sprintf(buffer, "[ - ] ");
    return;
  }
402

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

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

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

  /* Add CGS units */
441
  strncat(buffer, " [ ", 3);
442
443
444
445
446
447

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

459
460
  strncat(buffer, "]", 2);
}