units.c 17.7 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
#include "adiabatic_index.h"
41
42
#include "const.h"
#include "error.h"
43

44
45
46
47
48
49
50
51
52
53
54
55
56
57
/**
 * @brief Initialises the UnitSystem structure with CGS system
 *
 * @param us The UnitSystem to initialize
 */
void units_init_cgs(struct UnitSystem* us) {

  us->UnitMass_in_cgs = 1.;
  us->UnitLength_in_cgs = 1.;
  us->UnitTime_in_cgs = 1.;
  us->UnitCurrent_in_cgs = 1.;
  us->UnitTemperature_in_cgs = 1.;
}

58
/**
59
 * @brief Initialises the UnitSystem structure with the constants given in
60
61
62
63
 * rhe parameter file.
 *
 * @param us The UnitSystem to initialize.
 * @param params The parsed parameter file.
64
 * @param category The section of the parameter file to read from.
65
 */
66
67
void units_init(struct UnitSystem* us, const struct swift_params* params,
                const char* category) {
68

69
70
71
72
73
74
75
  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);
76
  us->UnitTime_in_cgs = us->UnitLength_in_cgs / unitVelocity;
77
78
79
80
  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);
81
82
}

83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
/**
 * @brief Initialises the UnitSystem structure with the constants given in
 * rhe parameter file. Uses a default if the values are not present in the file.
 *
 * @param us The UnitSystem to initialize.
 * @param params The parsed parameter file.
 * @param category The section of the parameter file to read from.
 * @param def The default unit system to copy from if required.
 */
void units_init_default(struct UnitSystem* us,
                        const struct swift_params* params, const char* category,
                        const struct UnitSystem* def) {

  if (!def) error("Default UnitSystem not allocated");

  char buffer[200];
  sprintf(buffer, "%s:UnitMass_in_cgs", category);
  us->UnitMass_in_cgs =
      parser_get_opt_param_double(params, buffer, def->UnitMass_in_cgs);
  sprintf(buffer, "%s:UnitLength_in_cgs", category);
  us->UnitLength_in_cgs =
      parser_get_opt_param_double(params, buffer, def->UnitLength_in_cgs);
  sprintf(buffer, "%s:UnitVelocity_in_cgs", category);
  const double defaultVelocity = def->UnitLength_in_cgs / def->UnitTime_in_cgs;
  const double unitVelocity =
      parser_get_opt_param_double(params, buffer, defaultVelocity);
  us->UnitTime_in_cgs = us->UnitLength_in_cgs / unitVelocity;
  sprintf(buffer, "%s:UnitCurrent_in_cgs", category);
  us->UnitCurrent_in_cgs =
      parser_get_opt_param_double(params, buffer, def->UnitCurrent_in_cgs);
  sprintf(buffer, "%s:UnitTemp_in_cgs", category);
  us->UnitTemperature_in_cgs =
      parser_get_opt_param_double(params, buffer, def->UnitTemperature_in_cgs);
}

118
119
120
/**
 * @brief Returns the base unit conversion factor for a given unit system
 * @param us The UnitSystem used
121
 * @param baseUnit The base unit
122
 */
123
124
double units_get_base_unit(const struct UnitSystem* us,
                           enum BaseUnits baseUnit) {
125
126
127
128
129
130
131
132
133
134
135
136
137
138
  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");
  }
139
  return 0.0;
140
141
142
}

/**
143
 * @brief Returns the base unit symbol used internally
144
 * @param baseUnit The base unit
145
 */
146
const char* units_get_base_unit_internal_symbol(enum BaseUnits baseUnit) {
147
148
149
150
151
152
153
154
155
156
157
158
159
160
  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");
  }
161
  return "";
162
163
164
165
}

/**
 * @brief Returns the base unit symbol in the cgs system
166
 * @param baseUnit The base unit
167
 */
168
const char* units_get_base_unit_cgs_symbol(enum BaseUnits baseUnit) {
169
170
171
172
173
174
175
176
177
178
179
180
181
182
  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");
  }
