ARTS 2.5.4 (git: 7d04b88e)
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
293{
294 const Index n = x.nelem();
295 ARTS_ASSERT(y.nelem() == n);
296 Numeric sum = 0.0;
297 for (Index i=1; i<n; ++i)
298 sum += (x[i]-x[i-1]) * (y[i]+y[i-1]);
299 return sum/2.0;
300}
301
303
314 ConstVectorView za_grid,
315 ConstVectorView aa_grid) {
316 Index n = za_grid.nelem();
317 Index m = aa_grid.nelem();
318 Vector res1(n);
319 ARTS_ASSERT(is_size(Integrand, n, m));
320
321 for (Index i = 0; i < n; ++i) {
322 res1[i] = 0.0;
323
324 for (Index j = 0; j < m - 1; ++j) {
325 res1[i] += 0.5 * DEG2RAD * (Integrand(i, j) + Integrand(i, j + 1)) *
326 (aa_grid[j + 1] - aa_grid[j]) * sin(za_grid[i] * DEG2RAD);
327 }
328 }
329 Numeric res = 0.0;
330 for (Index i = 0; i < n - 1; ++i) {
331 res +=
332 0.5 * DEG2RAD * (res1[i] + res1[i + 1]) * (za_grid[i + 1] - za_grid[i]);
333 }
334
335 return res;
336}
337
339
358 ConstVectorView za_grid,
359 ConstVectorView aa_grid,
360 ConstVectorView grid_stepsize) {
361 Numeric res = 0;
362 if ((grid_stepsize[0] > 0) && (grid_stepsize[1] > 0)) {
363 Index n = za_grid.nelem();
364 Index m = aa_grid.nelem();
365 Numeric stepsize_za = grid_stepsize[0];
366 Numeric stepsize_aa = grid_stepsize[1];
367 Vector res1(n);
368 ARTS_ASSERT(is_size(Integrand, n, m));
369
370 Numeric temp = 0.0;
371
372 for (Index i = 0; i < n; ++i) {
373 temp = Integrand(i, 0);
374 for (Index j = 1; j < m - 1; j++) {
375 temp += Integrand(i, j) * 2;
376 }
377 temp += Integrand(i, m - 1);
378 temp *= 0.5 * DEG2RAD * stepsize_aa * sin(za_grid[i] * DEG2RAD);
379 res1[i] = temp;
380 }
381
382 res = res1[0];
383 for (Index i = 1; i < n - 1; i++) {
384 res += res1[i] * 2;
385 }
386 res += res1[n - 1];
387 res *= 0.5 * DEG2RAD * stepsize_za;
388 } else {
389 res = AngIntegrate_trapezoid(Integrand, za_grid, aa_grid);
390 }
391
392 return res;
393}
394
396
411 ConstVectorView za_grid) {
412 Index n = za_grid.nelem();
413 ARTS_ASSERT(is_size(Integrand, n));
414
415 Numeric res = 0.0;
416 for (Index i = 0; i < n - 1; ++i) {
417 // in this place 0.5 * 2 * PI is calculated:
418 res += PI * DEG2RAD *
419 (Integrand[i] * sin(za_grid[i] * DEG2RAD) +
420 Integrand[i + 1] * sin(za_grid[i + 1] * DEG2RAD)) *
421 (za_grid[i + 1] - za_grid[i]);
422 }
423
424 return res;
425}
426
428
441 if (x < 0)
442 return -1.0;
443 else if (x == 0)
444 return 0.0;
445 else
446 return 1.0;
447}
448
469 const Vector& x,
470 const Numeric& n0,
471 const Numeric& mu,
472 const Numeric& la,
473 const Numeric& ga) {
474 const Index nx = x.nelem();
475
476 ARTS_ASSERT(psd.nelem() == nx);
477
478 if (ga == 1) {
479 if (mu == 0) {
480 // Exponential distribution
481 for (Index ix = 0; ix < nx; ix++) {
482 const Numeric eterm = exp(-la * x[ix]);
483 psd[ix] = n0 * eterm;
484 }
485 } else {
486 ARTS_USER_ERROR_IF (mu > 10,
487 "Given mu is ", mu, '\n',
488 "Seems unreasonable. Have you mixed up the inputs?")
489 // Gamma distribution
490 for (Index ix = 0; ix < nx; ix++) {
491 const Numeric eterm = exp(-la * x[ix]);
492 const Numeric xterm = pow(x[ix], mu);
493 psd[ix] = n0 * xterm * eterm;
494 psd[ix] = n0 * pow(x[ix], mu) * exp(-la * x[ix]);
495 }
496 }
497 } else {
498 // Complete MGD
499 ARTS_USER_ERROR_IF (mu > 10,
500 "Given mu is ", mu, '\n',
501 "Seems unreasonable. Have you mixed up the inputs?")
502 ARTS_USER_ERROR_IF (ga > 10,
503 "Given gamma is ", ga, '\n',
504 "Seems unreasonable. Have you mixed up the inputs?")
505 for (Index ix = 0; ix < nx; ix++) {
506 const Numeric pterm = pow(x[ix], ga);
507 const Numeric eterm = exp(-la * pterm);
508 const Numeric xterm = pow(x[ix], mu);
509 psd[ix] = n0 * xterm * eterm;
510 }
511 }
512}
513
539 MatrixView jac_data,
540 const Vector& x,
541 const Numeric& n0,
542 const Numeric& mu,
543 const Numeric& la,
544 const Numeric& ga,
545 const bool& do_n0_jac,
546 const bool& do_mu_jac,
547 const bool& do_la_jac,
548 const bool& do_ga_jac) {
549 const Index nx = x.nelem();
550
551 ARTS_ASSERT(psd.nelem() == nx);
552 ARTS_ASSERT(jac_data.nrows() == 4);
553 ARTS_ASSERT(jac_data.ncols() == nx);
554
555 if (ga == 1 && !do_ga_jac) {
556 if (mu == 0 && !do_mu_jac) {
557 // Exponential distribution
558 for (Index ix = 0; ix < nx; ix++) {
559 const Numeric eterm = exp(-la * x[ix]);
560 psd[ix] = n0 * eterm;
561 if (do_n0_jac) {
562 jac_data(0, ix) = eterm;
563 }
564 if (do_la_jac) {
565 jac_data(2, ix) = -x[ix] * psd[ix];
566 }
567 }
568 } else {
569 ARTS_USER_ERROR_IF (mu > 10,
570 "Given mu is ", mu, '\n',
571 "Seems unreasonable. Have you mixed up the inputs?")
572 // Gamma distribution
573 for (Index ix = 0; ix < nx; ix++) {
574 const Numeric eterm = exp(-la * x[ix]);
575 const Numeric xterm = pow(x[ix], mu);
576 psd[ix] = n0 * xterm * eterm;
577 if (do_n0_jac) {
578 jac_data(0, ix) = xterm * eterm;
579 }
580 if (do_mu_jac) {
581 jac_data(1, ix) = log(x[ix]) * psd[ix];
582 }
583 if (do_la_jac) {
584 jac_data(2, ix) = -x[ix] * psd[ix];
585 }
586 psd[ix] = n0 * pow(x[ix], mu) * exp(-la * x[ix]);
587 }
588 }
589 } else {
590 // Complete MGD
591 ARTS_USER_ERROR_IF (mu > 10,
592 "Given mu is ", mu, '\n',
593 "Seems unreasonable. Have you mixed up the inputs?")
594 ARTS_USER_ERROR_IF (ga > 10,
595 "Given gamma is ", ga, '\n',
596 "Seems unreasonable. Have you mixed up the inputs?")
597 for (Index ix = 0; ix < nx; ix++) {
598 const Numeric pterm = pow(x[ix], ga);
599 const Numeric eterm = exp(-la * pterm);
600 const Numeric xterm = pow(x[ix], mu);
601 psd[ix] = n0 * xterm * eterm;
602 if (do_n0_jac) {
603 jac_data(0, ix) = xterm * eterm;
604 }
605 if (do_mu_jac) {
606 jac_data(1, ix) = log(x[ix]) * psd[ix];
607 }
608 if (do_la_jac) {
609 jac_data(2, ix) = -pterm * psd[ix];
610 }
611 if (do_ga_jac) {
612 jac_data(3, ix) = -la * pterm * log(x[ix]) * psd[ix];
613 }
614 }
615 }
616}
617
619 MatrixView jac_data,
620 const Vector& x,
621 const Numeric& alpha,
622 const Numeric& beta) {
623 Numeric f_c = tgamma(4.0) / 256.0;
624 f_c *= pow(tgamma((alpha + 5.0) / beta), 4 + alpha);
625 f_c /= pow(tgamma((alpha + 4.0) / beta), 5 + alpha);
626
627 Numeric f_d = tgamma((alpha + 5.0) / beta);
628 f_d /= tgamma((alpha + 4.0) / beta);
629
630 for (Index i = 0; i < x.nelem(); ++i) {
631 Numeric xi = x[i];
632 psd[i] = beta * f_c * pow(xi, alpha) * exp(-pow(f_d * xi, beta));
633 jac_data(0, i) =
634 psd[i] * (alpha / xi - beta * f_d * pow(f_d * xi, beta - 1.0));
635 }
636}
637
639
653 Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma) {
654 Numeric dN;
655
656 if (x > 0. && N0 > 0. && Lambda > 0. && (mu + 1) / gamma > 0.) {
657 //Distribution function
658 dN = N0 * pow(x, mu) * exp(-Lambda * pow(x, gamma));
659
660 return dN;
661 } else {
663 "At least one argument is zero or negative.\n"
664 "Modified gamma distribution can not be calculated.\n"
665 "x = ", x, "\n"
666 "N0 = ", N0, "\n"
667 "lambda = ", Lambda, "\n"
668 "mu = ", mu, "\n"
669 "gamma = ", gamma, "\n")
670 }
671}
672
674
684void unitl(Vector& x) {
685 ARTS_ASSERT(x.nelem() > 0);
686
687 const Numeric l = sqrt(x * x);
688 for (Index i = 0; i < x.nelem(); i++) x[i] /= l;
689}
690
692
705 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols());
706
707 Index i = 0;
708
709 for (Index c = 0; c < X.ncols(); c++) {
710 for (Index r = 0; r < X.nrows(); r++) {
711 x[i] = X(r, c);
712 i += 1;
713 }
714 }
715}
716
718
731 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols() * X.npages());
732
733 Index i = 0;
734
735 for (Index c = 0; c < X.ncols(); c++) {
736 for (Index r = 0; r < X.nrows(); r++) {
737 for (Index p = 0; p < X.npages(); p++) {
738 x[i] = X(p, r, c);
739 i += 1;
740 }
741 }
742 }
743}
744
746
759 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols() * X.npages());
760
761 Index i = 0;
762
763 for (Index c = 0; c < X.ncols(); c++) {
764 for (Index r = 0; r < X.nrows(); r++) {
765 for (Index p = 0; p < X.npages(); p++) {
766 X(p, r, c) = x[i];
767 i += 1;
768 }
769 }
770 }
771}
772
774
787 ARTS_ASSERT(x.nelem() == X.nrows() * X.ncols());
788
789 Index i = 0;
790
791 for (Index c = 0; c < X.ncols(); c++) {
792 for (Index r = 0; r < X.nrows(); r++) {
793 X(r, c) = x[i];
794 i += 1;
795 }
796 }
797}
798
800
811 Index N = nph * 2;
812
813 //directions in total
814 nlinspace(x, -1, 1, N);
815
816 //allocate
817 w.resize(x.nelem());
818
819 // calculate weights
820 w[0] = (x[1] - x[0]) / 2.;
821
822 for (Index i = 1; i < nph * 2 - 1; i++) {
823 w[i] = (x[i + 1] - x[i - 1]) / 2.;
824 }
825 w[x.nelem() - 1] = (x[x.nelem() - 1] - x[x.nelem() - 2]) / 2.;
826}
This file contains the definition of Array.
Index nelem() const ARTS_NOEXCEPT
Definition: array.h:92
A constant view of a Matrix.
Definition: matpackI.h:1064
Index nrows() const noexcept
Definition: matpackI.h:1075
Index ncols() const noexcept
Definition: matpackI.h:1076
A constant view of a Tensor3.
Definition: matpackIII.h:130
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:140
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:143
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:146
A constant view of a Vector.
Definition: matpackI.h:530
Index nelem() const noexcept
Returns the number of elements.
Definition: matpackI.h:554
The MatrixView class.
Definition: matpackI.h:1182
The Tensor3View class.
Definition: matpackIII.h:246
The VectorView class.
Definition: matpackI.h:672
The Vector class.
Definition: matpackI.h:922
void resize(Index n)
Resize function.
Definition: matpackI.cc:411
#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: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: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:684
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:758
void mgd(VectorView psd, const Vector &x, const Numeric &n0, const Numeric &mu, const Numeric &la, const Numeric &ga)
Definition: math_funcs.cc:468
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:618
void calculate_weights_linear(Vector &x, Vector &w, const Index nph)
calculate_weights_linear
Definition: math_funcs.cc:810
Numeric trapz(ConstVectorView x, ConstVectorView y)
trapz
Definition: math_funcs.cc:291
Numeric mod_gamma_dist(Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma)
Generalized Modified Gamma Distribution.
Definition: math_funcs.cc:652
const Numeric PI
void flat(VectorView x, ConstMatrixView X)
flat
Definition: math_funcs.cc:704
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:313
const Numeric DEG2RAD
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:440
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:538
Numeric AngIntegrate_trapezoid_opti(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView grid_stepsize)
AngIntegrate_trapezoid_opti.
Definition: math_funcs.cc:357
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 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:80
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:431
#define w
#define a
#define c
#define b
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:736
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:728
#define N
Definition: rng.cc:164