41 lines
948 B
C
41 lines
948 B
C
#include "polifitgsl.h"
|
|
|
|
bool polynomialfit(int obs, int degree,
|
|
double *dx, double *dy, double *store) /* n, p */
|
|
{
|
|
gsl_multifit_linear_workspace *ws;
|
|
gsl_matrix *cov, *X;
|
|
gsl_vector *y, *c;
|
|
double chisq;
|
|
|
|
int i, j;
|
|
|
|
X = gsl_matrix_alloc(obs, degree);
|
|
y = gsl_vector_alloc(obs);
|
|
c = gsl_vector_alloc(degree);
|
|
cov = gsl_matrix_alloc(degree, degree);
|
|
|
|
for(i=0; i < obs; i++) {
|
|
for(j=0; j < degree; j++) {
|
|
gsl_matrix_set(X, i, j, pow(dx[i], j));
|
|
}
|
|
gsl_vector_set(y, i, dy[i]);
|
|
}
|
|
|
|
ws = gsl_multifit_linear_alloc(obs, degree);
|
|
gsl_multifit_linear(X, y, c, cov, &chisq, ws);
|
|
|
|
/* store result ... */
|
|
for(i=0; i < degree; i++)
|
|
{
|
|
store[i] = gsl_vector_get(c, i);
|
|
}
|
|
|
|
gsl_multifit_linear_free(ws);
|
|
gsl_matrix_free(X);
|
|
gsl_matrix_free(cov);
|
|
gsl_vector_free(y);
|
|
gsl_vector_free(c);
|
|
return true; /* we do not "analyse" the result (cov matrix mainly)
|
|
to know if the fit is "good" */
|
|
}
|