units.c 12.4 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
42
43
44
#include "const.h"
#include "error.h"
#include "units.h"

/**
45
46
 * @brief Initialises the UnitSystem structure with the constants given in
 * const.h
47
48
49
 * @param us The UnitSystem to initialize
 */

50
void initUnitSystem(struct UnitSystem* us) {
51
52
  us->UnitMass_in_cgs = const_unit_mass_in_cgs;
  us->UnitLength_in_cgs = const_unit_length_in_cgs;
53
54
  us->UnitTime_in_cgs = 1. / ((double)const_unit_velocity_in_cgs /
                              ((double)const_unit_length_in_cgs));
55
56
57
58
  us->UnitCurrent_in_cgs = 1.;
  us->UnitTemperature_in_cgs = 1.;
}

59
60
61
/**
 * @brief Returns the base unit conversion factor for a given unit system
 * @param us The UnitSystem used
62
 * @param baseUnit The base unit
63
 */
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
double getBaseUnit(struct UnitSystem* us, enum BaseUnits baseUnit) {
  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");
  }
79
  return 0.0;
80
81
82
83
}

/**
 * @brief Returns the base unit symbol
84
 * @param baseUnit The base unit
85
 */
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
const char* getBaseUnitSymbol(enum BaseUnits baseUnit) {
  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");
  }
101
  return "";
102
103
104
105
}

/**
 * @brief Returns the base unit symbol in the cgs system
106
 * @param baseUnit The base unit
107
 */
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
const char* getBaseUnitCGSSymbol(enum BaseUnits baseUnit) {
  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");
  }
123
  return "";
124
125
}

126
127
128
void getBaseUnitExponantsArray(float baseUnitsExp[5],
                               enum UnitConversionFactor unit) {
  switch (unit) {
129
130
131
    case UNIT_CONV_NO_UNITS:
      break;

132
133
134
    case UNIT_CONV_MASS:
      baseUnitsExp[UNIT_MASS] = 1.f;
      break;
135

136
137
138
    case UNIT_CONV_LENGTH:
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      break;
139

140
141
142
    case UNIT_CONV_TIME:
      baseUnitsExp[UNIT_TIME] = 1.f;
      break;
143

144
145
146
    case UNIT_CONV_FREQUENCY:
      baseUnitsExp[UNIT_TIME] = -1.f;
      break;
147

148
149
150
151
    case UNIT_CONV_DENSITY:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = -3.f;
      break;
152

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

158
159
160
161
    case UNIT_CONV_ACCELERATION:
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
162

163
164
165
166
167
    case UNIT_CONV_FORCE:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
168

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

175
176
177
178
    case UNIT_CONV_ENERGY_PER_UNIT_MASS:
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
179

180
181
182
183
184
    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;
185

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

192
193
194
195
196
    case UNIT_CONV_POWER:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -3.f;
      break;
197

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

204
    case UNIT_CONV_ELECTRIC_CHARGE:
205
206
207
      baseUnitsExp[UNIT_TIME] = 1.f;
      baseUnitsExp[UNIT_CURRENT] = 1.f;
      break;
208
209

    case UNIT_CONV_ELECTRIC_VOLTAGE:
210
211
212
213
214
215
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -3.f;
      baseUnitsExp[UNIT_CURRENT] = -1.f;
      break;

216
    case UNIT_CONV_ELECTRIC_CAPACITANCE:
217
218
219
220
221
      baseUnitsExp[UNIT_MASS] = -1.f;
      baseUnitsExp[UNIT_LENGTH] = -2.f;
      baseUnitsExp[UNIT_TIME] = 4;
      baseUnitsExp[UNIT_CURRENT] = 2.f;
      break;
222
223

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

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

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

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

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

    case UNIT_CONV_TEMPERATURE:
258
      baseUnitsExp[UNIT_TEMPERATURE] = 1.f;
259
  }
260
261
262
}

/**
263
264
 * @brief Returns the conversion factor for a given unit in the chosen unit
 * system
265
266
267
 * @param us The system of units in use
 * @param unit The unit to convert
 */
268
269
double conversionFactor(struct UnitSystem* us, enum UnitConversionFactor unit) {
  float baseUnitsExp[5] = {0.f};
270
271

  getBaseUnitExponantsArray(baseUnitsExp, unit);
272

273
274
275
276
277
278
279
280
  return generalConversionFactor(us, baseUnitsExp);
}

/**
 * @brief Returns the h factor exponentiation for a given unit
 * @param us The system of units in use
 * @param unit The unit to convert
 */
281
282
float hFactor(struct UnitSystem* us, enum UnitConversionFactor unit) {
  float baseUnitsExp[5] = {0.f};
283
284
285

  getBaseUnitExponantsArray(baseUnitsExp, unit);

286
  return generalhFactor(us, baseUnitsExp);
287
288
289
290
291
292
293
}

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

  getBaseUnitExponantsArray(baseUnitsExp, unit);

299
  return generalaFactor(us, baseUnitsExp);
300
301
}

302
/**
303
 * @brief Returns a string containing the exponents of the base units making up
304
 * the conversion factors
305
 */
