GRASS GIS 7 Programmer's Manual  7.9.dev(2021)-e5379bbd7
func2d.c
Go to the documentation of this file.
1 
2 /*!
3  * \file func2d.c
4  *
5  * \author
6  * Lubos Mitas (original program and various modifications)
7  *
8  * \author
9  * H. Mitasova,
10  * I. Kosinovsky, D. Gerdes,
11  * D. McCauley
12  * (GRASS4.1 version of the program and GRASS4.2 modifications)
13  *
14  * \author
15  * L. Mitas ,
16  * H. Mitasova ,
17  * I. Kosinovsky,
18  * D.Gerdes
19  * D. McCauley (1993, 1995)
20  *
21  * \author modified by McCauley in August 1995
22  * \author modified by Mitasova in August 1995, Nov. 1996
23  *
24  * \copyright
25  * (C) 1993-1999 by Lubos Mitas and the GRASS Development Team
26  *
27  * \copyright
28  * This program is free software under the
29  * GNU General Public License (>=v2).
30  * Read the file COPYING that comes with GRASS
31  * for details.
32  */
33 
34 
35 #include <stdio.h>
36 #include <math.h>
37 #include <grass/gis.h>
38 #include <grass/interpf.h>
39 
40 
41 /* parameter description from DESCRIPTION.INTERP */
42 /*!
43  * Radial basis function
44  *
45  * Radial basis function - completely regularized spline with tension (d=2)
46  *
47  */
48 double IL_crst(double r, /**< distance squared */
49  double fi /**< tension */
50  )
51 {
52  double rfsta2 = fi * fi * r / 4.;
53 
54  static double c[4] = { 8.5733287401, 18.0590169730, 8.6347608925,
55  0.2677737343
56  };
57  static double b[4] = { 9.5733223454, 25.6329561486, 21.0996530827,
58  3.9584969228
59  };
60  double ce = 0.57721566;
61 
62  static double u[10] = { 1.e+00, -.25e+00,
63  .055555555555556e+00, -.010416666666667e+00, /*fixed bug 415.. repl. by 416.. */
64  .166666666666667e-02, -2.31481481481482e-04,
65  2.83446712018141e-05, -3.10019841269841e-06,
66  3.06192435822065e-07, -2.75573192239859e-08
67  };
68  double x = rfsta2;
69  double res;
70 
71  double e1, ea, eb;
72 
73 
74  if (x < 1.e+00) {
75  res = x * (u[0] + x * (u[1] + x * (u[2] + x * (u[3] + x * (u[4] + x *
76  (u[5] +
77  x *
78  (u[6] +
79  x *
80  (u[7] +
81  x *
82  (u[8] +
83  x *
84  u
85  [9])))))))));
86  return (res);
87  }
88 
89  if (x > 25.e+00)
90  e1 = 0.00;
91  else {
92  ea = c[3] + x * (c[2] + x * (c[1] + x * (c[0] + x)));
93  eb = b[3] + x * (b[2] + x * (b[1] + x * (b[0] + x)));
94  e1 = (ea / eb) / (x * exp(x));
95  }
96  res = e1 + ce + log(x);
97  return (res);
98 }
99 
100 
101 /*!
102  * Function for calculating derivatives (d=2)
103  *
104  * Derivatives of radial basis function - regularized spline with tension(d=2)
105  */
106 int IL_crstg(double r, /**< distance squared */
107  double fi, /**< tension */
108  double *gd1, /**< G1(r) */
109  double *gd2 /**< G2(r) */
110  )
111 {
112  double r2 = r;
113  double rfsta2 = fi * fi * r / 4.;
114  double x, exm, oneme, hold;
115  double fsta2 = fi * fi / 2.;
116 
117  x = rfsta2;
118  if (x < 0.001) {
119  *gd1 = 1. - x / 2. + x * x / 6. - x * x * x / 24.;
120  *gd2 = fsta2 * (-.5 + x / 3. - x * x / 8. + x * x * x / 30.);
121  }
122  else {
123  if (x < 35.e+00) {
124  exm = exp(-x);
125  oneme = 1. - exm;
126  *gd1 = oneme / x;
127  hold = x * exm - oneme;
128  *gd2 = (hold + hold) / (r2 * x);
129  }
130  else {
131  *gd1 = 1. / x;
132  *gd2 = -2. / (x * r2);
133  }
134  }
135  return 1;
136 }
int IL_crstg(double r, double fi, double *gd1, double *gd2)
Definition: func2d.c:106
#define x
double IL_crst(double r, double fi)
Definition: func2d.c:48
double b
Definition: r_raster.c:39
double r
Definition: r_raster.c:39