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 }