GRASS GIS 8 Programmer's Manual 8.2.1(2023)-exported
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 */
48double 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 */
106int 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}
double b
double r
double IL_crst(double r, double fi)
Definition: func2d.c:48
int IL_crstg(double r, double fi, double *gd1, double *gd2)
Definition: func2d.c:106
#define x