ARTS 2.5.0 (git: 9ee3ac6c)
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 "logic.h"
42#include "mystring.h"
43
44extern const Numeric DEG2RAD;
45extern const Numeric PI;
46
47/*****************************************************************************
48 *** The functions (in alphabetical order)
49 *****************************************************************************/
50
52
63Numeric fac(const Index n) {
65
66 if (n == 0) return (1.0);
67
68 sum = 1.0;
69 for (Index i = 1; i <= n; i++) sum *= Numeric(i);
70
71 return (sum);
72}
73
75
87Index integer_div(const Index& x, const Index& y) {
89 return x / y;
90}
91
93
115 const Numeric a) {
116 // lowermost grid spacing on x-axis
117 const Numeric Dlimit = 1.00000e-15;
118
119 // Check that dimensions of x and y vector agree
120 const Index n_x = x.nelem();
121 const Index n_y = y.nelem();
122 ARTS_USER_ERROR_IF ((n_x != 4) || (n_y != 4),
123 "The vectors x and y must all have the same length of 4 elements!\n"
124 "Actual lengths:\n"
125 "x:", n_x, ", "
126 "y:", n_y, ".")
127
128 // assure that x1 =< a < x2 holds
129 ARTS_USER_ERROR_IF ((a < x[1]) || (a > x[2]),
130 "LagrangeInterpol4: the relation x[1] =< a < x[2] is not satisfied. "
131 "No interpolation can be calculated.\n")
132
133 // calculate the Lagrange polynomial coefficients for a polynomial of the order of 3
134 Numeric b[4];
135 for (Index i = 0; i < 4; ++i) {
136 b[i] = 1.000e0;
137 for (Index k = 0; k < 4; ++k) {
138 if ((k != i) && (fabs(x[i] - x[k]) > Dlimit))
139 b[i] = b[i] * ((a - x[k]) / (x[i] - x[k]));
140 };
141 };
142
143 Numeric ya = 0.000e0;
144 for (Index i = 0; i < n_x; ++i) ya = ya + b[i] * y[i];
145
146 return ya;
147}
148
150
160 ARTS_ASSERT(x.nelem() > 0);
161 return x[x.nelem() - 1];
162}
163
165
175 ARTS_ASSERT(x.nelem() > 0);
176 return x[x.nelem() - 1];
177}
178
180
199 const Numeric start,
200 const Numeric stop,
201 const Numeric step) {
202 Index n = (Index)floor((stop - start) / step) + 1;
203 if (n < 1) n = 1;
204 x.resize(n);
205 for (Index i = 0; i < n; i++) x[i] = start + (double)i * step;
206}
207
209
226 const Numeric start,
227 const Numeric stop,
228 const Index n) {
229 ARTS_ASSERT(1 < n); // Number of points must be greater 1.
230 x.resize(n);
231 Numeric step = (stop - start) / ((double)n - 1);
232 for (Index i = 0; i < n - 1; i++) x[i] = start + (double)i * step;
233 x[n - 1] = stop;
234}
236 const Numeric start,
237 const Numeric stop,
238 const Index n) {
239 Numeric step = (stop - start) / ((double)n - 1);
240 for (Index i = 0; i < n - 1; i++) x[i] = start + (double)i * step;
241 x[n - 1] = stop;
242}
243
245
262 const Numeric start,
263 const Numeric stop,
264 const Index n) {
265 // Number of points must be greater than 1:
266 ARTS_ASSERT(1 < n);
267 // Only positive numbers are allowed for start and stop:
268 ARTS_ASSERT(0 < start);
269 ARTS_ASSERT(0 < stop);
270
271 x.resize(n);
272 Numeric a = log(start);
273 Numeric step = (log(stop) - a) / ((double)n - 1);
274 x[0] = start;
275 for (Index i = 1; i < n - 1; i++) x[i] = exp(a + (double)i * step);
276 x[n - 1] = stop;
277}
278
280
291 ConstVectorView za_grid,
292 ConstVectorView aa_grid) {
293 Index n = za_grid.nelem();
294 Index m = aa_grid.nelem();
295 Vector res1(n);
296 ARTS_ASSERT(is_size(Integrand, n, m));
297
298 for (Index i = 0; i < n; ++i) {
299 res1[i] = 0.0;
300
301 for (Index j = 0; j < m - 1; ++j) {
302 res1[i] += 0.5 * DEG2RAD * (Integrand(i, j) + Integrand(i, j + 1)) *
303 (aa_grid[j + 1] - aa_grid[j]) * sin(za_grid[i] * DEG2RAD);
304 }
305 }
306 Numeric res = 0.0;
307 for (Index i = 0; i < n - 1; ++i) {
308 res +=
309 0.5 * DEG2RAD * (res1[i] + res1[i + 1]) * (za_grid[i + 1] - za_grid[i]);
310 }
311
312 return res;
313}
314
316
335 ConstVectorView za_grid,
336 ConstVectorView aa_grid,
337 ConstVectorView grid_stepsize) {
338 Numeric res = 0;
339 if ((grid_stepsize[0] > 0) && (grid_stepsize[1] > 0)) {
340 Index n = za_grid.nelem();
341 Index m = aa_grid.nelem();
342 Numeric stepsize_za = grid_stepsize[0];
343 Numeric stepsize_aa = grid_stepsize[1];
344 Vector res1(n);
345 ARTS_ASSERT(is_size(Integrand, n, m));
346
347 Numeric temp = 0.0;
348
349 for (Index i = 0; i < n; ++i) {
350 temp = Integrand(i, 0);
351 for (Index j = 1; j < m - 1; j++) {
352 temp += Integrand(i, j) * 2;
353 }
354 temp += Integrand(i, m - 1);
355 temp *= 0.5 * DEG2RAD * stepsize_aa * sin(za_grid[i] * DEG2RAD);
356 res1[i] = temp;
357 }
358
359 res = res1[0];
360 for (Index i = 1; i < n - 1; i++) {
361 res += res1[i] * 2;
362 }
363 res += res1[n - 1];
364 res *= 0.5 * DEG2RAD * stepsize_za;
365 } else {
366 res = AngIntegrate_trapezoid(Integrand, za_grid, aa_grid);
367 }
368
369 return res;
370}
371
373
388 ConstVectorView za_grid) {
389 Index n = za_grid.nelem();
390 ARTS_ASSERT(is_size(Integrand, n));
391
392 Numeric res = 0.0;
393 for (Index i = 0; i < n - 1; ++i) {
394 // in this place 0.5 * 2 * PI is calculated:
395 res += PI * DEG2RAD *
396 (Integrand[i] * sin(za_grid[i] * DEG2RAD) +
397 Integrand[i + 1] * sin(za_grid[i + 1] * DEG2RAD)) *
398 (za_grid[i + 1] - za_grid[i]);
399 }
400
401 return res;
402}
403
405
418 if (x < 0)
419 return -1.0;
420 else if (x == 0)
421 return 0.0;
422 else
423 return 1.0;
424}
425
446 const Vector& x,
447 const Numeric& n0,
448 const Numeric& mu,
449 const Numeric& la,
450 const Numeric& ga) {
451 const Index nx = x.nelem();
452
453 ARTS_ASSERT(psd.nelem() == nx);
454
455 if (ga == 1) {
456 if (mu == 0) {
457 // Exponential distribution
458 for (Index ix = 0; ix < nx; ix++) {
459 const Numeric eterm = exp(-la * x[ix]);
460 psd[ix] = n0 * eterm;
461 }
462 } else {
463 ARTS_USER_ERROR_IF (mu > 10,
464 "Given mu is ", mu, '\n',
465 "Seems unreasonable. Have you mixed up the inputs?")
466 // Gamma distribution
467 for (Index ix = 0; ix < nx; ix++) {
468 const Numeric eterm = exp(-la * x[ix]);
469 const Numeric xterm = pow(x[ix], mu);
470 psd[ix] = n0 * xterm * eterm;
471 psd[ix] = n0 * pow(x[ix], mu) * exp(-la * x[ix]);
472 }
473 }
474 } else {
475 // Complete MGD
476 ARTS_USER_ERROR_IF (mu > 10,
477 "Given mu is ", mu, '\n',
478 "Seems unreasonable. Have you mixed up the inputs?")
479 ARTS_USER_ERROR_IF (ga > 10,
480 "Given gamma is ", ga, '\n',
481 "Seems unreasonable. Have you mixed up the inputs?")
482 for (Index ix = 0; ix < nx; ix++) {
483 const Numeric pterm = pow(x[ix], ga);
484 const Numeric eterm = exp(-la * pterm);
485 const Numeric xterm = pow(x[ix], mu);
486 psd[ix] = n0 * xterm * eterm;
487 }
488 }
489}
490
516 MatrixView jac_data,
517 const Vector& x,
518 const Numeric& n0,
519 const Numeric& mu,
520 const Numeric& la,
521 const Numeric& ga,
522 const bool& do_n0_jac,
523 const bool& do_mu_jac,
524 const bool& do_la_jac,
525 const bool& do_ga_jac) {
526 const Index nx = x.nelem();
527
528 ARTS_ASSERT(psd.nelem() == nx);
529 ARTS_ASSERT(jac_data.nrows() == 4);
530 ARTS_ASSERT(jac_data.ncols() == nx);
531
532 if (ga == 1 && !do_ga_jac) {
533 if (mu == 0 && !do_mu_jac) {
534 // Exponential distribution
535 for (Index ix = 0; ix < nx; ix++) {
536 const Numeric eterm = exp(-la * x[ix]);
537 psd[ix] = n0 * eterm;
538 if (do_n0_jac) {
539 jac_data(0, ix) = eterm;
540 }
541 if (do_la_jac) {
542 jac_data(2, ix) = -x[ix] * psd[ix];
543 }
544 }
545 } else {
546 ARTS_USER_ERROR_IF (mu > 10,
547 "Given mu is ", mu, '\n',
548 "Seems unreasonable. Have you mixed up the inputs?")
549 // Gamma distribution
550 for (Index ix = 0; ix < nx; ix++) {
551 const Numeric eterm = exp(-la * x[ix]);
552 const Numeric xterm = pow(x[ix], mu);
553 psd[ix] = n0 * xterm * eterm;
554 if (do_n0_jac) {
555 jac_data(0, ix) = xterm * eterm;
556 }
557 if (do_mu_jac) {
558 jac_data(1, ix) = log(x[ix]) * psd[ix];
559 }
560 if (do_la_jac) {
561 jac_data(2, ix) = -x[ix] * psd[ix];
562 }
563 psd[ix] = n0 * pow(x[ix], mu) * exp(-la * x[ix]);
564 }
565 }
566 } else {
567 // Complete MGD
568 ARTS_USER_ERROR_IF (mu > 10,
569 "Given mu is ", mu, '\n',
570 "Seems unreasonable. Have you mixed up the inputs?")
571 ARTS_USER_ERROR_IF (ga > 10,
572 "Given gamma is ", ga, '\n',
573 "Seems unreasonable. Have you mixed up the inputs?")
574 for (Index ix = 0; ix < nx; ix++) {
575 const Numeric pterm = pow(x[ix], ga);
576 const Numeric eterm = exp(-la * pterm);
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) = -pterm * psd[ix];
587 }
588 if (do_ga_jac) {
589 jac_data(3, ix) = -la * pterm * log(x[ix]) * psd[ix];
590 }
591 }
592 }
593}
594
596 MatrixView jac_data,
597 const Vector& x,
598 const Numeric& alpha,
599 const Numeric& beta) {
600 Numeric f_c = tgamma(4.0) / 256.0;
601 f_c *= pow(tgamma((alpha + 5.0) / beta), 4 + alpha);
602 f_c /= pow(tgamma((alpha + 4.0) / beta), 5 + alpha);
603
604 Numeric f_d = tgamma((alpha + 5.0) / beta);
605 f_d /= tgamma((alpha + 4.0) / beta);
606
607 for (Index i = 0; i < x.nelem(); ++i) {
608 Numeric xi = x[i];
609 psd[i] = beta * f_c * pow(xi, alpha) * exp(-pow(f_d * xi, beta));
610 jac_data(0, i) =
611 psd[i] * (alpha / xi - beta * f_d * pow(f_d * xi, beta - 1.0));
612 }
613}
614
616
630 Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma) {
631 Numeric dN;
632
633 if (x > 0. && N0 > 0. && Lambda > 0. && (mu + 1) / gamma > 0.) {
634 //Distribution function
635 dN = N0 * pow(x, mu) * exp(-Lambda * pow(x, gamma));
636
637 return dN;
638 } else {
640 "At least one argument is zero or negative.\n"
641 "Modified gamma distribution can not be calculated.\n"
642 "x = ", x, "\n"
643 "N0 = ", N0, "\n"
644 "lambda = ", Lambda, "\n"
645 "mu = ", mu, "\n"
646 "gamma = ", gamma, "\n")
647 }
648}
649
651
661void unitl(Vector& x) {
662 ARTS_ASSERT(x.nelem() > 0);
663
664 const Numeric l = sqrt(x * x);
665 for (Index i = 0; i < x.nelem(); i++) x[i] /= l;
666}
667
669
682 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols());
683
684 Index i = 0;
685
686 for (Index c = 0; c < X.ncols(); c++) {
687 for (Index r = 0; r < X.nrows(); r++) {
688 x[i] = X(r, c);
689 i += 1;
690 }
691 }
692}
693
695
708 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols() * X.npages());
709
710 Index i = 0;
711
712 for (Index c = 0; c < X.ncols(); c++) {
713 for (Index r = 0; r < X.nrows(); r++) {
714 for (Index p = 0; p < X.npages(); p++) {
715 x[i] = X(p, r, c);
716 i += 1;
717 }
718 }
719 }
720}
721
723
736 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols() * X.npages());
737
738 Index i = 0;
739
740 for (Index c = 0; c < X.ncols(); c++) {
741 for (Index r = 0; r < X.nrows(); r++) {
742 for (Index p = 0; p < X.npages(); p++) {
743 X(p, r, c) = x[i];
744 i += 1;
745 }
746 }
747 }
748}
749
751
764 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols());
765
766 Index i = 0;
767
768 for (Index c = 0; c < X.ncols(); c++) {
769 for (Index r = 0; r < X.nrows(); r++) {
770 X(r, c) = x[i];
771 i += 1;
772 }
773 }
774}
775
777
788 Index N = nph * 2;
789
790 //directions in total
791 nlinspace(x, -1, 1, N);
792
793 //allocate
794 w.resize(x.nelem());
795
796 // calculate weights
797 w[0] = (x[1] - x[0]) / 2.;
798
799 for (Index i = 1; i < nph * 2 - 1; i++) {
800 w[i] = (x[i + 1] - x[i - 1]) / 2.;
801 }
802 w[x.nelem() - 1] = (x[x.nelem() - 1] - x[x.nelem() - 2]) / 2.;
803}
This file contains the definition of Array.
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Definition: array.h:195
A constant view of a Matrix.
Definition: matpackI.h:1014
Index nrows() const ARTS_NOEXCEPT
Returns the number of rows.
Definition: matpackI.cc:433
Index ncols() const ARTS_NOEXCEPT
Returns the number of columns.
Definition: matpackI.cc:436
A constant view of a Tensor3.
Definition: matpackIII.h:132
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
A constant view of a Vector.
Definition: matpackI.h:489
Index nelem() const ARTS_NOEXCEPT
Returns the number of elements.
Definition: matpackI.cc:48
The MatrixView class.
Definition: matpackI.h:1125
The Tensor3View class.
Definition: matpackIII.h:239
The VectorView class.
Definition: matpackI.h:626
The Vector class.
Definition: matpackI.h:876
void resize(Index n)
Resize function.
Definition: matpackI.cc:408
#define ARTS_ASSERT(condition,...)
Definition: debug.h:83
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define ARTS_USER_ERROR_IF(condition,...)
Definition: debug.h:134
#define beta
#define temp
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:79
bool is_multiple(const Index &x, const Index &y)
Checks if an integer is a multiple of another integer.
Definition: logic.cc:65
Header file for logic.cc.
void linspace(Vector &x, const Numeric start, const Numeric stop, const Numeric step)
linspace
Definition: math_funcs.cc:198
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:225
void unitl(Vector &x)
unitl
Definition: math_funcs.cc:661
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:63
void nlogspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlogspace
Definition: math_funcs.cc:261
void reshape(Tensor3View X, ConstVectorView x)
reshape
Definition: math_funcs.cc:735
void mgd(VectorView psd, const Vector &x, const Numeric &n0, const Numeric &mu, const Numeric &la, const Numeric &ga)
Definition: math_funcs.cc:445
Numeric LagrangeInterpol4(ConstVectorView x, ConstVectorView y, const Numeric a)
Lagrange Interpolation (internal function).
Definition: math_funcs.cc:113
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:595
void calculate_weights_linear(Vector &x, Vector &w, const Index nph)
calculate_weights_linear
Definition: math_funcs.cc:787
Numeric mod_gamma_dist(Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma)
Generalized Modified Gamma Distribution.
Definition: math_funcs.cc:629
const Numeric PI
void flat(VectorView x, ConstMatrixView X)
flat
Definition: math_funcs.cc:681
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:159
Index integer_div(const Index &x, const Index &y)
integer_div
Definition: math_funcs.cc:87
Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid)
AngIntegrate_trapezoid.
Definition: math_funcs.cc:290
const Numeric DEG2RAD
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:417
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:515
Numeric AngIntegrate_trapezoid_opti(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView grid_stepsize)
AngIntegrate_trapezoid_opti.
Definition: math_funcs.cc:334
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
This file contains the definition of String, the ARTS string class.
constexpr Numeric l(const Index p0, const Index n, const Numeric x, const SortedVectorType &xi, const Index j, const std::pair< Numeric, Numeric > cycle={ -180, 180}) noexcept
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:78
void sum(PropagationMatrix &pm, const ComplexVectorView &abs, const PolarizationVector &polvec, const bool do_phase) ARTS_NOEXCEPT
Sums the Zeeman components into a propagation matrix.
Definition: zeemandata.cc:375
#define w
#define a
#define c
#define b
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:747
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:739
#define N
Definition: rng.cc:164