← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to stdred/spec/libsrc/mutil.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
stdred/spec/libsrc/mutil.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1991-2009 European Southern Observatory (ESO)
 
3
 
 
4
  This program is free software; you can redistribute it and/or 
 
5
  modify it under the terms of the GNU General Public License as 
 
6
  published by the Free Software Foundation; either version 2 of 
 
7
  the License, or (at your option) any later version.
 
8
 
 
9
  This program is distributed in the hope that it will be useful,
 
10
  but WITHOUT ANY WARRANTY; without even the implied warranty of
 
11
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
12
  GNU General Public License for more details.
 
13
 
 
14
  You should have received a copy of the GNU General Public 
 
15
  License along with this program; if not, write to the Free 
 
16
  Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
 
17
  MA 02139, USA.
 
18
 
 
19
  Correspondence concerning ESO-MIDAS should be addressed as follows:
 
20
        Internet e-mail: midas@eso.org
 
21
        Postal address: European Southern Observatory
 
22
                        Data Management Division 
 
23
                        Karl-Schwarzschild-Strasse 2
 
24
                        D 85748 Garching bei Muenchen 
 
25
                        GERMANY
 
26
===========================================================================*/
 
27
 
 
28
 
 
29
/*+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
.IDENT        mutil.c
 
31
.MODULE       subroutines 
 
32
.ENVIRONMENT  UNIX
 
33
.LANGUAGE     C
 
34
.AUTHOR       Cristian Levin - ESO La Silla
 
35
.PURPOSE      Miscelanous mathematical algorithms:
 
36
              - linear square fitting.
 
37
              - fast polynomial evaluation.
 
38
              - search of the median.
 
39
              - polynomial interpolation.
 
40
.KEYWORDS     interpolation, fitting, polynomial evaluation.
 
41
.COMMENTS     Most of this routines were taken as they are from
 
42
              the book "Numerical Recipes in C" -- 1st edition.
 
43
.VERSION 1.0  1-Apr-1991   Implementation
 
44
 
 
45
 090723         last modif
 
46
------------------------------------------------------------*/
 
47
 
 
48
#include <stdio.h>
 
49
#include <math.h>
 
50
 
 
51
/******************************************************
 
52
 *
 
53
 * lfit(): general linear least squares solution.
 
54
 *
 
55
 ******************************************************/
 
56
void lfit( x, y, ndata, a, mfit, funcs )
 
57
int ndata, mfit;        /* IN: no. of values, no. of coefficients */
 
58
double x[], y[], a[];   /* IN: (Xi,Yi)   OUT: (Ai) */
 
59
void (*funcs)();
 
60
{
 
61
    int k, j, i;
 
62
    double **covar, **beta, *afunc;
 
63
    void spec_gaussj(), free_dmatrix(), free_dvector();
 
64
    double **dmatrix(), *dvector();
 
65
 
 
66
    beta  = dmatrix( 1, mfit, 1, 1 );
 
67
    covar = dmatrix( 1, mfit, 1, mfit );
 
68
    afunc = dvector( 1, mfit );
 
69
 
 
70
    for ( j = 1; j <= mfit; j++ ) {
 
71
        for ( k = 1; k <= mfit; k++ )
 
72
             covar[j][k] = 0.0;
 
73
        beta[j][1] = 0.0;
 
74
    }
 
75
 
 
76
    for ( i = 1; i <= ndata; i++ ) {
 
77
        (*funcs)( x[i], afunc, mfit );
 
78
        for ( j = 1; j <= mfit; j++ ) {
 
79
            for ( k = 1; k <= j; k++ )
 
80
                covar[j][k] += afunc[j] * afunc[k];
 
81
            beta[j][1] += y[i] * afunc[j];
 
82
        }
 
83
    }
 
84
 
 
85
    for ( j = 2; j <= mfit; j++ )
 
86
        for ( k = 1; k <= j-1; k++ )
 
87
            covar[k][j] = covar[j][k];
 
88
 
 
89
    spec_gaussj( covar, mfit, beta, 1 );        /* Matrix solution */
 
90
    
 
91
    for ( j = 1; j <= mfit; j++ )       /* solution to coefficients a */
 
92
        a[j] = beta[j][1];
 
93
 
 
94
    free_dvector( afunc, 1, mfit );
 
95
    free_dmatrix( beta, 1, mfit, 1, 1 );
 
96
    free_dmatrix( covar, 1, mfit, 1, mfit );
 
97
}
 
