units.c 17.8 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
      break;

    case UNIT_CONV_VOLUME:
322
      baseUnitsExp[UNIT_LENGTH] = 3.f;
323
      break;
324
325

    case UNIT_CONV_INV_VOLUME:
326
      baseUnitsExp[UNIT_LENGTH] = -3.f;
327
328
329
330
331
      break;

    default:
      error("Invalid choice of pre-defined units");
      break;
332
  }
333
334
335
}

/**
336
337
 * @brief Returns the conversion factor for a given unit in the chosen unit
 * system
338
339
340
 * @param us The system of units in use
 * @param unit The unit to convert
 */
341
342
double units_cgs_conversion_factor(const struct UnitSystem* us,
                                   enum UnitConversionFactor unit) {
343
  float baseUnitsExp[5] = {0.f};
344

345
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
346

347
  return units_general_cgs_conversion_factor(us, baseUnitsExp);
348
349
350
351
352
353
354
}

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

359
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
360

361
  return units_general_h_factor(us, baseUnitsExp);
362
363
364
365
366
367
368
}

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

373
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
374

375
  return units_general_a_factor(us, baseUnitsExp);
376
377
}

378
/**
379
 * @brief Returns a string containing the exponents of the base units making up
380
 * the conversion factors
381
 */
382
383
void units_cgs_conversion_string(char* buffer, const struct UnitSystem* us,
                                 enum UnitConversionFactor unit) {
384
  float baseUnitsExp[5] = {0.f};
385

386
  units_get_base_unit_exponants_array(baseUnitsExp, unit);
387

388
  units_general_cgs_conversion_string(buffer, us, baseUnitsExp);
389
390
}

391
/**
392
393
 * @brief Returns the conversion factor for a given unit (expressed in terms of
 * the 5 fundamental units) in the chosen unit system
394
 * @param us The unit system used
395
 * @param baseUnitsExponants The exponent of each base units required to form
396
 * the desired quantity. See conversionFactor() for a working example
397
 */
398
399
double units_general_cgs_conversion_factor(const struct UnitSystem* us,
                                           const float baseUnitsExponants[5]) {
400
401
  double factor = 1.;

402
  for (int i = 0; i < 5; ++i)
403
    if (baseUnitsExponants[i] != 0)
404
405
      factor *= pow(units_get_base_unit(us, (enum BaseUnits)i),
                    baseUnitsExponants[i]);
406
  return factor;
407
408
409
}

/**
410
411
 * @brief Returns the h factor exponentiation for a given unit (expressed in
 * terms of the 5 fundamental units)
412
 * @param us The unit system used
413
 * @param baseUnitsExponants The exponent of each base units required to form
414
 * the desired quantity. See conversionFactor() for a working example
415
 */
416
float units_general_h_factor(const struct UnitSystem* us,
417
                             const float baseUnitsExponants[5]) {
418
  float factor_exp = 0.f;
419

420
421
422
  factor_exp += -baseUnitsExponants[UNIT_MASS];
  factor_exp += -baseUnitsExponants[UNIT_LENGTH];
  factor_exp += -baseUnitsExponants[UNIT_TIME];
423

424
  return factor_exp;
425
426
427
}

/**
428
429
 * @brief Returns the scaling factor exponentiation for a given unit (expressed
 * in terms of the 5 fundamental units)
430
 * @param us The unit system used
431
 * @param baseUnitsExponants The exponent of each base units required to form
432
 * the desired quantity. See conversionFactor() for a working example
433
 */
434
float units_general_a_factor(const struct UnitSystem* us,
435
                             const float baseUnitsExponants[5]) {
436
  float factor_exp = 0.f;
437

438
  factor_exp += baseUnitsExponants[UNIT_LENGTH];
439
440

  return factor_exp;
441
}
442
443

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

459
460
  /* Check whether we are unitless or not */
  char isAllNonZero = 1;
461
  for (int i = 0; i < 5; ++i)
462
    if (baseUnitsExponants[i] != 0.) isAllNonZero = 0;
463

464
465
466
467
  if (isAllNonZero) {
    sprintf(buffer, "[ - ] ");
    return;
  }
468

469
  /* Add a-factor */
470
  if (a_exp == 0)
471
    sprintf(buffer, " ");
472
  else if (a_exp == 1)
473
    sprintf(buffer, "a ");
474
475
  else if (remainder(a_exp, 1.) == 0)
    sprintf(buffer, "a^%d ", (int)a_exp);
476
477
478
479
  else
    sprintf(buffer, "a^%7.4f ", a_exp);

  /* Add h-factor */
480
  if (h_exp == 0)
481
    sprintf(temp, " ");
482
  else if (h_exp == 1)
483
    sprintf(temp, "h ");
484
485
  else if (remainder(h_exp, 1.) == 0)
    sprintf(temp, "h^%d ", (int)h_exp);
486
487
488
489
490
  else
    sprintf(temp, "h^%7.4f ", h_exp);
  strncat(buffer, temp, 12);

  /* Add conversion units */
491
  for (int i = 0; i < 5; ++i)
492
493
494
495
    if (baseUnitsExponants[i] != 0) {
      if (baseUnitsExponants[i] == 0.)
        sprintf(temp, " ");
      else if (baseUnitsExponants[i] == 1.)
496
497
        sprintf(temp, "%s ",
                units_get_base_unit_internal_symbol((enum BaseUnits)i));
498
      else if (remainder(baseUnitsExponants[i], 1.) == 0)
499
500
        sprintf(temp, "%s^%d ",
                units_get_base_unit_internal_symbol((enum BaseUnits)i),
501
502
                (int)baseUnitsExponants[i]);
      else
503
504
        sprintf(temp, "%s^%7.4f ",
                units_get_base_unit_internal_symbol((enum BaseUnits)i),
505
                baseUnitsExponants[i]);
506
507
      strncat(buffer, temp, 12);
    }
508
509

  /* Add CGS units */
510
  strncat(buffer, " [ ", 3);
511

512
  for (int i = 0; i < 5; ++i) {
513
514
515
516
    if (baseUnitsExponants[i] != 0) {
      if (baseUnitsExponants[i] == 0.)
        continue;
      else if (baseUnitsExponants[i] == 1.)
517
        sprintf(temp, "%s ", units_get_base_unit_cgs_symbol((enum BaseUnits)i));
518
      else if (remainder(baseUnitsExponants[i], 1.) == 0)
519
520
        sprintf(temp, "%s^%d ",
                units_get_base_unit_cgs_symbol((enum BaseUnits)i),
521
522
                (int)baseUnitsExponants[i]);
      else
523
524
        sprintf(temp, "%s^%7.4f ",
                units_get_base_unit_cgs_symbol((enum BaseUnits)i),
525
526
                baseUnitsExponants[i]);
      strncat(buffer, temp, 12);
527
    }
528
529
  }

530
531
  strncat(buffer, "]", 2);
}
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
567
568
569
570
571
572
573
574
575

/**
 * @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
576
 * @param unit The unit we are converting
577
 *
Matthieu Schaller's avatar
Matthieu Schaller committed
578
 * @return The conversion factor
579
580
581
582
583
584
585
586
587
588
589
 */
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);
}