[system] / trunk / pg / lib / Hermite.pm Repository:
ViewVC logotype

Diff of /trunk/pg/lib/Hermite.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

Revision 1050 Revision 1079
1#!/usr/bin/perl -w 1
2 2
3use strict; 3use strict;
4use Carp; 4use Carp;
5 5
6=head1 NAME 6=head1 NAME
9 9
10=head1 SYNPOSIS 10=head1 SYNPOSIS
11 11
12 Usage: 12 Usage:
13 $obj = new Hermit(\@x_values, \y_valuses \@yp_values); 13 $obj = new Hermit(\@x_values, \y_valuses \@yp_values);
14 14
15 #get and set methods 15 #get and set methods
16 $ra_x_values = $obj -> ra_x(\@x_values); 16 $ra_x_values = $obj -> ra_x(\@x_values);
17 $ra_y_values = $obj -> ra_y; 17 $ra_y_values = $obj -> ra_y;
18 $ra_yp_values = $obj -> ra_yp; 18 $ra_yp_values = $obj -> ra_yp;
19 19
20 $obj -> initialize; # calculates the approximation 20 $obj -> initialize; # calculates the approximation
21 21
22 #get methods 22 #get methods
23 $rf_function = $obj -> rf_f; 23 $rf_function = $obj -> rf_f;
24 $rf_function_derivative = $obj -> rf_fp; 24 $rf_function_derivative = $obj -> rf_fp;
25 $rf_function_2nd_derivative = $obj -> rf_fpp; 25 $rf_function_2nd_derivative = $obj -> rf_fpp;
26 26
27 $rh_critical_points =$obj -> rh_critical_points 27 $rh_critical_points =$obj -> rh_critical_points
28 $rh_inflection_points =$obj -> rh_inflection_points 28 $rh_inflection_points =$obj -> rh_inflection_points
29 29
30 30
31 31
32 32
33=head1 DESCRIPTION 33=head1 DESCRIPTION
34 34
35This module defines an object containing a Hermite spline approximation to a function. 35This module defines an object containing a Hermite spline approximation to a function.
36 The approximation 36 The approximation
37consists of a piecewise cubic polynomial which agrees with the original 37consists of a piecewise cubic polynomial which agrees with the original
38function and its first derivative at 38function and its first derivative at
39the node points. 39the node points.
40 40
41This is useful for creating on the fly graphics. Care must be taken to use a small 41This is useful for creating on the fly graphics. Care must be taken to use a small
42number of points spaced reasonably far apart, preferably 42number of points spaced reasonably far apart, preferably
43points with alternating slopes, since this will minimize 43points with alternating slopes, since this will minimize
44the number of extraneous critical points 44the number of extraneous critical points
45introduced. Too many points will introduce many small variations and 45introduced. Too many points will introduce many small variations and
46a large number of extraneous critical points. 46a large number of extraneous critical points.
47 47
48There are even more extraneous inflections points. 48There are even more extraneous inflections points.
49This parameter is probably not very useful. A different 49This parameter is probably not very useful. A different
50approximation scheme needs to be used. 50approximation scheme needs to be used.
51 51
52=cut 52=cut
53 53
83 my $ra_y = shift; 83 my $ra_y = shift;
84 my $ra_yp = shift; 84 my $ra_yp = shift;
85 $self->ra_x($ra_x); 85 $self->ra_x($ra_x);
86 $self->ra_y($ra_y); 86 $self->ra_y($ra_y);
87 $self->ra_yp($ra_yp); 87 $self->ra_yp($ra_yp);
88 88
89 $self -> initialize(); 89 $self -> initialize();
90 ($self->{x_ref},$self->{y_ref}, $self->{yp_ref} ) 90 ($self->{x_ref},$self->{y_ref}, $self->{yp_ref} )
91} 91}
92 92
93sub initialize { 93sub initialize {
94 my $self = shift; 94 my $self = shift;
95 # define the functions rf_f 95 # define the functions rf_f
96 96
97 $self -> {rf_f} = hermite_spline($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} ); 97 $self -> {rf_f} = hermite_spline($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} );
98# # define the function rf_fp 98# # define the function rf_fp
99 $self -> {rf_fp} = hermite_spline_p($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} ); 99 $self -> {rf_fp} = hermite_spline_p($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} );
100# # define the function rf_fp 100# # define the function rf_fp
101 $self -> {rf_fpp} = hermite_spline_pp($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} ); 101 $self -> {rf_fpp} = hermite_spline_pp($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} );
106# # define the inflection_points 106# # define the inflection_points
107 %{$self -> {rh_inflection_points}} = inflection_points($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} , 107 %{$self -> {rh_inflection_points}} = inflection_points($self -> {ra_x}, $self -> {ra_y}, $self -> {ra_yp} ,
108 $self -> { rf_f }, 108 $self -> { rf_f },
109 ); 109 );
110# # define the maximum points 110# # define the maximum points
111 111
112 # define the minimum points 112 # define the minimum points
113 113
114} 114}
115 115
116# fix up these accesses to do more checking 116# fix up these accesses to do more checking
117sub ra_x { 117sub ra_x {
118 my $self = shift; 118 my $self = shift;
119 my $ra_x = shift; 119 my $ra_x = shift;
120 @{$self->{ra_x}} = @$ra_x if ref($ra_x) eq 'ARRAY'; 120 @{$self->{ra_x}} = @$ra_x if ref($ra_x) eq 'ARRAY';
121 $self->{ra_x}; 121 $self->{ra_x};
122 122
123} 123}
124 124
125sub ra_y { 125sub ra_y {
126 my $self = shift; 126 my $self = shift;
127 my $ra_y = shift; 127 my $ra_y = shift;
138 138
139sub rf_f { 139sub rf_f {
140 my $self = shift; 140 my $self = shift;
141 my $rf_f =shift; 141 my $rf_f =shift;
142 $self ->{rf_f} = $rf_f if defined( $rf_f ); 142 $self ->{rf_f} = $rf_f if defined( $rf_f );
143 $self ->{rf_f}; 143 $self ->{rf_f};
144} 144}
145sub rf_fp { 145sub rf_fp {
146 my $self = shift; 146 my $self = shift;
147 my $rf_fp =shift; 147 my $rf_fp =shift;
148 $self ->{rf_fp} = $rf_fp if defined( $rf_fp ); 148 $self ->{rf_fp} = $rf_fp if defined( $rf_fp );
149 $self ->{rf_fp}; 149 $self ->{rf_fp};
150} 150}
151sub rf_fpp { 151sub rf_fpp {
152 my $self = shift; 152 my $self = shift;
153 my $rf_fpp =shift; 153 my $rf_fpp =shift;
154 $self ->{rf_fpp} = $rf_fpp if defined( $rf_fpp ); 154 $self ->{rf_fpp} = $rf_fpp if defined( $rf_fpp );
155 $self ->{rf_fpp}; 155 $self ->{rf_fpp};
156} 156}
157 157
158sub rh_critical_points { 158sub rh_critical_points {
159 my $self = shift; 159 my $self = shift;
160 $self -> { rh_critical_points}; 160 $self -> { rh_critical_points};
198 %inflection_points; 198 %inflection_points;
199} 199}
200sub internal_critical_points{ 200sub internal_critical_points{
201 my ($x0,$l,$lp, $x1,$r,$rp, $rh_roots ,$rf_function) = @_; 201 my ($x0,$l,$lp, $x1,$r,$rp, $rh_roots ,$rf_function) = @_;
202 #data for one segment of the hermite spline 202 #data for one segment of the hermite spline
203 203
204 # coefficients for the approximating polynomial 204 # coefficients for the approximating polynomial
205 205
206 my @a = ( $l, 206 my @a = ( $l,
207 $lp, 207 $lp,
208 -$lp/2 + (3*(-$lp - 2*($l - $r) - $rp))/2 + $rp/2, 208 -$lp/2 + (3*(-$lp - 2*($l - $r) - $rp))/2 + $rp/2,
209 $lp + 2*($l - $r) + $rp 209 $lp + 2*($l - $r) + $rp
210 ); 210 );
214 if ( $a[2] == 0 ) { 214 if ( $a[2] == 0 ) {
215 } else { 215 } else {
216 $root1 = -$a[1]/( 2*$a[2] ); 216 $root1 = -$a[1]/( 2*$a[2] );
217 if ( 0 <= $root1 and $root1 < 1) { 217 if ( 0 <= $root1 and $root1 < 1) {
218 $z1 = $root1*($x1-$x0) + $x0; 218 $z1 = $root1*($x1-$x0) + $x0;
219 219
220 $rh_roots -> {$z1} = &$rf_function($z1); 220 $rh_roots -> {$z1} = &$rf_function($z1);
221 } 221 }
222 } 222 }
223 } else { 223 } else {
224 my $discriminent = (4*$a[2]**2 - 12*$a[1]*$a[3]); 224 my $discriminent = (4*$a[2]**2 - 12*$a[1]*$a[3]);
225 225
226 226
227 if ( $discriminent >= 0 ) { 227 if ( $discriminent >= 0 ) {
228 $discriminent = $discriminent**0.5; 228 $discriminent = $discriminent**0.5;
229 $root1 = (-2*$a[2] - $discriminent )/( 6*$a[3] ); 229 $root1 = (-2*$a[2] - $discriminent )/( 6*$a[3] );
230 $root2 = (-2*$a[2] + $discriminent )/( 6*$a[3] ); 230 $root2 = (-2*$a[2] + $discriminent )/( 6*$a[3] );
231 $z1 = $root1*($x1-$x0) + $x0; 231 $z1 = $root1*($x1-$x0) + $x0;
232 $z2 = $root2*($x1-$x0) + $x0; 232 $z2 = $root2*($x1-$x0) + $x0;
233 $rh_roots -> {$z1} = &$rf_function($z1) if 0 <= $root1 and $root1 < 1; 233 $rh_roots -> {$z1} = &$rf_function($z1) if 0 <= $root1 and $root1 < 1;
234 $rh_roots -> {$z2} = &$rf_function($z1) if 0 <= $root2 and $root2 < 1; 234 $rh_roots -> {$z2} = &$rf_function($z1) if 0 <= $root2 and $root2 < 1;
235 } 235 }
236 } 236 }
237} 237}
238 238
239sub internal_inflection_points { 239sub internal_inflection_points {
240 my ($x0,$l,$lp,,$x1,$r,$rp,$rh_roots,$rf_function) = @_; 240 my ($x0,$l,$lp,,$x1,$r,$rp,$rh_roots,$rf_function) = @_;
241 #data for one segment of the hermite spline 241 #data for one segment of the hermite spline
242 242
243 # coefficients for the approximating polynomial 243 # coefficients for the approximating polynomial
244 244
245 my @a = ( $l, 245 my @a = ( $l,
246 $lp, 246 $lp,
247 -$lp/2 + (3*(-$lp - 2*($l - $r) - $rp))/2 + $rp/2, 247 -$lp/2 + (3*(-$lp - 2*($l - $r) - $rp))/2 + $rp/2,
248 $lp + 2*($l - $r) + $rp 248 $lp + 2*($l - $r) + $rp
249 ); 249 );
250 if ($a[3] == 0 ) { 250 if ($a[3] == 0 ) {
251 } else { 251 } else {
252 my $root1 = -$a[2]/( 3*$a[3] ); 252 my $root1 = -$a[2]/( 3*$a[3] );
253 my $z1 = $root1*($x1-$x0) + $x0; 253 my $z1 = $root1*($x1-$x0) + $x0;
254 $rh_roots -> {$z1} = &$rf_function($z1) if 0 <= $root1 and $root1 < 1; 254 $rh_roots -> {$z1} = &$rf_function($z1) if 0 <= $root1 and $root1 < 1;
255 } 255 }
256 256
257} 257}
258 258
259 259
260 260
261sub cubic_hermite { 261sub cubic_hermite {
263 my @a; 263 my @a;
264 my $width = $x1 - $x0; 264 my $width = $x1 - $x0;
265 $yp0 = $yp0*$width; # normalize to unit width 265 $yp0 = $yp0*$width; # normalize to unit width
266 $yp1 = $yp1*$width; 266 $yp1 = $yp1*$width;
267 267
268 $a[0] = $y0; 268 $a[0] = $y0;
269 $a[1] = $yp0; 269 $a[1] = $yp0;
270 $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1; 270 $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1;
271 $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1; 271 $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1;
272 272
273 my $f = sub { 273 my $f = sub {
274 my $x = shift; 274 my $x = shift;
275 #normalize to unit width 275 #normalize to unit width
276 $x = ( $x - $x0 )/$width; 276 $x = ( $x - $x0 )/$width;
277 ( ($a[3]*$x + $a[2]) * $x + $a[1] )*$x + $a[0]; 277 ( ($a[3]*$x + $a[2]) * $x + $a[1] )*$x + $a[0];
278 278
279 }; 279 };
280 $f; 280 $f;
281} 281}
282 282
283sub cubic_hermite_p { 283sub cubic_hermite_p {
284 my ( $x0, $y0,$yp0, $x1, $y1, $yp1 )=@_; 284 my ( $x0, $y0,$yp0, $x1, $y1, $yp1 )=@_;
285 my @a; 285 my @a;
286 my $width = $x1 - $x0; 286 my $width = $x1 - $x0;
287 $yp0 = $yp0*$width; # normalize to unit width 287 $yp0 = $yp0*$width; # normalize to unit width
288 $yp1 = $yp1*$width; 288 $yp1 = $yp1*$width;
289 289
290 $a[0] = $y0; 290 $a[0] = $y0;
291 $a[1] = $yp0; 291 $a[1] = $yp0;
292 $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1; 292 $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1;
293 $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1; 293 $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1;
294 294
295 my $fp = sub { 295 my $fp = sub {
296 my $x = shift; 296 my $x = shift;
297 #normalize to unit width 297 #normalize to unit width
298 $x = ( $x - $x0 )/$width; 298 $x = ( $x - $x0 )/$width;
299 ( (3*$a[3]*$x + 2*$a[2]) * $x + $a[1] )/$width ; 299 ( (3*$a[3]*$x + 2*$a[2]) * $x + $a[1] )/$width ;
300 }; 300 };
301
302 301
303 302
303
304 $fp; 304 $fp;
305 305
306} 306}
307 307
308sub cubic_hermite_pp { 308sub cubic_hermite_pp {
309 my ( $x0, $y0,$yp0, $x1, $y1, $yp1 ) = @_; 309 my ( $x0, $y0,$yp0, $x1, $y1, $yp1 ) = @_;
310 my @a; 310 my @a;
311 my $width = $x1 - $x0; 311 my $width = $x1 - $x0;
312 $yp0 = $yp0*$width; # normalize to unit width 312 $yp0 = $yp0*$width; # normalize to unit width
313 $yp1 = $yp1*$width; 313 $yp1 = $yp1*$width;
314 314
315 $a[0] = $y0; 315 $a[0] = $y0;
316 $a[1] = $yp0; 316 $a[1] = $yp0;
317 $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1; 317 $a[2] = -3*$y0 - 2*$yp0 +3*$y1 -$yp1;
318 $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1; 318 $a[3] = 2*$y0 + $yp0 - 2*$y1 +$yp1;
319
320 319
321 320
321
322 my $fpp = sub { 322 my $fpp = sub {
323 my $x = shift; 323 my $x = shift;
324 #normalize to unit width 324 #normalize to unit width
325 $x = ( $x - $x0 )/$width; 325 $x = ( $x - $x0 )/$width;
326 (6*$a[3]*$x + 2*$a[2])/$width**2 ; 326 (6*$a[3]*$x + 2*$a[2])/$width**2 ;
327 }; 327 };
328 328
329 $fpp; 329 $fpp;
330} 330}
331 331
332 332
333 333
334sub hermite_spline { 334sub hermite_spline {
335 my ($xref, $yref, $ypref) = @_; 335 my ($xref, $yref, $ypref) = @_;
339 my $x0 = shift @xvals; 339 my $x0 = shift @xvals;
340 my $y0 = shift @yvals; 340 my $y0 = shift @yvals;
341 my $yp0 = shift @ypvals; 341 my $yp0 = shift @ypvals;
342 my ($x1,$y1,$yp1); 342 my ($x1,$y1,$yp1);
343 my @polys; #calculate a hermite polynomial evaluator for each region 343 my @polys; #calculate a hermite polynomial evaluator for each region
344 344
345 while (@xvals) { 345 while (@xvals) {
346 $x1 = shift @xvals; 346 $x1 = shift @xvals;
347 $y1 = shift @yvals; 347 $y1 = shift @yvals;
348 $yp1 = shift @ypvals; 348 $yp1 = shift @ypvals;
349 349
350 push @polys, cubic_hermite($x0, $y0, $yp0 , $x1, $y1 , $yp1); 350 push @polys, cubic_hermite($x0, $y0, $yp0 , $x1, $y1 , $yp1);
351 $x0 = $x1; 351 $x0 = $x1;
352 $y0 = $y1; 352 $y0 = $y1;
353 $yp0 = $yp1; 353 $yp0 = $yp1;
354 } 354 }
355 355
356 356
357 my $hermite_spline_function = sub { 357 my $hermite_spline_function = sub {
358 my $x = shift; 358 my $x = shift;
359 my $y; 359 my $y;
360 my $fun; 360 my $fun;
361 my @xvals = @$xref; 361 my @xvals = @$xref;
362 my @fns = @polys; 362 my @fns = @polys;
363 return &{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint 363 return &{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint
364 364
365 while (@xvals && $x > $xvals[0]) { # find the function for this range of x 365 while (@xvals && $x > $xvals[0]) { # find the function for this range of x
366 shift(@xvals); 366 shift(@xvals);
367 $fun = shift(@fns); 367 $fun = shift(@fns);
368 } 368 }
369 369
370 # now that we have the left hand of the input 370 # now that we have the left hand of the input
371 #check first that x isn't out of range to the left or right 371 #check first that x isn't out of range to the left or right
372 if (@xvals && defined($fun) ) { 372 if (@xvals && defined($fun) ) {
373 $y =&$fun($x); 373 $y =&$fun($x);
374 } 374 }
394 push @polys, cubic_hermite_p($x0, $y0, $yp0 , $x1, $y1 , $yp1); 394 push @polys, cubic_hermite_p($x0, $y0, $yp0 , $x1, $y1 , $yp1);
395 $x0 = $x1; 395 $x0 = $x1;
396 $y0 = $y1; 396 $y0 = $y1;
397 $yp0 = $yp1; 397 $yp0 = $yp1;
398 } 398 }
399 399
400 400
401 my $hermite_spline_function_p = sub { 401 my $hermite_spline_function_p = sub {
402 my $x = shift; 402 my $x = shift;
403 my $y; 403 my $y;
404 my $fun; 404 my $fun;
405 my @xvals = @$xref; 405 my @xvals = @$xref;
406 my @fns = @polys; 406 my @fns = @polys;
407 return $y=&{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint 407 return $y=&{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint
408 408
409 while (@xvals && $x > $xvals[0]) { # find the function for this range of x 409 while (@xvals && $x > $xvals[0]) { # find the function for this range of x
410 shift(@xvals); 410 shift(@xvals);
411 $fun = shift(@fns); 411 $fun = shift(@fns);
412 } 412 }
413 413
414 # now that we have the left hand of the input 414 # now that we have the left hand of the input
415 #check first that x isn't out of range to the left or right 415 #check first that x isn't out of range to the left or right
416 if (@xvals && defined($fun) ) { 416 if (@xvals && defined($fun) ) {
417 $y =&$fun($x); 417 $y =&$fun($x);
418 } 418 }
437 push @polys, cubic_hermite_pp($x0, $y0, $yp0 , $x1, $y1 , $yp1); 437 push @polys, cubic_hermite_pp($x0, $y0, $yp0 , $x1, $y1 , $yp1);
438 $x0 = $x1; 438 $x0 = $x1;
439 $y0 = $y1; 439 $y0 = $y1;
440 $yp0 = $yp1; 440 $yp0 = $yp1;
441 } 441 }
442 442
443 443
444 my $hermite_spline_function_pp = sub { 444 my $hermite_spline_function_pp = sub {
445 my $x = shift; 445 my $x = shift;
446 my $y; 446 my $y;
447 my $fun; 447 my $fun;
448 my @xvals = @$xref; 448 my @xvals = @$xref;
449 my @fns = @polys; 449 my @fns = @polys;
450 return $y=&{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint 450 return $y=&{$fns[0]} ($x) if $x == $xvals[0]; #handle left most endpoint
451 451
452 while (@xvals && $x > $xvals[0]) { # find the function for this range of x 452 while (@xvals && $x > $xvals[0]) { # find the function for this range of x
453 shift(@xvals); 453 shift(@xvals);
454 $fun = shift(@fns); 454 $fun = shift(@fns);
455 } 455 }
456 456
457 # now that we have the left hand of the input 457 # now that we have the left hand of the input
458 #check first that x isn't out of range to the left or right 458 #check first that x isn't out of range to the left or right
459 if (@xvals && defined($fun) ) { 459 if (@xvals && defined($fun) ) {
460 $y =&$fun($x); 460 $y =&$fun($x);
461 } 461 }

Legend:
Removed from v.1050  
changed lines
  Added in v.1079

aubreyja at gmail dot com
ViewVC Help
Powered by ViewVC 1.0.9