PLplot  5.10.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
dsplint.c
Go to the documentation of this file.
1 //
2 // Copyright (C) 2009 Alan W. Irwin
3 //
4 // This file is part of PLplot.
5 //
6 // PLplot is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU Library General Public License as published
8 // by the Free Software Foundation; either version 2 of the License, or
9 // (at your option) any later version.
10 //
11 // PLplot is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU Library General Public License for more details.
15 //
16 // You should have received a copy of the GNU Library General Public License
17 // along with PLplot; if not, write to the Free Software
18 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
19 //
20 // Provenance: This code was originally developed under the GPL as part of
21 // the FreeEOS project (revision 121). This code has been converted from
22 // Fortran to C with the aid of f2c and relicensed for PLplot under the LGPL
23 // with the permission of the FreeEOS copyright holder (Alan W. Irwin).
24 //
25 
26 #include "dsplint.h"
27 
28 # define MAX( a, b ) ( ( ( a ) > ( b ) ) ? ( a ) : ( b ) )
29 # define MIN( a, b ) ( ( ( a ) < ( b ) ) ? ( a ) : ( b ) )
30 
31 //int dsplint(double *xa, double *ya, double *y2a,
32 // int n, double x, double *y, double *dy, double *d2y)
33 int dsplint( double *xa, double *ya, double *y2a,
34  int n, double x, double *y )
35 {
36  // Initialized data
37 
38  static int nsave = 0, khi, klo;
39 
40  int i__1, i__2, k;
41  double a, b, h__;
42 
43 // evaluate spline = y and its derivatives dy and d2y at x given
44 // xa, ya, y2a from dspline.
45  // Parameter adjustments
46  --y2a;
47  --ya;
48  --xa;
49 
50  // Function Body
51  if ( n != nsave )
52  {
53 // if call with different n value, then redo range
54  nsave = n;
55  klo = 1;
56  khi = n;
57  if ( xa[klo] > x )
58  {
59  return 1;
60  }
61  if ( xa[khi] < x )
62  {
63  return 2;
64  }
65  }
66  else
67  {
68 // optimize range assuming continuous (ascending or
69 // descending x calls.
70  if ( xa[klo] > x )
71  {
72 // x is descending so try next range.
73  khi = MAX( 2, klo );
74  klo = khi - 1;
75 // if x smaller than next range try lower limit.
76  if ( xa[klo] > x )
77  {
78  klo = 1;
79  }
80  if ( xa[klo] > x )
81  {
82  return 1;
83  }
84  }
85  else if ( xa[khi] <= x )
86  {
87 // x is ascending so try next range.
88 // Computing MIN
89  i__1 = khi, i__2 = n - 1;
90  klo = MIN( i__1, i__2 );
91  khi = klo + 1;
92 // if x larger than next range try upper limit.
93  if ( xa[khi] <= x )
94  {
95  khi = n;
96  }
97  if ( xa[khi] < x )
98  {
99  return 2;
100  }
101  }
102  }
103  while ( khi - klo > 1 )
104  {
105  k = ( khi + klo ) / 2;
106  if ( xa[k] > x )
107  {
108  khi = k;
109  }
110  else
111  {
112  klo = k;
113  }
114  }
115  h__ = xa[khi] - xa[klo];
116  if ( h__ <= 0. )
117  {
118  return 3;
119  }
120  a = ( xa[khi] - x ) / h__;
121  b = ( x - xa[klo] ) / h__;
122  *y = a * ya[klo] + b * ya[khi] + ( a * ( a * a - 1. ) * y2a[klo] + b * ( b *
123  b - 1. ) * y2a[khi] ) * ( h__ * h__ ) / 6.;
124 // *dy = (-ya[klo] + ya[khi] + (-(a * 3. * a - 1.) * y2a[klo] + (b * 3. * b
125 // - 1.) * y2a[khi]) * (h__ * h__) / 6.) / h__;
126 //d2y = a * y2a[klo] + b * y2a[khi];
127  return 0;
128 }
129