5 #include <gsl/gsl_version.h>
6 #include <gsl/gsl_vector.h>
7 #include <gsl/gsl_blas.h>
8 #include <gsl/gsl_multifit_nlin.h>
19 : n(n), t(new double[n]), y(new double[n]), sigma(new double[n])
30 fprintf(stderr,
"i t y sigma\n");
31 for (
unsigned i = 0; i < n; ++i)
32 fprintf(stderr,
"%2d %.2f %.2f %.2f\n", i, t[i], y[i], sigma[i]);
37 expb_f (
const gsl_vector * x,
void *data,
40 const struct data& args = *(
struct data*)data;
41 double B = gsl_vector_get (x, 0);
42 double E = gsl_vector_get (x, 1);
43 double G = gsl_vector_get (x, 2);
44 double C = gsl_vector_get (x, 3);
45 double F = gsl_vector_get (x, 4);
46 for (
unsigned i = 0; i < args.n; i++)
50 double Yi = B * exp (- (args.t[i]-E)/G*(args.t[i]-E)/G ) + C+F*args.t[i];
51 gsl_vector_set (f, i, (Yi - args.y[i])/args.sigma[i]);
57 expb_df (
const gsl_vector * x,
void *data,
60 const struct data& args = *(
struct data*)data;
61 double B = gsl_vector_get (x, 0);
62 double E = gsl_vector_get (x, 1);
63 double G= gsl_vector_get (x, 2);
64 for (
unsigned i = 0; i < args.n; i++)
72 double arg =((args.t[i]-E)/G);
73 double ex=exp(-arg*arg);
74 double fac=B*ex*2.0*arg;
76 gsl_matrix_set (J, i, 0, ex);
77 gsl_matrix_set (J, i, 1, fac/G);
78 gsl_matrix_set (J, i, 2, fac*arg/G);
79 gsl_matrix_set (J, i, 3, 1.);
80 gsl_matrix_set (J, i, 4, args.t[i]);
87 expb_fdf (
const gsl_vector * x,
void *data,
88 gsl_vector * f, gsl_matrix * J)
97 print_state (
size_t iter, gsl_multifit_fdfsolver * s)
99 fprintf(stderr,
"iter: %3zu x = % 15.8f % 15.8f % 15.8f % 15.8f % 15.8f "
102 gsl_vector_get (s->x, 0),
103 gsl_vector_get (s->x, 1),
104 gsl_vector_get (s->x, 2),
105 gsl_vector_get (s->x, 3),
106 gsl_vector_get (s->x, 4),
107 gsl_blas_dnrm2 (s->f));
110 int testfit(
double a[])
112 if (a[0]<0. || a[0] >15.)
return 1;
113 if (a[1] >10.)
return 1;
114 if (a[2]<0.2 || a[2] > 0.6 )
return 1;
115 if (a[3]<0. )
return 1;
116 if (a[4]>0 )
return 1;
120 float lineargauss(
double xint ,
double a[])
122 return a[0] * exp (- (xint-a[1])/a[2]*(xint-a[1])/a[2] ) + a[3]+a[4]*xint;
131 int InterpolaVPR_GSL::interpola_VPR(
const float* vpr,
int hvprmax,
int livmin)
133 LOG_CATEGORY(
"radar.vpr");
134 static const unsigned N = 10;
135 const gsl_multifit_fdfsolver_type *T;
136 gsl_multifit_fdfsolver *s;
141 char file_vprint[512];
142 gsl_matrix *covar = gsl_matrix_alloc (p, p);
145 gsl_multifit_function_fdf f;
146 double x_init[5] = { 4, 0.2, 3. , 1.4, -0.4 };
147 gsl_vector_view x = gsl_vector_view_array (x_init, p);
153 int in1=(int)((hvprmax-TCK_VPR/2)/TCK_VPR);
154 int in2=(int)((hvprmax+HALF_BB)/TCK_VPR);
157 if (in4 > NMAXLAYER-1) {
166 F=vpr[in4]<vpr[in3]?(vpr[in4]-vpr[in3])/((in4-in3)*TCK_VPR/1000.):0.;
191 for (i = 0; i < n; i++)
193 d.t[i]= ((hvprmax-1000.)>livmin)? (i*TCK_VPR+(hvprmax-800)-TCK_VPR)/1000. : (livmin+i*TCK_VPR)/1000.;
194 d.y[i]= ((hvprmax-1000.)>livmin)? vpr[i+(int)(((hvprmax-800)-TCK_VPR)/TCK_VPR)] : vpr[i+(int)(livmin/TCK_VPR)];
198 T = gsl_multifit_fdfsolver_lmsder;
199 s = gsl_multifit_fdfsolver_alloc (T, n, p);
200 gsl_multifit_fdfsolver_set (s, &f, &x.vector);
204 for (
unsigned iter = 0; !found && iter < 500; ++iter)
208 int status = gsl_multifit_fdfsolver_iterate (s);
211 LOG_ERROR(
"gsl_multifit_fdfsolver_iterate: %s", gsl_strerror(status));
217 status = gsl_multifit_test_delta (s->dx, s->x,
221 case GSL_SUCCESS: found =
true;
break;
222 case GSL_CONTINUE:
break;
224 LOG_ERROR(
"gsl_multifit_test_delta: %s", gsl_strerror(status));
229 #if GSL_MAJOR_VERSION == 2
231 gsl_matrix* J = gsl_matrix_alloc(s->fdf->n, s->fdf->p);
232 gsl_multifit_fdfsolver_jac(s, J);
233 gsl_multifit_covar(J, 0.0, covar);
235 gsl_multifit_covar(s->J, 0.0, covar);
238 #define FIT(i) gsl_vector_get(s->x, i)
239 #define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
242 double chi = gsl_blas_dnrm2(s->f);
244 double c = GSL_MAX_DBL(1, chi / sqrt(dof));
261 gsl_multifit_fdfsolver_free (s);
262 gsl_matrix_free (covar);
271 xint=(i*TCK_VPR-TCK_VPR/2)/1000.;
272 yint= lineargauss(xint, a);