PLplot  5.10.0
dsplint.c
Go to the documentation of this file.
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Defines