arrays.c (6437B)
1 #include <stdlib.h> 2 #include <stdio.h> 3 #include <math.h> 4 #include <err.h> 5 6 #define DEG2RAD(x) (x*M_PI/180.0) 7 8 void 9 check_magnitude(const char *func_name, int limit, int value) 10 { 11 if (value < limit) 12 errx(1, "%s: input size %d is less than %d\n", 13 func_name, value, limit); 14 } 15 16 /* Translate a i,j,k index in grid with dimensions nx, ny, nz into a 17 * linear index */ 18 unsigned int 19 idx3(const unsigned int i, 20 const unsigned int j, 21 const unsigned int k, 22 const unsigned int nx, 23 const unsigned int ny) 24 { 25 return i + nx * j + nx * ny * k; 26 } 27 28 /* Translate a i,j,k index in grid with dimensions nx, ny, nz and a 29 * padding of single ghost nodes into a linear index */ 30 unsigned int 31 idx3g(const unsigned int i, 32 const unsigned int j, 33 const unsigned int k, 34 const unsigned int nx, 35 const unsigned int ny) 36 { 37 return i + 1 + (nx + 2) * (j + 1) + (nx + 2) * (ny + 2) * (k + 1); 38 } 39 40 /* Translate a i,j index in grid with dimensions nx, ny into a linear 41 * index */ 42 unsigned int 43 idx2(const unsigned int i, const unsigned int j, const unsigned int nx) 44 { 45 return i + nx * j; 46 } 47 48 /* Translate a i,j index in grid with dimensions nx, ny and a padding of 49 * single ghost nodes into a linear index */ 50 unsigned int 51 idx2g(const unsigned int i, const unsigned int j, const unsigned int nx) 52 { 53 return i + 1 + (nx + 2) * (j + 1); 54 } 55 56 /* Translate a i index in grid with a padding of single into a linear 57 * index */ 58 unsigned int 59 idx1g(const unsigned int i) 60 { 61 return i + 1; 62 } 63 64 /* Return an array of `n` linearly spaced values in the range [lower; 65 * upper] */ 66 double * 67 linspace(const double lower, const double upper, const int n) 68 { 69 int i; 70 double *x; 71 double dx; 72 73 check_magnitude(__func__, 1, n); 74 if (!(x = calloc(n, sizeof(double)))) 75 err(1, "%s: x calloc", __func__); 76 dx = (upper - lower) / (double) (n - 1); 77 for (i = 0; i < n; ++i) 78 x[i] = lower + dx * i; 79 80 return x; 81 } 82 83 /* Return an array of `n-1` values with the intervals between `x` values */ 84 double * 85 spacing(const double *x, const int n) 86 { 87 int i; 88 double *dx; 89 90 check_magnitude(__func__, 2, n); 91 if (!(dx = calloc((n - 1), sizeof(double)))) 92 err(1, "%s: dx calloc", __func__); 93 for (i = 0; i < n - 1; ++i) 94 dx[i] = x[i + 1] - x[i]; 95 96 return dx; 97 } 98 99 /* Return an array of `n` values with the value 0.0 */ 100 double * 101 zeros(const int n) 102 { 103 int i; 104 double *x; 105 106 check_magnitude(__func__, 1, n); 107 if (!(x = calloc(n, sizeof(double)))) 108 err(1, "%s: x calloc", __func__); 109 for (i = 0; i < n; ++i) 110 x[i] = 0.0; 111 112 return x; 113 } 114 115 /* Return an array of `n` values with the value 1.0 */ 116 double * 117 ones(const int n) 118 { 119 int i; 120 double *x; 121 122 check_magnitude(__func__, 1, n); 123 if (!(x = calloc(n, sizeof(double)))) 124 err(1, "%s: x calloc", __func__); 125 for (i = 0; i < n; ++i) 126 x[i] = 1.0; 127 128 return x; 129 } 130 131 /* Return an array of `n` values with a specified value */ 132 double * 133 initval(const double value, const int n) 134 { 135 int i; 136 double *x; 137 138 check_magnitude(__func__, 1, n); 139 if (!(x = calloc(n, sizeof(double)))) 140 err(1, "%s: x calloc", __func__); 141 for (i = 0; i < n; ++i) 142 x[i] = value; 143 144 return x; 145 } 146 147 /* Return an array of `n` uninitialized values */ 148 double * 149 empty(const int n) 150 { 151 double *out; 152 153 check_magnitude(__func__, 1, n); 154 155 if (!(out = calloc(n, sizeof(double)))) 156 err(1, "%s: calloc", __func__); 157 158 return out; 159 } 160 161 /* Return largest value in array `a` with size `n` */ 162 double 163 max(const double *a, const int n) 164 { 165 int i; 166 double maxval; 167 168 check_magnitude(__func__, 1, n); 169 maxval = -INFINITY; 170 for (i = 0; i < n; ++i) 171 if (a[i] > maxval) 172 maxval = a[i]; 173 174 return maxval; 175 } 176 177 /* Return smallest value in array `a` with size `n` */ 178 double 179 min(const double *a, const int n) 180 { 181 int i; 182 double minval; 183 184 check_magnitude(__func__, 1, n); 185 minval = +INFINITY; 186 for (i = 0; i < n; ++i) 187 if (a[i] < minval) 188 minval = a[i]; 189 190 return minval; 191 } 192 193 void 194 print_array(const double *a, const int n) 195 { 196 int i; 197 198 check_magnitude(__func__, 1, n); 199 for (i = 0; i < n; ++i) 200 printf("%.17g\n", a[i]); 201 } 202 203 void 204 print_arrays(const double *a, const double *b, const int n) 205 { 206 int i; 207 208 check_magnitude(__func__, 1, n); 209 for (i = 0; i < n; ++i) 210 printf("%.17g\t%.17g\n", a[i], b[i]); 211 } 212 213 void 214 print_arrays_2nd_normalized(const double *a, const double *b, const int n) 215 { 216 int i; 217 double max_b; 218 219 check_magnitude(__func__, 1, n); 220 max_b = max(b, n); 221 for (i = 0; i < n; ++i) 222 printf("%.17g\t%.17g\n", a[i], b[i] / max_b); 223 } 224 225 void 226 print_three_arrays(const double *a, 227 const double *b, 228 const double *c, 229 const int n) 230 { 231 int i; 232 233 check_magnitude(__func__, 1, n); 234 for (i = 0; i < n; ++i) 235 printf("%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]); 236 } 237 238 void 239 fprint_arrays(FILE *fp, const double *a, const double *b, const int n) 240 { 241 int i; 242 243 check_magnitude(__func__, 1, n); 244 for (i = 0; i < n; ++i) 245 fprintf(fp, "%.17g\t%.17g\n", a[i], b[i]); 246 } 247 248 void 249 fprint_three_arrays(FILE *fp, 250 const double *a, 251 const double *b, 252 const double *c, 253 const int n) 254 { 255 int i; 256 257 check_magnitude(__func__, 1, n); 258 for (i = 0; i < n; ++i) 259 fprintf(fp, "%.17g\t%.17g\t%.17g\n", a[i], b[i], c[i]); 260 } 261 262 void 263 copy_values(const double *in, double *out, const int n) 264 { 265 int i; 266 267 check_magnitude(__func__, 1, n); 268 for (i = 0; i < n; ++i) 269 out[i] = in[i]; 270 } 271 272 double * 273 copy(const double *in, const int n) 274 { 275 double *out; 276 277 check_magnitude(__func__, 1, n); 278 out = empty(n); 279 copy_values(in, out, n); 280 return out; 281 } 282 283 double * 284 normalize(const double *in, const int n) 285 { 286 int i; 287 double max_val; 288 double *out; 289 290 check_magnitude(__func__, 1, n); 291 out = empty(n); 292 copy_values(in, out, n); 293 max_val = max(out, n); 294 295 if (max_val == 0.0) 296 max_val = 1.0; 297 298 for (i = 0; i < n; ++i) 299 out[i] /= max_val; 300 301 return out; 302 } 303 304 double 305 euclidean_norm(const double *a, const int n) 306 { 307 int i; 308 double out = 0.0; 309 310 for (i = 0; i < n; i++) 311 out += pow(a[i], 2.0); 312 313 return sqrt(out); 314 } 315 316 double 317 euclidean_distance(const double *a, const double *b, const int n) 318 { 319 int i; 320 double out = 0.0; 321 322 for (i = 0; i < n; i++) 323 out += pow(b[i] - a[i], 2.0); 324 325 return sqrt(out); 326 } 327 328 double 329 dot(const double *a, const double *b, const int n) 330 { 331 int i; 332 double out = 0.0; 333 334 for (i = 0; i < n; i++) 335 out += a[i] * b[i]; 336 337 return out; 338 } 339 340 double * 341 cross(const double a[3], const double b[3]) 342 { 343 double *out = empty(3); 344 345 out[0] = a[1]*b[2] - a[2]*b[1]; 346 out[1] = a[2]*b[0] - a[0]*b[2]; 347 out[2] = a[0]*b[1] - a[1]*b[0]; 348 349 return out; 350 }