← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/spec/libsrc/fitnol.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/fitnol.c
 
1
/*===========================================================================
 
2
  Copyright (C) 1995-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        fitnol.c
 
31
.MODULE       subroutines 
 
32
.LANGUAGE     C
 
33
.AUTHOR       Cristian Levin - ESO La Silla
 
34
.PURPOSE      gaussian fitting.
 
35
.KEYWORDS     gaussian fitting.
 
36
.COMMENTS     Most of this routines were taken as they are from
 
37
              the book "Numerical Recipes in C" -- 1st edition.
 
38
.ENVIRONMENT  UNIX
 
39
.VERSION 1.0  1-Apr-1991   Implementation
 
40
 
 
41
 090723         last modif
 
42
------------------------------------------------------------*/
 
43
 
 
44
#include <math.h>
 
45
 
 
46
/*
 
47
  the root version of mrqmin is in /midas/test/libsrc/math ...
 
48
  this is essentially the same as function mrqmin of 
 
49
  /midas/test/stdred/echelle/libsrc/mrqmin.c
 
50
  using double instead of float and not returning an int 
 
51
*/
 
52
 
 
53
void mmrqmin(x,y,sig,ndata,a,ma,lista,mfit,covar,alpha,chisq,funcs,alamda)
 
54
double x[],y[],sig[],a[],**covar,**alpha,*chisq,*alamda;
 
55
int ndata,ma,lista[],mfit;
 
56
void (*funcs)();
 
57
{
 
58
        int k,kk,j,ihit;
 
59
        static double *da,*atry,**oneda,*beta,ochisq;
 
60
        double *dvector(),**dmatrix();
 
61
        void mmrqcof(),spec_gaussj(),spec_covsrt(),nrerror();
 
62
        void free_dmatrix(),free_dvector();
 
63
 
 
64
        if (*alamda < 0.0) {
 
65
                oneda=dmatrix(1,mfit,1,1);
 
66
                atry=dvector(1,ma);
 
67
                da=dvector(1,ma);
 
68
                beta=dvector(1,ma);
 
69
                kk=mfit+1;
 
70
                for (j=1;j<=ma;j++) {
 
71
                        ihit=0;
 
72
                        for (k=1;k<=mfit;k++)
 
73
                                if (lista[k] == j) ihit++;
 
74
                        if (ihit == 0)
 
75
                                lista[kk++]=j;
 
76
                        else if (ihit > 1) nrerror("Error in non linear fitting");
 
77
                }
 
78
                if (kk != ma+1) nrerror("Error in non linear fitting");
 
79
                *alamda=0.001;
 
80
                mmrqcof(x,y,sig,ndata,a,ma,lista,mfit,alpha,beta,chisq,funcs);
 
81
                ochisq=(*chisq);
 
82
        }
 
83
        for (j=1;j<=mfit;j++) {
 
84
                for (k=1;k<=mfit;k++) covar[j][k]=alpha[j][k];
 
85
                covar[j][j]=alpha[j][j]*(1.0+(*alamda));
 
86
                oneda[j][1]=beta[j];
 
87
        }
 
88
        spec_gaussj(covar,mfit,oneda,1);
 
89
        for (j=1;j<=mfit;j++)
 
90
                da[j]=oneda[j][1];
 
91
        if (*alamda == 0.0) {
 
92
                spec_covsrt(covar,ma,lista,mfit);
 
93
                free_dvector(beta,1,ma);
 
94
                free_dvector(da,1,ma);
 
95
                free_dvector(atry,1,ma);
 
96
                free_dmatrix(oneda,1,mfit,1,1);
 
97
                return;
 
98
        }
 
99
        for (j=1;j<=ma;j++) atry[j]=a[j];
 
100
        for (j=1;j<=mfit;j++)
 
101
                atry[lista[j]] = a[lista[j]]+da[j];
 
102
        mmrqcof(x,y,sig,ndata,atry,ma,lista,mfit,covar,da,chisq,funcs);
 
103
        if (*chisq < ochisq) {
 
104
                *alamda *= 0.1;
 
105
                ochisq=(*chisq);
 
106
                for (j=1;j<=mfit;j++) {
 
107
                        for (k=1;k<=mfit;k++) alpha[j][k]=covar[j][k];
 
108
                        beta[j]=da[j];
 
109
                        a[lista[j]]=atry[lista[j]];
 
110
                }
 
111
        } else {
 
112
                *alamda *= 10.0;
 
113
                *chisq=ochisq;
 
114
        }
 
115
        return;
 
116
}
 
117
 
 
118
void mmrqcof(x,y,sig,ndata,a,ma,lista,mfit,alpha,beta,chisq,funcs)
 
119
double x[],y[],sig[],a[],**alpha,beta[],*chisq;
 
120
int ndata,ma,lista[],mfit;
 
121
void (*funcs)();        /* ANSI: void (*funcs)(double,double *,double *,double *,int); */
 
