PLplot 5.15.0
Loading...
Searching...
No Matches
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)
33int 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
#define MIN(a, b)
Definition: dsplint.c:29
int dsplint(double *xa, double *ya, double *y2a, int n, double x, double *y)
Definition: dsplint.c:33
#define MAX(a, b)
Definition: dsplint.c:28