1
+ # Specify Cython compiler flags.
1
2
# cython: boundscheck=False
2
3
# cython: cdivision=True
3
4
# cython: language_level=3
4
5
# distutils: language=c++
5
6
6
- # Python imports
7
+ # Python imports.
7
8
import collections
8
9
import itertools
9
10
import numpy as np
10
11
11
- # Cython imports
12
+ # Cython imports.
12
13
cimport numpy as np
13
14
cimport libc.math
14
15
from libcpp.vector cimport vector as cpp_vector
@@ -20,6 +21,7 @@ ctypedef np.uint16_t _UINT_t
20
21
DTYPE_REAL = np.float64
21
22
DTYPE_UINT = np.uint16
22
23
24
+ # Define a value to signal errors.
23
25
DEF _ERROR_REAL = - 999999999999.
24
26
ERROR_REAL = DTYPE_REAL(_ERROR_REAL)
25
27
@@ -34,54 +36,84 @@ class OutOfBoundsError(Exception):
34
36
35
37
36
38
class EikonalSolver (object ):
39
+ '''
40
+ A class to solver the Eikonal equation in 3D Cartesian or
41
+ spherical coordinates.
42
+
43
+ Properties
44
+ **********
45
+ .. autoattribute:: close
46
+ .. autoattribute:: iax_null
47
+ .. autoattribute:: is_alive
48
+ .. autoattribute:: is_far
49
+
50
+ Methods
51
+ *******
52
+ .. automethod:: trace_ray(self, end, step_size=None)
53
+ '''
54
+
37
55
def __init__ (self , coord_sys = ' cartesian' ):
38
- '''
39
- Solves the Eikonal equation in 3D cartesian coordinates.
40
- '''
41
56
self ._ndim = 3
42
57
self ._class = str (self .__class__).strip(' >\' ' ).split(' .' )[- 1 ]
43
58
self ._coord_sys = coord_sys
44
59
self ._vgrid = GridND(ndim = self ._ndim, coord_sys = self .coord_sys)
45
60
46
61
@property
47
62
def close (self ):
63
+ '''
64
+ A list of node indices in the *Trial* region.
65
+ '''
48
66
if not hasattr (self , ' _close' ):
49
67
self ._close = Heap(self .uu)
50
68
return (self ._close)
51
69
52
70
@property
53
71
def iax_null (self ):
72
+ '''
73
+ '''
54
74
return (self .pgrid.iax_null)
55
75
56
76
@property
57
77
def is_alive (self ):
78
+ '''
79
+ '''
58
80
if not hasattr (self , ' _is_alive' ):
59
81
self ._is_alive = np.full(self .pgrid.npts, fill_value = False , dtype = np.bool)
60
82
return (self ._is_alive)
61
83
62
84
@property
63
85
def is_far (self ):
86
+ '''
87
+ '''
64
88
if not hasattr (self , ' _is_far' ):
65
89
self ._is_far = np.full(self .pgrid.npts, fill_value = True , dtype = np.bool)
66
90
return (self ._is_far)
67
91
68
92
@property
69
93
def coord_sys (self ):
94
+ '''
95
+ '''
70
96
return (self ._coord_sys)
71
97
72
98
@coord_sys.setter
73
99
def coord_sys (self , value ):
100
+ '''
101
+ '''
74
102
value = value.lower()
75
103
if value not in (' cartesian' , ' spherical' ):
76
104
raise (ValueError (f' Invalid coord_sys specification: {value}' ))
77
105
self ._coord_sys = value
78
106
79
107
@property
80
108
def ndim (self ):
109
+ '''
110
+ '''
81
111
return (self ._ndim)
82
112
83
113
@property
84
114
def norm (self ):
115
+ '''
116
+ '''
85
117
if not hasattr (self , ' _norm' ):
86
118
self ._norm = np.tile(
87
119
self .pgrid.node_intervals,
@@ -95,6 +127,8 @@ class EikonalSolver(object):
95
127
96
128
@property
97
129
def pgrid (self ):
130
+ '''
131
+ '''
98
132
if not hasattr (self , ' _pgrid' ):
99
133
self ._pgrid = GridND(ndim = self ._ndim, coord_sys = self .vgrid.coord_sys)
100
134
for attr in (' min_coords' , ' node_intervals' , ' npts' ):
@@ -103,10 +137,14 @@ class EikonalSolver(object):
103
137
104
138
@property
105
139
def src_loc (self ):
140
+ '''
141
+ '''
106
142
return (self ._src_loc)
107
143
108
144
@src_loc.setter
109
145
def src_loc (self , value ):
146
+ '''
147
+ '''
110
148
if not isinstance (value, collections.Iterable):
111
149
raise (TypeError (f' {self._class}.src_loc value must be <Iterable> type' ))
112
150
if len (value) != self ._ndim:
@@ -118,6 +156,8 @@ class EikonalSolver(object):
118
156
119
157
@property
120
158
def src_rtp (self ):
159
+ '''
160
+ '''
121
161
if self .coord_sys == ' spherical' :
122
162
return (self .src_loc)
123
163
else :
@@ -128,6 +168,8 @@ class EikonalSolver(object):
128
168
129
169
@property
130
170
def src_xyz (self ):
171
+ '''
172
+ '''
131
173
if self .coord_sys == ' cartesian' :
132
174
return (self .src_loc)
133
175
else :
@@ -139,16 +181,22 @@ class EikonalSolver(object):
139
181
140
182
@property
141
183
def uu (self ):
184
+ '''
185
+ '''
142
186
if not hasattr (self , ' _uu' ):
143
187
self ._uu = np.full(self .pgrid.npts, fill_value = np.inf, dtype = DTYPE_REAL)
144
188
return (self ._uu)
145
189
146
190
@property
147
191
def vgrid (self ):
192
+ '''
193
+ '''
148
194
return (self ._vgrid)
149
195
150
196
@property
151
197
def vv (self ):
198
+ '''
199
+ '''
152
200
return (self ._vv)
153
201
154
202
@vv.setter
@@ -159,6 +207,8 @@ class EikonalSolver(object):
159
207
160
208
@property
161
209
def vvp (self ):
210
+ '''
211
+ '''
162
212
cdef Py_ssize_t i1, i2, i3
163
213
cdef np.ndarray[_REAL_t, ndim= 3 ] vvp
164
214
@@ -186,10 +236,32 @@ class EikonalSolver(object):
186
236
187
237
188
238
def solve (self ):
239
+ '''
240
+ '''
189
241
self ._update()
190
242
191
243
192
- def trace_ray (self , start , step_size = None ):
244
+ def trace_ray (self , end , step_size = None ):
245
+ '''
246
+ Trace the ray ending at *end*.
247
+
248
+ This method traces the ray that ends at *end* in reverse
249
+ direction by taking small steps along the path of steepest
250
+ travel-time descent. The resulting ray is reversed before being
251
+ returned, so it is in the normal forward-time orientation.
252
+
253
+ :param end: Coordinates of the ray's end point.
254
+ :type end: tuple, list, np.ndarray
255
+
256
+ :param step_size: The distance between points on the ray.
257
+ The smaller this value is the more accurate the resulting
258
+ ray will be. By default, this parameter is assigned the
259
+ smallest node interval of the propagation grid.
260
+ :type step_size: float, optional
261
+
262
+ :return: The ray path ending at *end*.
263
+ :rtype: np.ndarray(Nx3)
264
+ '''
193
265
cdef cpp_vector[_REAL_t * ] ray
194
266
cdef _REAL_t g0, g1, g2, norm
195
267
cdef _REAL_t * point_new
@@ -198,7 +270,7 @@ class EikonalSolver(object):
198
270
cdef np.ndarray[_REAL_t, ndim= 2 ] ray_np
199
271
200
272
point_new = < _REAL_t * > malloc(3 * sizeof(_REAL_t))
201
- point_new[0 ], point_new[1 ], point_new[2 ] = start
273
+ point_new[0 ], point_new[1 ], point_new[2 ] = end
202
274
ray.push_back(point_new)
203
275
if step_size is None :
204
276
step_size = self ._get_step_size()
@@ -243,14 +315,26 @@ class EikonalSolver(object):
243
315
ray_np[i, 1 ] = ray[i][1 ]
244
316
ray_np[i, 2 ] = ray[i][2 ]
245
317
free(ray[i])
246
- return (ray_np)
318
+ return (np.flipud( ray_np) )
247
319
248
320
249
321
def transfer_travel_times_from (self , old , origin , rotate = False , set_alive = False ):
250
322
'''
251
323
Transfer the velocity model from old EikonalSolver to self
252
- :param pykonal.EikonalSolver old: The old EikonalSolver to transfer from.
253
- :param tuple old_origin: The coordinates of the origin of old w.r.t. to the self frame of reference.
324
+
325
+ :param old: The old EikonalSolver to transfer from.
326
+ :type old: pykonal.EikonalSolver
327
+
328
+ :param origin: The coordinates of the origin of old w.r.t. to
329
+ the self frame of reference.
330
+ :type origin: tuple, list, np.ndarray
331
+
332
+ :param rotate: Rotate the coordinates?
333
+ :type rotate: bool
334
+
335
+ :return: None
336
+ :rtype: NoneType
337
+
254
338
'''
255
339
256
340
pgrid_new = self .pgrid.map_to(old.coord_sys, origin, rotate = rotate)
0 commit comments