6#define DTOR (M_PI/180.)
13Livmin::Livmin(
const VPR& vpr)
15 for (
unsigned ilay = 0; ilay < vpr.size(); ++ilay)
17 if (vpr.val[ilay] <= NODATAVPR)
continue;
18 livmin = ilay * TCK_VPR + TCK_VPR / 2;
26InstantaneousVPR::InstantaneousVPR(
const Volume<double>& volume,
const Volume<unsigned char>& qual, Volume<unsigned char>& flag_vpr,
int az_min,
int az_max)
27 : volume(volume), qual(qual), flag_vpr(flag_vpr), az_min(az_min), az_max(az_max), dbz(volume)
62void InstantaneousVPR::compute()
64 LOG_CATEGORY(
"radar.vpr");
65 int i,iA,ilay,il,ilast;
66 long int dist,dist_plain,vol_rain;
67 float area,quota_true_st;
79 LOG_DEBUG(
" Iaz_min %d iaz_max %d",iaz_min,iaz_max);
81 for (
unsigned l=0; l < volume.size(); ++l)
83 const PolarScan<double>& scan = volume[l];
85 for (
unsigned k=0; k < scan.beam_size; k++)
89 dist=k*(
long int)(scan.cell_size)+(int)(scan.cell_size)/2.;
92 for (iA=iaz_min; iA<iaz_max; iA++)
95 i = (iA + volume.beam_count) % volume.beam_count;
98 quota_true_st = scan.sample_height_real(i, k);
101 ilay=floor(quota_true_st/TCK_VPR);
103 if (ilay <0 || ilay >= NMAXLAYER) {
114 const double elevaz = scan.elevations_real(i) * M_PI / 180.;
115 dist_plain=(
long int)(dist*cos(elevaz));
116 if (dist_plain <RMIN_VPR || dist_plain > RMAX_VPR )
117 flag_vpr.scan(l).set(i, k, 0);
120 if (qual.scan(l).get(i, k) < QMIN_VPR) flag_vpr.scan(l).set(i, k, 0);
124 if(cum_bac.do_class){
125 if(conv(i,k)>= CONV_VAL){
126 flag_vpr.scan(l).set(i, k, 0);
132 area = scan.cell_size * dist_plain * AMPLITUDE * DTOR * flag_vpr.scan(l).get(i, k)/1000.;
136 cv = cv + (
long int)(area);
140 double sample = scan.get(i, k);
141 if (sample > THR_VPR && flag_vpr.scan(l).get(i, k) > 0 )
144 vol_rain=(
long int)(dbz.DBZ_to_mp_func(sample)*area);
147 if (sample > THR_PDF)
148 ct = ct + (
long int)(area);
151 if (vpr.area[ilay]> 0) vpr.val[ilay]=vpr.val[ilay]+(float)(vol_rain);
152 else vpr.val[ilay]=(float)(vol_rain);
155 vpr.area[ilay]=vpr.area[ilay]+area;
161 LOG_INFO(
"calcolati ct e cv ct= %li cv= %li", ct, cv);
167 if (ct > CT_MIN * cv) {
176 for (ilay=0; ilay<NMAXLAYER; ilay++){
179 LOG_INFO(
" ilay %d area_vpr= %ld ct= %ld cv= %ld", ilay, vpr.area[ilay], ct, cv);
181 if (vpr.area[ilay]>=MIN_AREA) {
182 vert_ext=vert_ext+TCK_VPR;
183 vpr.val[ilay]=vpr.val[ilay]/(float)(vpr.area[ilay]);
188 vpr.val[ilay]=NODATAVPR;
195 if (ilay*TCK_VPR+TCK_VPR/2 > MINPOINTBASE || ilay== NMAXLAYER -1){
196 LOG_INFO(
"raggiunta cima profilo");
200 if (vert_ext<VEXTMIN_VPR ){
201 LOG_INFO(
"estensione profilo verticale troppo bassa");
204 for (il=0; il<ilast; il++) vpr.val[il]=NODATAVPR;
216 LOG_INFO(
"volume precipitante troppo piccolo");
219 for (il=0; il<NMAXLAYER; il++) vpr.val[il]=NODATAVPR;
226 grad=((vpr.val[ilast]-vpr.val[ilast-1]) + (vpr.val[ilast-1]-vpr.val[ilast-2])/(2.)+ (vpr.val[ilast-2]-vpr.val[ilast-3])/(3.) ) /3.;
229 LOG_INFO(
" %f", grad);
233 for (ilay=ilast+1; ilay<NMAXLAYER; ilay++) {
234 if (vpr.val[ilay-1] + grad > 0.002)
235 vpr.val[ilay]= vpr.val[ilay-1]+grad;
248float comp_levels(
float v0,
float v1,
float nodata,
float peso)
254 if ((v0>nodata) && (v1>nodata) ) result=((1.-peso)*v0+peso*v1);
261VPR combine_profiles(
const VPR& _vpr0,
const VPR& _vpr1,
long int cv,
long int ct)
266 vpr.area = vpr1.area;
268 long int c0 = 2 * cv;
273 for (
unsigned ilay=0; ilay<NMAXLAYER; ilay++){
274 if ( vpr0.val[ilay]> NODATAVPR && vpr1.val[ilay]>NODATAVPR ){
275 diff=diff + vpr0.val[ilay]-vpr1.val[ilay];
281 Livmin livmin1(vpr0);
282 Livmin livmin2(vpr1);
284 for (
unsigned ilay=0; ilay<max(livmin1.livmin/TCK_VPR,livmin2.livmin/TCK_VPR); ilay++){
285 if (vpr0.val[ilay]<= NODATAVPR && vpr1.val[ilay] > NODATAVPR)
286 vpr0.val[ilay]=vpr1.val[ilay]-diff;
287 vpr.area[ilay]=vpr1.area[ilay];
288 if (vpr1.val[ilay]<= NODATAVPR && vpr0.val[ilay] > NODATAVPR)
289 vpr1.val[ilay]=vpr0.val[ilay]+diff;
290 vpr.area[ilay]=vpr0.area[ilay];
295 float alfat = (float)ct / (c0 + ct);
296 for (
unsigned ilay=0; ilay<NMAXLAYER; ilay++){
297 if (vpr0.val[ilay] > NODATAVPR && vpr1.val[ilay] > NODATAVPR)
298 vpr.val[ilay]=
comp_levels(vpr0.val[ilay],vpr1.val[ilay],NODATAVPR,alfat);
299 vpr.area[ilay]=vpr1.area[ilay];
float comp_levels(float v0, float v1, float nodata, float peso)
combina livelli