[system] / trunk / pg / lib / Hermite.pm Repository: Repository Listing bbplugincoursesdistsnplrochestersystemwww

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

Revision 1079 - (download) (as text) (annotate)
Mon Jun 9 17:36:12 2003 UTC (16 years, 8 months ago) by apizer
File size: 11526 byte(s)
```removed unneccesary shebang lines
```

```    1
2
3 use strict;
4 use Carp;
5
7
8   Hermite.pm
9
11
12       Usage:
13                 \$obj = new Hermit(\@x_values, \y_valuses \@yp_values);
14
15     #get and set methods
16                 \$ra_x_values = \$obj -> ra_x(\@x_values);
17     \$ra_y_values = \$obj -> ra_y;
18     \$ra_yp_values = \$obj -> ra_yp;
19
20     \$obj -> initialize;           # calculates the approximation
21
22     #get methods
23     \$rf_function                  = \$obj -> rf_f;
24     \$rf_function_derivative       = \$obj -> rf_fp;
25     \$rf_function_2nd_derivative   = \$obj -> rf_fpp;
26
27     \$rh_critical_points           =\$obj -> rh_critical_points
28     \$rh_inflection_points         =\$obj -> rh_inflection_points
29
30
31
32
34
35 This module defines an object containing a Hermite spline approximation to a function.
36   The approximation
37 consists of a piecewise cubic polynomial which agrees with the original
38 function and its first derivative at
39 the node points.
40
41 This is useful for creating on the fly graphics.  Care must be taken to use a small
42 number of points spaced reasonably far apart, preferably
43 points with alternating slopes, since this will minimize
44 the number of extraneous critical points
45 introduced.  Too many points will introduce many small variations and
46 a large number of extraneous critical points.
47
48 There are even more extraneous inflections points.
49 This parameter is probably not very useful.   A different
50 approximation scheme needs to be used.
51
52 =cut
53
54
55 # an object that contains a cubic hermite spline
56
57 package Hermite;
58
59
60 sub new {
61     my \$class = shift;
62     my (\$ra_x,\$ra_y,\$ra_yp) = @_;
63     my \$self = {
64            ra_x       =>  [],
65      ra_y       =>  [],
66      ra_yp        =>  [],
67      rf_f       =>  0,   # refers to a subroutine for the function
68      rf_fp        =>  0,   # refers to a subroutine for the derivative of the function
69      rf_fpp       =>  0,   # refers to a subroutine for the second derivative of the function
70      rh_critical_points   =>  {},
71      rh_inflection_points   =>  {},
72      rh_maximum_points    =>  {},
73      rh_minimum_points    =>  {},
74     };
75     bless \$self, \$class;
76     \$self->define(\$ra_x,\$ra_y,\$ra_yp) if ref(\$ra_x) eq 'ARRAY' and ref(\$ra_y) eq 'ARRAY' and ref(\$ra_yp) eq 'ARRAY';
77     \$self;
78 }
79
80 sub define{
81     my \$self   = shift;
82     my \$ra_x  = shift;
83     my \$ra_y  = shift;
84     my \$ra_yp = shift;
85     \$self->ra_x(\$ra_x);
86     \$self->ra_y(\$ra_y);
87     \$self->ra_yp(\$ra_yp);
88
89     \$self -> initialize();
90     (\$self->{x_ref},\$self->{y_ref}, \$self->{yp_ref} )
91 }
92
93 sub initialize {
94     my \$self = shift;
95     # define the functions rf_f
96
97     \$self -> {rf_f} = hermite_spline(\$self -> {ra_x}, \$self -> {ra_y}, \$self -> {ra_yp} );
98 #     # define the function  rf_fp
99      \$self -> {rf_fp} = hermite_spline_p(\$self -> {ra_x}, \$self -> {ra_y}, \$self -> {ra_yp} );
100 #     # define the function  rf_fp
101      \$self -> {rf_fpp} = hermite_spline_pp(\$self -> {ra_x}, \$self -> {ra_y}, \$self -> {ra_yp} );
102 #     # define the critical points
103      %{\$self -> {rh_critical_points}} = critical_points(\$self -> {ra_x}, \$self -> {ra_y}, \$self -> {ra_yp},
104                                                      \$self -> { rf_f },
105                  );
106 #     # define the inflection_points
107      %{\$self -> {rh_inflection_points}} = inflection_points(\$self -> {ra_x}, \$self -> {ra_y}, \$self -> {ra_yp} ,
108                                                      \$self -> { rf_f },
109                  );
110 #     # define the maximum points
111
112     # define the minimum points
113
114 }
115
116 # fix up these accesses to do more checking
117 sub ra_x {
118     my \$self = shift;
119     my \$ra_x  = shift;
120     @{\$self->{ra_x}} = @\$ra_x if ref(\$ra_x) eq 'ARRAY';
121     \$self->{ra_x};
122
123 }
124
125 sub ra_y {
126     my \$self = shift;
127     my \$ra_y  = shift;
128     @{\$self->{ra_y}} = @\$ra_y if ref(\$ra_y) eq 'ARRAY';
129     \$self->{ra_y};
130 }
131
132 sub ra_yp {
133     my \$self = shift;
134     my \$ra_yp  = shift;
135     @{\$self->{ra_yp}} = @\$ra_yp if ref(\$ra_yp) eq 'ARRAY';
136     \$self->{ra_yp};
137 }
138
139 sub rf_f {
140     my \$self = shift;
141     my \$rf_f =shift;
142     \$self ->{rf_f} = \$rf_f if defined( \$rf_f );
143     \$self ->{rf_f};
144 }
145 sub rf_fp {
146     my \$self = shift;
147     my \$rf_fp =shift;
148     \$self ->{rf_fp} = \$rf_fp if defined( \$rf_fp );
149     \$self ->{rf_fp};
150 }
151 sub rf_fpp {
152     my \$self = shift;
153     my \$rf_fpp =shift;
154     \$self ->{rf_fpp} = \$rf_fpp if defined( \$rf_fpp );
155     \$self ->{rf_fpp};
156 }
157
158 sub rh_critical_points {
159     my \$self = shift;
160     \$self -> { rh_critical_points};
161 }
162
163 sub rh_inflection_points {
164     my \$self = shift;
165     \$self -> { rh_inflection_points};
166 }
167
168 ##Internal subroutines
169
170
171 sub critical_points{
172     my \$ra_x = shift;
173     my \$ra_y =shift;
174     my \$ra_yp = shift;
175     my \$rf_hermite_fun =shift;
176     my %critical_points = ();
177     my \$last_index = @\$ra_x -2;
178     foreach my \$i (0 .. \$last_index ) {
179         internal_critical_points(\$ra_x->[\$i],   \$ra_y->[\$i],   \$ra_yp->[\$i],
180                            \$ra_x->[\$i+1], \$ra_y->[\$i+1], \$ra_yp->[\$i+1],
181          \%critical_points,\$rf_hermite_fun);
182     }
183     %critical_points;
184 }
185
186 sub inflection_points{
187     my \$ra_x = shift;
188     my \$ra_y =shift;
189     my \$ra_yp = shift;
190     my \$rf_hermite_fun =shift;
191     my %inflection_points = ();
192     my \$last_index = @\$ra_x -2;
193     foreach my \$i (0 .. \$last_index ) {
194         internal_inflection_points(\$ra_x->[\$i],   \$ra_y->[\$i],   \$ra_yp->[\$i],
195                            \$ra_x->[\$i+1], \$ra_y->[\$i+1], \$ra_yp->[\$i+1],
196          \%inflection_points,\$rf_hermite_fun);
197     }
198     %inflection_points;
199 }
200 sub internal_critical_points{
201     my (\$x0,\$l,\$lp, \$x1,\$r,\$rp, \$rh_roots ,\$rf_function) = @_;
202      #data for one segment of the hermite spline
203
204      # coefficients for the approximating polynomial
205
206      my @a = (   \$l,
207                   \$lp,
208       -\$lp/2 + (3*(-\$lp - 2*(\$l - \$r) - \$rp))/2 + \$rp/2,
209       \$lp + 2*(\$l - \$r) + \$rp
210                );
211
212     my (\$root1, \$root2, \$z1,\$z2);
213     if ( \$a[3] == 0 ) {
214   if ( \$a[2] == 0 ) {
215   } else {
216       \$root1 = -\$a[1]/( 2*\$a[2] );
217       if ( 0 <= \$root1 and \$root1 < 1) {
218     \$z1 = \$root1*(\$x1-\$x0) + \$x0;
219
220           \$rh_roots -> {\$z1} = &\$rf_function(\$z1);
221       }
222   }
223    } else {
224   my \$discriminent = (4*\$a[2]**2 - 12*\$a[1]*\$a[3]);
225
226
227   if ( \$discriminent >= 0 ) {
228       \$discriminent = \$discriminent**0.5;
229       \$root1 = (-2*\$a[2] - \$discriminent )/( 6*\$a[3] );
230       \$root2 = (-2*\$a[2] + \$discriminent )/( 6*\$a[3] );
231       \$z1 = \$root1*(\$x1-\$x0) + \$x0;
232       \$z2 = \$root2*(\$x1-\$x0) + \$x0;
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;
235   }
236     }
237 }
238
239 sub internal_inflection_points {
240     my (\$x0,\$l,\$lp,,\$x1,\$r,\$rp,\$rh_roots,\$rf_function) = @_;
241      #data for one segment of the hermite spline
242
243      # coefficients for the approximating polynomial
244
245      my @a = (   \$l,
246                   \$lp,
247       -\$lp/2 + (3*(-\$lp - 2*(\$l - \$r) - \$rp))/2 + \$rp/2,
248       \$lp + 2*(\$l - \$r) + \$rp
249                );
250     if (\$a[3] == 0 ) {
251     } else {
252   my \$root1 = -\$a[2]/( 3*\$a[3] );
253   my \$z1 = \$root1*(\$x1-\$x0) + \$x0;
254   \$rh_roots -> {\$z1} = &\$rf_function(\$z1) if  0 <= \$root1 and \$root1 < 1;
255     }
256
257 }
258
259
260
261 sub cubic_hermite {
262     my ( \$x0, \$y0,\$yp0, \$x1, \$y1, \$yp1 ) = @_;
263     my @a;
264     my \$width = \$x1 - \$x0;
265     \$yp0 = \$yp0*\$width;  # normalize to unit width
266     \$yp1 = \$yp1*\$width;
267
268     \$a[0] = \$y0;
269     \$a[1] = \$yp0;
270     \$a[2] = -3*\$y0 - 2*\$yp0 +3*\$y1 -\$yp1;
271     \$a[3] = 2*\$y0 + \$yp0 - 2*\$y1 +\$yp1;
272
273     my \$f = sub {
274                         my \$x = shift;
275                         #normalize to unit width
276       \$x = ( \$x - \$x0 )/\$width;
277       ( (\$a[3]*\$x + \$a[2]) * \$x + \$a[1] )*\$x + \$a[0];
278
279                 };
280     \$f;
281 }
282
283 sub cubic_hermite_p {
284     my ( \$x0, \$y0,\$yp0, \$x1, \$y1, \$yp1 )=@_;
285     my @a;
286     my \$width = \$x1 - \$x0;
287     \$yp0 = \$yp0*\$width;  # normalize to unit width
288     \$yp1 = \$yp1*\$width;
289
290     \$a[0] = \$y0;
291     \$a[1] = \$yp0;
292     \$a[2] = -3*\$y0 - 2*\$yp0 +3*\$y1 -\$yp1;
293     \$a[3] = 2*\$y0 + \$yp0 - 2*\$y1 +\$yp1;
294
295     my \$fp = sub {
296                         my \$x = shift;
297                         #normalize to unit width
298       \$x = ( \$x - \$x0 )/\$width;
299       ( (3*\$a[3]*\$x + 2*\$a[2]) * \$x + \$a[1] )/\$width ;
300                 };
301
302
303
304      \$fp;
305
306 }
307
308 sub cubic_hermite_pp {
309     my ( \$x0, \$y0,\$yp0, \$x1, \$y1, \$yp1 ) = @_;
310     my @a;
311     my \$width = \$x1 - \$x0;
312     \$yp0 = \$yp0*\$width;  # normalize to unit width
313     \$yp1 = \$yp1*\$width;
314
315     \$a[0] = \$y0;
316     \$a[1] = \$yp0;
317     \$a[2] = -3*\$y0 - 2*\$yp0 +3*\$y1 -\$yp1;
318     \$a[3] = 2*\$y0 + \$yp0 - 2*\$y1 +\$yp1;
319
320
321
322     my \$fpp = sub {
323                         my \$x = shift;
324                         #normalize to unit width
325       \$x = ( \$x - \$x0 )/\$width;
326        (6*\$a[3]*\$x + 2*\$a[2])/\$width**2 ;
327                 };
328
329     \$fpp;
330 }
331
332
333
334 sub hermite_spline {
335   my (\$xref, \$yref, \$ypref) = @_;
336   my @xvals  = @\$xref;
337   my @yvals  = @\$yref;
338   my @ypvals = @\$ypref;
339   my \$x0 = shift @xvals;
340   my \$y0 = shift @yvals;
341   my \$yp0 = shift @ypvals;
342   my (\$x1,\$y1,\$yp1);
343   my @polys;  #calculate a hermite polynomial evaluator for each region
344
345   while (@xvals) {
346     \$x1 = shift @xvals;
347     \$y1 = shift @yvals;
348     \$yp1 = shift @ypvals;
349
350     push @polys, cubic_hermite(\$x0, \$y0, \$yp0 , \$x1, \$y1 , \$yp1);
351     \$x0  = \$x1;
352     \$y0  = \$y1;
353     \$yp0 = \$yp1;
354   }
355
356
357   my \$hermite_spline_function = sub {
358     my \$x = shift;
359     my \$y;
360     my \$fun;
361     my @xvals = @\$xref;
362     my @fns = @polys;
363     return &{\$fns[0]} (\$x) if \$x == \$xvals[0]; #handle left most endpoint
364
365     while (@xvals && \$x > \$xvals[0]) {  # find the function for this range of x
366       shift(@xvals);
367       \$fun = shift(@fns);
368     }
369
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
372     if (@xvals  && defined(\$fun) )  {
373       \$y =&\$fun(\$x);
374     }
375     \$y;
376   };
377   \$hermite_spline_function;
378 }
379
380 sub hermite_spline_p {
381   my (\$xref, \$yref, \$ypref) = @_;
382   my @xvals  = @\$xref;
383   my @yvals  = @\$yref;
384   my @ypvals = @\$ypref;
385   my \$x0 = shift @xvals;
386   my \$y0 = shift @yvals;
387   my \$yp0 = shift @ypvals;
388   my (\$x1,\$y1,\$yp1);
389   my @polys;  #calculate a hermite polynomial evaluator for each region
390   while (@xvals) {
391     \$x1 = shift @xvals;
392     \$y1 = shift @yvals;
393     \$yp1 = shift @ypvals;
394     push @polys, cubic_hermite_p(\$x0, \$y0, \$yp0 , \$x1, \$y1 , \$yp1);
395     \$x0  = \$x1;
396     \$y0  = \$y1;
397     \$yp0 = \$yp1;
398   }
399
400
401   my \$hermite_spline_function_p = sub {
402     my \$x = shift;
403     my \$y;
404     my \$fun;
405     my @xvals = @\$xref;
406     my @fns = @polys;
407     return \$y=&{\$fns[0]} (\$x) if \$x == \$xvals[0]; #handle left most endpoint
408
409     while (@xvals && \$x > \$xvals[0]) {  # find the function for this range of x
410       shift(@xvals);
411       \$fun = shift(@fns);
412     }
413
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
416     if (@xvals  && defined(\$fun) )  {
417       \$y =&\$fun(\$x);
418     }
419     \$y;
420   };
421   \$hermite_spline_function_p;
422 }
423 sub hermite_spline_pp {
424   my (\$xref, \$yref, \$ypref) = @_;
425   my @xvals  = @\$xref;
426   my @yvals  = @\$yref;
427   my @ypvals = @\$ypref;
428   my \$x0 = shift @xvals;
429   my \$y0 = shift @yvals;
430   my \$yp0 = shift @ypvals;
431   my (\$x1,\$y1,\$yp1);
432   my @polys;  #calculate a hermite polynomial evaluator for each region
433   while (@xvals) {
434     \$x1 = shift @xvals;
435     \$y1 = shift @yvals;
436     \$yp1 = shift @ypvals;
437     push @polys, cubic_hermite_pp(\$x0, \$y0, \$yp0 , \$x1, \$y1 , \$yp1);
438     \$x0  = \$x1;
439     \$y0  = \$y1;
440     \$yp0 = \$yp1;
441   }
442
443
444   my \$hermite_spline_function_pp = sub {
445     my \$x = shift;
446     my \$y;
447     my \$fun;
448     my @xvals = @\$xref;
449     my @fns = @polys;
450     return \$y=&{\$fns[0]} (\$x) if \$x == \$xvals[0]; #handle left most endpoint
451
452     while (@xvals && \$x > \$xvals[0]) {  # find the function for this range of x
453       shift(@xvals);
454       \$fun = shift(@fns);
455     }
456
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
459     if (@xvals  && defined(\$fun) )  {
460       \$y =&\$fun(\$x);
461     }
462     \$y;
463   };
464   \$hermite_spline_function_pp;
465 }
466
467
468 1;
```

 aubreyja at gmail dot com ViewVC Help Powered by ViewVC 1.0.9