ARTS  2.4.0(git:4fb77825)
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 
44 extern const Numeric DEG2RAD;
45 extern const Numeric PI;
46 
47 /*****************************************************************************
48  *** The functions (in alphabetical order)
49  *****************************************************************************/
50 
52 
63 Numeric fac(const Index n) {
64  Numeric sum;
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 
87 Index integer_div(const Index& x, const Index& y) {
88  assert(is_multiple(x, 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  if ((n_x != 4) || (n_y != 4)) {
123  ostringstream os;
124  os << "The vectors x and y must all have the same length of 4 elements!\n"
125  << "Actual lengths:\n"
126  << "x:" << n_x << ", "
127  << "y:" << n_y << ".";
128  throw runtime_error(os.str());
129  }
130 
131  // assure that x1 =< a < x2 holds
132  if ((a < x[1]) || (a > x[2])) {
133  ostringstream os;
134  os << "LagrangeInterpol4: the relation x[1] =< a < x[2] is not satisfied. "
135  << "No interpolation can be calculated.\n";
136  throw runtime_error(os.str());
137  };
138 
139  // calculate the Lagrange polynomial coefficients for a polynomial of the order of 3
140  Numeric b[4];
141  for (Index i = 0; i < 4; ++i) {
142  b[i] = 1.000e0;
143  for (Index k = 0; k < 4; ++k) {
144  if ((k != i) && (fabs(x[i] - x[k]) > Dlimit))
145  b[i] = b[i] * ((a - x[k]) / (x[i] - x[k]));
146  };
147  };
148 
149  Numeric ya = 0.000e0;
150  for (Index i = 0; i < n_x; ++i) ya = ya + b[i] * y[i];
151 
152  return ya;
153 }
154 
156 
166  assert(x.nelem() > 0);
167  return x[x.nelem() - 1];
168 }
169 
171 
181  assert(x.nelem() > 0);
182  return x[x.nelem() - 1];
183 }
184 
186 
205  const Numeric start,
206  const Numeric stop,
207  const Numeric step) {
208  Index n = (Index)floor((stop - start) / step) + 1;
209  if (n < 1) n = 1;
210  x.resize(n);
211  for (Index i = 0; i < n; i++) x[i] = start + (double)i * step;
212 }
213 
215 
232  const Numeric start,
233  const Numeric stop,
234  const Index n) {
235  assert(1 < n); // Number of points must be greater 1.
236  x.resize(n);
237  Numeric step = (stop - start) / ((double)n - 1);
238  for (Index i = 0; i < n - 1; i++) x[i] = start + (double)i * step;
239  x[n - 1] = stop;
240 }
242  const Numeric start,
243  const Numeric stop,
244  const Index n) {
245  Numeric step = (stop - start) / ((double)n - 1);
246  for (Index i = 0; i < n - 1; i++) x[i] = start + (double)i * step;
247  x[n - 1] = stop;
248 }
249 
251 
268  const Numeric start,
269  const Numeric stop,
270  const Index n) {
271  // Number of points must be greater than 1:
272  assert(1 < n);
273  // Only positive numbers are allowed for start and stop:
274  assert(0 < start);
275  assert(0 < stop);
276 
277  x.resize(n);
278  Numeric a = log(start);
279  Numeric step = (log(stop) - a) / ((double)n - 1);
280  x[0] = start;
281  for (Index i = 1; i < n - 1; i++) x[i] = exp(a + (double)i * step);
282  x[n - 1] = stop;
283 }
284 
286 
299  Index n = za_grid.nelem();
300  Index m = aa_grid.nelem();
301  Vector res1(n);
302  assert(is_size(Integrand, n, m));
303 
304  for (Index i = 0; i < n; ++i) {
305  res1[i] = 0.0;
306 
307  for (Index j = 0; j < m - 1; ++j) {
308  res1[i] += 0.5 * DEG2RAD * (Integrand(i, j) + Integrand(i, j + 1)) *
309  (aa_grid[j + 1] - aa_grid[j]) * sin(za_grid[i] * DEG2RAD);
310  }
311  }
312  Numeric res = 0.0;
313  for (Index i = 0; i < n - 1; ++i) {
314  res +=
315  0.5 * DEG2RAD * (res1[i] + res1[i + 1]) * (za_grid[i + 1] - za_grid[i]);
316  }
317 
318  return res;
319 }
320 
322 
343  ConstVectorView grid_stepsize) {
344  Numeric res = 0;
345  if ((grid_stepsize[0] > 0) && (grid_stepsize[1] > 0)) {
346  Index n = za_grid.nelem();
347  Index m = aa_grid.nelem();
348  Numeric stepsize_za = grid_stepsize[0];
349  Numeric stepsize_aa = grid_stepsize[1];
350  Vector res1(n);
351  assert(is_size(Integrand, n, m));
352 
353  Numeric temp = 0.0;
354 
355  for (Index i = 0; i < n; ++i) {
356  temp = Integrand(i, 0);
357  for (Index j = 1; j < m - 1; j++) {
358  temp += Integrand(i, j) * 2;
359  }
360  temp += Integrand(i, m - 1);
361  temp *= 0.5 * DEG2RAD * stepsize_aa * sin(za_grid[i] * DEG2RAD);
362  res1[i] = temp;
363  }
364 
365  res = res1[0];
366  for (Index i = 1; i < n - 1; i++) {
367  res += res1[i] * 2;
368  }
369  res += res1[n - 1];
370  res *= 0.5 * DEG2RAD * stepsize_za;
371  } else {
372  res = AngIntegrate_trapezoid(Integrand, za_grid, aa_grid);
373  }
374 
375  return res;
376 }
377 
379 
395  Index n = za_grid.nelem();
396  assert(is_size(Integrand, n));
397 
398  Numeric res = 0.0;
399  for (Index i = 0; i < n - 1; ++i) {
400  // in this place 0.5 * 2 * PI is calculated:
401  res += PI * DEG2RAD *
402  (Integrand[i] * sin(za_grid[i] * DEG2RAD) +
403  Integrand[i + 1] * sin(za_grid[i + 1] * DEG2RAD)) *
404  (za_grid[i + 1] - za_grid[i]);
405  }
406 
407  return res;
408 }
409 
411 
424  if (x < 0)
425  return -1.0;
426  else if (x == 0)
427  return 0.0;
428  else
429  return 1.0;
430 }
431 
451 void mgd(VectorView psd,
452  const Vector& x,
453  const Numeric& n0,
454  const Numeric& mu,
455  const Numeric& la,
456  const Numeric& ga) {
457  const Index nx = x.nelem();
458 
459  assert(psd.nelem() == nx);
460 
461  if (ga == 1) {
462  if (mu == 0) {
463  // Exponential distribution
464  for (Index ix = 0; ix < nx; ix++) {
465  const Numeric eterm = exp(-la * x[ix]);
466  psd[ix] = n0 * eterm;
467  }
468  } else {
469  if (mu > 10) {
470  ostringstream os;
471  os << "Given mu is " << mu << endl
472  << "Seems unreasonable. Have you mixed up the inputs?";
473  throw runtime_error(os.str());
474  }
475  // Gamma distribution
476  for (Index ix = 0; ix < nx; ix++) {
477  const Numeric eterm = exp(-la * x[ix]);
478  const Numeric xterm = pow(x[ix], mu);
479  psd[ix] = n0 * xterm * eterm;
480  psd[ix] = n0 * pow(x[ix], mu) * exp(-la * x[ix]);
481  }
482  }
483  } else {
484  // Complete MGD
485  if (mu > 10) {
486  ostringstream os;
487  os << "Given mu is " << mu << endl
488  << "Seems unreasonable. Have you mixed up the inputs?";
489  throw runtime_error(os.str());
490  }
491  if (ga > 10) {
492  ostringstream os;
493  os << "Given gamma is " << ga << endl
494  << "Seems unreasonable. Have you mixed up the inputs?";
495  throw runtime_error(os.str());
496  }
497  for (Index ix = 0; ix < nx; ix++) {
498  const Numeric pterm = pow(x[ix], ga);
499  const Numeric eterm = exp(-la * pterm);
500  const Numeric xterm = pow(x[ix], mu);
501  psd[ix] = n0 * xterm * eterm;
502  }
503  }
504 }
505 
531  MatrixView jac_data,
532  const Vector& x,
533  const Numeric& n0,
534  const Numeric& mu,
535  const Numeric& la,
536  const Numeric& ga,
537  const bool& do_n0_jac,
538  const bool& do_mu_jac,
539  const bool& do_la_jac,
540  const bool& do_ga_jac) {
541  const Index nx = x.nelem();
542 
543  assert(psd.nelem() == nx);
544  assert(jac_data.nrows() == 4);
545  assert(jac_data.ncols() == nx);
546 
547  if (ga == 1 && !do_ga_jac) {
548  if (mu == 0 && !do_mu_jac) {
549  // Exponential distribution
550  for (Index ix = 0; ix < nx; ix++) {
551  const Numeric eterm = exp(-la * x[ix]);
552  psd[ix] = n0 * eterm;
553  if (do_n0_jac) {
554  jac_data(0, ix) = eterm;
555  }
556  if (do_la_jac) {
557  jac_data(2, ix) = -x[ix] * psd[ix];
558  }
559  }
560  } else {
561  if (mu > 10) {
562  ostringstream os;
563  os << "Given mu is " << mu << endl
564  << "Seems unreasonable. Have you mixed up the inputs?";
565  throw runtime_error(os.str());
566  }
567  // Gamma distribution
568  for (Index ix = 0; ix < nx; ix++) {
569  const Numeric eterm = exp(-la * x[ix]);
570  const Numeric xterm = pow(x[ix], mu);
571  psd[ix] = n0 * xterm * eterm;
572  if (do_n0_jac) {
573  jac_data(0, ix) = xterm * eterm;
574  }
575  if (do_mu_jac) {
576  jac_data(1, ix) = log(x[ix]) * psd[ix];
577  }
578  if (do_la_jac) {
579  jac_data(2, ix) = -x[ix] * psd[ix];
580  }
581  psd[ix] = n0 * pow(x[ix], mu) * exp(-la * x[ix]);
582  }
583  }
584  } else {
585  // Complete MGD
586  if (mu > 10) {
587  ostringstream os;
588  os << "Given mu is " << mu << endl
589  << "Seems unreasonable. Have you mixed up the inputs?";
590  throw runtime_error(os.str());
591  }
592  if (ga > 10) {
593  ostringstream os;
594  os << "Given gamma is " << ga << endl
595  << "Seems unreasonable. Have you mixed up the inputs?";
596  throw runtime_error(os.str());
597  }
598  for (Index ix = 0; ix < nx; ix++) {
599  const Numeric pterm = pow(x[ix], ga);
600  const Numeric eterm = exp(-la * pterm);
601  const Numeric xterm = pow(x[ix], mu);
602  psd[ix] = n0 * xterm * eterm;
603  if (do_n0_jac) {
604  jac_data(0, ix) = xterm * eterm;
605  }
606  if (do_mu_jac) {
607  jac_data(1, ix) = log(x[ix]) * psd[ix];
608  }
609  if (do_la_jac) {
610  jac_data(2, ix) = -pterm * psd[ix];
611  }
612  if (do_ga_jac) {
613  jac_data(3, ix) = -la * pterm * log(x[ix]) * psd[ix];
614  }
615  }
616  }
617 }
618 
620  MatrixView jac_data,
621  const Vector& x,
622  const Numeric& alpha,
623  const Numeric& beta) {
624  Numeric f_c = tgamma(4.0) / 256.0;
625  f_c *= pow(tgamma((alpha + 5.0) / beta), 4 + alpha);
626  f_c /= pow(tgamma((alpha + 4.0) / beta), 5 + alpha);
627 
628  Numeric f_d = tgamma((alpha + 5.0) / beta);
629  f_d /= tgamma((alpha + 4.0) / beta);
630 
631  for (Index i = 0; i < x.nelem(); ++i) {
632  Numeric xi = x[i];
633  psd[i] = beta * f_c * pow(xi, alpha) * exp(-pow(f_d * xi, beta));
634  jac_data(0, i) =
635  psd[i] * (alpha / xi - beta * f_d * pow(f_d * xi, beta - 1.0));
636  }
637 }
638 
640 
654  Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma) {
655  Numeric dN;
656 
657  if (x > 0. && N0 > 0. && Lambda > 0. && (mu + 1) / gamma > 0.) {
658  //Distribution function
659  dN = N0 * pow(x, mu) * exp(-Lambda * pow(x, gamma));
660 
661  return dN;
662  } else {
663  ostringstream os;
664  os << "At least one argument is zero or negative.\n"
665  << "Modified gamma distribution can not be calculated.\n"
666  << "x = " << x << "\n"
667  << "N0 = " << N0 << "\n"
668  << "lambda = " << Lambda << "\n"
669  << "mu = " << mu << "\n"
670  << "gamma = " << gamma << "\n";
671 
672  throw runtime_error(os.str());
673  }
674 }
675 
677 
687 void unitl(Vector& x) {
688  assert(x.nelem() > 0);
689 
690  const Numeric l = sqrt(x * x);
691  for (Index i = 0; i < x.nelem(); i++) x[i] /= l;
692 }
693 
695 
708  assert(x.nelem() == X.nrows() * X.ncols());
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  x[i] = X(r, c);
715  i += 1;
716  }
717  }
718 }
719 
721 
734  assert(x.nelem() == X.nrows() * X.ncols() * X.npages());
735 
736  Index i = 0;
737 
738  for (Index c = 0; c < X.ncols(); c++) {
739  for (Index r = 0; r < X.nrows(); r++) {
740  for (Index p = 0; p < X.npages(); p++) {
741  x[i] = X(p, r, c);
742  i += 1;
743  }
744  }
745  }
746 }
747 
749 
762  assert(x.nelem() == X.nrows() * X.ncols() * X.npages());
763 
764  Index i = 0;
765 
766  for (Index c = 0; c < X.ncols(); c++) {
767  for (Index r = 0; r < X.nrows(); r++) {
768  for (Index p = 0; p < X.npages(); p++) {
769  X(p, r, c) = x[i];
770  i += 1;
771  }
772  }
773  }
774 }
775 
777 
790  assert(x.nelem() == X.nrows() * X.ncols());
791 
792  Index i = 0;
793 
794  for (Index c = 0; c < X.ncols(); c++) {
795  for (Index r = 0; r < X.nrows(); r++) {
796  X(r, c) = x[i];
797  i += 1;
798  }
799  }
800 }
801 
803 
814  Index N = nph * 2;
815 
816  //directions in total
817  nlinspace(x, -1, 1, N);
818 
819  //allocate
820  w.resize(x.nelem());
821 
822  // calculate weights
823  w[0] = (x[1] - x[0]) / 2.;
824 
825  for (Index i = 1; i < nph * 2 - 1; i++) {
826  w[i] = (x[i + 1] - x[i - 1]) / 2.;
827  }
828  w[x.nelem() - 1] = (x[x.nelem() - 1] - x[x.nelem() - 2]) / 2.;
829 }
mgd
void mgd(VectorView psd, const Vector &x, const Numeric &n0, const Numeric &mu, const Numeric &la, const Numeric &ga)
Definition: math_funcs.cc:451
fac
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:63
MatrixView
The MatrixView class.
Definition: matpackI.h:1093
LagrangeInterpol4
Numeric LagrangeInterpol4(ConstVectorView x, ConstVectorView y, const Numeric a)
Lagrange Interpolation (internal function).
Definition: math_funcs.cc:113
w
Complex w(Complex z) noexcept
The Faddeeva function.
Definition: linefunctions.cc:42
sign
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:423
temp
#define temp
Definition: legacy_continua.cc:20951
last
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:165
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
AngIntegrate_trapezoid_opti
Numeric AngIntegrate_trapezoid_opti(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView grid_stepsize)
AngIntegrate_trapezoid_opti.
Definition: math_funcs.cc:340
PI
const Numeric PI
array.h
This file contains the definition of Array.
is_size
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:79
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620
ConstTensor3View::npages
Index npages() const
Returns the number of pages.
Definition: matpackIII.h:144
nlogspace
void nlogspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlogspace
Definition: math_funcs.cc:267
Tensor3View
The Tensor3View class.
Definition: matpackIII.h:239
Array
This can be used to make arrays out of anything.
Definition: array.h:108
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
flat
void flat(VectorView x, ConstMatrixView X)
flat
Definition: math_funcs.cc:707
unitl
void unitl(Vector &x)
unitl
Definition: math_funcs.cc:687
DEG2RAD
const Numeric DEG2RAD
reshape
void reshape(Tensor3View X, ConstVectorView x)
reshape
Definition: math_funcs.cc:761
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
mod_gamma_dist
Numeric mod_gamma_dist(Numeric x, Numeric N0, Numeric Lambda, Numeric mu, Numeric gamma)
Generalized Modified Gamma Distribution.
Definition: math_funcs.cc:653
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
VectorView
The VectorView class.
Definition: matpackI.h:610
mgd_with_derivatives
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:530
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
math_funcs.h
delanoe_shape_with_derivative
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:619
AngIntegrate_trapezoid
Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid)
AngIntegrate_trapezoid.
Definition: math_funcs.cc:296
ARTS::Var::za_grid
Vector za_grid(Workspace &ws) noexcept
Definition: autoarts.h:7771
ConstTensor3View::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackIII.h:147
integer_div
Index integer_div(const Index &x, const Index &y)
integer_div
Definition: math_funcs.cc:87
nlinspace
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:231
ConstMatrixView
A constant view of a Matrix.
Definition: matpackI.h:982
N
#define N
Definition: rng.cc:164
logic.h
Header file for logic.cc.
Zeeman::start
constexpr Rational start(Rational Ju, Rational Jl, Polarization type) noexcept
Gives the lowest M for a polarization type of this transition.
Definition: zeemandata.h:77
is_multiple
bool is_multiple(const Index &x, const Index &y)
Checks if an integer is a multiple of another integer.
Definition: logic.cc:65
ConstTensor3View::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackIII.h:150
ConstTensor3View
A constant view of a Tensor3.
Definition: matpackIII.h:132
calculate_weights_linear
void calculate_weights_linear(Vector &x, Vector &w, const Index nph)
calculate_weights_linear
Definition: math_funcs.cc:813
ARTS::Var::x
Vector x(Workspace &ws) noexcept
Definition: autoarts.h:7346
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
linspace
void linspace(Vector &x, const Numeric start, const Numeric stop, const Numeric step)
linspace
Definition: math_funcs.cc:204
Vector
The Vector class.
Definition: matpackI.h:860
beta
#define beta
Definition: legacy_continua.cc:21435
ARTS::Var::aa_grid
Vector aa_grid(Workspace &ws) noexcept
Definition: autoarts.h:1717
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
mystring.h
This file contains the definition of String, the ARTS string class.