← Back to branch summary

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

« back to all changes in this revision

Viewing changes to libsrc/math/mutil.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2015-06-10 14:20:37 UTC
  • mfrom: (1.2.1) (6.1.9 experimental)
  • Revision ID: package-import@ubuntu.com-20150610142037-6iowpbtyjrpou36o
Tags: 15.02pl1.3-1
* New upstream version
* Add CI tests
* Move back to unstable

Show diffs side-by-side

added added

removed removed

Lines of Context:
libsrc/math/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
.VERSION 1.0  1-Apr-1991   Implementation
 
42
 
 
43
 090723         last modif
 
44
------------------------------------------------------------*/
 
45
 
 
46
#include "cpl_matrix.h"
 
47
#include "mutil.h"
 
48
#include "nrutil.h"
 
49
#include <stdio.h>
 
50
#include <math.h>
 
51
 
 
52
/******************************************************
 
53
 *
 
54
 * lfit(): general linear least squares solution.
 
55
 *
 
56
 ******************************************************/
 
57
void lfit( x, y, ndata, a, mfit, funcs )
 
58
int ndata, mfit;        /* IN: no. of values, no. of coefficients */
 
59
double x[], y[], a[];   /* IN: (Xi,Yi)   OUT: (Ai) */
 
60
void (*funcs)();
 
61
{
 
62
    lsqfit_nr(x, y, NULL, ndata, a, mfit, funcs);
 
63
}
 
64
 
 
65
void fit_poly (float inimage[], float outimage[], double start,
 
66
          double step, int npix, int fit_deg)
 
67
{
 
68
  float xout;
 
69
  double *af;
 
70
  double *ad, *xin, *yin;
 
71
  int i;
 
72
 
 
73
  xin = dvector (0, npix - 1);
 
74
  yin = dvector (0, npix - 1);
 
75
 
 
76
  for (i = 0; i < npix; i++)
 
77
    {
 
78
      xin[i] = start + i * step;
 
79
      yin[i] = (double) inimage[i];
 
80
    }
 
81
 
 
82
  ad = dvector (1, fit_deg);
 
83
  af = dvector (1, fit_deg);
 
84
 
 
85
  lfit (xin, yin, npix, ad, fit_deg, dpoly);
 
86
  for (i = 1; i <= fit_deg; i++)
 
87
    {
 
88
      af[i] = ad[i];
 
89
    }
 
90
 
 
91
  for (i = 0; i < npix; i++)
 
92
    {
 
93
      xout = start + step * i;
 
94
      outimage[i] = (float) eval_dpoly (xout, ad, fit_deg);
 
95
    }
 
96
 
 
97
  free_dvector(xin, 0, npix - 1), free_dvector(yin, 0, npix - 1);
 
98
  free_dvector(ad, 1, fit_deg),  free_dvector(af, 1, fit_deg);
 
99
 
 
100
}
 
101
 
 
102
 
 
103
 
 
104
void dpoly(double x, double p[], int np)
 
105
{
 
106
  int j;
 
107
 
 
108
  p[1] = 1.0;
 
109
  for (j = 2; j <= np; j++)
 
110
    p[j] = p[j - 1] * x;
 
111
}
 
112
 
 
113
void dleg(double x, double pl[], int nl)
 
114
{
 
115
  int j;
 
116
  double twox, f2, f1, d;
 
117
 
 
118
  pl[1] = 1.0;
 
119
  pl[2] = x;
 
120
  if (nl > 2)
 
121
    {
 
122
      twox = 2.0 * x;
 
123
      f2 = x;
 
124
      d = 1.0;
 
125
      for (j = 3; j <= nl; j++)
 
126
        {
 
127
          f1 = d;
 
128
          d += 1.0;
 
129
          f2 += twox;
 
130
          pl[j] = (f2 * pl[j - 1] - f1 * pl[j - 2]) / d;
 
131
        }
 
132
    }
 
133
}
 
134
 
 
135
void dcheb(double x, double pl[], int nl)
 
136
{
 
137
  int j;
 
138
  double twox;
 
139
 
 
140
  pl[1] = 1.0;
 
141
  pl[2] = x;
 
142
  if (nl > 2)
 
143
    {
 
144
      twox = 2.0 * x;
 
145
      for (j = 3; j <= nl; j++)
 
146
          pl[j] = twox * pl[j - 1] -  pl[j - 2];
 
147
    }
 
148
}
 
149
 
 
150
/************************************************************
 
151
 *
 
152
 * eval_fpoly(), eval_dpoly(): fast polynomial evaluation.
 
153
 *
 
154
 ************************************************************/
 
155
 
 
156
float eval_fpoly(x, a, n)
 
157
float x;
 
158
float *a;       /* coefficients (range: 1..n) */
 
159
int n;
 
