ARTS 2.5.9 (git: 825fa5f2)
math_funcs.cc
Go to the documentation of this file.
1/* Copyright (C) 2002-2012
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
6 under the terms of the GNU General Public License as published by the
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
316 ConstVectorView za_grid,
317 ConstVectorView aa_grid) {
318 Index n = za_grid.nelem();
319 Index m = aa_grid.nelem();
320 Vector res1(n);
321 ARTS_ASSERT(is_size(Integrand, n, m));
322
323 for (Index i = 0; i < n; ++i) {
324 res1[i] = 0.0;
325
326 for (Index j = 0; j < m - 1; ++j) {
327 res1[i] += 0.5 * DEG2RAD * (Integrand(i, j) + Integrand(i, j + 1)) *
328 (aa_grid[j + 1] - aa_grid[j]) * sin(za_grid[i] * DEG2RAD);
329 }
330 }
331 Numeric res = 0.0;
332 for (Index i = 0; i < n - 1; ++i) {
333 res +=
334 0.5 * DEG2RAD * (res1[i] + res1[i + 1]) * (za_grid[i + 1] - za_grid[i]);
335 }
336
337 return res;
338}
339
341
360 ConstVectorView za_grid,
361 ConstVectorView aa_grid,
362 ConstVectorView grid_stepsize) {
363 Numeric res = 0;
364 if ((grid_stepsize[0] > 0) && (grid_stepsize[1] > 0)) {
365 Index n = za_grid.nelem();
366 Index m = aa_grid.nelem();
367 Numeric stepsize_za = grid_stepsize[0];
368 Numeric stepsize_aa = grid_stepsize[1];
369 Vector res1(n);
370 ARTS_ASSERT(is_size(Integrand, n, m));
371
372 Numeric temp = 0.0;
373
374 for (Index i = 0; i < n; ++i) {
375 temp = Integrand(i, 0);
376 for (Index j = 1; j < m - 1; j++) {
377 temp += Integrand(i, j) * 2;
378 }
379 temp += Integrand(i, m - 1);
380 temp *= 0.5 * DEG2RAD * stepsize_aa * sin(za_grid[i] * DEG2RAD);
381 res1[i] = temp;
382 }
383
384 res = res1[0];
385 for (Index i = 1; i < n - 1; i++) {
386 res += res1[i] * 2;
387 }
388 res += res1[n - 1];
389 res *= 0.5 * DEG2RAD * stepsize_za;
390 } else {
391 res = AngIntegrate_trapezoid(Integrand, za_grid, aa_grid);
392 }
393
394 return res;
395}
396
398
413 ConstVectorView za_grid) {
414 Index n = za_grid.nelem();
415 ARTS_ASSERT(is_size(Integrand, n));
416
417 Numeric res = 0.0;
418 for (Index i = 0; i < n - 1; ++i) {
419 // in this place 0.5 * 2 * PI is calculated:
420 res += PI * DEG2RAD *
421 (Integrand[i] * sin(za_grid[i] * DEG2RAD) +
422 Integrand[i + 1] * sin(za_grid[i + 1] * DEG2RAD)) *
423 (za_grid[i + 1] - za_grid[i]);
424 }
425
426 return res;
427}
428
430
443 if (x < 0)
444 return -1.0;
445 else if (x == 0)
446 return 0.0;
447 else
448 return 1.0;
449}
450
471 const Vector& x,
472 const Numeric& n0,
473 const Numeric& mu,
474 const Numeric& la,
475 const Numeric& ga) {
476 const Index nx = x.nelem();
477
478 ARTS_ASSERT(psd.nelem() == nx);
479
480 if (ga == 1) {
481 if (mu == 0) {
482 // Exponential distribution
483 for (Index ix = 0; ix < nx; ix++) {
484 const Numeric eterm = exp(-la * x[ix]);
485 psd[ix] = n0 * eterm;
486 }
487 } else {
488 ARTS_USER_ERROR_IF (mu > 10,
489 "Given mu is ", mu, '\n',
490 "Seems unreasonable. Have you mixed up the inputs?")
491 // Gamma distribution
492 for (Index ix = 0; ix < nx; ix++) {
493 const Numeric eterm = exp(-la * x[ix]);
494 const Numeric xterm = pow(x[ix], mu);
495 psd[ix] = n0 * xterm * eterm;
496 psd[ix] = n0 * pow(x[ix], mu) * exp(-la * x[ix]);
497 }
498 }
499 } else {
500 // Complete MGD
501 ARTS_USER_ERROR_IF (mu > 10,
502 "Given mu is ", mu, '\n',
503 "Seems unreasonable. Have you mixed up the inputs?")
504 ARTS_USER_ERROR_IF (ga > 10,
505 "Given gamma is ", ga, '\n',
506 "Seems unreasonable. Have you mixed up the inputs?")
507 for (Index ix = 0; ix < nx; ix++) {
508 const Numeric pterm = pow(x[ix], ga);
509 const Numeric eterm = exp(-la * pterm);
510 const Numeric xterm = pow(x[ix], mu);
511 psd[ix] = n0 * xterm * eterm;
512 }
513 }
514}
515
541 MatrixView jac_data,
542 const Vector& x,
543 const Numeric& n0,
544 const Numeric& mu,
545 const Numeric& la,
546 const Numeric& ga,
547 const bool& do_n0_jac,
548 const bool& do_mu_jac,
549 const bool& do_la_jac,
550 const bool& do_ga_jac) {
551 const Index nx = x.nelem();
552
553 ARTS_ASSERT(psd.nelem() == nx);
554 ARTS_ASSERT(jac_data.nrows() == 4);
555 ARTS_ASSERT(jac_data.ncols() == nx);
556
557 if (ga == 1 && !do_ga_jac) {
558 if (mu == 0 && !do_mu_jac) {
559 // Exponential distribution
560 for (Index ix = 0; ix < nx; ix++) {
561 const Numeric eterm = exp(-la * x[ix]);
562 psd[ix] = n0 * eterm;
563 if (do_n0_jac) {
564 jac_data(0, ix) = eterm;
565 }
566 if (do_la_jac) {
567 jac_data(2, ix) = -x[ix] * psd[ix];
568 }
569 }
570 } else {
571 ARTS_USER_ERROR_IF (mu > 10,
572 "Given mu is ", mu, '\n',
573 "Seems unreasonable. Have you mixed up the inputs?")
574 // Gamma distribution
575 for (Index ix = 0; ix < nx; ix++) {
576 const Numeric eterm = exp(-la * x[ix]);
577 const Numeric xterm = pow(x[ix], mu);
578 psd[ix] = n0 * xterm * eterm;
579 if (do_n0_jac) {
580 jac_data(0, ix) = xterm * eterm;
581 }
582 if (do_mu_jac) {
583 jac_data(1, ix) = log(x[ix]) * psd[ix];
584 }
585 if (do_la_jac) {
586 jac_data(2, ix) = -x[ix] * psd[ix];
587 }
588 psd[ix] = n0 * pow(x[ix], mu) * exp(-la * x[ix]);
589 }
590 }
591 } else {
592 // Complete MGD
593 ARTS_USER_ERROR_IF (mu > 10,
594 "Given mu is ", mu, '\n',
595 "Seems unreasonable. Have you mixed up the inputs?")
596 ARTS_USER_ERROR_IF (ga > 10,
597 "Given gamma is ", ga, '\n',
598 "Seems unreasonable. Have you mixed up the inputs?")
599 for (Index ix = 0; ix < nx; ix++) {
600 const Numeric pterm = pow(x[ix], ga);
601 const Numeric eterm = exp(-la * pterm);
602 const Numeric xterm = pow(x[ix], mu);
603 psd[ix] = n0 * xterm * eterm;
604 if (do_n0_jac) {
605 jac_data(0, ix) = xterm * eterm;
606 }
607 if (do_mu_jac) {
608 jac_data(1, ix) = log(x[ix]) * psd[ix];
609 }
610 if (do_la_jac) {
611 jac_data(2, ix) = -pterm * psd[ix];
612 }
613 if (do_ga_jac) {
614 jac_data(3, ix) = -la * pterm * log(x[ix]) * psd[ix];
615 }
616 }
617 }
618}
619
621 MatrixView jac_data,
622 const Vector& x,
623 const Numeric& alpha,
624 const Numeric& beta) {
625 Numeric f_c = tgamma(4.0) / 256.0;
626 f_c *= pow(tgamma((alpha + 5.0) / beta), 4 + alpha);
627 f_c /= pow(tgamma((alpha + 4.0) / beta), 5 + alpha);
628
629 Numeric f_d = tgamma((alpha + 5.0) / beta);
630 f_d /= tgamma((alpha + 4.0) / beta);
631
632 for (Index i = 0; i < x.nelem(); ++i) {
633 Numeric xi = x[i];
634 psd[i] = beta * f_c * pow(xi, alpha) * exp(-pow(f_d * xi, beta));
635 jac_data(0, i) =
636 psd[i] * (alpha / xi - beta * f_d * pow(f_d * xi, beta - 1.0));
637 }
638}
639
641
655 Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma) {
656 Numeric dN;
657
658 if (x > 0. && N0 > 0. && Lambda > 0. && (mu + 1) / gamma > 0.) {
659 //Distribution function
660 dN = N0 * pow(x, mu) * exp(-Lambda * pow(x, gamma));
661
662 return dN;
663 } else {
665 "At least one argument is zero or negative.\n"
666 "Modified gamma distribution can not be calculated.\n"
667 "x = ", x, "\n"
668 "N0 = ", N0, "\n"
669 "lambda = ", Lambda, "\n"
670 "mu = ", mu, "\n"
671 "gamma = ", gamma, "\n")
672 }
673}
674
676
686void unitl(Vector& x) {
687 ARTS_ASSERT(x.nelem() > 0);
688
689 const Numeric l = sqrt(x * x);
690 for (Index i = 0; i < x.nelem(); i++) x[i] /= l;
691}
692
694
707 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols());
708
709 Index i = 0;
710
711 for (Index c = 0; c < X.ncols(); c++) {
712 for (Index r = 0; r < X.nrows(); r++) {
713 x[i] = X(r, c);
714 i += 1;
715 }
716 }
717}
718
720
733 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols() * X.npages());
734
735 Index i = 0;
736
737 for (Index c = 0; c < X.ncols(); c++) {
738 for (Index r = 0; r < X.nrows(); r++) {
739 for (Index p = 0; p < X.npages(); p++) {
740 x[i] = X(p, r, c);
741 i += 1;
742 }
743 }
744 }
745}
746
748
761 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols() * X.npages());
762
763 Index i = 0;
764
765 for (Index c = 0; c < X.ncols(); c++) {
766 for (Index r = 0; r < X.nrows(); r++) {
767 for (Index p = 0; p < X.npages(); p++) {
768 X(p, r, c) = x[i];
769 i += 1;
770 }
771 }
772 }
773}
774
776
789 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols());
790
791 Index i = 0;
792
793 for (Index c = 0; c < X.ncols(); c++) {
794 for (Index r = 0; r < X.nrows(); r++) {
795 X(r, c) = x[i];
796 i += 1;
797 }
798 }
799}
800
802
813 Index N = nph * 2;
814
815 //directions in total
816 nlinspace(x, -1, 1, N);
817
818 //allocate
819 w.resize(x.nelem());
820
821 // calculate weights
822 w[0] = (x[1] - x[0]) / 2.;
823
824 for (Index i = 1; i < nph * 2 - 1; i++) {
825 w[i] = (x[i + 1] - x[i - 1]) / 2.;
826 }
827 w[x.nelem() - 1] = (x[x.nelem() - 1] - x[x.nelem() - 2]) / 2.;
828}
829
831
832 ARTS_USER_ERROR_IF(x.nelem() <1, "Grid needs at least 2 points." );
833
834 Index N = x.nelem();
835
836 w.resize(N);
837
838 w[0] = (x[1] - x[0])/2;
839 if (N>2){
840
841 for (Index i = 1; i < N-1; i++ ){
842 w[i] = (x[i+1] - x[i-1])/2;
843 }
844
845 }
846 w[N-1] = (x[N-1] - x[N-2])/2;
847}
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:82
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
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
Header file for logic.cc.
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 unitl(Vector &x)
unitl
Definition: math_funcs.cc:686
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:65
constexpr Numeric DEG2RAD
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:760
void mgd(VectorView psd, const Vector &x, const Numeric &n0, const Numeric &mu, const Numeric &la, const Numeric &ga)
Definition: math_funcs.cc:470
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:620
void calculate_weights_linear(Vector &x, Vector &w, const Index nph)
calculate_weights_linear
Definition: math_funcs.cc:812
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:654
void flat(VectorView x, ConstMatrixView X)
flat
Definition: math_funcs.cc:706
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:315
constexpr Numeric PI
Definition: math_funcs.cc:47
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:442
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:540
void calculate_int_weights_arbitrary_grid(Vector &w, const Vector &x)
Calculates trapezoidal integration weights for arbitray grid.
Definition: math_funcs.cc:830
Numeric AngIntegrate_trapezoid_opti(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView grid_stepsize)
AngIntegrate_trapezoid_opti.
Definition: math_funcs.cc:359
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:
constexpr auto deg2rad(auto x) noexcept
Converts degrees to radians.
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:713
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:705
#define N
Definition: rng.cc:164
#define w
#define a
#define c
#define b