1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
|
/*-
*
* Original program and various modifications:
* Lubos Mitas
*
* GRASS4.1 version of the program and GRASS4.2 modifications:
* H. Mitasova,
* I. Kosinovsky, D. Gerdes
* D. McCauley
*
* Copyright 1993, 1995:
* L. Mitas ,
* H. Mitasova ,
* I. Kosinovsky,
* D.Gerdes
* D. McCauley
*
* modified by McCauley in August 1995
* modified by Mitasova in August 1995, Nov. 1996
*
*/
#include <stdio.h>
#include <math.h>
#include <grass/gis.h>
#include <grass/interpf.h>
double IL_crst(double r, double fi)
/*
* Radial basis function - completely regularized spline with
* tension (d=2)
*/
{
double rfsta2 = fi * fi * r / 4.;
static double c[4] = { 8.5733287401, 18.0590169730, 8.6347608925,
0.2677737343
};
static double b[4] = { 9.5733223454, 25.6329561486, 21.0996530827,
3.9584969228
};
double ce = 0.57721566;
static double u[10] = { 1.e+00, -.25e+00,
.055555555555556e+00, -.010416666666667e+00, /*fixed bug 415.. repl. by 416.. */
.166666666666667e-02, -2.31481481481482e-04,
2.83446712018141e-05, -3.10019841269841e-06,
3.06192435822065e-07, -2.75573192239859e-08
};
double x = rfsta2;
double res;
double e1, ea, eb;
if (x < 1.e+00) {
res = x * (u[0] + x * (u[1] + x * (u[2] + x * (u[3] + x * (u[4] + x *
(u[5] +
x *
(u[6] +
x *
(u[7] +
x *
(u[8] +
x *
u
[9])))))))));
return (res);
}
if (x > 25.e+00)
e1 = 0.00;
else {
ea = c[3] + x * (c[2] + x * (c[1] + x * (c[0] + x)));
eb = b[3] + x * (b[2] + x * (b[1] + x * (b[0] + x)));
e1 = (ea / eb) / (x * exp(x));
}
res = e1 + ce + log(x);
return (res);
}
int IL_crstg(double r, double fi, double *gd1, /* G1(r) */
double *gd2 /* G2(r) */
)
/*
* Function for calculating derivatives (d=2)
*/
{
double r2 = r;
double rfsta2 = fi * fi * r / 4.;
double x, exm, oneme, hold;
double fsta2 = fi * fi / 2.;
x = rfsta2;
if (x < 0.001) {
*gd1 = 1. - x / 2. + x * x / 6. - x * x * x / 24.;
*gd2 = fsta2 * (-.5 + x / 3. - x * x / 8. + x * x * x / 30.);
}
else {
if (x < 35.e+00) {
exm = exp(-x);
oneme = 1. - exm;
*gd1 = oneme / x;
hold = x * exm - oneme;
*gd2 = (hold + hold) / (r2 * x);
}
else {
*gd1 = 1. / x;
*gd2 = -2. / (x * r2);
}
}
return 1;
}
|