183
  return "";
184
185
}

186
187
void units_get_base_unit_exponants_array(float baseUnitsExp[5],
                                         enum UnitConversionFactor unit) {
188
  switch (unit) {
189
190
191
    case UNIT_CONV_NO_UNITS:
      break;

192
193
194
    case UNIT_CONV_MASS:
      baseUnitsExp[UNIT_MASS] = 1.f;
      break;
195

196
197
198
    case UNIT_CONV_LENGTH:
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      break;
199

200
201
202
    case UNIT_CONV_TIME:
      baseUnitsExp[UNIT_TIME] = 1.f;
      break;
203

204
205
206
    case UNIT_CONV_FREQUENCY:
      baseUnitsExp[UNIT_TIME] = -1.f;
      break;
207

208
209
210
211
    case UNIT_CONV_DENSITY:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = -3.f;
      break;
212

213
214
215
216
    case UNIT_CONV_SPEED:
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      baseUnitsExp[UNIT_TIME] = -1.f;
      break;
217

218
219
220
221
    case UNIT_CONV_ACCELERATION:
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
222

223
224
225
226
227
    case UNIT_CONV_FORCE:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
228

229
230
231
232
233
    case UNIT_CONV_ENERGY:
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
234

235
236
237
238
    case UNIT_CONV_ENERGY_PER_UNIT_MASS:
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
239

240
    case UNIT_CONV_ENTROPY:
241
242
      baseUnitsExp[UNIT_MASS] = 1.f - hydro_gamma;
      baseUnitsExp[UNIT_LENGTH] = 3.f * hydro_gamma - 1.f;
243
244
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
245

246
    case UNIT_CONV_ENTROPY_PER_UNIT_MASS:
247
248
      baseUnitsExp[UNIT_MASS] = -hydro_gamma;
      baseUnitsExp[UNIT_LENGTH] = 3.f * hydro_gamma - 1.f;
249
250
      baseUnitsExp[UNIT_TIME] = -2.f;
      break;
251

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

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

264
    case UNIT_CONV_ELECTRIC_CHARGE:
265
266
267
      baseUnitsExp[UNIT_TIME] = 1.f;
      baseUnitsExp[UNIT_CURRENT] = 1.f;
      break;
268
269

    case UNIT_CONV_ELECTRIC_VOLTAGE:
270
271
272
273
274
275
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -3.f;
      baseUnitsExp[UNIT_CURRENT] = -1.f;
      break;

276
    case UNIT_CONV_ELECTRIC_CAPACITANCE:
277
278
279
280
281
      baseUnitsExp[UNIT_MASS] = -1.f;
      baseUnitsExp[UNIT_LENGTH] = -2.f;
      baseUnitsExp[UNIT_TIME] = 4;
      baseUnitsExp[UNIT_CURRENT] = 2.f;
      break;
282
283

    case UNIT_CONV_ELECTRIC_RESISTANCE:
284
285
286
287
288
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -3.f;
      baseUnitsExp[UNIT_CURRENT] = -2.f;
      break;
289
290

    case UNIT_CONV_ELECTRIC_CONDUCTANCE:
291
292
293
294
295
296
      baseUnitsExp[UNIT_MASS] = -1.f;
      baseUnitsExp[UNIT_LENGTH] = -2.f;
      baseUnitsExp[UNIT_TIME] = 3.f;
      baseUnitsExp[UNIT_CURRENT] = 2.f;
      break;

297
    case UNIT_CONV_MAGNETIC_FLUX:
298
299
300
301
302
303
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      baseUnitsExp[UNIT_CURRENT] = -1.f;
      break;

304
    case UNIT_CONV_MAGNETIC_FIELD:
305
306
307
308
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      baseUnitsExp[UNIT_CURRENT] = -1.f;
      break;
309
310

    case UNIT_CONV_MAGNETIC_INDUCTANCE:
311
312
313
314
315
      baseUnitsExp[UNIT_MASS] = 1.f;
      baseUnitsExp[UNIT_LENGTH] = 2.f;
      baseUnitsExp[UNIT_TIME] = -2.f;
      baseUnitsExp[UNIT_CURRENT] = -2.f;
      break;
316
317

    case UNIT_CONV_TEMPERATURE:
318
      baseUnitsExp[UNIT_TEMPERATURE] = 1.f;
319
320
321
322
      break;

    case UNIT_CONV_VOLUME:
      baseUnitsExp[UNIT_LENGTH] = -3.f;
323
  }
324
325
326
}

