GRASS GIS 8 Programmer's Manual 8.4.1(2025)-45ca3179ab
Loading...
Searching...
No Matches
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
46double 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
106int 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 b
double r
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
#define x