ARTS  2.2.66
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 
22 /*****************************************************************************
23  *** File description
24  *****************************************************************************/
25 
36 /*****************************************************************************
37  *** External declarations
38  *****************************************************************************/
39 
40 #include <iostream>
41 #include <cmath>
42 #include <stdexcept>
43 #include "array.h"
44 #include "math_funcs.h"
45 #include "logic.h"
46 #include "mystring.h"
47 extern const Numeric DEG2RAD;
48 extern const Numeric PI;
49 
50 
51 
52 /*****************************************************************************
53  *** The functions (in alphabetical order)
54  *****************************************************************************/
55 
57 
68 Numeric fac(const Index n)
69 {
70  Numeric sum;
71 
72  if (n == 0) return (1.0);
73 
74  sum = 1.0;
75  for (Index i = 1; i <= n; i++)
76  sum *= Numeric(i);
77 
78  return(sum);
79 }
80 
81 
83 
95 Index integer_div( const Index& x, const Index& y )
96 {
97  assert( is_multiple( x, y ) );
98  return x/y;
99 }
100 
101 
102 
104 
125  ConstVectorView y,
126  const Numeric a)
127 {
128  // lowermost grid spacing on x-axis
129  const Numeric Dlimit = 1.00000e-15;
130 
131  // Check that dimensions of x and y vector agree
132  const Index n_x = x.nelem();
133  const Index n_y = y.nelem();
134  if ( (n_x != 4) || (n_y != 4) )
135  {
136  ostringstream os;
137  os << "The vectors x and y must all have the same length of 4 elements!\n"
138  << "Actual lengths:\n"
139  << "x:" << n_x << ", " << "y:" << n_y << ".";
140  throw runtime_error(os.str());
141  }
142 
143  // assure that x1 =< a < x2 holds
144  if ( (a < x[1]) || (a > x[2]) )
145  {
146  ostringstream os;
147  os << "LagrangeInterpol4: the relation x[1] =< a < x[2] is not satisfied. "
148  << "No interpolation can be calculated.\n";
149  throw runtime_error(os.str());
150  };
151 
152  // calculate the Lagrange polynomial coefficients for a polynomial of the order of 3
153  Numeric b[4];
154  for (Index i=0 ; i < 4 ; ++i)
155  {
156  b[i] = 1.000e0;
157  for (Index k=0 ; k < 4 ; ++k)
158  {
159  if ( (k != i) && (fabs(x[i]-x[k]) > Dlimit) )
160  b[i] = b[i] * ( (a-x[k]) / (x[i]-x[k]) );
161  };
162  };
163 
164  Numeric ya = 0.000e0;
165  for (Index i=0 ; i < n_x ; ++i) ya = ya + b[i]*y[i];
166 
167  return ya;
168 }
169 
170 
171 
172 
174 
184 {
185  assert( x.nelem() > 0 );
186  return x[x.nelem()-1];
187 }
188 
189 
190 
192 
202 {
203  assert( x.nelem() > 0 );
204  return x[x.nelem()-1];
205 }
206 
207 
208 
210 
228 void linspace(
229  Vector& x,
230  const Numeric start,
231  const Numeric stop,
232  const Numeric step )
233 {
234  Index n = (Index) floor( (stop-start)/step ) + 1;
235  if ( n<1 )
236  n=1;
237  x.resize(n);
238  for ( Index i=0; i<n; i++ )
239  x[i] = start + (double)i*step;
240 }
241 
242 
243 
245 
262  Vector& x,
263  const Numeric start,
264  const Numeric stop,
265  const Index n )
266 {
267  assert( 1<n ); // Number of points must be greatere 1.
268  x.resize(n);
269  Numeric step = (stop-start)/((double)n-1) ;
270  for ( Index i=0; i<n-1; i++ )
271  x[i] = start + (double)i*step;
272  x[n-1] = stop;
273 }
274 
275 
276 
278 
294 void nlogspace(
295  Vector& x,
296  const Numeric start,
297  const Numeric stop,
298  const Index n )
299 {
300  // Number of points must be greater than 1:
301  assert( 1<n );
302  // Only positive numbers are allowed for start and stop:
303  assert( 0<start );
304  assert( 0<stop );
305 
306  x.resize(n);
307  Numeric a = log(start);
308  Numeric step = (log(stop)-a)/((double)n-1);
309  x[0] = start;
310  for ( Index i=1; i<n-1; i++ )
311  x[i] = exp(a + (double)i*step);
312  x[n-1] = stop;
313 }
314 
315 
317 
328  ConstVectorView za_grid,
329  ConstVectorView aa_grid)
330 {
331 
332  Index n = za_grid.nelem();
333  Index m = aa_grid.nelem();
334  Vector res1(n);
335  assert (is_size(Integrand, n, m));
336 
337  for (Index i = 0; i < n ; ++i)
338  {
339  res1[i] = 0.0;
340 
341  for (Index j = 0; j < m - 1; ++j)
342  {
343  res1[i] += 0.5 * DEG2RAD * (Integrand(i, j) + Integrand(i, j + 1)) *
344  (aa_grid[j + 1] - aa_grid[j]) * sin(za_grid[i] * DEG2RAD);
345  }
346  }
347  Numeric res = 0.0;
348  for (Index i = 0; i < n - 1; ++i)
349  {
350  res += 0.5 * DEG2RAD * (res1[i] + res1[i + 1]) *
351  (za_grid[i + 1] - za_grid[i]);
352  }
353 
354  return res;
355 }
356 
357 
359 
378  ConstVectorView za_grid,
379  ConstVectorView aa_grid,
380  ConstVectorView grid_stepsize)
381 {
382  Numeric res = 0;
383  if ((grid_stepsize[0] > 0) && (grid_stepsize[1] > 0))
384  {
385  Index n = za_grid.nelem();
386  Index m = aa_grid.nelem();
387  Numeric stepsize_za = grid_stepsize[0];
388  Numeric stepsize_aa = grid_stepsize[1];
389  Vector res1(n);
390  assert (is_size(Integrand, n, m));
391 
392  Numeric temp = 0.0;
393 
394  for (Index i = 0; i < n ; ++i)
395  {
396  temp = Integrand(i, 0);
397  for (Index j = 1; j < m - 1; j++)
398  {
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  {
409  res += res1[i] * 2;
410  }
411  res += res1[n-1];
412  res *= 0.5 * DEG2RAD * stepsize_za;
413  }
414  else
415  {
416  res = AngIntegrate_trapezoid(Integrand, za_grid, aa_grid);
417  }
418 
419  return res;
420 }
421 
422 
424 
439  ConstVectorView za_grid)
440 {
441 
442  Index n = za_grid.nelem();
443  assert (is_size(Integrand, n));
444 
445  Numeric res = 0.0;
446  for (Index i = 0; i < n - 1; ++i)
447  {
448  // in this place 0.5 * 2 * PI is calculated:
449  res += PI * DEG2RAD * (Integrand[i]* sin(za_grid[i] * DEG2RAD)
450  + Integrand[i + 1] * sin(za_grid[i + 1] * DEG2RAD))
451  * (za_grid[i + 1] - za_grid[i]);
452  }
453 
454  return res;
455 }
456 
457 
458 
459 
461 
473 Numeric sign( const Numeric& x )
474 {
475  if( x < 0 )
476  return -1.0;
477  else if( x == 0 )
478  return 0.0;
479  else
480  return 1.0;
481 }
482 
483 
485 
498 {
499  //double lgamma(double xx);
500 
501  Numeric gam;
502  Index i;
503 
504 
505  if (xx > 0.0) {
506  if (xx == (int)xx) {
507  gam = 1.0; // use factorial
508  for (i=2;i<xx;i++) {
509  gam *= (Numeric)i;
510  }
511  }
512  else {
513  return exp(lgamma_func(xx));
514  }
515  } else {
516  ostringstream os;
517  os << "Argument is zero or negative."
518  << "Gamma function can not be calculated.\n";
519  throw runtime_error(os.str());
520  }
521 
522 
523  return gam;
524 }
525 
526 
528 
539 {
540 
541  Numeric x,y,tmp,ser;
542  static const Numeric cof[6] = {
543  76.18009172947146, -86.50532032941677, 24.01409824083091,
544  -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5
545  };
546 
547  if (xx > 0.0)
548  {
549  y=x=xx;
550  tmp = x+5.5;
551  tmp -= (x+0.5)*log(tmp);
552  ser = 1.000000000190015;
553  for (Index j=0;j<=5;j++) ser += cof[j]/++y;
554  return -tmp+log(2.5066282746310005*ser/x);
555  }
556  else
557  {
558  ostringstream os;
559  os << "Argument is zero or negative.\n"
560  << "log Gamma function can not be calculated.\n";
561  throw runtime_error(os.str());
562  }
563 }
564 
565 
566 
568 
578 void unitl( Vector& x )
579 {
580  assert( x.nelem() > 0 );
581 
582  const Numeric l = sqrt(x*x);
583  for(Index i=0; i<x.nelem(); i++ )
584  x[i] /= l;
585 }
586 
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
void unitl(Vector &x)
lunit
Definition: math_funcs.cc:578
const Numeric PI
Index nelem() const
Number of elements.
Definition: array.h:176
Numeric AngIntegrate_trapezoid_opti(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid, ConstVectorView grid_stepsize)
AngIntegrate_trapezoid_opti.
Definition: math_funcs.cc:377
The Vector class.
Definition: matpackI.h:556
const Numeric DEG2RAD
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:68
Numeric LagrangeInterpol4(ConstVectorView x, ConstVectorView y, const Numeric a)
Lagrange Interpolation (internal function).
Definition: math_funcs.cc:124
Numeric last(ConstVectorView x)
last
Definition: math_funcs.cc:183
bool is_multiple(const Index &x, const Index &y)
Checks if an integer is a multiple of another integer.
Definition: logic.cc:74
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
Numeric sign(const Numeric &x)
sign
Definition: math_funcs.cc:473
This file contains the definition of Array.
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
Definition: math_funcs.cc:261
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
Index integer_div(const Index &x, const Index &y)
integer_div
Definition: math_funcs.cc:95
void nlogspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlogspace
Definition: math_funcs.cc:294
Header file for logic.cc.
This can be used to make arrays out of anything.
Definition: array.h:40
void resize(Index n)
Assignment operator from VectorView.
Definition: matpackI.cc:798
A constant view of a Vector.
Definition: matpackI.h:292
#define temp
Definition: continua.cc:20773
A constant view of a Matrix.
Definition: matpackI.h:596
bool is_size(ConstVectorView x, const Index &n)
Verifies that the size of x is l.
Definition: logic.cc:91
Numeric gamma_func(Numeric xx)
Gamma Function.
Definition: math_funcs.cc:497
void linspace(Vector &x, const Numeric start, const Numeric stop, const Numeric step)
linspace
Definition: math_funcs.cc:228
Numeric lgamma_func(Numeric xx)
ln Gamma Function
Definition: math_funcs.cc:538
Numeric AngIntegrate_trapezoid(ConstMatrixView Integrand, ConstVectorView za_grid, ConstVectorView aa_grid)
AngIntegrate_trapezoid.
Definition: math_funcs.cc:327
This file contains the definition of String, the ARTS string class.