PLplot
5.10.0
|
00001 // 00002 // Copyright (C) 2009 Alan W. Irwin 00003 // 00004 // This file is part of PLplot. 00005 // 00006 // PLplot is free software; you can redistribute it and/or modify 00007 // it under the terms of the GNU Library General Public License as published 00008 // by the Free Software Foundation; either version 2 of the License, or 00009 // (at your option) any later version. 00010 // 00011 // PLplot is distributed in the hope that it will be useful, 00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 // GNU Library General Public License for more details. 00015 // 00016 // You should have received a copy of the GNU Library General Public License 00017 // along with PLplot; if not, write to the Free Software 00018 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 00019 // 00020 // Provenance: This code was originally developed under the GPL as part of 00021 // the FreeEOS project (revision 121). This code has been converted from 00022 // Fortran to C with the aid of f2c and relicensed for PLplot under the LGPL 00023 // with the permission of the FreeEOS copyright holder (Alan W. Irwin). 00024 // 00025 00026 #include "dsplint.h" 00027 00028 # define MAX( a, b ) ( ( ( a ) > ( b ) ) ? ( a ) : ( b ) ) 00029 # define MIN( a, b ) ( ( ( a ) < ( b ) ) ? ( a ) : ( b ) ) 00030 00031 //int dsplint(double *xa, double *ya, double *y2a, 00032 // int n, double x, double *y, double *dy, double *d2y) 00033 int dsplint( double *xa, double *ya, double *y2a, 00034 int n, double x, double *y ) 00035 { 00036 // Initialized data 00037 00038 static int nsave = 0, khi, klo; 00039 00040 int i__1, i__2, k; 00041 double a, b, h__; 00042 00043 // evaluate spline = y and its derivatives dy and d2y at x given 00044 // xa, ya, y2a from dspline. 00045 // Parameter adjustments 00046 --y2a; 00047 --ya; 00048 --xa; 00049 00050 // Function Body 00051 if ( n != nsave ) 00052 { 00053 // if call with different n value, then redo range 00054 nsave = n; 00055 klo = 1; 00056 khi = n; 00057 if ( xa[klo] > x ) 00058 { 00059 return 1; 00060 } 00061 if ( xa[khi] < x ) 00062 { 00063 return 2; 00064 } 00065 } 00066 else 00067 { 00068 // optimize range assuming continuous (ascending or 00069 // descending x calls. 00070 if ( xa[klo] > x ) 00071 { 00072 // x is descending so try next range. 00073 khi = MAX( 2, klo ); 00074 klo = khi - 1; 00075 // if x smaller than next range try lower limit. 00076 if ( xa[klo] > x ) 00077 { 00078 klo = 1; 00079 } 00080 if ( xa[klo] > x ) 00081 { 00082 return 1; 00083 } 00084 } 00085 else if ( xa[khi] <= x ) 00086 { 00087 // x is ascending so try next range. 00088 // Computing MIN 00089 i__1 = khi, i__2 = n - 1; 00090 klo = MIN( i__1, i__2 ); 00091 khi = klo + 1; 00092 // if x larger than next range try upper limit. 00093 if ( xa[khi] <= x ) 00094 { 00095 khi = n; 00096 } 00097 if ( xa[khi] < x ) 00098 { 00099 return 2; 00100 } 00101 } 00102 } 00103 while ( khi - klo > 1 ) 00104 { 00105 k = ( khi + klo ) / 2; 00106 if ( xa[k] > x ) 00107 { 00108 khi = k; 00109 } 00110 else 00111 { 00112 klo = k; 00113 } 00114 } 00115 h__ = xa[khi] - xa[klo]; 00116 if ( h__ <= 0. ) 00117 { 00118 return 3; 00119 } 00120 a = ( xa[khi] - x ) / h__; 00121 b = ( x - xa[klo] ) / h__; 00122 *y = a * ya[klo] + b * ya[khi] + ( a * ( a * a - 1. ) * y2a[klo] + b * ( b * 00123 b - 1. ) * y2a[khi] ) * ( h__ * h__ ) / 6.; 00124 // *dy = (-ya[klo] + ya[khi] + (-(a * 3. * a - 1.) * y2a[klo] + (b * 3. * b 00125 // - 1.) * y2a[khi]) * (h__ * h__) / 6.) / h__; 00126 //d2y = a * y2a[klo] + b * y2a[khi]; 00127 return 0; 00128 } 00129