160
{
 
161
    int i;
 
162
    float y = 0.0;
 
163
 
 
164
    for ( i = n; i >= 1; i-- )
 
165
        y = x*y + a[i];
 
166
 
 
167
    return(y);
 
168
}
 
169
 
 
170
double eval_dpoly(x, a, n)
 
171
double x;
 
172
double *a;      /* coefficients (range: 1..n) */
 
173
int n;
 
174
{
 
175
    int i;
 
176
    double y = 0.0;
 
177
 
 
178
    for ( i = n; i >= 1; i-- )
 
179
        y = x*y + a[i];
 
180
 
 
181
    return(y);
 
182
}
 
183
 
 
184
#ifdef __STDC__
 
185
double select_pos (unsigned long k, unsigned long n, double arr[])
 
186
#else
 
187
double select_pos (k, n, arr)
 
188
     unsigned long k, n;
 
189
     double arr[];
 
190
#endif
 
191
 
 
192
/*      select the kth smallest value in arr[1...n], NR, pg. 342 */
 
193
 
 
194
#define SWAPD(a, b)     {double temp = (a); (a) = (b); (b) = temp; }
 
195
 
 
196
{
 
197
 
 
198
  unsigned long i, ir, j, l, mid;
 
199
/*  double a, temp; */
 
200
  double a; 
 
201
 
 
202
  l = 1;
 
203
  ir = n;
 
204
  for (;;)
 
205
    {
 
206
      if (ir <= l + 1)
 
207
        {
 
208
          if (ir == l + 1 && arr[ir] < arr[l])
 
209
            {
 
210
              SWAPD (arr[l], arr[ir])
 
211
            }
 
212
          return arr[k];
 
213
        }
 
214
      else
 
215
        {
 
216
          mid = (l + ir) >> 1;
 
217
          SWAPD (arr[mid], arr[l + 1])
 
218
            if (arr[l + 1] > arr[ir])
 
219
            {
 
220
              SWAPD (arr[l + 1], arr[ir])
 
221
            }
 
222
          if (arr[l] > arr[ir])
 
223
            {
 
224
              SWAPD (arr[l], arr[ir])
 
225
            }
 
226
          if (arr[l + 1] > arr[l])
 
227
            {
 
228
              SWAPD (arr[l + 1], arr[l])
 
229
            }
 
230
          i = l + 1;
 
231
          j = ir;
 
232
          a = arr[l];
 
233
          for (;;)
 
234
            {
 
235
              do
 
236
                i++;
 
237
              while (arr[i] < a);
 
238
              do
 
239
                j--;
 
240
              while (arr[j] > a);
 
241
              if (j < i)
 
242
                break;
 
243
              SWAPD (arr[i], arr[j])
 
244
            }
 
245
          arr[l] = arr[j];
 
246
          arr[j] = a;
 
247
          if (j >= k)
 
248
            ir = j - 1;
 
249
          if (j <= k)
 
250
            l = i;
 
251
        }
 
252
    }
 
253
}
 
254
#undef SWAPD
 
255
#undef SWAP
 
256
 
 
257
#ifdef __STDC__
 
258
void heap_copy(int n, float arr[], float out[])
 
259
#else
 
260
void heap_copy (n, arr, out)
 
261
  float arr[];
 
262
  float out[];
 
263
  int   n;
 
264
#endif
 
265
/*              copy arrays         */
 
266
{
 
267
  int    j; 
 
268
  float  diff; 
 
269
 
 
270
  for (j=0; j<n; j++)  out[j] = arr[j];
 
271
 
 
272
}
 
273
 
 
274
 
 
275
#ifdef __STDC__
 
276
int heap_compare(int n, float arr1[], float arr2[])
 
277
#else
 
278
int heap_compare (n, arr1, arr2)
 
279
  float arr1[];
 
280
  float arr2[]; 
 
281
  int   n;
 
282
#endif
 
283
/*              sort distribution (heapsort)           */
 
284
{
 
285
  int j;
 
286
  int check; 
 
287
 
 
288
  printf("Comparing arrays of size %d\n",n); 
 
289
 
 
290
  check = 0;
 
291
 
 
292
 
 
293
  for (j = 0; j < 4; j++) 
 
294
       printf("HEAPSORT: Array elements [%d] = %f %f\n",j,arr1[j],arr2[j]); 
 
295
  for (j = n-4; j < n; j++) 
 
296
       printf("HEAPSORT: Array elements [%d] = %f %f\n",j,arr1[j],arr2[j]); 
 
297
 
 
298
 
 
299
  for (j = 0; j < n; j++)
 
300
    {
 
301
      // printf("Index : %d\n",j);      
 
302
      if (arr1[j] != arr2[j]) { 
 
303
          check=1;
 
304
          printf("HEAPSORT: Array difference at index %d (%f, %f)\n",j,arr1[j],arr2[j]); 
 
305
      }
 
306
    }
 
307
  printf ("Comparison flag = %d\n",check); 
 
308
 
 
309
  return(check); 
 
310
}
 
