tutil.c - granular - granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/granular
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
tutil.c (3304B)
---
1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <math.h>
4 #include <err.h>
5 #include "granular.h"
6
7 void
8 warn_parameter_value(const char message[],
9 const double value,
10 int *status)
11 {
12 fprintf(stderr,
13 "%s: %s (%.17g)\n",
14 __func__, message, value);
15 *status = 1;
16 }
17
18 void
19 check_float(const char name[], const double value, int *status)
20 {
21 int ret;
22 char message[100];
23
24 if (isnan(value)) {
25 ret = snprintf(message, sizeof(message), "%s is NaN", name);
26 if (ret < 0 || (size_t)ret >= sizeof(message))
27 err(1, "%s: message parsing", __func__);
28 warn_parameter_value(message, value, status);
29 *status = 1;
30 } else if (isinf(value)) {
31 ret = snprintf(message, sizeof(message), "%s is infinite", name);
32 if (ret < 0 || (size_t)ret >= sizeof(message))
33 err(1, "%s: message parsing", __func__);
34 warn_parameter_value(message, value, status);
35 *status = 1;
36 }
37 }
38
39 void
40 check_float_non_negative(const char name[], const double value, int *status)
41 {
42 int ret;
43 char message[100];
44
45 check_float(name, value, status);
46 if (value < 0.0) {
47 ret = snprintf(message, sizeof(message), "%s is negative", name);
48 if (ret < 0 || (size_t)ret >= sizeof(message))
49 err(1, "%s: message parsing", __func__);
50 warn_parameter_value(message, value, status);
51 *status = 1;
52 }
53 }
54
55 void
56 check_float_positive(const char name[], const double value, int *status)
57 {
58 int ret;
59 char message[100];
60
61 check_float(name, value, status);
62 if (value <= 0.0) {
63 ret = snprintf(message, sizeof(message), "%s is not positive", name);
64 if (ret < 0 || (size_t)ret >= sizeof(message))
65 err(1, "%s: message parsing", __func__);
66 warn_parameter_value(message, value, status);
67 *status = 1;
68 }
69 }
70
71 void
72 check_int_bool(const char name[], const int value, int *status)
73 {
74 int ret;
75 char message[100];
76
77 if (value < 0 || value > 1) {
78 ret = snprintf(message, sizeof(message), "%s is not 0 or 1", name);
79 if (ret < 0 || (size_t)ret >= sizeof(message))
80 err(1, "%s: message parsing", __func__);
81 warn_parameter_value(message, (double)value, status);
82 *status = 1;
83 }
84 }
85
86 void
87 check_int_non_negative(const char name[], const int value, int *status)
88 {
89 int ret;
90 char message[100];
91
92 if (value < 0) {
93 ret = snprintf(message, sizeof(message), "%s is negative", name);
94 if (ret < 0 || (size_t)ret >= sizeof(message))
95 err(1, "%s: message parsing", __func__);
96 warn_parameter_value(message, (double)value, status);
97 *status = 1;
98 }
99 }
100
101 double
102 residual(const double new_val, const double old_val)
103 {
104 return (new_val - old_val) / (old_val + 1e-16);
105 }
106
107 double
108 randomval(void)
109 {
110 return (double)rand() / RAND_MAX;
111 }
112
113 double
114 random_value_uniform(double min, double max)
115 {
116 return randomval() * (max - min) + min;
117 }
118
119 double
120 random_value_powerlaw(double min, double max)
121 {
122 double power = -1.8;
123 return pow((pow(max, power + 1.0) - pow(min, power + 1.0)) * randomval()
124 + pow(min, power + 1.0), 1.0 / (power + 1.0));
125 }
126
127 void *
128 xreallocarray(void *m, size_t n, size_t s)
129 {
130 void *nm;
131
132 if (n == 0 || s == 0) {
133 free(m);
134 return NULL;
135 }
136 if (s && n > (size_t) - 1 / s)
137 err(1, "realloc: overflow");
138 if (!(nm = realloc(m, n * s)))
139 err(1, "realloc");
140
141 return nm;
142 }
143
144 double
145 harmonic_mean(double a, double b)
146 {
147 if (a < 1e-16 || b < 1e-16)
148 return 0.0;
149 else
150 return 2.0 * a * b / (a + b);
151 }