dsurf.algoritm source code

1 module dsurf.algoritm;
2 
3 
4 import std.math;
5 import std.range;
6 import std.parallelism;
7 
8 import dsurf.cartesian;
9 
10 /***********************************
11 * Set of algoritms to work with cartesian surfaces
12 * Class CartesianSurface represents regular rectangular-cell surface
13 * Authors: Dmitriy Linev
14 * License: MIT
15 * Date: 2019-2020
16 */
17 
18 /// Returns surface gradient along X axis in the node with index i and j
19 double gX(CartesianSurface surface, int i, int j) {
20     if (i == 0)
21         return (surface.z[i + 1][j] - surface.z[i][j]) / surface.dx;
22     if (i == surface.nx - 1)
23         return (surface.z[i][j] - surface.z[i - 1][j]) / surface.dx;
24     return (surface.z[i + 1][j] - surface.z[i - 1][j]) / (2 * surface.dx);
25 }
26 
27 /// Returns surface gradient along Y axis in the node with index i and j
28 double gY(CartesianSurface surface, int i, int j) {
29     if (j == 0)
30         return (surface.z[i][j + 1] - surface.z[i][j]) / surface.dy;
31     if (j == surface.ny-1)
32         return (surface.z[i][j] - surface.z[i][j - 1]) / surface.dy;
33     return (surface.z[i][j + 1] - surface.z[i][j - 1]) / (2 * surface.dy);
34 }
35 
36 /// Returns CartesianSurface containing gradient map along X axis based on a given surface
37 CartesianSurface buildGradientXMap(CartesianSurface surface) {
38     CartesianSurface gradient = new CartesianSurface(surface);
39     foreach (i; taskPool.parallel(iota(surface.nx))) {
40         foreach (j; 0..surface.ny) {
41             gradient.z[i][j] = surface.gX(i, j);
42         }
43     }
44     return gradient;
45 }
46 
47 /// Returns CartesianSurface containing gradient map along Y axis based on a given surface
48 CartesianSurface buildGradientYMap(CartesianSurface surface) {
49     CartesianSurface gradient = new CartesianSurface(surface);
50     foreach (i; taskPool.parallel(iota(surface.nx))) {
51         foreach (j; 0..surface.ny) {
52             gradient.z[i][j] = surface.gY(i, j);
53         }
54     }
55     return gradient;
56 }
57 
58 /// Returns dip angle in radians of the surface in the node with index i and j
59 double dipAngle(CartesianSurface surface, int i, int j) {
60     immutable double gx = surface.gX(i, j);
61     immutable double gy = surface.gY(i, j);
62     if (isNaN(gx) || isNaN(gy))
63         return double.nan;
64     return atan(sqrt(gx * gx + gy * gy));
65 }
66 
67 /// Returns CartesianSurface containing dip angle map in radians based on a given surface
68 CartesianSurface buildDipAngle(CartesianSurface surface) {
69     CartesianSurface result = new CartesianSurface(surface);
70     foreach (i; taskPool.parallel(iota(surface.nx))) {
71         foreach (j; 0..surface.ny) {
72             result.z[i][j] = surface.dipAngle(i, j);
73         }
74     }
75     return result;
76 }