311
 
 
312
#ifdef __STDC__
 
313
void heap_sort (int n, float arr[])
 
314
#else
 
315
void heap_sort (n, arr)
 
316
  float arr[];
 
317
  int   n;
 
318
#endif
 
319
/*              sort distribution (heapsort)           */
 
320
{
 
321
  int i, j, index;
 
322
  float a, *b;
 
323
 
 
324
  for (j = 1; j < n; j++)
 
325
    {
 
326
      a = arr[j];
 
327
      i = j - 1;
 
328
      while (i >= 0 && arr[i] > a)
 
329
        {
 
330
          arr[i + 1] = arr[i];
 
331
          i--;
 
332
        }
 
333
      arr[i + 1] = a;
 
334
    }
 
335
 
 
336
  return;
 
337
}
 
338
 
 
339
#ifdef __STDC__
 
340
float heap_median (int n, float arr[])
 
341
#else
 
342
float heap_median (n, arr)
 
343
  float arr[];
 
344
  int   n;
 
345
#endif
 
346
/*              find median of distribution (heapsort)           */
 
347
{
 
348
  int i, j, index;
 
349
  float a, *b, hmed;
 
350
  int   n2, n2p; 
 
351
 
 
352
  b = (float*) malloc ((size_t)((n)*sizeof (float)));
 
353
 
 
354
  heap_copy(n, arr, b); 
 
355
  // printf("heap_sort: Comparing copies\n"); 
 
356
  // heap_compare (n, arr, b); 
 
357
  heap_sort(n,b); 
 
358
 
 
359
  index = (n - 1) / 2;
 
360
  hmed = b[index]; 
 
361
 
 
362
  // n2p=(n2=n/2)+1;
 
363
  // median = (n % 2 ? b[n2p] : 0.5*(b[n2]+b[n2p])); // b[(n-1)/2];
 
364
 
 
365
  free(b); 
 
366
  // printf("heap_sort: Median value %f\n",hmed); 
 
367
  return(hmed);
 
368
}
 
369
 
 
370
 
 
371
/************************************************************
 
372
 *
 
373
 * median(): efficient search of the median.
 
374
 *
 
375
 ************************************************************/
 
376
 
 
377
float median(x, n)
 
378
float x[];
 
379
int n;
 
380
{
 
381
    int n2 ,n2p;
 
382
    void piksrt();
 
383
    
 
384
    piksrt(n, x);
 
385
    n2p = (n2 = n/2) + 1;
 
386
    return((n % 2) ? x[n2p] : 0.5 * (x[n2] + x[n2p]));
 
387
}
 
388
 
 
389
/************************************************************
 
390
 *
 
391
 * piksrt(): sorting by straight insertion.
 
392
 *
 
393
 ************************************************************/
 
394
 
 
395
void piksrt(n,arr)
 
396
int n;
 
397
float arr[];
 
398
{
 
399
        int i,j;
 
400
        float a;
 
401
 
 
402
        for (j=2;j<=n;j++) {
 
403
                a=arr[j];
 
404
                i=j-1;
 
405
                while (i > 0 && arr[i] > a) {
 
406
                        arr[i+1]=arr[i];
 
407
                        i--;
 
408
                }
 
409
                arr[i+1]=a;
 
410
        }
 
411
}
 
412
 
 
413
#ifdef __STDC__
 
414
  float pik_median (int n, float arr[])
 
415
#else
 
416
  float pik_median ( n, arr )
 
417
  float arr[];
 
418
  int   n;
 
419
#endif
 
420
/*              find median of distribution (use piksort)  
 
421
                NR page 330  modified for mos by Sabine ??  
 
422
                input array arr[] is conserved by buffer b[] -- here we start
 
423
                with arr[0,..,n-1] both in contrary to select_pos     TS-0896  
 
424
*/
 
425
 
 
426
 
 
427
 
 
428
{
 
429
  int i, j, index;
 
430
  float a, b[101];
 
431
 
 
432
 
 
433
  for (j = 0; j < n; j++)
 
434
    b[j] = arr[j];
 
435
 
 
436
  for (j = 1; j < n; j++)
 
437
    {
 
438
      a = b[j];
 
439
      i = j - 1;
 
440
      while (i >= 0 && b[i] > a)
 
441
        {
 
442
          b[i + 1] = b[i];
 
443
          i--;
 
444
        }
 
445
      b[i + 1] = a;
 
446
    }
 
447
  index = (n - 1) / 2;
 
448
  return(b[index]);
 
449
}
 
450