RSTT 3.2.1 Regional Seismic Travel Time

Loading...
Searching...
No Matches
IntegrateFunction.h
Go to the documentation of this file.
1//- ****************************************************************************
2//-
3//- Copyright 2009 National Technology & Engineering Solutions of Sandia, LLC
4//- (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the U.S.
5//- Government retains certain rights in this software.
6//-
7//- BSD Open Source License
8//- All rights reserved.
9//-
10//- Redistribution and use in source and binary forms, with or without
11//- modification, are permitted provided that the following conditions are met:
12//-
13//- 1. Redistributions of source code must retain the above copyright notice,
14//- this list of conditions and the following disclaimer.
15//-
16//- 2. Redistributions in binary form must reproduce the above copyright
17//- notice, this list of conditions and the following disclaimer in the
18//- documentation and/or other materials provided with the distribution.
19//-
20//- 3. Neither the name of the copyright holder nor the names of its
21//- contributors may be used to endorse or promote products derived from
22//- this software without specific prior written permission.
23//-
24//- THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
25//- AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
26//- IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
27//- ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
28//- LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
29//- CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
30//- SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
31//- INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
32//- CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
33//- ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
34//- POSSIBILITY OF SUCH DAMAGE.
35//-
36//- ****************************************************************************
37
38#ifndef INTEGRATEFUNCTION_H
39#define INTEGRATEFUNCTION_H
40
41// **** _SYSTEM INCLUDES_ ******************************************************
42
43#include <float.h>
44#include <vector>
45#include <iostream>
46#include <limits>
47
48 using namespace std;
49// use standard library objects
50
51// **** _LOCAL INCLUDES_ *******************************************************
52
53#include "UtilGlobals.h"
54
55// **** _BEGIN UTIL NAMESPACE_ *************************************************
56
57 namespace util {
58
59// **** _CLASS CONSTANTS_ ******************************************************
60// **** _FORWARD REFERENCES_ ***************************************************
61
62// **** _CLASS DEFINITION_ *****************************************************
63//
69//
110// *****************************************************************************
111 template<class F>
112 class UTIL_EXP IntegrateFunction
113{
114 public:
115
116 // **** _PUBLIC LIFCYCLES_ *************************************************
117
122 IntegrateFunction(F& f, double tol);
123
125 IntegrateFunction(const IntegrateFunction<F>& ifctn);
126
128 virtual ~IntegrateFunction();
129
130 // **** _PUBLIC OPERATORS_ *************************************************
131
133 IntegrateFunction<F>& operator=(const IntegrateFunction<F>& ifctn);
134
135 // **** _PUBLIC METHODS_ ***************************************************
136
141 double integrateClosed(double a, double b);
142
148 double integrateAOpenS(double a, double b);
149
151 F& getF();
152
154 const F& getF() const;
155
157 void setTolerance(double tol);
158
160 double getTolerance() const;
161
163 //static double getTest(int f, double ftrn, double tol);
164
165 private:
166
167 // **** _PRIVATE METHODS_ **************************************************
168
173 //
177 double integrateClosedRcrsv(double a, double b, double* fr);
178
180 //
201 double simpson(double a, double b, double* fa, double& s);
202
203 // **** _PRIVATE DATA_ *****************************************************
204
206 double ifTol;
207
209 F& ifF;
210};
211
212// **** _INLINE FUNCTION IMPLEMENTATIONS_ **************************************
213
214// Returns a reference to the integrated function.
215 template<class F>
216 inline F& IntegrateFunction<F>::getF()
217{
218 return ifF;
219}
220
221// Returns a const reference to the integrated function.
222 template<class F>
223 inline const F& IntegrateFunction<F>::getF() const
224{
225 return ifF;
226}
227
228// Sets the integration tolerance value.
229 template<class F>
230 inline void IntegrateFunction<F>::setTolerance(double tol)
231{
232 ifTol = tol;
233}
234
235// Returns the integration tolerance value.
236 template<class F>
237 inline double IntegrateFunction<F>::getTolerance() const
238{
239 return ifTol;
240}
241
242// Adaptive closed form numerical integration of the input function
243// over the closed interval a to b. The function is assumed to be
244// defined over the entire interval [a, b]. The function uses the 3 point
245// Simpsons rule to evaluate all sub-intervals.
246 template<class F>
247 inline double IntegrateFunction<F>::integrateClosed(double a, double b)
248{
249 double fr[3];
250
251 // evaluate the function at a, b, and the midpoint and save into fr.
252
253 fr[0] = ifF(a);
254 fr[1] = ifF(0.5 * (a + b));
255 fr[2] = ifF(b);
256
257 // integrate and return result
258
259 return integrateClosedRcrsv(a, b, fr);
260}
261
262 template<class F>
263 inline double IntegrateFunction<F>::integrateAOpenS(double a, double b)
264{
265 double fs, fa, ae, be, senew;
266
267 // set the integration result to 0.0 and calculate a small value near a
268 // but slightly larger (e). Set the smallest legal value for e (minstep)
269 // and account for precision loss in the step size due to significant
270 // digits in front of the decimal place (of a).
271
272 double s = 0.0;
273 double e = (b - a) * ifTol;
274 fa = fabs(a);
275 double minstep = 10.0 * DBL_EPSILON;
276 if (fa > 1.0) minstep *= fa;
277
278 // calculate integral from just shy of the asymptote to b and
279 // loop over successive steps getting closer to a each step
280 // (e to e/10 each loop). Loop while the stepsize is decremented by
281 // an order of magnitude at each step.
282
283 ae = a + e;
284 s = integrateClosed(ae, b);
285 be = ae; ae = a + 0.1 * e;
286 do
287 {
288 // get contribution from a+e/10 to a+e ... sum contribution to s
289
290 senew = integrateClosed(ae, be);
291 s += senew;
292
293 // see if error criteria has been met ... if so return s, otherwise
294 // decrement e by 10 and try again
295
296 fs = fabs(s);
297 if ((fabs(senew) < ifTol * fs) || (fs < ifTol))
298 return s;
299 else
300 e /= 10.0;
301
302 // continue while e is larger than minimum step size
303
304 be = ae; ae = a + 0.1 * e;
305 } while ((e > minstep) && (be > ae) && (ae > a));
306
307 // output error tolerance condtion and return best estimate
308
309 if (fabs(senew) > ifTol)
310 {
311 // should throw here ... need to resume after the fact
312 // should output a, b, e, ae, be, s, senew, iftol
313 cout << " Error:: Function Error Tolerance Exceeded ... " << endl
314 << " Tolerance Condition Was Not Met." << endl;
315 }
316 return s;
317}
318
319// Adaptive closed form numerical integration of the input function
320// over the closed interval a to b. The function is assumed to be
321// defined over the entire interval [a, b]. The function uses the 3 point
322// Simpsons rule to evaluate all sub-intervals.
323//
324// The array fr contains previous function evaluations at ifF(a),
325// ifF((a + b)/2), and ifF(b). This private function is only called by the
326// public cover function integrateClosed(a, b).
327 template<class F>
328 inline double IntegrateFunction<F>::integrateClosedRcrsv(double a, double b,
329 double* fr)
330{
331 // set the integration result to 0.0 and save the input previously
332 // calculated function evalutions in fr into fa. fa[1] and fa[3] are
333 // calculated in function simpson(...) below.
334
335 double s = 0.0;
336 double fa[5];
337 fa[0] = fr[0];
338 fa[2] = fr[1];
339 fa[4] = fr[2];
340
341 // evaluate Simpsons rule and return result if tolerance is met.
342
343 double se = simpson(a, b, fa, s);
344 double fs = fabs(s);
345 if ((se < ifTol * fs) || (fs < ifTol))
346 return s;
347 else
348 {
349 // not accurate enough yet ... calculate interval mid-point and verify
350 // machine precison will allow evaluation.
351
352 double m = 0.5 * (a + b);
353 if ((m <= a) | (b <= m))
354 {
355 if (se > ifTol)
356 {
357 // write out warning
358 // should throw error here ... write out a, b, m, s, se, iftol
359 cout << " Error:: Function Error Tolerance Exceeded ... " << endl
360 << " Tolerance Condition Was Not Met." << endl;
361 }
362 return s;
363 }
364
365 // integrate over left and right sub-intervals ... return the summed
366 // result
367
368 return integrateClosedRcrsv(a, m, &fa[0]) +
369 integrateClosedRcrsv(m, b, &fa[2]);
370 }
371}
372
373// A 3 point 4th order Simpsons integration rule.
374//
375// (h/3)(f0 + 4*f1 + f2)
376//
377// where h is half the interval width, f0 is the function to be integrated
378// defined at the left most range, f2 is defined at the right most range,
379// and f1 is defined at the mid-point of the range.
380//
381// The rule evaluates the integral across a sub-interval range using the
382// rule 3 times, once across the sub-interval, once across the 1st half of
383// the sub-interval and once across the second half of the sub-interval.
384// The split halves are summed to form one result while the entire
385// interval evaluation forms the other. The former has an accuracy that is
386// about 4 times greater than the entire interval evaluation. Convergence
387// is attained when the difference between the two is less than the input
388// tolerance.
389//
390// Five function evaluations are required for the closed form solution.
391// However, 3 evaluations are passed from the previous call so only 2 new
392// function evalutions are required per call to this function. The 3
393// previous functions evaluations are provided in the array fa at
394// locations fa[0], fa[2], and fa[4]. Evaltions for fa[1] and fa[3] are
395// performed in this function and saved into the array.
396 template<class F>
397 inline double IntegrateFunction<F>::simpson(double a, double b,
398 double* fa, double& s)
399{
400 // calculate interval width and the two remaining function evalutions
401
402 double h = b - a;
403 fa[1] = ifF(a + 0.25 * h);
404 fa[3] = ifF(a + 0.75 * h);
405 h *= 0.5;
406
407 // calculate results over entire interval (s1) and the interval split in
408 // half (s). Return s and the difference between both evaluations.
409
410 double dh = h / 6.0;
411 double s1 = 2.0 * dh * (fa[0] + 4.0 * fa[2] + fa[4]);
412 s = dh * (fa[0] + 4.0 * (fa[1] + fa[3]) + 2.0 * fa[2] + fa[4]);
413 return fabs(s - s1);
414}
415
416// **** _TEST FUNCTION IMPLEMENTATIONS_ **************************************
417
420//
425
426// f = velocity function
427// ftrn = is fractional radius where rtrn = ftrn * (rtop - rbot) + rbot
428// fa = is fractional radius where ra = fa * (rtop - rtrn) + rtrn above the turn radius
429// fb = is fractional radius where rb = fb * (rtop - rtrn) + rtrn above the turn radius
430// (Note: fb > fa, fb == 1 gives rb == rtop)
431// tif.setF(f);
432// tif.setRTurn(ftrn);
433// double a = tif.getA(fa);
434// double b = tif.getB(fb);
435//class TestIntegrateFunction
436//{
437// public:
438//
439// //! Default Constructor.
440// TestIntegrateFunction() {f=1;};
441//
442// //! Copy Constructor.
443// TestIntegrateFunction(const TestIntegrateFunction& tif) {f=1;};
444//
445// //! Destructor.
446// ~TestIntegrateFunction() {f=1;};
447//
448// //! Assignment Operator.
449// TestIntegrateFunction& operator=(const TestIntegrateFunction& tif)
450// {return *this;};
451//
452// //! \brief Function parentheses operator(). Returns a double given a
453// //! single double argument. This function is the integrand.
454// double operator()(double r)
455// {
456// ++fccnt;
457// double vr = v(r);
458// double pv = p * vr;
459//
460// //return sqrt((r - pv) * (r + pv)) / r / vr; // Tau Function
461// return 1.0 / sqrt((r - pv) * (r + pv)); // Turning Asymptote
462// //return pv / sqrt((r - pv) * (r + pv)) / r; // Distance
463// //return r / sqrt((r - pv) * (r + pv)) / vr; // Time
464// }
465//
466// //! Zeros the function call count.
467// void resetFCCount() {fccnt = 0;};
468//
469// //! Returns the function call count.
470// int getFCCount() const {return fccnt;};
471//
472// //! Set velocity model
473// void setF(int fin) {f = fin;};
474//
475// //! Set turning radius
476// void setRTurn(double ftrn)
477// {
478// rturn = ftrn * (rbot() - rtop()) + rtop();
479// p = rturn / v(rturn);
480// };
481//
482// //! Get a limit.
483// double getLimit(double fin) const
484// {
485// return fin * (rtop() - rturn) + rturn;
486// };
487//
488// // mantle1 p vel from r = 6336 to 6251
489// double v1(double r) const
490// {
491// double rn = r / 6371.0;
492// return 8.78541 - 0.74953 * rn;
493// };
494//
495// // mantle2 p vel from r = 6251 to 6161
496// double v2(double r) const
497// {
498// double rn = r / 6371.0;
499// return 25.41389 - 17.69722 * rn;
500// };
501//
502// // mantle3 p vel from r = 6161 to 5961
503// double v3(double r) const
504// {
505// double rn = r / 6371.0;
506// return 30.78765 - 23.25415 * rn;
507// };
508//
509// // mantle4 p vel from r = 5961 to 5711
510// double v4(double r) const
511// {
512// double rn = r / 6371.0;
513// return 29.38896 - 21.40656 * rn;
514// };
515//
516// // mantle5 p vel from r = 5711 to 5611
517// double v5(double r) const
518// {
519// double rn = r / 6371.0;
520// return 25.96984 - 16.93412 * rn;
521// };
522//
523// // mantle p vel from r = 5611 to 3631
524// double v6(double r) const
525// {
526// double rn = r / 6371.0;
527// return 25.1486 + rn * (-41.1538 + rn * (51.9932 + rn * -26.6083));
528// };
529//
530// // c-m trns p vel from r = 3631 to 3482
531// double v7(double r) const
532// {
533// double rn = r / 6371.0;
534// return 14.49470 - 1.47089 * rn;
535// };
536//
537// // outer core p vel from r = 3482 to 1217.1
538// double v8(double r) const
539// {
540// double rn = r / 6371.0;
541// return 10.03904 + rn * (3.75665 + rn * -13.67046);
542// };
543//
544// // inner core p vel from r = 1217.1 to 0
545// double v9(double r) const
546// {
547// double rn = r / 6371.0;
548// return 11.24094 + rn * rn * -4.09689;
549// };
550//
551// double rtop() const
552// {
553// if (f < 5)
554// {
555// if (f < 3)
556// {
557// if (f == 1)
558// return 6336.0;
559// else
560// return 6251.0;
561// }
562// else if (f == 3)
563// return 6161.0;
564// else
565// return 5961.0;
566// }
567// else if (f < 7)
568// {
569// if (f == 5)
570// return 5711.0;
571// else
572// return 5611.0;
573// }
574// else if (f < 9)
575// {
576// if (f == 7)
577// return 3631.0;
578// else
579// return 3482.0;
580// }
581// else
582// return 1217.1;
583//
584// return 0.0;
585// };
586//
587// double rbot() const
588// {
589// if (f < 5)
590// {
591// if (f < 3)
592// {
593// if (f == 1)
594// return 6251.0;
595// else
596// return 6161.0;
597// }
598// else if (f == 3)
599// return 5961.0;
600// else
601// return 5711.0;
602// }
603// else if (f < 7)
604// {
605// if (f == 5)
606// return 5611.0;
607// else
608// return 3631.0;
609// }
610// else if (f < 9)
611// {
612// if (f == 7)
613// return 3482.0;
614// else
615// return 1217.1;
616// }
617// else
618// return 0.0;
619//
620// return 0.0;
621// };
622//
623// double v(double r) const
624// {
625// if (f < 5)
626// {
627// if (f < 3)
628// {
629// if (f == 1)
630// return v1(r);
631// else
632// return v2(r);
633// }
634// else if (f == 3)
635// return v3(r);
636// else
637// return v4(r);
638// }
639// else if (f < 7)
640// {
641// if (f == 5)
642// return v5(r);
643// else
644// return v6(r);
645// }
646// else if (f < 9)
647// {
648// if (f == 7)
649// return v7(r);
650// else
651// return v8(r);
652// }
653// else
654// return v9(r);
655//
656// return 0.0;
657// };
658//
659// private:
660//
661// //! \brief The integrand function call count. Used to examine the
662// //! number of times the operator()(double x) function is called.
663// int fccnt;
664//
665// int f;
666// double p;
667//
668// double rturn;
669//};
670
671} // end namespace util
672
673#endif // INTEGRATEFUNCTION_H
UtilGlobals.h
UTIL_EXP
#define UTIL_EXP
Definition UtilGlobals.h:66
util::IntegrateFunction::setTolerance
void setTolerance(double tol)
Sets the integration tolerance value.
Definition IntegrateFunction.h:230
util::IntegrateFunction::integrateAOpenS
double integrateAOpenS(double a, double b)
Adaptive closed form numerical integration of the input function over the open interval a to b,...
Definition IntegrateFunction.h:263
util::IntegrateFunction::IntegrateFunction
IntegrateFunction(F &f, double tol)
Standard Constructor.
util::IntegrateFunction::operator=
IntegrateFunction< F > & operator=(const IntegrateFunction< F > &ifctn)
Assignment Operator.
util::IntegrateFunction::IntegrateFunction
IntegrateFunction(const IntegrateFunction< F > &ifctn)
Copy Constructor.
util::IntegrateFunction::~IntegrateFunction
virtual ~IntegrateFunction()
Destructor.
util::IntegrateFunction::getTolerance
double getTolerance() const
Returns the integration tolerance value.
Definition IntegrateFunction.h:237
util::IntegrateFunction::integrateClosed
double integrateClosed(double a, double b)
Adaptive closed form numerical integration of the input function over the closed interval a to b....
Definition IntegrateFunction.h:247
util::IntegrateFunction::getF
F & getF()
Returns a reference to the integrated function.
Definition IntegrateFunction.h:216
util
Definition Brents.h:54