98
 
 
99
 
 
100
/******************************************************************
 
101
 *
 
102
 * spec_gaussj(): linear equation solution by Gauss-Jordan elimination
 
103
 *
 
104
 ******************************************************************/
 
105
 
 
106
#define SWAP(a, b)      {float temp = (a); (a) = (b); (b) = temp; }
 
107
 
 
108
void spec_gaussj( a, n, b, m )
 
109
double **a, /* IN    : a[1..n][1..n] */
 
110
       **b; /* IN/OUT: b[1..n][1..m] */
 
111
int n, m;
 
112
{
 
113
    int *indxc, *indxr, *ipiv; /* for bookeeping in the pivoting */
 
114
    int i, icol, irow, j, k;
 
115
    int *ivector();
 
116
    double big, dum, pivinv;
 
117
    void nrerror(), free_ivector();
 
118
 
 
119
 
 
120
    irow = icol = 0;
 
121
    indxc = ivector( 1, n );
 
122
    indxr = ivector( 1, n );
 
123
    ipiv  = ivector( 1, n );
 
124
 
 
125
    for ( i = 1; i <= n; i++ )
 
126
        ipiv[i] = 0.0;
 
127
 
 
128
    for ( i = 1; i <= n; i++ ) { /* main loop over columns to be reduced */
 
129
        big = 0.0;
 
130
        for ( j = 1; j <= n; j++ )
 
131
            if ( ipiv[j] != 1 )
 
132
                for ( k = 1; k <= n; k++ ) {
 
133
                    if ( ipiv[k] == 0 ) {
 
134
                        if ( fabs(a[j][k]) >= big ) {
 
135
                            big = fabs(a[j][k]);
 
136
                            irow = j;
 
137
                            icol = k;
 
138
                        }
 
139
                    } 
 
140
                    else if ( ipiv[k] > 1 ) {
 
141
                        nrerror( "Tolerance too small for fitting" );
 
142
                        return;
 
143
                    }
 
144
                }
 
145
        ++(ipiv[icol]);
 
146
 
 
147
        if ( irow != icol ) {
 
148
            for ( j = 1; j <= n; j++ ) SWAP( a[irow][j], a[icol][j] );
 
149
            for ( j = 1; j <= m; j++ ) SWAP( b[irow][j], b[icol][j] );
 
150
        }
 
151
 
 
152
        indxr[i] = irow;
 
153
        indxc[i] = icol;
 
154
        if ( a[icol][icol] == 0.0 ) {
 
155
            nrerror( "Tolerance too small for fitting" );
 
156
            return;
 
157
        }
 
158
 
 
159
        pivinv = 1.0 / a[icol][icol];
 
160
        a[icol][icol] = 1.0;
 
161
        for ( j = 1; j <= n; j++ ) a[icol][j] *= pivinv;
 
162
        for ( j = 1; j <= m; j++ ) b[icol][j] *= pivinv;
 
163
 
 
164
        for ( j = 1; j <= n; j++ )
 
165
            if ( j != icol ) {
 
166
                dum = a[j][icol];
 
167
                a[j][icol] = 0.0;
 
168
                for ( k = 1; k <= n; k++ ) a[j][k] -= a[icol][k] * dum;
 
169
                for ( k = 1; k <= m; k++ ) b[j][k] -= b[icol][k] * dum;
 
170
            }
 
171
    } /* main loop */
 
172
 
 
173
    for ( i = n; i >= 1; i-- )
 
174
        if ( indxr[i] != indxc[i] )
 
175
            for ( k = 1; k <= n; k++ )
 
176
                SWAP( a[k][indxr[i]], a[k][indxc[i]] );
 
177
 
 
178
    free_ivector( ipiv, 1, n );
 
179
    free_ivector( indxc, 1, n );
 
180
    free_ivector( indxr, 1, n );
 
181
}
 