122
{
 
123
        int k,j,i;
 
124
        double ymod,wt,sig2i,dy,*dyda,*dvector();
 
125
        void free_dvector();
 
126
 
 
127
        dyda=dvector(1,ma);
 
128
        for (j=1;j<=mfit;j++) {
 
129
                for (k=1;k<=j;k++) alpha[j][k]=0.0;
 
130
                beta[j]=0.0;
 
131
        }
 
132
        *chisq=0.0;
 
133
        for (i=1;i<=ndata;i++) {
 
134
                (*funcs)(x[i],a,&ymod,dyda,ma);
 
135
                sig2i=1.0/(sig[i]*sig[i]);
 
136
                dy=y[i]-ymod;
 
137
                for (j=1;j<=mfit;j++) {
 
138
                        wt=dyda[lista[j]]*sig2i;
 
139
                        for (k=1;k<=j;k++)
 
140
                                alpha[j][k] += wt*dyda[lista[k]];
 
141
                        beta[j] += dy*wt;
 
142
                }
 
143
                (*chisq) += dy*dy*sig2i;
 
144
        }
 
145
        for (j=2;j<=mfit;j++)
 
146
                for (k=1;k<=j-1;k++) alpha[k][j]=alpha[j][k];
 
147
        free_dvector(dyda,1,ma);
 
148
}
 
149
 
 
150
void spec_covsrt(covar,ma,lista,mfit)
 
151
double **covar;
 
152
int ma,lista[],mfit;
 
153
{
 
154
        int i,j;
 
155
        double swap;
 
156
 
 
157
        for (j=1;j<ma;j++)
 
158
                for (i=j+1;i<=ma;i++) covar[i][j]=0.0;
 
159
        for (i=1;i<mfit;i++)
 
160
                for (j=i+1;j<=mfit;j++) {
 
161
                        if (lista[j] > lista[i])
 
162
                                covar[lista[j]][lista[i]]=covar[i][j];
 
163
                        else
 
164
                                covar[lista[i]][lista[j]]=covar[i][j];
 
165
                }
 
166
        swap=covar[1][1];
 
167
        for (j=1;j<=ma;j++) {
 
168
                covar[1][j]=covar[j][j];
 
169
                covar[j][j]=0.0;
 
170
        }
 
171
        covar[lista[1]][lista[1]]=swap;
 
172
        for (j=2;j<=mfit;j++) covar[lista[j]][lista[j]]=covar[1][j];
 
173
        for (j=2;j<=ma;j++)
 
174
                for (i=1;i<=j-1;i++) covar[i][j]=covar[j][i];
 
175
}
 
176
 
 
177
/************************************************************
 
178
  fgauss(): optimized adding fac1, fac2. (C.Levin)
 
179
            optimized using only 3 coefs. (1 gaussian) (C.Levin).
 
180
*/
 
181
void fgauss(x, a, y, dyda, na)
 
182
double x, a[], *y, dyda[];
 
183
int na;
 
184
{
 
185
        double fac1, fac2, ex, arg;
 
186
 
 
187
        *y = 0.0;
 
188
        
 
189
        arg = (x - a[2]) / a[3];
 
190
        dyda[1] = ex = exp(-0.5 * arg * arg);
 
191
        *y = fac1 = a[1] * ex;
 
192
        dyda[2] = fac2 = fac1 * 2.0 * arg / a[3];
 
193
        dyda[3] = fac2 * arg;
 
194
}
 
195
 
 
196
/************************************************************
 
197
 *
 
198
 * fit_gauss(): Gaussian fitting. 
 
199
 * 
 
200
 * calls   : fitnol.c{mmrqmin} 
 
201
 * modified: Criterium of stopping is more relaxed (C.Levin).
 
202
 *
 
203
 ************************************************************/
 
204
 
 
205
#define EPS     0.001
 
206
 
 
207
void fit_gauss( x, y, n, a )
 
208
double *x, *y;  /* coordinates */
 
209
int n;          /*number of points */
 
210
double *a;      /* parameters of Gauss function: a[1], a[2], a[3] */
 
211
{
 
212
    int *lista;
 
213
    int nfit = 3, ncoefs = 3;
 
214
    int i, iter = 1;
 
215
    double **covar, **alpha;
 
216
    double *sig, chisq, ochisq, alamda = -1;
 
217
    void fgauss();
 
218
    double **dmatrix(), *dvector();
 
219
    int *ivector();
 
220
    void free_dmatrix(), free_dvector(), free_ivector();
 
221
 
 
222
    sig = dvector( 1, n );
 
223
    lista = ivector( 1, 3 );
 
224
    covar = dmatrix( 1, 3, 1, 3 );
 
225
    alpha = dmatrix( 1, 3, 1, 3 );
 
226
 
 
227
    for ( i = 1; i <= n; i++ )
 
228
        sig[i] = 1.0;
 
229
 
 
230
    for ( i = 1; i <= 3; i++ )
 
231
        lista[i] = i;
 
232
 
 
233
    mmrqmin( x, y, sig, n, a, ncoefs, lista, nfit, covar, alpha,
 
234
            &chisq, fgauss, &alamda );
 
235
 
 
236
    do {
 
237
        iter++;
 
238
        ochisq = chisq;
 
239
        mmrqmin( x, y, sig, n, a, ncoefs, lista, nfit, covar, alpha,
 
240
                &chisq, fgauss, &alamda );
 
241
    } while ( (ochisq - chisq) / chisq > EPS );
 
242
 
 
243
    free_dvector( sig, 1, n );
 
244
    free_ivector( lista, 1, 3 );
 
245
    free_dmatrix( covar, 1, 3, 1, 3 );
 
246
    free_dmatrix( alpha, 1, 3, 1, 3 );
 
247
}
 
248
 
 
249