/**
327
328
 * @brief Returns the conversion factor for a given unit in the chosen unit
 * system
329
330
331
 * @param us The system of units in use
 * @param unit The unit to convert
 */
332
333
double units_cgs_conversion_factor(const struct UnitSystem* us,
                                   enum UnitConversionFactor unit) {
334
  float baseUnitsExp[5] = {0.f};
335

336
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
337

338
  return units_general_cgs_conversion_factor(us, baseUnitsExp);
339
340
341
342
343
344
345
}

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

350
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
351

352
  return units_general_h_factor(us, baseUnitsExp);
353
354
355
356
357
358
359
}

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

364
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
365

366
  return units_general_a_factor(us, baseUnitsExp);
367
368
}

369
/**
370
 * @brief Returns a string containing the exponents of the base units making up
371
 * the conversion factors
372
 */
373
374
void units_cgs_conversion_string(char* buffer, const struct UnitSystem* us,
                                 enum UnitConversionFactor unit) {
375
  float baseUnitsExp[5] = {0.f};
376

377
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
378

379
  units_general_cgs_conversion_string(buffer, us, baseUnitsExp);
380
381
}

382
/**
383
384
 * @brief Returns the conversion factor for a given unit (expressed in terms of
 * the 5 fundamental units) in the chosen unit system
385
 * @param us The unit system used
386
 * @param baseUnitsExponants The exponent of each base units required to form
387
 * the desired quantity. See conversionFactor() for a working example
388
 */
389
390
double units_general_cgs_conversion_factor(const struct UnitSystem* us,
                                           const float baseUnitsExponants[5]) {
391
392
  double factor = 1.;

393
  for (int i = 0; i < 5; ++i)
394
    if (baseUnitsExponants[i] != 0)
395
396
      factor *= pow(units_get_base_unit(us, (enum BaseUnits)i),
                    baseUnitsExponants[i]);
397
  return factor;
398
399
400
}

/**
401
402
 * @brief Returns the h factor exponentiation for a given unit (expressed in
 * terms of the 5 fundamental units)
403
 * @param us The unit system used
404
 * @param baseUnitsExponants The exponent of each base units required to form
405
 * the desired quantity. See conversionFactor() for a working example
406
 */
407
float units_general_h_factor(const struct UnitSystem* us,
408
                             const float baseUnitsExponants[5]) {
409
  float factor_exp = 0.f;
410

411
412
413
  factor_exp += -baseUnitsExponants[UNIT_MASS];
  factor_exp += -baseUnitsExponants[UNIT_LENGTH];
  factor_exp += -baseUnitsExponants[UNIT_TIME];
414

415
  return factor_exp;
416
417
418
}

/**
419
420
 * @brief Returns the scaling factor exponentiation for a given unit (expressed
 * in terms of the 5 fundamental units)
421
 * @param us The unit system used
422
 * @param baseUnitsExponants The exponent of each base units required to form
423
 * the desired quantity. See conversionFactor() for a working example
424
 */
425
float units_general_a_factor(const struct UnitSystem* us,
426
                             const float baseUnitsExponants[5]) {
427
  float factor_exp = 0.f;
428

429
  factor_exp += baseUnitsExponants[UNIT_LENGTH];
430
431

  return factor_exp;
432
}
433
434