306
307
308
void conversionString(char* buffer, struct UnitSystem* us,
                      enum UnitConversionFactor unit) {
  float baseUnitsExp[5] = {0.f};
309
310

  getBaseUnitExponantsArray(baseUnitsExp, unit);
311

312
313
314
  generalConversionString(buffer, us, baseUnitsExp);
}

315
/**
316
317
 * @brief Returns the conversion factor for a given unit (expressed in terms of
 * the 5 fundamental units) in the chosen unit system
318
 * @param us The unit system used
319
 * @param baseUnitsExponants The exponent of each base units required to form
320
 * the desired quantity. See conversionFactor() for a working example
321
 */
322
323
double generalConversionFactor(struct UnitSystem* us,
                               float baseUnitsExponants[5]) {
324
325
326
  double factor = 1.;
  int i;

327
328
329
330
  for (i = 0; i < 5; ++i)
    if (baseUnitsExponants[i] != 0)
      factor *= pow(getBaseUnit(us, i), baseUnitsExponants[i]);
  return factor;
331
332
333
}

/**
334
335
 * @brief Returns the h factor exponentiation for a given unit (expressed in
 * terms of the 5 fundamental units)
336
 * @param us The unit system used
337
 * @param baseUnitsExponants The exponent of each base units required to form
338
 * the desired quantity. See conversionFactor() for a working example
339
 */
340
float generalhFactor(struct UnitSystem* us, float baseUnitsExponants[5]) {
341
  float factor_exp = 0.f;
342

343
344
345
  factor_exp += -baseUnitsExponants[UNIT_MASS];
  factor_exp += -baseUnitsExponants[UNIT_LENGTH];
  factor_exp += -baseUnitsExponants[UNIT_TIME];
346

347
  return factor_exp;
348
349
350
}

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

360
  factor_exp += baseUnitsExponants[UNIT_LENGTH];
361
362

  return factor_exp;
363
}
364
365

/**
366
 * @brief Returns a string containing the exponents of the base units making up
367
368
369
 * 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)
370
371
 * @param us The UnitsSystem in use.
 * @param baseUnitsExponants The exponent of each base units required to form
372
 * the desired quantity. See conversionFactor() for a working example
373
 */
374
375
void generalConversionString(char* buffer, struct UnitSystem* us,
                             float baseUnitsExponants[5]) {
376
377
378
379
380
  char temp[14];
  double a_exp = generalaFactor(us, baseUnitsExponants);
  double h_exp = generalhFactor(us, baseUnitsExponants);
  int i;

381
382
  /* Check whether we are unitless or not */
  char isAllNonZero = 1;
383
384
  for (i = 0; i < 5; ++i)
    if (baseUnitsExponants[i] != 0.) isAllNonZero = 0;
385

386
387
388
389
  if (isAllNonZero) {
    sprintf(buffer, "[ - ] ");
    return;
  }
390

391
  /* Add a-factor */
392
  if (a_exp == 0)
393
    sprintf(buffer, " ");
394
  else if (a_exp == 1)
395
    sprintf(buffer, "a ");
396
397
  else if (remainder(a_exp, 1.) == 0)
    sprintf(buffer, "a^%d ", (int)a_exp);
398
399
400
401
  else
    sprintf(buffer, "a^%7.4f ", a_exp);

  /* Add h-factor */
402
  if (h_exp == 0)
403
    sprintf(temp, " ");
404
  else if (h_exp == 1)
405
    sprintf(temp, "h ");
406
407
  else if (remainder(h_exp, 1.) == 0)
    sprintf(temp, "h^%d ", (int)h_exp);
408
409
410
411
412
  else
    sprintf(temp, "h^%7.4f ", h_exp);
  strncat(buffer, temp, 12);

  /* Add conversion units */
413
414
415
416
417
418
419
420
421
422
423
424
425
  for (i = 0; i < 5; ++i)
    if (baseUnitsExponants[i] != 0) {
      if (baseUnitsExponants[i] == 0.)
        sprintf(temp, " ");
      else if (baseUnitsExponants[i] == 1.)
        sprintf(temp, "%s ", getBaseUnitSymbol(i));
      else if (remainder(baseUnitsExponants[i], 1.) == 0)
        sprintf(temp, "%s^%d ", getBaseUnitSymbol(i),
                (int)baseUnitsExponants[i]);
      else
        sprintf(temp, "%s^%7.4f ", getBaseUnitSymbol(i), baseUnitsExponants[i]);
      strncat(buffer, temp, 12);
    }
426
427

  /* Add CGS units */
428
  strncat(buffer, " [ ", 3);
429
430
431
432
433
434
435
436
437
438
439
440
441
442

  for (i = 0; i < 5; ++i) {
    if (baseUnitsExponants[i] != 0) {
      if (baseUnitsExponants[i] == 0.)
        continue;
      else if (baseUnitsExponants[i] == 1.)
        sprintf(temp, "%s ", getBaseUnitCGSSymbol(i));
      else if (remainder(baseUnitsExponants[i], 1.) == 0)
        sprintf(temp, "%s^%d ", getBaseUnitCGSSymbol(i),
                (int)baseUnitsExponants[i]);
      else
        sprintf(temp, "%s^%7.4f ", getBaseUnitCGSSymbol(i),
                baseUnitsExponants[i]);
      strncat(buffer, temp, 12);
443
    }
444
445
  }

446
447
  strncat(buffer, "]", 2);
}