ARTS 2.5.10 (git: 2f1c442c)
math_funcs.cc
Go to the documentation of this file.
2 Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3 Stefan Buehler <sbuehler@ltu.se>
4
5 This program is free software; you can redistribute it and/or modify it
7 Free Software Foundation; either version 2, or (at your option) any
8 later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18 USA. */
19
20/*****************************************************************************
21 *** File description
22 *****************************************************************************/
23
32/*****************************************************************************
33 *** External declarations
34 *****************************************************************************/
35
36#include "math_funcs.h"
37#include <cmath>
38#include <iostream>
39#include <stdexcept>
40#include "array.h"
41#include "arts_constants.h"
42#include "arts_conversions.h"
43#include "logic.h"
44#include "mystring.h"
45
47inline constexpr Numeric PI=Constant::pi;
48
49/*****************************************************************************
50 *** The functions (in alphabetical order)
51 *****************************************************************************/
52
54
65Numeric fac(const Index n) {
66 Numeric sum;
67
68 if (n == 0) return (1.0);
69
70 sum = 1.0;
71 for (Index i = 1; i <= n; i++) sum *= Numeric(i);
72
73 return (sum);
74}
75
77
89Index integer_div(const Index& x, const Index& y) {
91 return x / y;
92}
93
95
117 const Numeric a) {
118 // lowermost grid spacing on x-axis
119 const Numeric Dlimit = 1.00000e-15;
120
121 // Check that dimensions of x and y vector agree
122 const Index n_x = x.nelem();
123 const Index n_y = y.nelem();
124 ARTS_USER_ERROR_IF ((n_x != 4) || (n_y != 4),
125 "The vectors x and y must all have the same length of 4 elements!\n"
126 "Actual lengths:\n"
127 "x:", n_x, ", "
128 "y:", n_y, ".")
129
130 // assure that x1 =< a < x2 holds
131 ARTS_USER_ERROR_IF ((a < x[1]) || (a > x[2]),
132 "LagrangeInterpol4: the relation x[1] =< a < x[2] is not satisfied. "
133 "No interpolation can be calculated.\n")
134
135 // calculate the Lagrange polynomial coefficients for a polynomial of the order of 3
136 Numeric b[4];
137 for (Index i = 0; i < 4; ++i) {
138 b[i] = 1.000e0;
139 for (Index k = 0; k < 4; ++k) {
140 if ((k != i) && (fabs(x[i] - x[k]) > Dlimit))
141 b[i] = b[i] * ((a - x[k]) / (x[i] - x[k]));
142 };
143 };
144
145 Numeric ya = 0.000e0;
146 for (Index i = 0; i < n_x; ++i) ya = ya + b[i] * y[i];
147
148 return ya;
149}
150
152
162 ARTS_ASSERT(x.nelem() > 0);
163 return x[x.nelem() - 1];
164}
165
167
177 ARTS_ASSERT(x.nelem() > 0);
178 return x[x.nelem() - 1];
179}
180
182
201 const Numeric start,
202 const Numeric stop,
203 const Numeric step) {
204 Index n = (Index)floor((stop - start) / step) + 1;
205 if (n < 1) n = 1;
206 x.resize(n);
207 for (Index i = 0; i < n; i++) x[i] = start + (double)i * step;
208}
209
211
228 const Numeric start,
229 const Numeric stop,
230 const Index n) {
231 ARTS_ASSERT(1 < n); // Number of points must be greater 1.
232 x.resize(n);
233 Numeric step = (stop - start) / ((double)n - 1);
234 for (Index i = 0; i < n - 1; i++) x[i] = start + (double)i * step;
235 x[n - 1] = stop;
236}
238 const Numeric start,
239 const Numeric stop,
240 const Index n) {
241 Numeric step = (stop - start) / ((double)n - 1);
242 for (Index i = 0; i < n - 1; i++) x[i] = start + (double)i * step;
243 x[n - 1] = stop;
244}
245
247
264 const Numeric start,
265 const Numeric stop,
266 const Index n) {
267 // Number of points must be greater than 1:
268 ARTS_ASSERT(1 < n);
269 // Only positive numbers are allowed for start and stop:
270 ARTS_ASSERT(0 < start);
271 ARTS_ASSERT(0 < stop);
272
273 x.resize(n);
274 Numeric a = log(start);
275 Numeric step = (log(stop) - a) / ((double)n - 1);
276 x[0] = start;
277 for (Index i = 1; i < n - 1; i++) x[i] = exp(a + (double)i * step);
278 x[n - 1] = stop;
279}
280
282
295{
296 const Index n = x.nelem();
297 ARTS_ASSERT(y.nelem() == n);
298 Numeric sum = 0.0;
299 for (Index i=1; i<n; ++i)
300 sum += (x[i]-x[i-1]) * (y[i]+y[i-1]);
301 return sum/2.0;
302}
303
305
318{
319 const Index n = x.nelem();
320 ARTS_ASSERT(csum.nelem() == n);
321 csum[0] = x[0];
322 for (Index i=1; i<n; ++i)
323 csum[i] = csum[i-1] + x[i];
324}
325
327
338 ConstVectorView za_grid,
339 ConstVectorView aa_grid) {
340 Index n = za_grid.nelem();
341 Index m = aa_grid.nelem();
342 Vector res1(n);
343 ARTS_ASSERT(is_size(Integrand, n, m));
344
345 for (Index i = 0; i < n; ++i) {
346 res1[i] = 0.0;
347
348 for (Index j = 0; j < m - 1; ++j) {
349 res1[i] += 0.5 * DEG2RAD * (Integrand(i, j) + Integrand(i, j + 1)) *
350 (aa_grid[j + 1] - aa_grid[j]) * sin(za_grid[i] * DEG2RAD);
351 }
352 }
353 Numeric res = 0.0;
354 for (Index i = 0; i < n - 1; ++i) {
355 res +=
356 0.5 * DEG2RAD * (res1[i] + res1[i + 1]) * (za_grid[i + 1] - za_grid[i]);
357 }
358
359 return res;
360}
361
363
382 ConstVectorView za_grid,
383 ConstVectorView aa_grid,
384 ConstVectorView grid_stepsize) {
385 Numeric res = 0;
386 if ((grid_stepsize[0] > 0) && (grid_stepsize[1] > 0)) {
387 Index n = za_grid.nelem();
388 Index m = aa_grid.nelem();
389 Numeric stepsize_za = grid_stepsize[0];
390 Numeric stepsize_aa = grid_stepsize[1];
391 Vector res1(n);
392 ARTS_ASSERT(is_size(Integrand, n, m));
393
394 Numeric temp = 0.0;
395
396 for (Index i = 0; i < n; ++i) {
397 temp = Integrand(i, 0);
398 for (Index j = 1; j < m - 1; j++) {
399 temp += Integrand(i, j) * 2;
400 }
401 temp += Integrand(i, m - 1);
402 temp *= 0.5 * DEG2RAD * stepsize_aa * sin(za_grid[i] * DEG2RAD);
403 res1[i] = temp;
404 }
405
406 res = res1[0];
407 for (Index i = 1; i < n - 1; i++) {
408 res += res1[i] * 2;
409 }
410 res += res1[n - 1];
411 res *= 0.5 * DEG2RAD * stepsize_za;
412 } else {
413 res = AngIntegrate_trapezoid(Integrand, za_grid, aa_grid);
414 }
415
416 return res;
417}
418
420
435 ConstVectorView za_grid) {
436 Index n = za_grid.nelem();
437 ARTS_ASSERT(is_size(Integrand, n));
438
439 Numeric res = 0.0;
440 for (Index i = 0; i < n - 1; ++i) {
441 // in this place 0.5 * 2 * PI is calculated:
442 res += PI * DEG2RAD *
443 (Integrand[i] * sin(za_grid[i] * DEG2RAD) +
444 Integrand[i + 1] * sin(za_grid[i + 1] * DEG2RAD)) *
445 (za_grid[i + 1] - za_grid[i]);
446 }
447
448 return res;
449}
450
452
465 if (x < 0)
466 return -1.0;
467 else if (x == 0)
468 return 0.0;
469 else
470 return 1.0;
471}
472
493 const Vector& x,
494 const Numeric& n0,
495 const Numeric& mu,
496 const Numeric& la,
497 const Numeric& ga) {
498 const Index nx = x.nelem();
499
500 ARTS_ASSERT(psd.nelem() == nx);
501
502 if (ga == 1) {
503 if (mu == 0) {
504 // Exponential distribution
505 for (Index ix = 0; ix < nx; ix++) {
506 const Numeric eterm = exp(-la * x[ix]);
507 psd[ix] = n0 * eterm;
508 }
509 } else {
510 ARTS_USER_ERROR_IF (mu > 10,
511 "Given mu is ", mu, '\n',
512 "Seems unreasonable. Have you mixed up the inputs?")
513 // Gamma distribution
514 for (Index ix = 0; ix < nx; ix++) {
515 const Numeric eterm = exp(-la * x[ix]);
516 const Numeric xterm = pow(x[ix], mu);
517 psd[ix] = n0 * xterm * eterm;
518 psd[ix] = n0 * pow(x[ix], mu) * exp(-la * x[ix]);
519 }
520 }
521 } else {
522 // Complete MGD
523 ARTS_USER_ERROR_IF (mu > 10,
524 "Given mu is ", mu, '\n',
525 "Seems unreasonable. Have you mixed up the inputs?")
526 ARTS_USER_ERROR_IF (ga > 10,
527 "Given gamma is ", ga, '\n',
528 "Seems unreasonable. Have you mixed up the inputs?")
529 for (Index ix = 0; ix < nx; ix++) {
530 const Numeric pterm = pow(x[ix], ga);
531 const Numeric eterm = exp(-la * pterm);
532 const Numeric xterm = pow(x[ix], mu);
533 psd[ix] = n0 * xterm * eterm;
534 }
535 }
536}
537
563 MatrixView jac_data,
564 const Vector& x,
565 const Numeric& n0,
566 const Numeric& mu,
567 const Numeric& la,
568 const Numeric& ga,
569 const bool& do_n0_jac,
570 const bool& do_mu_jac,
571 const bool& do_la_jac,
572 const bool& do_ga_jac) {
573 const Index nx = x.nelem();
574
575 ARTS_ASSERT(psd.nelem() == nx);
576 ARTS_ASSERT(jac_data.nrows() == 4);
577 ARTS_ASSERT(jac_data.ncols() == nx);
578
579 if (ga == 1 && !do_ga_jac) {
580 if (mu == 0 && !do_mu_jac) {
581 // Exponential distribution
582 for (Index ix = 0; ix < nx; ix++) {
583 const Numeric eterm = exp(-la * x[ix]);
584 psd[ix] = n0 * eterm;
585 if (do_n0_jac) {
586 jac_data(0, ix) = eterm;
587 }
588 if (do_la_jac) {
589 jac_data(2, ix) = -x[ix] * psd[ix];
590 }
591 }
592 } else {
593 ARTS_USER_ERROR_IF (mu > 10,
594 "Given mu is ", mu, '\n',
595 "Seems unreasonable. Have you mixed up the inputs?")
596 // Gamma distribution
597 for (Index ix = 0; ix < nx; ix++) {
598 const Numeric eterm = exp(-la * x[ix]);
599 const Numeric xterm = pow(x[ix], mu);
600 psd[ix] = n0 * xterm * eterm;
601 if (do_n0_jac) {
602 jac_data(0, ix) = xterm * eterm;
603 }
604 if (do_mu_jac) {
605 jac_data(1, ix) = log(x[ix]) * psd[ix];
606 }
607 if (do_la_jac) {
608 jac_data(2, ix) = -x[ix] * psd[ix];
609 }
610 psd[ix] = n0 * pow(x[ix], mu) * exp(-la * x[ix]);
611 }
612 }
613 } else {
614 // Complete MGD
615 ARTS_USER_ERROR_IF (mu > 10,
616 "Given mu is ", mu, '\n',
617 "Seems unreasonable. Have you mixed up the inputs?")
618 ARTS_USER_ERROR_IF (ga > 10,
619 "Given gamma is ", ga, '\n',
620 "Seems unreasonable. Have you mixed up the inputs?")
621 for (Index ix = 0; ix < nx; ix++) {
622 const Numeric pterm = pow(x[ix], ga);
623 const Numeric eterm = exp(-la * pterm);
624 const Numeric xterm = pow(x[ix], mu);
625 psd[ix] = n0 * xterm * eterm;
626 if (do_n0_jac) {
627 jac_data(0, ix) = xterm * eterm;
628 }
629 if (do_mu_jac) {
630 jac_data(1, ix) = log(x[ix]) * psd[ix];
631 }
632 if (do_la_jac) {
633 jac_data(2, ix) = -pterm * psd[ix];
634 }
635 if (do_ga_jac) {
636 jac_data(3, ix) = -la * pterm * log(x[ix]) * psd[ix];
637 }
638 }
639 }
640}
641
643 MatrixView jac_data,
644 const Vector& x,
645 const Numeric& alpha,
646 const Numeric& beta) {
647 Numeric f_c = tgamma(4.0) / 256.0;
648 f_c *= pow(tgamma((alpha + 5.0) / beta), 4 + alpha);
649 f_c /= pow(tgamma((alpha + 4.0) / beta), 5 + alpha);
650
651 Numeric f_d = tgamma((alpha + 5.0) / beta);
652 f_d /= tgamma((alpha + 4.0) / beta);
653
654 for (Index i = 0; i < x.nelem(); ++i) {
655 Numeric xi = x[i];
656 psd[i] = beta * f_c * pow(xi, alpha) * exp(-pow(f_d * xi, beta));
657 jac_data(0, i) =
658 psd[i] * (alpha / xi - beta * f_d * pow(f_d * xi, beta - 1.0));
659 }
660}
661
663
677 Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma) {
678 Numeric dN;
679
680 if (x > 0. && N0 > 0. && Lambda > 0. && (mu + 1) / gamma > 0.) {
681 //Distribution function
682 dN = N0 * pow(x, mu) * exp(-Lambda * pow(x, gamma));
683
684 return dN;
685 } else {
687 "At least one argument is zero or negative.\n"
688 "Modified gamma distribution can not be calculated.\n"
689 "x = ", x, "\n"
690 "N0 = ", N0, "\n"
691 "lambda = ", Lambda, "\n"
692 "mu = ", mu, "\n"
693 "gamma = ", gamma, "\n")
694 }
695}
696
698
708void unitl(Vector& x) {
709 ARTS_ASSERT(x.nelem() > 0);
710
711 const Numeric l = sqrt(x * x);
712 for (Index i = 0; i < x.nelem(); i++) x[i] /= l;
713}
714
716
729 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols());
730
731 Index i = 0;
732
733 for (Index c = 0; c < X.ncols(); c++) {
734 for (Index r = 0; r < X.nrows(); r++) {
735 x[i] = X(r, c);
736 i += 1;
737 }
738 }
739}
740
742
755 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols() * X.npages());
756
757 Index i = 0;
758
759 for (Index c = 0; c < X.ncols(); c++) {
760 for (Index r = 0; r < X.nrows(); r++) {
761 for (Index p = 0; p < X.npages(); p++) {
762 x[i] = X(p, r, c);
763 i += 1;
764 }
765 }
766 }
767}
768
770
783 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols() * X.npages());
784
785 Index i = 0;
786
787 for (Index c = 0; c < X.ncols(); c++) {
788 for (Index r = 0; r < X.nrows(); r++) {
789 for (Index p = 0; p < X.npages(); p++) {
790 X(p, r, c) = x[i];
791 i += 1;
792 }
793 }
794 }
795}
796
798
811 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols());
812
813 Index i = 0;
814
815 for (Index c = 0; c < X.ncols(); c++) {
816 for (Index r = 0; r < X.nrows(); r++) {
817 X(r, c) = x[i];
818 i += 1;
819 }
820 }
821}
822
824
835 Index N = nph * 2;
836
837 //directions in total
838 nlinspace(x, -1, 1, N);
839
840 //allocate
841 w.resize(x.nelem());
842
843 // calculate weights
844 w[0] = (x[1] - x[0]) / 2.;
845
846 for (Index i = 1; i < nph * 2 - 1; i++) {
847 w[i] = (x[i + 1] - x[i - 1]) / 2.;
848 }
849 w[x.nelem() - 1] = (x[x.nelem() - 1] - x[x.nelem() - 2]) / 2.;
850}
851
853
854 ARTS_USER_ERROR_IF(x.nelem() <1, "Grid needs at least 2 points." );
855
856 Index N = x.nelem();
857
858 w.resize(N);
859
860 w[0] = (x[1] - x[0])/2;
861 if (N>2){
862
863 for (Index i = 1; i < N-1; i++ ){
864 w[i] = (x[i+1] - x[i-1])/2;
865 }
866
867 }
868 w[N-1] = (x[N-1] - x[N-2])/2;
869}
This file contains the definition of Array.
Constants of physical expressions as constexpr.
Common ARTS conversions.
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Matrix.
Definition: matpackI.h:1065
Index nrows() const noexcept
Definition: matpackI.h:1079
Index ncols() const noexcept
Definition: matpackI.h:1080
A constant view of a Tensor3.
Definition: matpackIII.h:133
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:145
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:148
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:151
A constant view of a Vector.
Definition: matpackI.h:521
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:547
The MatrixView class.
Definition: matpackI.h:1188
The Tensor3View class.
Definition: matpackIII.h:251
The VectorView class.
Definition: matpackI.h:674
The Vector class.
Definition: matpackI.h:910
void resize(Index n)
Resize function.
Definition: matpackI.cc:390
#define ARTS_ASSERT(condition,...)
Definition: debug.h:102
#define ARTS_USER_ERROR(...)
Definition: debug.h:169
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:153
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:82
bool is_multiple(const Index &x, const Index &y)
Checks if an integer is a multiple of another integer.
Definition: logic.cc:68
void linspace(Vector &x, const Numeric start, const Numeric stop, const Numeric step)
linspace
Definition: math_funcs.cc:200
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:227
void cumsum(VectorView csum, ConstVectorView x)
cumsum
Definition: math_funcs.cc:316
void unitl(Vector &x)
unitl
Definition: math_funcs.cc:708
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:65
Definition: math_funcs.cc:46
void nlogspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlogspace
Definition: math_funcs.cc:263
void reshape(Tensor3View X, ConstVectorView x)
reshape
Definition: math_funcs.cc:782
void mgd(VectorView psd, const Vector &x, const Numeric &n0, const Numeric &mu, const Numeric &la, const Numeric &ga)
Definition: math_funcs.cc:492
Numeric LagrangeInterpol4(ConstVectorView x, ConstVectorView y, const Numeric a)
Lagrange Interpolation (internal function).
Definition: math_funcs.cc:115
void delanoe_shape_with_derivative(VectorView psd, MatrixView jac_data, const Vector &x, const Numeric &alpha, const Numeric &beta)
! Shape functions for normalized PSD.
Definition: math_funcs.cc:642
void calculate_weights_linear(Vector &x, Vector &w, const Index nph)
calculate_weights_linear
Definition: math_funcs.cc:834
Numeric trapz(ConstVectorView x, ConstVectorView y)
trapz
Definition: math_funcs.cc:293
Numeric mod_gamma_dist(Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma)
Generalized Modified Gamma Distribution.
Definition: math_funcs.cc:676
void flat(VectorView x, ConstMatrixView X)
flat
Definition: math_funcs.cc:728
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:161
Index integer_div(const Index &x, const Index &y)
integer_div
Definition: math_funcs.cc:89
Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid)
AngIntegrate_trapezoid.
Definition: math_funcs.cc:337
constexpr Numeric PI
Definition: math_funcs.cc:47
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:464
void mgd_with_derivatives(VectorView psd, MatrixView jac_data, const Vector &x, const Numeric &n0, const Numeric &mu, const Numeric &la, const Numeric &ga, const bool &do_n0_jac, const bool &do_mu_jac, const bool &do_la_jac, const bool &do_ga_jac)
Definition: math_funcs.cc:562
void calculate_int_weights_arbitrary_grid(Vector &w, const Vector &x)
Calculates trapezoidal integration weights for arbitray grid.
Definition: math_funcs.cc:852
Numeric AngIntegrate_trapezoid_opti(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView grid_stepsize)
AngIntegrate_trapezoid_opti.
Definition: math_funcs.cc:381
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
This file contains the definition of String, the ARTS string class.
constexpr Numeric pi
The following mathematical constants are generated in python Decimal package by the code: