Elaboradar 0.1
Caricamento in corso...
Ricerca in corso...
Nessun risultato
interpola_vpr_nr.cpp
2#ifdef DO_INTERPOLA_VPR_NR
3//#include "cum_bac.h"
4#include <radarelab/logging.h>
5
6#ifdef __cplusplus
7extern "C" {
8#endif
9// nr
10#include <nrutil.h>
11#include <nr.h>
12#ifdef __cplusplus
13}
14#endif
15
16//#include <vpr_par.h> presente in interpola_vpr.h
17
18namespace {
19
31void lineargauss(float x, float a[], float *y, float dyda[], int na)
32{
33 float fac, ex, arg;
34
35 *y=0.0;
36 arg=(x-a[2])/a[3];
37 ex=exp(-arg*arg);
38 fac=a[1]*ex*2.0*arg;
39 *y+=a[1]*ex+a[4]+a[5]*x;
40 dyda[1]=ex;
41 dyda[2]=fac/a[3];
42 dyda[3]=fac*arg/a[3];
43 dyda[4]=1.;
44 dyda[5]=x;
45}
46
57int testfit(float a[], float chisq, float chisqin)
58{
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; //da analisi set dati
62 if (a[4]<0. ) return 1;
63 if (a[5]>0 ) return 1;
64 if (chisq>chisqin ) return 1;
65 return 0;
66}
67
68}
69
70namespace radarelab {
71
72/*
73 comstart interpola_VPR
74 idx interpola il profilo verticale tramite una funzione lingauss
75 interpola il profilo verticale tramite una funzione gaussiana + lineare del tipo
76
77 y= B*exp(-((x-E)/G)^2)+C+Fx
78
79 usa la funzione mrqmin delle numerical recipes in C: tutti i vettori passati a mrqmin devono essere allocati e deallcocati usando le funzioni di NR (vector, matrix, free_vector, free_matrix.. etc) che definiscono vettori con indice a partire da 1.
80 NB gli ndata dati considerati partono da 1000 m sotto il massimo (in caso il massimo sia più basso di 1000 m partono da 0 m)
81 A ogni iterazione si esegue un test sui parametri. Se ritorna 1 si torna ai valori dell'iterazione precedente.
82 A fine interpolazione si verifica che il chisquare non superi una soglia prefissata, in tal caso ritorna 1 e interpol. fallisce.
83
84 INIZIALIZZAZIONE PARAMETRI:
85 a[1]=B=vpr(liv del massimo)-vpr(liv. del massimo+500m);
86 a[2]=E= quota liv. massimo vpr (in KM);
87 a[3]=G=semiampiezza BB (quota liv .massimo- quota massimo decremento vpr nei 600 m sopra il massimo ) in KM;
88 a[4]=C=vpr(liv massimo + 700 m);
89 a[5]=F=coeff. angolare segmento con estremi nel vpr ai livelli max+900m e max+1700m, se negativo =0.;
90
91
92 float a[ma], int ma: vettore parametri e n0 parametri
93 float *x, *y: quote in KM e valori del vpr usati per l'interpolazione
94 float *sig,alamda : vettore dev st. e variabile che decrementa al convergere delle iterazioni
95 float *dyda: vettore derivate rispetto ai parametri
96 float B,E,C,G,F: parametri da ottimizzare, first guess
97 float chisq; scarto quadratico
98 int i,in1,in2,in3,in4,*ia,ifit,ii,ndati_ok,k;
99 int ndata=15; numero di dati considerati
100 float **covar,**alpha; matrice ovarianze, matrice alpha
101
102 comend
103*/
104int InterpolaVPR_NR::interpola_VPR(const float* vpr, int hvprmax, int livmin)
105{
106 LOG_CATEGORY("radar.vpr");
107 static const unsigned npar=5;
108 float *x, *y,*sig,alamda,y1=0,*dyda,xint,qdist,*abest;
109 float chisq=100.;
110 float chisqold=0.0;
111 float chisqin=0.0;
112 int i,in1,in2,in3,in4,*ia,ifit,ii,ndati_nok,k,ier_int;
113 //int ma=5;
114 int ndata=10;
115 float **covar;
116 float **alpha;
117 FILE *file;
118
119 float *a=vector(1,npar);
120 for (i=1;i<=npar;i++){
121 a[i]=NODATAVPR;
122 }
123
124 LOG_INFO("sono in interpola_vpr");
125 ier_int=0;
126
127 in1=(hvprmax-TCK_VPR/2)/TCK_VPR; //indice del massimo
128 in2=(hvprmax+HALF_BB)/TCK_VPR; //indice del massimo + 500 m
129 in3=in2+1;
130 in4=in2+5; //indice del massimo + 1000 m
131 LOG_INFO("in1 in2 %i %i %f %f",in1,in2,vpr[in1],vpr[in2]);
132
133 if (in4 > NMAXLAYER-1) {
134 ier_int=1;
135 return ier_int;
136 }
137
138 /* inizializzazione vettore parametri */
139 abest=vector(1,npar);
140 x=vector(1,ndata);
141 y=vector(1,ndata);
142 sig=vector(1,ndata);
143 ia=ivector(1,npar);
144 covar=matrix(1,npar,1,npar);
145 alpha=matrix(1,npar,1,npar);
146
147 for (k=in1+2; k<=in3; k++)
148 {
149 ier_int=0;
150
151 dyda=vector(1,npar);
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.;
155 // a[3]= G=0.25;
156 a[4]=C=vpr[in2];
157 a[5]=F=vpr[in4]<vpr[in3]?(vpr[in4]-vpr[in3])/((in4-in3)*TCK_VPR/1000.):0.;
158 //fprintf(stderr, "k:%d, a1:%f a2:%f a3:%f a4:%f a5:%f\n", k, a[1], a[2], a[3], a[4], a[5]);
159
160 alamda=-0.01;
161
162 for (i=1;i<=npar;i++) ia[i]=1;
163 qdist=0;
164 ii=1;
165 ndati_nok=0;
166
167 for (i=1; i<=ndata; i++)
168 {
169 sig[ii]=0.5;
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];
172 // x[ii]= ((hvprmax-800.)>livmin)? (i*TCK_VPR+(hvprmax-600)-TCK_VPR)/1000. : (livmin+(i-1)*TCK_VPR)/1000.;
173 //y[ii]= ((hvprmax-800.)>livmin)? vpr[i+((hvprmax-600)-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]);
176 //fprintf(stderr, "i:%d, ii:%d, xii:%f, yii:%f, y1:%f, qdist:%f, chisqin:%f\n", i, ii, x[ii], y[ii], y1, qdist, chisqin);
177 if (sqrt(qdist) < DIST_MAX)
178 {
179 ii+=1;
180 chisqin=qdist+chisqin;
181 }
182 else
183 ndati_nok=ndati_nok+1;
184
185 if ( ndati_nok > 2 )
186 {
187 LOG_WARN(" first guess troppo lontano dai punti , interpolazione fallisce");
188 ier_int=1;
189 break;
190 }
191 }
192
193 if (!ier_int)
194 {
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]);
196 //fprintf(stderr, "i t y sigma\n");
197 //for (unsigned i = 1; i <= ndata; ++i)
198 //fprintf(stderr, "%2d %.2f %.2f %.2f\n", i, x[i], y[i], sig[i]);
199 ifit=0;
200 while (fabs(chisq-chisqold) > DCHISQTHR && ifit < 1)
201 {
202 chisqold=chisq;
203 B=a[1];
204 E=a[2];
205 G=a[3];
206 C=a[4];
207 F=a[5];
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); /*test sul risultato del fit */
211 if (ifit)
212 {/*test sul risultato del fit */
213 a[1]=B; /*test sul risultato del fit */
214 a[2]=E; /*test sul risultato del fit */
215 a[3]=G; /*test sul risultato del fit */
216 a[4]=C; /*test sul risultato del fit */
217 a[5]=F; /*test sul risultato del fit */
218 chisq=chisqin;
219 }
220 }
221
222
223 if (chisq < chisqfin)
224 {
225 chisqfin=chisq;
226 for (i=1;i<=npar;i++) abest[i]=a[i];
227 }
228 }
229 }
230
231 for (i=1; i<=ndata-ndati_nok; i++)
232 {
233 lineargauss(x[i], abest, &y1, dyda, ndata);
234 rmsefin=rmsefin+ (y[i]-y1)*(y[i]-y1) ;
235 }
236 rmsefin=sqrt(rmsefin/(float)((ndata-ndati_nok)*(ndata-ndati_nok)));
237 LOG_INFO("RMSEFIN %f", rmsefin );
238
239
240 if (chisqfin>CHISQ_MAX)
241 {
242 ier_int=1;
243 }
244 else {
245 // Calcola il profilo interpolato
246 for (i=1;i<=npar;i++) a[i]=abest[i];
247 for (i=1; i<=NMAXLAYER; i++)
248 {
249 xint=(i*TCK_VPR-TCK_VPR/2)/1000.;
250 lineargauss(xint, a, &y1, dyda, ndata);
251 vpr_int[i-1] = y1;
252 }
253 }
254 B=a[1];
255 E=a[2];
256 G=a[3];
257 C=a[4];
258 F=a[5];
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);
268
269 return ier_int;
270}
271
272}
273#endif
String functions.
Definition: cart.cpp:4