2#ifdef DO_INTERPOLA_VPR_NR
31void lineargauss(
float x,
float a[],
float *y,
float dyda[],
int na)
39 *y+=a[1]*ex+a[4]+a[5]*x;
57int testfit(
float a[],
float chisq,
float chisqin)
59 if (a[1]<0. || a[1] >15.)
return 1;
60 if (a[2] >10.)
return 1;
61 if (a[3]<0.2 || a[3] > 0.6 )
return 1;
62 if (a[4]<0. )
return 1;
63 if (a[5]>0 )
return 1;
64 if (chisq>chisqin )
return 1;
104int InterpolaVPR_NR::interpola_VPR(
const float* vpr,
int hvprmax,
int livmin)
106 LOG_CATEGORY(
"radar.vpr");
107 static const unsigned npar=5;
108 float *x, *y,*sig,alamda,y1=0,*dyda,xint,qdist,*abest;
112 int i,in1,in2,in3,in4,*ia,ifit,ii,ndati_nok,k,ier_int;
119 float *a=vector(1,npar);
120 for (i=1;i<=npar;i++){
124 LOG_INFO(
"sono in interpola_vpr");
127 in1=(hvprmax-TCK_VPR/2)/TCK_VPR;
128 in2=(hvprmax+HALF_BB)/TCK_VPR;
131 LOG_INFO(
"in1 in2 %i %i %f %f",in1,in2,vpr[in1],vpr[in2]);
133 if (in4 > NMAXLAYER-1) {
139 abest=vector(1,npar);
144 covar=matrix(1,npar,1,npar);
145 alpha=matrix(1,npar,1,npar);
147 for (k=in1+2; k<=in3; k++)
152 a[1]=B=vpr[in1]-vpr[in2];
153 a[2]=E=hvprmax/1000.;
154 a[3]=G=(k-in1-0.5)*TCK_VPR/1000.;
157 a[5]=F=vpr[in4]<vpr[in3]?(vpr[in4]-vpr[in3])/((in4-in3)*TCK_VPR/1000.):0.;
162 for (i=1;i<=npar;i++) ia[i]=1;
167 for (i=1; i<=ndata; i++)
170 x[ii]= ((hvprmax-1000.)>livmin)? (i*TCK_VPR+(hvprmax-800)-TCK_VPR)/1000. : (livmin+(i-1)*TCK_VPR)/1000.;
171 y[ii]= ((hvprmax-1000.)>livmin)? vpr[i+((hvprmax-800)-TCK_VPR)/TCK_VPR] : vpr[i-1+livmin/TCK_VPR];
174 lineargauss(x[ii], a, &y1, dyda, ndata);
175 qdist=(y1-y[ii])*(y1-y[ii]);
177 if (sqrt(qdist) < DIST_MAX)
180 chisqin=qdist+chisqin;
183 ndati_nok=ndati_nok+1;
187 LOG_WARN(
" first guess troppo lontano dai punti , interpolazione fallisce");
195 LOG_INFO(
"\n alamda %f chisqin % f a[1] % f a[2] % f a[3] % f a[4] % f a[5] % f", alamda,chisqin,a[1], a[2], a[3],a[4],a[5]);
200 while (fabs(chisq-chisqold) > DCHISQTHR && ifit < 1)
208 mrqmin(x, y, sig, ndata, a, ia, npar, covar, alpha, &chisq, &lineargauss, &alamda);
209 LOG_INFO(
"alamda %f chisq % f a[1] % f a[2] % f a[3] % f a[4] % f a[5] % f", alamda,chisq,a[1], a[2], a[3],a[4],a[5]);
210 ifit=testfit(a,chisq,chisqin);
223 if (chisq < chisqfin)
226 for (i=1;i<=npar;i++) abest[i]=a[i];
231 for (i=1; i<=ndata-ndati_nok; i++)
233 lineargauss(x[i], abest, &y1, dyda, ndata);
234 rmsefin=rmsefin+ (y[i]-y1)*(y[i]-y1) ;
236 rmsefin=sqrt(rmsefin/(
float)((ndata-ndati_nok)*(ndata-ndati_nok)));
237 LOG_INFO(
"RMSEFIN %f", rmsefin );
240 if (chisqfin>CHISQ_MAX)
246 for (i=1;i<=npar;i++) a[i]=abest[i];
247 for (i=1; i<=NMAXLAYER; i++)
249 xint=(i*TCK_VPR-TCK_VPR/2)/1000.;
250 lineargauss(xint, a, &y1, dyda, ndata);
259 free_vector(dyda,1,npar);
260 free_vector(abest,1,npar);
261 free_ivector(ia,1,npar);
262 free_vector(x,1,ndata);
263 free_vector(y,1,ndata);
264 free_vector(sig,1,ndata);
265 free_matrix(alpha,1,ndata,1,ndata);
266 free_matrix(covar,1,ndata,1,ndata);
267 free_vector(a,1,npar);