39
|
1 |
#======================================================================
|
|
2 |
# . N M I N T E R P . S P L I N E
|
|
3 |
# doc: Wed Nov 22 21:01:09 2000
|
|
4 |
# dlm: Fri Sep 23 16:52:31 2011
|
|
5 |
# (c) 2000 A.M. Thurnherr
|
|
6 |
# uE-Info: 19 50 NIL 0 0 72 2 2 4 NIL ofnI
|
|
7 |
#======================================================================
|
|
8 |
|
|
9 |
# spline interpolation for [fillgaps]
|
|
10 |
|
|
11 |
# HISTORY:
|
|
12 |
# Apr 3, 2006: - created from [.interp.spline]
|
|
13 |
# Jul 28, 2006: - made to work
|
|
14 |
# Jul 30, 2006: - added max-gap support
|
|
15 |
# - BUG: extended interpolation into last interval
|
|
16 |
# - BUG: assertions made it incompatible with decreasing data
|
|
17 |
# Aug 8, 2008: - renamed
|
|
18 |
# Sep 23, 2011: - removed $xv from interpolate() args
|
|
19 |
# - added error when xf<0 (%RECNO)
|
|
20 |
|
|
21 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
22 |
#
|
|
23 |
# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
|
|
24 |
#
|
|
25 |
# $ISOpts string of allowed options
|
|
26 |
# $ISOptsUsage usage information string for options
|
|
27 |
#
|
|
28 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
29 |
|
|
30 |
$ISOpts = '';
|
|
31 |
$ISOptsUsage = '';
|
|
32 |
|
|
33 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
34 |
#
|
|
35 |
# &ISUsage() mangle parameters (options, really)
|
|
36 |
#
|
|
37 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
38 |
|
|
39 |
sub ISUsage() {}
|
|
40 |
|
|
41 |
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
42 |
#
|
|
43 |
# &ISInit(f,xf) init interpolation of field f
|
|
44 |
# f field number
|
|
45 |
# xf x field number
|
|
46 |
# <ret val> none
|
|
47 |
#
|
|
48 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
49 |
|
|
50 |
my(@y2);
|
|
51 |
|
|
52 |
sub ISInit($$) {
|
|
53 |
my($f,$xf) = @_;
|
|
54 |
my($k,$i,$p,$qn,$sig,$un,@u);
|
|
55 |
|
|
56 |
croak(".nminterp.spline: cannot use %RECNO as x-field (implementation restriction)\n")
|
|
57 |
unless ($xf >= 0);
|
|
58 |
|
|
59 |
my($pi,$ppi,$li); # find idx of first 3 & last valid recs
|
|
60 |
for (my($r)=0; $r<=$#ants_; $r++) {
|
|
61 |
if (numberp($ants_[$r][$f])) {
|
|
62 |
unless (defined($ppi)) {
|
|
63 |
$ppi = $pi; $pi = $i; $i = $r;
|
|
64 |
}
|
|
65 |
$li = $r;
|
|
66 |
}
|
|
67 |
}
|
|
68 |
unless (defined($ppi)) {
|
|
69 |
&antsInfo("WARNING: not enough <$antsLayout[$f]> data for spline interpolation");
|
|
70 |
return;
|
|
71 |
}
|
|
72 |
|
|
73 |
|
|
74 |
my($fpi) = $pi;
|
|
75 |
$y2[$pi][$f] = $u[$pi] = 0;
|
|
76 |
|
|
77 |
for (; $i<=$li; $i++) {
|
|
78 |
next unless numberp($ants_[$i][$f]);
|
|
79 |
|
|
80 |
$sig = ($ants_[$pi][$xf]-$ants_[$ppi][$xf]) /
|
|
81 |
($ants_[$i][$xf]-$ants_[$ppi][$xf]);
|
|
82 |
$p = $sig*$y2[$pi][$f] + 2;
|
|
83 |
$y2[$i][$f] = ($sig-1)/$p;
|
|
84 |
$u[$i] = ($ants_[$i][$f]-$ants_[$pi][$f]) /
|
|
85 |
($ants_[$i][$xf]-$ants_[$pi][$xf])
|
|
86 |
- ($ants_[$pi][$f]-$ants_[$ppi][$f]) /
|
|
87 |
($ants_[$pi][$xf]-$ants_[$ppi][$xf]);
|
|
88 |
$u[$i] = (6*$u[$i]/($ants_[$i][$xf]-$ants_[$ppi][$xf]) -
|
|
89 |
$sig * $u[$pi]) / $p;
|
|
90 |
|
|
91 |
$ppi = $pi; $pi = $i;
|
|
92 |
}
|
|
93 |
|
|
94 |
$qn = $un = 0;
|
|
95 |
|
|
96 |
$y2[$li][$f] = ($un-$qn*$u[$pi]) / ($qn*$y2[$pi][$f]+1);
|
|
97 |
my($pk) = $li;
|
|
98 |
for ($k=$pi; $k>=$fpi; $k--) {
|
|
99 |
next unless numberp($ants_[$k][$f]);
|
|
100 |
$y2[$k][$f] = $y2[$k][$f] * $y2[$pk][$f] + $u[$k];
|
|
101 |
$pk = $k;
|
|
102 |
}
|
|
103 |
}
|
|
104 |
|
|
105 |
#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
106 |
#
|
|
107 |
# &interpolate(xf,mg,f,lvr,cr) interpolate field f
|
|
108 |
# xf x field
|
|
109 |
# mg max allowed x gap
|
|
110 |
# f field number to interpolate
|
|
111 |
# lvr last record with valid y val
|
|
112 |
# cr current record
|
|
113 |
# <ret val> interpolated value
|
|
114 |
#
|
|
115 |
# NB:
|
|
116 |
# - handle f == xf
|
|
117 |
# - return undef if interpolation cannot be carried out
|
|
118 |
# - x VALUES MAY NOT BE MONOTONIC
|
|
119 |
#
|
|
120 |
#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
|
|
121 |
|
|
122 |
sub interpolate($$$$$$)
|
|
123 |
{
|
|
124 |
my($xf,$mg,$f,$lvr,$cr) = @_;
|
|
125 |
|
|
126 |
return $ants_[$cr][$xf] if ($xf == $f);
|
|
127 |
|
|
128 |
return $ants_[$lvr][$f] # edge values are ok
|
|
129 |
if equal($ants_[$cr][$xf],$ants_[$lvr][$xf]);
|
|
130 |
|
|
131 |
my($nvr1) = $cr + 1; # find next
|
|
132 |
while ($nvr1 <= $#y2 && !numberp($y2[$nvr1][$f])) { $nvr1++; }
|
|
133 |
return undef # none or gap too large
|
|
134 |
if ($nvr1>$#ants_ || abs($ants_[$nvr1][$xf]-$ants_[$lvr][$xf])>$mg);
|
|
135 |
return $ants_[$nvr1][$f] # edge values are ok
|
|
136 |
if equal($ants_[$cr][$xf],$ants_[$nvr1][$xf]);
|
|
137 |
|
|
138 |
my($nvr2) = $nvr1 + 1;
|
|
139 |
while ($nvr2 <= $#y2 && !numberp($y2[$nvr2][$f])) { $nvr2++; }
|
|
140 |
|
|
141 |
my($h) = $ants_[$nvr1][$xf] - $ants_[$lvr][$xf];
|
|
142 |
my($a) = ($ants_[$nvr1][$xf] - $ants_[$cr][$xf]) / $h;
|
|
143 |
my($b) = ($ants_[$cr][$xf] - $ants_[$lvr][$xf]) / $h;
|
|
144 |
|
|
145 |
return $a*$ants_[$lvr][$f] + $b*$ants_[$nvr1][$f] +
|
|
146 |
(($a*$a*$a-$a) * $y2[$nvr1][$f] +
|
|
147 |
($b*$b*$b-$b) * $y2[$nvr2][$f]) *
|
|
148 |
($h*$h)/6;
|
|
149 |
}
|
|
150 |
|
|
151 |
#----------------------------------------------------------------------
|
|
152 |
|
|
153 |
1;
|