-
Notifications
You must be signed in to change notification settings - Fork 2
/
splinesc.c
73 lines (63 loc) · 1.78 KB
/
splinesc.c
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
/* part of the shansyn spherical harmonics package, see COPYRIGHT for license */
/* $Id: splinesc.c,v 1.3 2001/03/11 18:20:04 becker Exp becker $ */
#include "shansyn.h"
/*
spline interpolation routine, calls FORTRAN code
due to Ritsema as distributed with the s20rts model
spline knots are equally spaced from -1 .. 1 in m steps
when first called
x1, x2: left and right boundaries of integral, don't change
once initialized
c[0...m-1]: array with weights for splines 0 ... m-1
m: number of splines
x: location for interpolation
$Id: splinesc.c,v 1.3 2001/03/11 18:20:04 becker Exp becker $
*/
#define MAX_SPLINE_N 50
COMP_PRECISION spline_base(COMP_PRECISION x1,COMP_PRECISION x2,
COMP_PRECISION *c,int m,
COMP_PRECISION x)
{
static int init=0;
int i;
// static arrays for spline setup
static COMP_PRECISION
qq0[MAX_SPLINE_N*MAX_SPLINE_N],
qq[3*MAX_SPLINE_N*MAX_SPLINE_N],
spknt[MAX_SPLINE_N],
qqwk[3*MAX_SPLINE_N],dx;
COMP_PRECISION xp,val;
// check range
if ((x-x1)*(x-x2) > 0.0){
fprintf(stderr,"spline_base: x (%g) not in range (%g - %g)\n",
x,x1,x2);
exit(-1);
}
if(!init){
if(m > MAX_SPLINE_N){
fprintf(stderr,"spline_base: works only for order <= %i (m: %i)\n",
MAX_SPLINE_N,m);
exit(-1);
}
splhsetup(spknt,qq0,qq,qqwk,&m);
dx=(x2 - x1)/2.0;
init=1;
}
// local coordinate from -1 to 1
xp=-1.0+(x-x1)/dx;
val=0.0;
for(i=0;i<m;i++)// add up contribution from all base functions
val += c[i]*splh(&m,spknt,qq0,qq,&i,&xp);
return(val);
}
#undef MAX_SPLINE_N
// norm damping for splines
COMP_PRECISION spline_norm_damping(int i,int j,int n)
{
if(i!=j)
return 0.0;
else{
// all damped the same way
return 1.0;
}
}