dsurf.algoritm source code

1 module dsurf.algoritm;
2 
3 import std.math;
4 import std.range;
5 import std.parallelism;
6 
7 import dsurf.cartesian;
8 
9 /***********************************
10 * Set of algoritms to work with surfaces
11 * Authors: Dmitriy Linev
12 * License: MIT
13 * Date: 2019-2020
14 */
15 
16 /// Returns surface gradient along X axis in the node with index i and j
17 double gX(CartesianSurface surface, int i, int j) {
18     if (i == 0)
19         return (surface.z[i + 1][j] - surface.z[i][j]) / surface.dx;
20     if (i == surface.nx - 1)
21         return (surface.z[i][j] - surface.z[i - 1][j]) / surface.dx;
22     return (surface.z[i + 1][j] - surface.z[i - 1][j]) / (2 * surface.dx);
23 }
24 
25 /// Returns surface gradient along Y axis in the node with index i and j
26 double gY(CartesianSurface surface, int i, int j) {
27     if (j == 0)
28         return (surface.z[i][j + 1] - surface.z[i][j]) / surface.dy;
29     if (j == surface.ny-1)
30         return (surface.z[i][j] - surface.z[i][j - 1]) / surface.dy;
31     return (surface.z[i][j + 1] - surface.z[i][j - 1]) / (2 * surface.dy);
32 }
33 
34 /// Returns CartesianSurface containing gradient map along X axis based on a given surface
35 CartesianSurface buildGradientXMap(CartesianSurface surface) {
36     CartesianSurface gradient = new CartesianSurface(surface);
37     foreach (i; taskPool.parallel(iota(surface.nx))) {
38         foreach (j; 0..surface.ny) {
39             gradient.z[i][j] = surface.gX(i, j);
40         }
41     }
42     return gradient;
43 }
44 
45 /// Returns CartesianSurface containing gradient map along Y axis based on a given surface
46 CartesianSurface buildGradientYMap(CartesianSurface surface) {
47     CartesianSurface gradient = new CartesianSurface(surface);
48     foreach (i; taskPool.parallel(iota(surface.nx))) {
49         foreach (j; 0..surface.ny) {
50             gradient.z[i][j] = surface.gY(i, j);
51         }
52     }
53     return gradient;
54 }
55 
56 /// Returns dip angle in radians of the surface in the node with index i and j
57 double dipAngle(CartesianSurface surface, int i, int j) {
58     immutable double gx = surface.gX(i, j);
59     immutable double gy = surface.gY(i, j);
60     if (isNaN(gx) || isNaN(gy))
61         return double.nan;
62     return atan(sqrt(gx * gx + gy * gy));
63 }
64 
65 /// Returns CartesianSurface containing dip angle map in radians based on a given surface
66 CartesianSurface buildDipAngle(CartesianSurface surface) {
67     CartesianSurface result = new CartesianSurface(surface);
68     foreach (i; taskPool.parallel(iota(surface.nx))) {
69         foreach (j; 0..surface.ny) {
70             result.z[i][j] = surface.dipAngle(i, j);
71         }
72     }
73     return result;
74 }
75 
76 /** 
77  * Scales the given surface around its origin point
78  * Params:
79  *   surface = surface to scale
80  *   xf = scale factor along X direction
81  *   xf = scale factor value along Y direction
82  * Returns: translated surface for more convenient call chaining
83  */
84 CartesianSurface scale(CartesianSurface surface, double xf, double yf) {  //scales around origin point
85     //TODO filter negative factors
86     surface.m_dx *= xf;
87     surface.m_dy *= yf;
88     return surface;
89 }
90 
91 /** 
92  * Normalizes the given surface (minimum value will be 0, maximum value will be 1)
93  * Doesn`t alter surface limits and origin point
94  * Case of equal min and max is not controlled (division by zero will occur)
95  * Params:
96  *   surface = surface to normalize
97  * Returns: Normalized surface for more convenient call chaining
98  */
99 CartesianSurface normalize(CartesianSurface surface) {
100     immutable double zmax = surface.max();
101     immutable double zmin = surface.min();
102     surface -= zmin;
103     surface /= (zmax - zmin);
104     return surface;
105 }