← Back to branch summary

~uhh-ssd/+junk/humidity_readout

« back to all changes in this revision

Viewing changes to plplot/plplot-5.9.9/lib/qsastime/dsplint.c

  • Committer: Joachim Erfle
  • Date: 2013-07-24 13:53:41 UTC
  • Revision ID: joachim.erfle@desy.de-20130724135341-1qojpp701zsn009p
initial commit

Show diffs side-by-side

added added

removed removed

Lines of Context:
plplot/plplot-5.9.9/lib/qsastime/dsplint.c
 
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
# define MAX( a, b )    ( ( ( a ) > ( b ) ) ? ( a ) : ( b ) )
 
26
# define MIN( a, b )    ( ( ( a ) < ( b ) ) ? ( a ) : ( b ) )
 
27
 
 
28
//int dsplint(double *xa, double *ya, double *y2a,
 
29
//          int n, double x, double *y, double *dy, double *d2y)
 
30
int dsplint( double *xa, double *ya, double *y2a,
 
31
             int n, double x, double *y )
 
32
{
 
33
    // Initialized data
 
34
 
 
35
    static int nsave = 0, khi, klo;
 
36
 
 
37
    int        i__1, i__2, k;
 
38
    double     a, b, h__;
 
39
 
 
40
//      evaluate spline = y and its derivatives dy and d2y at x given
 
41
//      xa, ya, y2a from dspline.
 
42
    // Parameter adjustments
 
43
    --y2a;
 
44
    --ya;
 
45
    --xa;
 
46
 
 
47
    // Function Body
 
48
    if ( n != nsave )
 
49
    {
 
50
//        if call with different n value, then redo range
 
51
        nsave = n;
 
52
        klo   = 1;
 
53
        khi   = n;
 
54
        if ( xa[klo] > x )
 
55
        {
 
56
            return 1;
 
57
        }
 
58
        if ( xa[khi] < x )
 
59
        {
 
60
            return 2;
 
61
        }
 
62
    }
 
63
    else
 
64
    {
 
65
//        optimize range assuming continuous (ascending or
 
66
//        descending x calls.
 
67
        if ( xa[klo] > x )
 
68
        {
 
69
//          x is descending so try next range.
 
70
            khi = MAX( 2, klo );
 
71
            klo = khi - 1;
 
72
//          if x smaller than next range try lower limit.
 
73
            if ( xa[klo] > x )
 
74
            {
 
75
                klo = 1;
 
76
            }
 
77
            if ( xa[klo] > x )
 
78
            {
 
79
                return 1;
 
80
            }
 
81
        }
 
82
        else if ( xa[khi] <= x )
 
83
        {
 
84
//          x is ascending so try next range.
 
85
// Computing MIN
 
86
            i__1 = khi, i__2 = n - 1;
 
87
            klo  = MIN( i__1, i__2 );
 
88
            khi  = klo + 1;
 
89
//          if x larger than next range try upper limit.
 
90
            if ( xa[khi] <= x )
 
91
            {
 
92
                khi = n;
 
93
            }
 
94
            if ( xa[khi] < x )
 
95
            {
 
96
                return 2;
 
97
            }
 
98
        }
 
99
    }
 
100
    while ( khi - klo > 1 )
 
101
    {
 
102
        k = ( khi + klo ) / 2;
 
103
        if ( xa[k] > x )
 
104
        {
 
105
            khi = k;
 
106
        }
 
107
        else
 
108
        {
 
109
            klo = k;
 
110
        }
 
111
    }
 
112
    h__ = xa[khi] - xa[klo];
 
113
    if ( h__ <= 0. )
 
114
    {
 
115
        return 3;
 
116
    }
 
117
    a  = ( xa[khi] - x ) / h__;
 
118
    b  = ( x - xa[klo] ) / h__;
 
119
    *y = a * ya[klo] + b * ya[khi] + ( a * ( a * a - 1. ) * y2a[klo] + b * ( b *
 
120
                                                                             b - 1. ) * y2a[khi] ) * ( h__ * h__ ) / 6.;
 
121
//    *dy = (-ya[klo] + ya[khi] + (-(a * 3. * a - 1.) * y2a[klo] + (b * 3. * b
 
122
//          - 1.) * y2a[khi]) * (h__ * h__) / 6.) / h__;
 
123
//d2y = a * y2a[klo] + b * y2a[khi];
 
124
    return 0;
 
125
}
 
126