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