h5geo 0.4.0
C++17 and python API to work with geo-data (seismic, wells, maps, other in process) based on HDF5. Aimed at geoscientists and developers.
Loading...
Searching...
No Matches
h5deviation.h
1#ifndef H5DEVIATION_H
2#define H5DEVIATION_H
3
4#include "../private/h5enum.h"
5
6#define _USE_MATH_DEFINES // should be before <cmath>, include 'pi' val
7
8#include <math.h>
9#include <Eigen/Dense>
10#include <units/units.hpp>
11
12namespace h5geo
13{
14/* */
15template<typename D>
16Eigen::MatrixX<typename D::Scalar> MdAzIncl2ALL(
17 const Eigen::DenseBase<D> &M,
18 const double& x0, const double& y0, const double& kb,
19 const std::string& angularUnitsFrom,
20 const std::string& angularUnitsTo,
21 const bool& XNorth);
22
23template<typename D>
24Eigen::MatrixX<typename D::Scalar> TvdXY2ALL(
25 const Eigen::DenseBase<D> &M,
26 const double& x0, const double& y0, const double& kb,
27 const std::string& angularUnits,
28 const bool& XNorth);
29
30template<typename D>
31Eigen::MatrixX<typename D::Scalar> TvdDxDy2ALL(
32 const Eigen::DenseBase<D> &M,
33 const double& x0, const double& y0, const double& kb,
34 const std::string& angularUnits,
35 const bool& XNorth);
36
37template<typename D>
38Eigen::MatrixX<typename D::Scalar> TvdssXY2ALL(
39 const Eigen::DenseBase<D> &M,
40 const double& x0, const double& y0, const double& kb,
41 const std::string& angularUnits,
42 const bool& XNorth);
43
44template<typename D>
45Eigen::MatrixX<typename D::Scalar> TvdssDxDy2ALL(
46 const Eigen::DenseBase<D> &M,
47 const double& x0, const double& y0, const double& kb,
48 const std::string& angularUnits,
49 const bool& XNorth);
50
51template<typename D>
52Eigen::MatrixX<typename D::Scalar> traj2ALL(
53 const Eigen::DenseBase<D> &M,
54 const double& x0, const double& y0, const double& kb,
55 const std::string& angularUnitsFrom,
56 const std::string& angularUnitsTo,
57 const h5geo::TrajectoryFormat& trajFormat,
58 const bool& XNorth);
59
60/* */
61template<typename D>
62Eigen::MatrixX<typename D::Scalar> MdAzIncl2MdXYTvd(
63 const Eigen::DenseBase<D> &M,
64 const double& x0, const double& y0,
65 const std::string& angularUnits,
66 const bool& XNorth);
67
68template<typename D>
69Eigen::MatrixX<typename D::Scalar> TvdXY2MdXYTvd(
70 const Eigen::DenseBase<D> &M,
71 const double& x0, const double& y0,
72 const bool& XNorth);
73
74template<typename D>
75Eigen::MatrixX<typename D::Scalar> TvdDxDy2MdXYTvd(
76 const Eigen::DenseBase<D> &M,
77 const double& x0, const double& y0,
78 const bool& XNorth);
79
80template<typename D>
81Eigen::MatrixX<typename D::Scalar> TvdssXY2MdXYTvd(
82 const Eigen::DenseBase<D> &M,
83 const double& x0, const double& y0, const double& kb,
84 const bool& XNorth);
85
86template<typename D>
87Eigen::MatrixX<typename D::Scalar> TvdssDxDy2MdXYTvd(
88 const Eigen::DenseBase<D> &M,
89 const double& x0, const double& y0, const double& kb,
90 const bool& XNorth);
91
92template<typename D>
93Eigen::MatrixX<typename D::Scalar> TvdXY2MdAzIncl(
94 const Eigen::DenseBase<D> &M,
95 const double& x0, const double& y0, const std::string& angularUnits, const bool& XNorth);
96
97template<typename D>
98Eigen::MatrixX<typename D::Scalar> TvdDxDy2MdAzIncl(
99 const Eigen::DenseBase<D> &M, const std::string& angularUnits, const bool& XNorth);
100
101template<typename D>
102Eigen::MatrixX<typename D::Scalar> TvdssXY2MdAzIncl(
103 const Eigen::DenseBase<D> &M,
104 const double& x0, const double& y0, const double& kb, const std::string& angularUnits, const bool& XNorth);
105
106template<typename D>
107Eigen::MatrixX<typename D::Scalar> TvdssDxDy2MdAzIncl(
108 const Eigen::DenseBase<D> &M,
109 const double& kb, const std::string& angularUnits, const bool& XNorth);
110
111
118inline Eigen::VectorX<double> angleAndRadius(
119 const double& arcLen,
120 const double& chordLen,
121 const double& eps,
122 const size_t& max_iter);
123
124
125inline double _betaAngle(
126 double I1,
127 double I2,
128 double A1,
129 double A2,
130 const std::string& angularUnits)
131{
132 if (!angularUnits.empty()){
133 double coef = units::convert(
134 units::unit_from_string(angularUnits),
135 units::unit_from_string("radian"));
136 I1 *= coef;
137 I2 *= coef;
138 A1 *= coef;
139 A2 *= coef;
140 }
141
142 return acos(cos(I2-I1) - (sin(I1)*sin(I2)*(1 - cos(A2-A1))));
143}
144
145
146inline double _ratioFactor(double B, const std::string& angularUnits)
147{
148 if (!angularUnits.empty()){
149 double coef = units::convert(
150 units::unit_from_string(angularUnits),
151 units::unit_from_string("radian"));
152 B *= coef;
153 }
154
155 if (B == 0)
156 return 1 + pow(B,2)/12 + pow(B,4)/120 + pow(B,6)/20160; // or it is possible to simply return 1
157 else
158 return 2 / B * tan(B / 2);
159}
160
161
162inline double _deltaEast(
163 const double& dMD,
164 double I1,
165 double I2,
166 double A1,
167 double A2,
168 const double& RF,
169 const std::string& angularUnits)
170{
171 if (!angularUnits.empty()){
172 double coef = units::convert(
173 units::unit_from_string(angularUnits),
174 units::unit_from_string("radian"));
175 I1 *= coef;
176 I2 *= coef;
177 A1 *= coef;
178 A2 *= coef;
179 }
180
181 return dMD/2 * (sin(I1)*sin(A1) + sin(I2)*sin(A2))*RF;
182}
183
184
185inline double _deltaNorth(
186 const double& dMD,
187 double I1,
188 double I2,
189 double A1,
190 double A2,
191 const double& RF,
192 const std::string& angularUnits)
193{
194 if (!angularUnits.empty()){
195 double coef = units::convert(
196 units::unit_from_string(angularUnits),
197 units::unit_from_string("radian"));
198 I1 *= coef;
199 I2 *= coef;
200 A1 *= coef;
201 A2 *= coef;
202 }
203
204 return dMD/2 * (sin(I1)*cos(A1) + sin(I2)*cos(A2))*RF;
205}
206
207
208inline double _deltaZ(
209 const double& dMD,
210 double I1,
211 double I2,
212 const double& RF,
213 const std::string& angularUnits)
214{
215 if (!angularUnits.empty()){
216 double coef = units::convert(
217 units::unit_from_string(angularUnits),
218 units::unit_from_string("radian"));
219 I1 *= coef;
220 I2 *= coef;
221 }
222
223 return dMD/2 * (cos(I1) + cos(I2))*RF;
224}
225
226
227inline double _dMDFromDeltaEast(
228 const double& dE,
229 double I1,
230 double I2,
231 double A1,
232 double A2,
233 const double& RF,
234 const std::string& angularUnits)
235{
236 if (!angularUnits.empty()){
237 double coef = units::convert(
238 units::unit_from_string(angularUnits),
239 units::unit_from_string("radian"));
240 I1 *= coef;
241 I2 *= coef;
242 A1 *= coef;
243 A2 *= coef;
244 }
245
246 return 2*dE / ((sin(I1)*sin(A1) + sin(I2)*sin(A2))*RF);
247}
248
249
250inline double _dMDFromDeltaNorth(
251 const double& dN,
252 double I1,
253 double I2,
254 double A1,
255 double A2,
256 const double& RF,
257 const std::string& angularUnits)
258{
259 if (!angularUnits.empty()){
260 double coef = units::convert(
261 units::unit_from_string(angularUnits),
262 units::unit_from_string("radian"));
263 I1 *= coef;
264 I2 *= coef;
265 A1 *= coef;
266 A2 *= coef;
267 }
268
269 return 2*dN / ((sin(I1)*cos(A1) + sin(I2)*cos(A2))*RF);
270}
271
272
273inline double _dMDFromDeltaZ(
274 const double& dZ,
275 double I1,
276 double I2,
277 const double& RF,
278 const std::string& angularUnits)
279{
280 if (!angularUnits.empty()){
281 double coef = units::convert(
282 units::unit_from_string(angularUnits),
283 units::unit_from_string("radian"));
284 I1 *= coef;
285 I2 *= coef;
286 }
287
288 return 2*dZ / ((cos(I1) + cos(I2))*RF);
289}
290
291
292inline double _fa(
293 const double& arc,
294 const double& chord,
295 double a,
296 const std::string& angularUnits)
297{
298 if (!angularUnits.empty()){
299 double coef = units::convert(
300 units::unit_from_string(angularUnits),
301 units::unit_from_string("radian"));
302 a *= coef;
303 }
304
305 return chord*a-2*arc*sin(a/2);
306}
307
308
309inline double _dfa(
310 const double& arc,
311 const double& chord,
312 double a,
313 const std::string& angularUnits)
314{
315 if (!angularUnits.empty()){
316 double coef = units::convert(
317 units::unit_from_string(angularUnits),
318 units::unit_from_string("radian"));
319 a *= coef;
320 }
321
322 return chord-arc*cos(a/2);
323}
324
325
326inline double _ddfa(
327 const double& arc,
328 double a,
329 const std::string& angularUnits)
330{
331 if (!angularUnits.empty()){
332 double coef = units::convert(
333 units::unit_from_string(angularUnits),
334 units::unit_from_string("radian"));
335 a *= coef;
336 }
337
338 return arc*sin(a/2)/2;
339}
340
341};
342
343
344
345
346
347template<typename D>
348Eigen::MatrixX<typename D::Scalar> h5geo::MdAzIncl2ALL(
349 const Eigen::DenseBase<D> &M,
350 const double& x0, const double& y0, const double& kb,
351 const std::string& angularUnitsFrom,
352 const std::string& angularUnitsTo,
353 const bool& XNorth)
354{
355 if (M.cols() != 3)
356 return Eigen::MatrixX<typename D::Scalar>();
357
358 Eigen::MatrixX<typename D::Scalar> M_MdXYTvd =
359 MdAzIncl2MdXYTvd(M, x0, y0, angularUnitsFrom, XNorth);
360
361 Eigen::MatrixX<typename D::Scalar> M_OUT(M.rows(), 9);
362 M_OUT.col(0) = M_MdXYTvd.col(0); // MD
363 M_OUT.col(1) = M_MdXYTvd.col(1); // X
364 M_OUT.col(2) = M_MdXYTvd.col(2); // Y
365 M_OUT.col(3) = M_MdXYTvd.col(3).array() - kb; // TVDSS
366 M_OUT.col(4) = M_MdXYTvd.col(3); // TVD
367 M_OUT.col(5) = M_MdXYTvd.col(1).array() - x0; // DX
368 M_OUT.col(6) = M_MdXYTvd.col(2).array() - y0; // DY
369
370 if (!angularUnitsTo.empty()){
371 double coef = units::convert(
372 units::unit_from_string(angularUnitsFrom),
373 units::unit_from_string(angularUnitsTo));
374 M_OUT.col(7) = M.col(1)*coef; // AZ
375 M_OUT.col(8) = M.col(2)*coef; // INCL
376 } else {
377 M_OUT.col(7) = M.col(1); // AZ
378 M_OUT.col(8) = M.col(2); // INCL
379 }
380
381 return M_OUT;
382}
383
384template<typename D>
385Eigen::MatrixX<typename D::Scalar> h5geo::TvdXY2ALL(
386 const Eigen::DenseBase<D> &M,
387 const double& x0, const double& y0, const double& kb,
388 const std::string& angularUnits,
389 const bool& XNorth)
390{
391 if (M.cols() != 3)
392 return Eigen::MatrixX<typename D::Scalar>();
393
394 Eigen::MatrixX<typename D::Scalar> M_MdAzIncl =
395 TvdXY2MdAzIncl(M, x0, y0, angularUnits, XNorth);
396
397 Eigen::MatrixX<typename D::Scalar> M_OUT(M.rows(), 9);
398 M_OUT.col(0) = M_MdAzIncl.col(0); // MD
399 M_OUT.col(1) = M.col(1); // X
400 M_OUT.col(2) = M.col(2); // Y
401 M_OUT.col(3) = M.col(0).array() - kb; // TVDSS
402 M_OUT.col(4) = M.col(0); // TVD
403 M_OUT.col(5) = M.col(1).array() - x0; // DX
404 M_OUT.col(6) = M.col(2).array() - y0; // DY
405 M_OUT.col(7) = M_MdAzIncl.col(1); // AZ
406 M_OUT.col(8) = M_MdAzIncl.col(2); // INCL
407 return M_OUT;
408}
409
410template<typename D>
411Eigen::MatrixX<typename D::Scalar> h5geo::TvdDxDy2ALL(
412 const Eigen::DenseBase<D> &M,
413 const double& x0, const double& y0, const double& kb,
414 const std::string& angularUnits,
415 const bool& XNorth)
416{
417 if (M.cols() != 3)
418 return Eigen::MatrixX<typename D::Scalar>();
419
420 Eigen::MatrixX<typename D::Scalar> M_MdAzIncl =
421 TvdDxDy2MdAzIncl(M, angularUnits, XNorth);
422
423 Eigen::MatrixX<typename D::Scalar> M_OUT(M.rows(), 9);
424 M_OUT.col(0) = M_MdAzIncl.col(0); // MD
425 M_OUT.col(1) = M.col(1).array() + x0; // X
426 M_OUT.col(2) = M.col(2).array() + y0; // Y
427 M_OUT.col(3) = M.col(0).array() - kb; // TVDSS
428 M_OUT.col(4) = M.col(0); // TVD
429 M_OUT.col(5) = M.col(1); // DX
430 M_OUT.col(6) = M.col(2); // DY
431 M_OUT.col(7) = M_MdAzIncl.col(1); // AZ
432 M_OUT.col(8) = M_MdAzIncl.col(2); // INCL
433 return M_OUT;
434}
435
436template<typename D>
437Eigen::MatrixX<typename D::Scalar> h5geo::TvdssXY2ALL(
438 const Eigen::DenseBase<D> &M,
439 const double& x0, const double& y0, const double& kb,
440 const std::string& angularUnits,
441 const bool& XNorth)
442{
443 if (M.cols() != 3)
444 return Eigen::MatrixX<typename D::Scalar>();
445
446 Eigen::MatrixX<typename D::Scalar> M_MdAzIncl =
447 TvdssXY2MdAzIncl(M, x0, y0, kb, angularUnits, XNorth);
448
449 Eigen::MatrixX<typename D::Scalar> M_OUT(M.rows(), 9);
450 M_OUT.col(0) = M_MdAzIncl.col(0); // MD
451 M_OUT.col(1) = M.col(1); // X
452 M_OUT.col(2) = M.col(2); // Y
453 M_OUT.col(3) = M.col(0); // TVDSS
454 M_OUT.col(4) = M.col(0).array() + kb; // TVD
455 M_OUT.col(5) = M.col(1).array() - x0; // DX
456 M_OUT.col(6) = M.col(2).array() - y0; // DY
457 M_OUT.col(7) = M_MdAzIncl.col(1); // AZ
458 M_OUT.col(8) = M_MdAzIncl.col(2); // INCL
459 return M_OUT;
460}
461
462template<typename D>
463Eigen::MatrixX<typename D::Scalar> h5geo::TvdssDxDy2ALL(
464 const Eigen::DenseBase<D> &M,
465 const double& x0, const double& y0, const double& kb,
466 const std::string& angularUnits,
467 const bool& XNorth)
468{
469 if (M.cols() != 3)
470 return Eigen::MatrixX<typename D::Scalar>();
471
472 Eigen::MatrixX<typename D::Scalar> M_MdAzIncl =
473 TvdssDxDy2MdAzIncl(M, kb, angularUnits, XNorth);
474
475 Eigen::MatrixX<typename D::Scalar> M_OUT(M.rows(), 9);
476 M_OUT.col(0) = M_MdAzIncl.col(0); // MD
477 M_OUT.col(1) = M.col(1).array() + x0; // X
478 M_OUT.col(2) = M.col(2).array() + y0; // Y
479 M_OUT.col(3) = M.col(0); // TVDSS
480 M_OUT.col(4) = M.col(0).array() + kb; // TVD
481 M_OUT.col(5) = M.col(1); // DX
482 M_OUT.col(6) = M.col(2); // DY
483 M_OUT.col(7) = M_MdAzIncl.col(1); // AZ
484 M_OUT.col(8) = M_MdAzIncl.col(2); // INCL
485 return M_OUT;
486}
487
488template<typename D>
489Eigen::MatrixX<typename D::Scalar> h5geo::traj2ALL(
490 const Eigen::DenseBase<D> &M,
491 const double& x0, const double& y0, const double& kb,
492 const std::string& angularUnitsFrom,
493 const std::string& angularUnitsTo,
494 const h5geo::TrajectoryFormat& trajFormat,
495 const bool& XNorth)
496{
497 switch (trajFormat) {
498 case h5geo::TrajectoryFormat::MD_AZIM_INCL:
499 return h5geo::MdAzIncl2ALL(M, x0, y0, kb, angularUnitsFrom, angularUnitsTo, XNorth);
500 case h5geo::TrajectoryFormat::TVDSS_DX_DY:
501 return h5geo::TvdssDxDy2ALL(M, x0, y0, kb, angularUnitsTo, XNorth);
502 case h5geo::TrajectoryFormat::TVDSS_X_Y:
503 return h5geo::TvdssXY2ALL(M, x0, y0, kb, angularUnitsTo, XNorth);
504 case h5geo::TrajectoryFormat::TVD_DX_DY:
505 return h5geo::TvdDxDy2ALL(M, x0, y0, kb, angularUnitsTo, XNorth);
506 case h5geo::TrajectoryFormat::TVD_X_Y:
507 return h5geo::TvdXY2ALL(M, x0, y0, kb, angularUnitsTo, XNorth);
508 default:
509 return Eigen::MatrixX<typename D::Scalar>();
510 }
511}
512
513template<typename D>
514Eigen::MatrixX<typename D::Scalar> h5geo::MdAzIncl2MdXYTvd(
515 const Eigen::DenseBase<D> &M,
516 const double& x0, const double& y0,
517 const std::string& angularUnits,
518 const bool& XNorth)
519{
520 if (M.cols() != 3)
521 return Eigen::MatrixX<typename D::Scalar>();
522
523 double coef = 1;
524 if (!angularUnits.empty()){
525 coef = units::convert(
526 units::unit_from_string(angularUnits),
527 units::unit_from_string("radian"));
528 }
529
530 Eigen::MatrixX<typename D::Scalar> M_OUT(M.rows(), 4);
531 M_OUT.col(0) = M.col(0);
532 typename D::Scalar I1, I2, A1, A2, dMD, B, RF, dx, dy, dz;
533
534 for (size_t i = 0; i < M.rows(); i++){
535 if (i == 0){
536 A1 = M(0, 1) * coef;
537 A2 = M(0, 1) * coef;
538 I1 = M(0, 2) * coef;
539 I2 = M(0, 2) * coef;
540 dMD = M(0, 0);
541 } else {
542 A1 = M(i, 1) * coef;
543 A2 = M(i-1, 1) * coef;
544 I1 = M(i, 2) * coef;
545 I2 = M(i-1, 2) * coef;
546 dMD = M(i, 0) - M(i-1, 0);
547 }
548
549 B = _betaAngle(I1, I2, A1, A2, "");
550 RF = _ratioFactor(B, "");
551
552 if (XNorth){
553 dx = _deltaNorth(dMD, I1, I2, A1, A2, RF, "");
554 dy = _deltaEast(dMD, I1, I2, A1, A2, RF, "");
555 } else {
556 dx = _deltaEast(dMD, I1, I2, A1, A2, RF, "");
557 dy = _deltaNorth(dMD, I1, I2, A1, A2, RF, "");
558 }
559
560 dz = _deltaZ(dMD, I1, I2, RF, "");
561
562 if (i == 0){
563 M_OUT(0, 1) = x0 + dx;
564 M_OUT(0, 2) = y0 + dy;
565 M_OUT(0, 3) = dz;
566 } else {
567 M_OUT(i, 1) = M_OUT(i-1, 1) + dx;
568 M_OUT(i, 2) = M_OUT(i-1, 2) + dy;
569 M_OUT(i, 3) = M_OUT(i-1, 3) + dz;
570 }
571 }
572 return M_OUT;
573}
574
575template<typename D>
576Eigen::MatrixX<typename D::Scalar> h5geo::TvdXY2MdXYTvd(
577 const Eigen::DenseBase<D> &M,
578 const double& x0, const double& y0,
579 const bool& XNorth)
580{
581 if (M.cols() != 3)
582 return Eigen::MatrixX<typename D::Scalar>();
583
584 Eigen::MatrixX<typename D::Scalar> M_MdAzIncl =
585 TvdXY2MdAzIncl(M, x0, y0, "", XNorth);
586
587 Eigen::MatrixX<typename D::Scalar> M_MdXYTvd(M.rows(), 4);
588
589 M_MdXYTvd.col(0) = M_MdAzIncl.col(0);
590 M_MdXYTvd.col(1) = M.col(1);
591 M_MdXYTvd.col(2) = M.col(2);
592 M_MdXYTvd.col(3) = M.col(0);
593
594 return M_MdXYTvd;
595}
596
597template<typename D>
598Eigen::MatrixX<typename D::Scalar> h5geo::TvdDxDy2MdXYTvd(
599 const Eigen::DenseBase<D> &M,
600 const double& x0, const double& y0,
601 const bool& XNorth)
602{
603 if (M.cols() != 3)
604 return Eigen::MatrixX<typename D::Scalar>();
605
606 Eigen::MatrixX<typename D::Scalar> M_MdAzIncl =
607 TvdDxDy2MdAzIncl(M, "", XNorth);
608
609 Eigen::MatrixX<typename D::Scalar> M_MdXYTvd(M.rows(), 4);
610
611 M_MdXYTvd.col(0) = M_MdAzIncl.col(0);
612 M_MdXYTvd.col(1) = M.col(1).array() + x0;
613 M_MdXYTvd.col(2) = M.col(2).array() + y0;
614 M_MdXYTvd.col(3) = M.col(0);
615
616 return M_MdXYTvd;
617}
618
619template<typename D>
620Eigen::MatrixX<typename D::Scalar> h5geo::TvdssXY2MdXYTvd(
621 const Eigen::DenseBase<D> &M,
622 const double& x0, const double& y0, const double& kb,
623 const bool& XNorth)
624{
625 if (M.cols() != 3)
626 return Eigen::MatrixX<typename D::Scalar>();
627
628 Eigen::MatrixX<typename D::Scalar> M_MdAzIncl =
629 TvdssXY2MdAzIncl(M, x0, y0, kb, "", XNorth);
630
631 Eigen::MatrixX<typename D::Scalar> M_MdXYTvd(M.rows(), 4);
632
633 M_MdXYTvd.col(0) = M_MdAzIncl.col(0);
634 M_MdXYTvd.col(1) = M.col(1);
635 M_MdXYTvd.col(2) = M.col(2);
636 M_MdXYTvd.col(3) = M.col(0).array() + kb;
637
638 return M_MdXYTvd;
639}
640
641template<typename D>
642Eigen::MatrixX<typename D::Scalar> h5geo::TvdssDxDy2MdXYTvd(
643 const Eigen::DenseBase<D> &M,
644 const double& x0, const double& y0, const double& kb,
645 const bool& XNorth)
646{
647 if (M.cols() != 3)
648 return Eigen::MatrixX<typename D::Scalar>();
649
650 Eigen::MatrixX<typename D::Scalar> M_MdAzIncl =
651 TvdssDxDy2MdAzIncl(M, kb, "", XNorth);
652
653 Eigen::MatrixX<typename D::Scalar> M_MdXYTvd(M.rows(), 4);
654
655 M_MdXYTvd.col(0) = M_MdAzIncl.col(0);
656 M_MdXYTvd.col(1) = M.col(1).array() + x0;
657 M_MdXYTvd.col(2) = M.col(2).array() + y0;
658 M_MdXYTvd.col(3) = M.col(0).array() + kb;
659
660 return M_MdXYTvd;
661}
662
663template<typename D>
664Eigen::MatrixX<typename D::Scalar> h5geo::TvdXY2MdAzIncl(
665 const Eigen::DenseBase<D> &M,
666 const double& x0, const double& y0,
667 const std::string& angularUnits,
668 const bool& XNorth)
669{
670 if (M.cols() != 3)
671 return Eigen::MatrixX<typename D::Scalar>();
672
673 Eigen::MatrixX<typename D::Scalar> MM = M;
674 MM.col(1) = MM.col(1).array() - x0;
675 MM.col(2) = MM.col(2).array() - y0;
676
677 return TvdDxDy2MdAzIncl(MM, angularUnits, XNorth);
678}
679
680template<typename D>
681Eigen::MatrixX<typename D::Scalar> h5geo::TvdDxDy2MdAzIncl(
682 const Eigen::DenseBase<D> &M, const std::string& angularUnits, const bool& XNorth)
683{
684 if (M.cols() != 3)
685 return Eigen::MatrixX<typename D::Scalar>();
686
687 Eigen::MatrixX<typename D::Scalar> M_OUT(M.rows(), 3);
688 typename D::Scalar I1, I2, A1, A2, dMD, B, RF, dx, dy, dz;
689
690 for (size_t i = 0; i < M.rows(); i++){
691 if (i == 0){
692 dx = M(0, 1);
693 dy = M(0, 2);
694 dz = M(0, 0);
695
696 if (XNorth){
697 A1 = atan2(dy, dx);
698 } else {
699 A1 = atan2(dx, dy);
700 }
701
702 A2 = A1;
703 I1 = atan2(sqrt(pow(dx,2)+pow(dy,2)), dz);
704 I2 = I1;
705 } else {
706 dx = M(i, 1) - M(i-1, 1);
707 dy = M(i, 2) - M(i-1, 2);
708 dz = M(i, 0) - M(i-1, 0);
709 A1 = M_OUT(i-1, 1);
710
711 if (XNorth){
712 A2 = atan2(dy, dx);
713 } else {
714 A2 = atan2(dx, dy);
715 }
716
717 I1 = M_OUT(i-1, 2);
718 I2 = atan2(sqrt(pow(dx,2)+pow(dy,2)), dz);
719 }
720
721 B = _betaAngle(I1, I2, A1, A2, "");
722 RF = _ratioFactor(B, "");
723
724 dMD = _dMDFromDeltaZ(dz, I1, I2, RF, "");
725
726 if (i == 0)
727 M_OUT(0, 0) = dMD;
728 else
729 M_OUT(i, 0) = M_OUT(i-1, 0) + dMD;
730
731 M_OUT(i, 1) = A2;
732 M_OUT(i, 2) = I2;
733 }
734
735 if (!angularUnits.empty()){
736 double coef = units::convert(
737 units::unit_from_string("radian"),
738 units::unit_from_string(angularUnits));
739 M_OUT.col(1) = M_OUT.col(1).array() * coef;
740 M_OUT.col(2) = M_OUT.col(2).array() * coef;
741 }
742
743 return M_OUT;
744}
745
746template<typename D>
747Eigen::MatrixX<typename D::Scalar> h5geo::TvdssXY2MdAzIncl(
748 const Eigen::DenseBase<D> &M,
749 const double& x0, const double& y0, const double& kb,
750 const std::string& angularUnits,
751 const bool& XNorth)
752{
753 if (M.cols() != 3)
754 return Eigen::MatrixX<typename D::Scalar>();
755
756 Eigen::MatrixX<typename D::Scalar> MM = M;
757 MM.col(0) = MM.col(0).array() + kb;
758 MM.col(1) = MM.col(1).array() - x0;
759 MM.col(2) = MM.col(2).array() - y0;
760
761 return TvdDxDy2MdAzIncl(MM, angularUnits, XNorth);
762}
763
764template<typename D>
765Eigen::MatrixX<typename D::Scalar> h5geo::TvdssDxDy2MdAzIncl(
766 const Eigen::DenseBase<D> &M,
767 const double& kb, const std::string& angularUnits, const bool& XNorth)
768{
769 if (M.cols() != 3)
770 return Eigen::MatrixX<typename D::Scalar>();
771
772 Eigen::MatrixX<typename D::Scalar> MM = M;
773 MM.col(0) = MM.col(0).array() + kb;
774
775 return TvdDxDy2MdAzIncl(MM, angularUnits, XNorth);
776}
777
778
779Eigen::VectorX<double> h5geo::angleAndRadius(
780 const double& arcLen,
781 const double& chordLen,
782 const double& eps,
783 const size_t& max_iter)
784{
785 if (arcLen <= chordLen ||
786 arcLen == 0 ||
787 chordLen == 0)
788 return Eigen::VectorX<double>();
789
790 double a_min = 0;
791 double a_max = 2*M_PI;
792 double a0, a1, f, df, ddf;
793 for (a0 = a_min; a0 < a_max; a0 = a0+0.1){
794 f = _fa(arcLen, chordLen, a0, "");
795 ddf = _ddfa(arcLen, a0, "");
796 if (f*ddf > 0)
797 break;
798 }
799
800 df = _dfa(arcLen, chordLen, a0, "");
801 if (df == 0)
802 return Eigen::VectorX<double>();
803
804 a1 = a0 - f/df;
805 size_t i = 1, ii = 0;
806 while (abs(f) > eps && i < max_iter){
807 f = _fa(arcLen, chordLen, a0, "");
808 df = _dfa(arcLen, chordLen, a0, "");
809 if (df == 0)
810 return Eigen::VectorX<double>();
811
812 a0 = a1;
813 a1 = a0 - f/df;
814 if (a1 < a_min || a1 > a_max){
815 while (a1 < a_min || a1 > a_max){
816 a1 = (a0+a1)/2;
817 ii++;
818 if (ii == max_iter)
819 return Eigen::VectorX<double>();
820 }
821 ii = 0;
822 }
823 i++;
824 }
825 double radius = arcLen/a1;
826
827 Eigen::Vector2<double> v;
828 v(0) = a1;
829 v(1) = radius;
830 return v;
831}
832
833#endif // H5DEVIATION_H
h5geo
Basic namespace.
Definition h5base.h:29
h5geo::angleAndRadius
Eigen::VectorX< double > angleAndRadius(const double &arcLen, const double &chordLen, const double &eps, const size_t &max_iter)
angleAndRadius Calculate angle and radius from arc and chord
Definition h5deviation.h:779
Generated by doxygen 1.11.0