/**
435
 * @brief Returns a string containing the exponents of the base units making up
436
437
438
 * 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)
439
440
 * @param us The UnitsSystem in use.
 * @param baseUnitsExponants The exponent of each base units required to form
441
 * the desired quantity. See conversionFactor() for a working example
442
 */
443
444
445
void units_general_cgs_conversion_string(char* buffer,
                                         const struct UnitSystem* us,
                                         const float baseUnitsExponants[5]) {
446
  char temp[14];
447
448
  const double a_exp = units_general_a_factor(us, baseUnitsExponants);
  const double h_exp = units_general_h_factor(us, baseUnitsExponants);
449

450
451
  /* Check whether we are unitless or not */
  char isAllNonZero = 1;
452
  for (int i = 0; i < 5; ++i)
453
    if (baseUnitsExponants[i] != 0.) isAllNonZero = 0;
454

455
456
457
458
  if (isAllNonZero) {
    sprintf(buffer, "[ - ] ");
    return;
  }
459

460
  /* Add a-factor */
461
  if (a_exp == 0)
462
    sprintf(buffer, " ");
463
  else if (a_exp == 1)
464
    sprintf(buffer, "a ");
465
466
  else if (remainder(a_exp, 1.) == 0)
    sprintf(buffer, "a^%d ", (int)a_exp);
467
468
469
470
  else
    sprintf(buffer, "a^%7.4f ", a_exp);

  /* Add h-factor */
471
  if (h_exp == 0)
472
    sprintf(temp, " ");
473
  else if (h_exp == 1)
474
    sprintf(temp, "h ");
475
476
  else if (remainder(h_exp, 1.) == 0)
    sprintf(temp, "h^%d ", (int)h_exp);
477
478
479
480
481
  else
    sprintf(temp, "h^%7.4f ", h_exp);
  strncat(buffer, temp, 12);

  /* Add conversion units */
482
  for (int i = 0; i < 5; ++i)
483
484
485
486
    if (baseUnitsExponants[i] != 0) {
      if (baseUnitsExponants[i] == 0.)
        sprintf(temp, " ");
      else if (baseUnitsExponants[i] == 1.)
487
488
        sprintf(temp, "%s ",
                units_get_base_unit_internal_symbol((enum BaseUnits)i));
489
      else if (remainder(baseUnitsExponants[i], 1.) == 0)
490
491
        sprintf(temp, "%s^%d ",
                units_get_base_unit_internal_symbol((enum BaseUnits)i),
492
493
                (int)baseUnitsExponants[i]);
      else
494
495
        sprintf(temp, "%s^%7.4f ",
                units_get_base_unit_internal_symbol((enum BaseUnits)i),
496
                baseUnitsExponants[i]);
497
498
      strncat(buffer, temp, 12);
    }
499
500

  /* Add CGS units */
501
  strncat(buffer, " [ ", 3);
502

503
  for (int i = 0; i < 5; ++i) {
504
505
506
507
    if (baseUnitsExponants[i] != 0) {
      if (baseUnitsExponants[i] == 0.)
        continue;
      else if (baseUnitsExponants[i] == 1.)
508
        sprintf(temp, "%s ", units_get_base_unit_cgs_symbol((enum BaseUnits)i));
509
      else if (remainder(baseUnitsExponants[i], 1.) == 0)
510
511
        sprintf(temp, "%s^%d ",
                units_get_base_unit_cgs_symbol((enum BaseUnits)i),
512
513
                (int)baseUnitsExponants[i]);
      else
514
515
        sprintf(temp, "%s^%7.4f ",
                units_get_base_unit_cgs_symbol((enum BaseUnits)i),
516
517
                baseUnitsExponants[i]);
      strncat(buffer, temp, 12);
518
    }
519
520
  }

521
522
  strncat(buffer, "]", 2);
}
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566

/**
 * @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
567
 * @param unit The unit we are converting
568
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
569
 * @return The conversion factor
570
571
572
573
574
575
576
577
578
579
580
 */
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);
}