Elaboradar  0.1
interpola_vpr_nr.cpp
2 #ifdef DO_INTERPOLA_VPR_NR
3 //#include "cum_bac.h"
4 #include <radarelab/logging.h>
5 
6 #ifdef __cplusplus
7 extern "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 
18 namespace {
19 
31 void 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 
57 int 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 
70 namespace 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 */
104 int 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