182
 
 
183
/********************************************************************
 
184
 *
 
185
 * polint(): polynomial interpolation or extrapolation from N input
 
186
 *           points, using the Neville's algorithm.
 
187
 *
 
188
 ********************************************************************/
 
189
 
 
190
void polint( xa, ya, n, x, y, dy )
 
191
float xa[], ya[], x, *y, *dy;
 
192
int n;
 
193
{
 
194
    int i, m, ns = 1;
 
195
    float den, dif, dift, ho, hp, w;
 
196
    float *c, *d, *fvector();
 
197
    void nrerror(), free_vector();
 
198
 
 
199
    dif = fabs(x - xa[1]);
 
200
    c = fvector(1, n);
 
201
    d = fvector(1, n);
 
202
    for ( i = 1; i <= n; i++ ) {        /* find the index ns of the closest */
 
203
        if ( (dift = fabs(x - xa[i])) < dif ) { /* table entry. */
 
204
            ns = i;
 
205
            dif = dift;
 
206
        }
 
207
        c[i] = ya[i];   /* initialize the tableau of c's and d's */
 
208
        d[i] = ya[i];
 
209
    }
 
210
    *y = ya[ns--];      /* initial aprox. of y */
 
211
    for ( m = 1; m < n; m++ ) {
 
212
        for ( i = 1; i <= n-m; i++ ) {
 
213
            ho = xa[i] - x;
 
214
            hp = xa[i+m] - x;
 
215
            w = c[i+1] - d[i];
 
216
            if ( (den = ho - hp) == 0.0 ) {
 
217
                nrerror( "Error in polynomial interpolation" );
 
218
                return;
 
219
            }
 
220
            den = w / den;
 
221
            d[i] = hp * den;
 
222
            c[i] = ho * den;
 
223
        }
 
224
        *y += (*dy = (2*ns < (n-m) ? c[ns+1] : d[ns--]));
 
225
    }
 
226
    fvector(c, 1, n);
 
227
    fvector(d, 1, n);
 
228
}
 
229
 
 
230
/************************************************************
 
231
 *
 
232
 * eval_fpoly(), eval_dpoly(): fast polynomial evaluation.
 
233
 *
 
234
 ************************************************************/
 
235
 
 
236
float eval_fpoly(x, a, n)
 
237
float x;
 
238
float *a;       /* coefficients (range: 1..n) */
 
239
int n;
 
240
{
 
241
    int i;
 
242
    float y = 0.0;
 
243
 
 
244
    for ( i = n; i >= 1; i-- )
 
245
        y = x*y + a[i];
 
246
 
 
247
    return(y);
 
248
}
 
249
 
 
250
double eval_dpoly(x, a, n)
 
251
double x;
 
252
double *a;      /* coefficients (range: 1..n) */
 
253
int n;
 
254
{
 
255
    int i;
 
256
    double y = 0.0;
 
257
 
 
258
    for ( i = n; i >= 1; i-- )
 
259
        y = x*y + a[i];
 
260
 
 
261
    return(y);
 
262
}
 
263
 
 
264
/************************************************************
 
265
 *
 
266
 * median(): efficient search of the median.
 
267
 *
 
268
 ************************************************************/
 
269
 
 
270
float median(x, n)
 
271
float x[];
 
272
int n;
 
273
{
 
274
    int n2 ,n2p;
 
275
    void piksrt();
 
276
    
 
277
    piksrt(n, x);
 
278
    n2p = (n2 = n/2) + 1;
 
279
    return((n % 2) ? x[n2p] : 0.5 * (x[n2] + x[n2p]));
 
280
}
 
281
 
 
282
/************************************************************
 
283
 *
 
284
 * piksrt(): sorting by straight insertion.
 
285
 *
 
286
 ************************************************************/
 
287
 
 
288
void piksrt(n,arr)
 
289
int n;
 
290
float arr[];
 
291
{
 
292
        int i,j;
 
293
        float a;
 
294
 
 
295
        for (j=2;j<=n;j++) {
 
296
                a=arr[j];
 
297
                i=j-1;
 
298
                while (i > 0 && arr[i] > a) {
 
299
                        arr[i+1]=arr[i];
 
300
                        i--;
 
301
                }
 
302
                arr[i+1]=a;
 
303
        }
 
304
}