Hi,
I have found this function that will calculate the best fit polynomial for a set of dx's and dy's and then save the coefficients in the 1D store[] array using the GSL GNU Scientific Library. I'm not familiar with GSL so I'm not sure exactly what the gsl_* functions are doing, but I've gotten it to work for my purposes.
However, now I would like to fit my data with a polynomial and use an initial fixed point. Is there a way I can modify this to use an initial fixed point? If so, do you have any suggestions (or resources you could point me towards) about how I would go about doing this?
Thanks!
Code:
int polynomialfit(int obs, int degree, double *dx, double *dy, double *store, int degree_max) /* n, p */
{
gsl_multifit_linear_workspace *ws;
gsl_matrix *cov, *X;
gsl_vector *y, *c;
double chisq;
int i, j;
int coeff_num = degree+1; //number of coeff is one more than the degree of the polynomial
X = gsl_matrix_alloc(obs, coeff_num);
y = gsl_vector_alloc(obs);
c = gsl_vector_alloc(coeff_num);
cov = gsl_matrix_alloc(coeff_num, coeff_num);
for(i=0; i < obs; i++) {
gsl_matrix_set(X, i, 0, 1.0);
for(j=0; j < coeff_num; j++) {
gsl_matrix_set(X, i, j, pow(dx[i], j));
}
gsl_vector_set(y, i, dy[i]);
}
ws = gsl_multifit_linear_alloc(obs, coeff_num);
gsl_multifit_linear(X, y, c, cov, &chisq, ws);
/* store result ... */
for(i=0; i < degree_max; i++) {
if (i<coeff_num){
store[i] = gsl_vector_get(c, i); //coeff[consnt, co*x, c1*x^2....etc]
}
if (i>=coeff_num){
store[i] = 0.0;
}
}
gsl_multifit_linear_free(ws);
gsl_matrix_free(X);
gsl_matrix_free(cov);
gsl_vector_free(y);
gsl_vector_free(c);
return 0;
}