Updated with options and some other things to get a usable binary
[baltrad-wrwp.git] / lib / wrwp.c
1 /* --------------------------------------------------------------------
2 Copyright (C) 2011 Swedish Meteorological and Hydrological Institute, SMHI
3
4 This is free software: you can redistribute it and/or modify
5 it under the terms of the GNU Lesser General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
8
9 This software is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12 GNU Lesser General Public License for more details.
13
14 You should have received a copy of the GNU Lesser General Public License
15 along with HLHDF.  If not, see <http://www.gnu.org/licenses/>.
16 ------------------------------------------------------------------------*/
17
18 /** Function for deriving weather radar wind and reflectivity profiles
19  *
20  * @file
21  * @author Gunther Haase, SMHI
22  * @date 2013-02-06
23  *
24  * @co-author Ulf E. Nordh, SMHI
25  * @date 2017-02-23, started overhaul of the code to achieve better
26  * resemblance with N2 and requirements from customer E-profile
27  * 
28  * @date 2018-02-08, fixed a bug that resulted in a process hang if the
29  * incoming polar volume file only contained scans with elevation angle
30  * smaller than the one stated in the call to the pgf. These kind of files
31  * are rather common when using the central software EDGE coming with the
32  * modernized Swedish radars. If suck a file is encountered, the return is just
33  * NULL and the exception following is treated in the pgf.
34  */
35
36 #include "wrwp.h"
37 #include "vertical_profile.h"
38 #include "rave_debug.h"
39 #include "rave_alloc.h"
40 #include <string.h>
41 #include <stdio.h>
42 #include <stdlib.h>
43 #include "rave_attribute.h"
44 #include "rave_utilities.h"
45 #include "rave_datetime.h"
46
47 /**
48  * Represents one wrwp generator
49  */
50 struct _Wrwp_t
51 {
52   RAVE_OBJECT_HEAD /** Always on top */
53   int dz; /**< Height interval for deriving a profile [m] */
54   int hmax; /**< Maximum height of the profile [m] */
55   int dmin; /**< Minimum distance for deriving a profile [m] */
56   int dmax; /**< Maximum distance for deriving a profile [m]*/
57   double emin; /**< Minimum elevation angle [deg] */
58   double vmin; /**< Radial velocity threshold [m/s] */
59   double nodata_VP; /**< Nodata value for vertical profile */
60   double gain_VP; /**< Gain for VP fields */
61   double offset_VP; /**< Offset for VP fields */
62   double undetect_VP; /**<Undetect for VP fields */
63 };
64
65 /*@{ Private functions */
66 /**
67  * Constructor
68  */
69 static int Wrwp_constructor(RaveCoreObject* obj)
70 {
71   Wrwp_t* wrwp = (Wrwp_t*)obj;
72   wrwp->hmax = HMAX;
73   wrwp->dz = DZ;
74   wrwp->dmin = DMIN;
75   wrwp->dmax = DMAX;
76   wrwp->emin = EMIN;
77   wrwp->vmin = VMIN;
78   wrwp->nodata_VP = NODATA_VP;
79   wrwp->gain_VP = GAIN_VP;
80   wrwp->offset_VP = OFFSET_VP;
81   wrwp->undetect_VP = UNDETECT_VP;
82   return 1;
83 }
84
85 /**
86  * Destructor
87  */
88 static void Wrwp_destructor(RaveCoreObject* obj)
89 {
90 }
91
92 static int WrwpInternal_findAndAddAttribute(VerticalProfile_t* vp, PolarVolume_t* pvol, const char* name, double minSelAng)
93 {
94   int nscans = PolarVolume_getNumberOfScans(pvol);
95   int i = 0;
96   int found = 0;
97   double elangleForThisScan = 0.0;
98   
99   for (i = 0; i < nscans && found == 0; i++) {
100     PolarScan_t* scan = PolarVolume_getScan(pvol, i);
101     if (scan != NULL && PolarScan_hasAttribute(scan, name)) {
102         elangleForThisScan = PolarScan_getElangle(scan);
103         
104         /* Filter with respect to the selected min elangle
105            that is given in the web GUI route for WRWP generation */
106         if ((elangleForThisScan * RAD2DEG) >= minSelAng) {
107           RaveAttribute_t* attr = PolarScan_getAttribute(scan, name);
108           VerticalProfile_addAttribute(vp, attr);
109           found = 1;
110           RAVE_OBJECT_RELEASE(attr);
111         }
112     }
113     RAVE_OBJECT_RELEASE(scan);
114   }
115   return found;
116 }
117
118 static int WrwpInternal_addIntAttribute(VerticalProfile_t* vp, const char* name, int value)
119 {
120   RaveAttribute_t* attr = RaveAttributeHelp_createLong(name, value);
121   int result = 0;
122   if (attr != NULL) {
123     result = VerticalProfile_addAttribute(vp, attr);
124   }
125   RAVE_OBJECT_RELEASE(attr);
126   return result;
127 }
128
129 /* Function that adds various quantities under a field's what in order to
130    better resemble the function existing in vertical profiles from N2 */
131 static int WrwpInternal_addDoubleAttr2Field(RaveField_t* field, const char* name, double quantity)
132 {
133   int result = 0;
134   RaveAttribute_t* attr = RaveAttributeHelp_createDouble(name, quantity);
135   result = RaveField_addAttribute(field, attr);
136   if (attr == NULL || !result) {
137     RAVE_ERROR1("Failed to add %s attribute to field", name);
138     goto done;
139   }
140   RAVE_OBJECT_RELEASE(attr);
141 done:
142   return result;
143 }
144
145 /**
146  * Returns a tokenized list (or empty) of fields that is wanted
147  * @param[in] fieldsToGenerate - the fields to generate
148  * @returns the list of wanted fields
149  */
150 static RaveList_t* WrwpInternal_createFieldsList(const char* fieldsToGenerate)
151 {
152   RaveList_t* result = RaveUtilities_getTrimmedTokens(fieldsToGenerate, ',');
153   if (RaveList_size(result) == 0) {
154     RaveList_add(result, RAVE_STRDUP("ff"));
155     RaveList_add(result, RAVE_STRDUP("ff_dev"));
156     RaveList_add(result, RAVE_STRDUP("dd"));
157     RaveList_add(result, RAVE_STRDUP("dbzh"));
158     RaveList_add(result, RAVE_STRDUP("dbzh_dev"));
159     RaveList_add(result, RAVE_STRDUP("nz"));
160   }
161   return result;
162 }
163
164 /**
165  * Returns if the specified id exist as an id in the fieldIds.
166  * @param[in] fieldIds - the list of ids
167  * @param[in] id - the id to query for
168  * @returns if the list of ids contains the specified id
169  */
170 static int WrwpInternal_containsField(RaveList_t* fieldIds, const char* id)
171 {
172   int i = 0, n = 0;
173   n = RaveList_size(fieldIds);
174   for (i = 0; i < n; i++) {
175     if (RaveList_get(fieldIds, i) != NULL && strcmp((char*)RaveList_get(fieldIds, i), id) == 0) {
176       return 1;
177     }
178   }
179   return 0;
180 }
181
182 static int WrwpInternal_addNodataUndetectGainOffset(RaveField_t* field, double nodata, double undetect, double gain, double offset)
183 {
184   if (!WrwpInternal_addDoubleAttr2Field(field, "what/nodata", nodata) ||
185       !WrwpInternal_addDoubleAttr2Field(field, "what/undetect", undetect) ||
186       !WrwpInternal_addDoubleAttr2Field(field, "what/gain", gain) ||
187       !WrwpInternal_addDoubleAttr2Field(field, "what/offset", offset)) {
188     return 0;
189   }
190   return 1;
191 }
192
193 static RaveDateTime_t* WrwpInternal_getStartDateTimeFromScan(PolarScan_t* scan)
194 {
195   RaveDateTime_t* dt = RAVE_OBJECT_NEW(&RaveDateTime_TYPE);
196   RaveDateTime_t* result = NULL;
197   if (dt != NULL) {
198     if (!RaveDateTime_setDate(dt, PolarScan_getStartDate(scan)) ||
199         !RaveDateTime_setTime(dt, PolarScan_getStartTime(scan))) {
200       RAVE_WARNING0("Failed to initialize datetime object with start date/time");
201       goto done;
202     }
203   }
204   result = RAVE_OBJECT_COPY(dt);
205 done:
206   RAVE_OBJECT_RELEASE(dt);
207   return result;
208 }
209
210 static RaveDateTime_t* WrwpInternal_getEndDateTimeFromScan(PolarScan_t* scan)
211 {
212   RaveDateTime_t* dt = RAVE_OBJECT_NEW(&RaveDateTime_TYPE);
213   RaveDateTime_t* result = NULL;
214   if (dt != NULL) {
215     if (!RaveDateTime_setDate(dt, PolarScan_getEndDate(scan)) ||
216         !RaveDateTime_setTime(dt, PolarScan_getEndTime(scan))) {
217       RAVE_WARNING0("Failed to initialize datetime object with end date/time");
218       goto done;
219     }
220   }
221   result = RAVE_OBJECT_COPY(dt);
222 done:
223   RAVE_OBJECT_RELEASE(dt);
224   return result;
225 }
226
227 /*@} End of Private functions */
228
229 /*@{ Interface functions */
230 void Wrwp_setDZ(Wrwp_t* self, int dz)
231 {
232   RAVE_ASSERT((self != NULL), "self == NULL");
233   self->dz = dz;
234 }
235
236 int Wrwp_getDZ(Wrwp_t* self)
237 {
238   RAVE_ASSERT((self != NULL), "self == NULL");
239   return self->dz;
240 }
241
242 void Wrwp_setNODATA_VP(Wrwp_t* self, int nodata_VP)
243 {
244   RAVE_ASSERT((self != NULL), "self == NULL");
245   self->nodata_VP = nodata_VP;
246 }
247
248 int Wrwp_getNODATA_VP(Wrwp_t* self)
249 {
250   RAVE_ASSERT((self != NULL), "self == NULL");
251   return self->nodata_VP;
252 }
253
254 void Wrwp_setUNDETECT_VP(Wrwp_t* self, int undetect_VP)
255 {
256   RAVE_ASSERT((self != NULL), "self == NULL");
257   self->undetect_VP = undetect_VP;
258 }
259
260 int Wrwp_getUNDETECT_VP(Wrwp_t* self)
261 {
262   RAVE_ASSERT((self != NULL), "self == NULL");
263   return self->undetect_VP;
264 }
265
266 void Wrwp_setGAIN_VP(Wrwp_t* self, double gain_VP)
267 {
268   RAVE_ASSERT((self != NULL), "self == NULL");
269   self->gain_VP = gain_VP;
270 }
271
272 double Wrwp_getGAIN_VP(Wrwp_t* self)
273 {
274   RAVE_ASSERT((self != NULL), "self == NULL");
275   return self->gain_VP;
276 }
277
278 void Wrwp_setOFFSET_VP(Wrwp_t* self, double offset_VP)
279 {
280   RAVE_ASSERT((self != NULL), "self == NULL");
281   self->offset_VP = offset_VP;
282 }
283
284 double Wrwp_getOFFSET_VP(Wrwp_t* self)
285 {
286   RAVE_ASSERT((self != NULL), "self == NULL");
287   return self->offset_VP;
288 }
289
290 void Wrwp_setHMAX(Wrwp_t* self, int hmax)
291 {
292   RAVE_ASSERT((self != NULL), "self == NULL");
293   self->hmax = hmax;
294 }
295
296 int Wrwp_getHMAX(Wrwp_t* self)
297 {
298   RAVE_ASSERT((self != NULL), "self == NULL");
299   return self->hmax;
300 }
301
302 void Wrwp_setDMIN(Wrwp_t* self, int dmin)
303 {
304   RAVE_ASSERT((self != NULL), "self == NULL");
305   self->dmin = dmin;
306 }
307
308 int Wrwp_getDMIN(Wrwp_t* self)
309 {
310   RAVE_ASSERT((self != NULL), "self == NULL");
311   return self->dmin;
312 }
313
314 void Wrwp_setDMAX(Wrwp_t* self, int dmax)
315 {
316   RAVE_ASSERT((self != NULL), "self == NULL");
317   self->dmax = dmax;
318 }
319
320 int Wrwp_getDMAX(Wrwp_t* self)
321 {
322   RAVE_ASSERT((self != NULL), "self == NULL");
323   return self->dmax;
324 }
325
326 void Wrwp_setEMIN(Wrwp_t* self, double emin)
327 {
328   RAVE_ASSERT((self != NULL), "self == NULL");
329   self->emin = emin;
330 }
331
332 double Wrwp_getEMIN(Wrwp_t* self)
333 {
334   RAVE_ASSERT((self != NULL), "self == NULL");
335   return self->emin;
336 }
337
338 void Wrwp_setVMIN(Wrwp_t* self, double vmin)
339 {
340   RAVE_ASSERT((self != NULL), "self == NULL");
341   self->vmin = vmin;
342 }
343
344 double Wrwp_getVMIN(Wrwp_t* self)
345 {
346   RAVE_ASSERT((self != NULL), "self == NULL");
347   return self->vmin;
348 }
349
350 VerticalProfile_t* Wrwp_generate(Wrwp_t* self, PolarVolume_t* inobj, const char* fieldsToGenerate)
351 {
352   VerticalProfile_t* result = NULL;
353   PolarNavigator_t* polnav = NULL;
354   int nrhs = NRHS, lda = LDA, ldb = LDB;
355   int nscans = 0, nv, nz, i, iz, is, ib, ir;
356   int firstInit;
357
358   double gain, offset, nodata, undetect, val;
359   double d, h;
360   double alpha, beta, gamma, vvel, vdir, vstd, zsum, zmean, zstd;
361   double centerOfLayer, u_wnd_comp, v_wnd_comp, vdir_rad;
362   int ysize = 0, yindex = 0;
363   int countAcceptedScans = 0; /* counter for accepted scans i.e. scans with elangle >= selected
364                                  minimum elevatiuon angle and not being set as malfunc */
365
366   const char* product = "VP";
367   RaveDateTime_t *firstStartDT = NULL, *lastEndDT = NULL;
368
369   /* Field defs */
370   RaveField_t *nv_field = NULL, *hght_field = NULL;
371   RaveField_t *uwnd_field = NULL, *vwnd_field = NULL;
372   RaveField_t *ff_field = NULL, *ff_dev_field = NULL, *dd_field = NULL;
373   RaveField_t *dbzh_field = NULL, *dbzh_dev_field = NULL, *nz_field = NULL;
374
375   RaveList_t* wantedFields = NULL;
376
377   RAVE_ASSERT((self != NULL), "self == NULL");
378
379   wantedFields = WrwpInternal_createFieldsList(fieldsToGenerate);
380
381   if (WrwpInternal_containsField(wantedFields, "NV")) nv_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
382   if (WrwpInternal_containsField(wantedFields, "HGHT")) hght_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
383   if (WrwpInternal_containsField(wantedFields, "UWND")) uwnd_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
384   if (WrwpInternal_containsField(wantedFields, "VWND")) vwnd_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
385   if (WrwpInternal_containsField(wantedFields, "ff")) ff_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
386   if (WrwpInternal_containsField(wantedFields, "ff_dev")) ff_dev_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
387   if (WrwpInternal_containsField(wantedFields, "dd")) dd_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
388   if (WrwpInternal_containsField(wantedFields, "dbzh")) dbzh_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
389   if (WrwpInternal_containsField(wantedFields, "dbzh_dev")) dbzh_dev_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
390   if (WrwpInternal_containsField(wantedFields, "nz")) nz_field = RAVE_OBJECT_NEW(&RaveField_TYPE);
391
392   if ((WrwpInternal_containsField(wantedFields, "NV") && nv_field == NULL) ||
393       (WrwpInternal_containsField(wantedFields, "HGHT") && hght_field == NULL) ||
394       (WrwpInternal_containsField(wantedFields, "UWND") && uwnd_field == NULL) ||
395       (WrwpInternal_containsField(wantedFields, "VWND") && vwnd_field == NULL) ||
396       (WrwpInternal_containsField(wantedFields, "ff") && ff_field == NULL) ||
397       (WrwpInternal_containsField(wantedFields, "ff_dev") && ff_dev_field == NULL) ||
398       (WrwpInternal_containsField(wantedFields, "dd") && dd_field == NULL) ||
399       (WrwpInternal_containsField(wantedFields, "dbzh") && dbzh_field == NULL) ||
400       (WrwpInternal_containsField(wantedFields, "dbzh_dev") && dbzh_dev_field == NULL) ||
401       (WrwpInternal_containsField(wantedFields, "nz") && nz_field == NULL))
402   {
403     RAVE_ERROR0("Failed to allocate memory for the resulting vp fields");
404     goto done;
405   }
406
407   ysize = self->hmax / self->dz;
408
409   if ((nv_field != NULL && !RaveField_createData(nv_field, 1, ysize, RaveDataType_INT)) ||
410       (hght_field != NULL && !RaveField_createData(hght_field, 1, ysize, RaveDataType_DOUBLE)) ||
411       (uwnd_field != NULL && !RaveField_createData(uwnd_field, 1, ysize, RaveDataType_DOUBLE)) ||
412       (vwnd_field != NULL && !RaveField_createData(vwnd_field, 1, ysize, RaveDataType_DOUBLE)) ||
413       (ff_field != NULL && !RaveField_createData(ff_field, 1, ysize, RaveDataType_DOUBLE)) ||
414       (ff_dev_field != NULL && !RaveField_createData(ff_dev_field, 1, ysize, RaveDataType_DOUBLE)) ||
415       (dd_field != NULL && !RaveField_createData(dd_field, 1, ysize, RaveDataType_DOUBLE)) ||
416       (dbzh_field != NULL && !RaveField_createData(dbzh_field, 1, ysize, RaveDataType_DOUBLE)) ||
417       (dbzh_dev_field != NULL && !RaveField_createData(dbzh_dev_field, 1, ysize, RaveDataType_DOUBLE)) ||
418       (nz_field != NULL && !RaveField_createData(nz_field, 1, ysize, RaveDataType_DOUBLE))) {
419     RAVE_ERROR0("Failed to allocate arrays for the resulting vp fields");
420     goto done;
421   }
422
423   polnav = RAVE_OBJECT_NEW(&PolarNavigator_TYPE);
424   PolarNavigator_setLat0(polnav, PolarVolume_getLatitude(inobj));
425   PolarNavigator_setLon0(polnav, PolarVolume_getLongitude(inobj));
426   PolarNavigator_setAlt0(polnav, PolarVolume_getHeight(inobj));
427             
428   nscans = PolarVolume_getNumberOfScans (inobj);
429
430   // We use yindex for filling in the arrays even though we loop to hmax...
431   yindex = 0;
432
433   // loop over atmospheric layers
434   for (iz = 0; iz < self->hmax; iz += self->dz) {
435     /* allocate memory and initialize with zeros */
436     double *A = RAVE_CALLOC((size_t)(NOR*NOC), sizeof (double));
437     double *b = RAVE_CALLOC((size_t)(NOR), sizeof (double));
438     double *v = RAVE_CALLOC((size_t)(NOR), sizeof (double));
439     double *z = RAVE_CALLOC((size_t)(NOR), sizeof (double));
440     double *az = RAVE_CALLOC((size_t)(NOR), sizeof (double));
441
442     vdir = -9999.0;
443     vvel = -9999.0;
444     vstd = 0.0;
445     zsum = 0.0;
446     zmean = -9999.0;
447     zstd = 0.0;
448     nv = 0;
449     nz = 0;
450     
451     /* Define the center height of each verical layer, this will later
452        become the HGHT array */
453     centerOfLayer = iz + (self->dz / 2.0);
454     firstInit = 0;
455
456     // loop over scans
457     for (is = 0; is < nscans; is++) {
458       char* malfuncString = NULL;
459       PolarScan_t* scan = PolarVolume_getScan(inobj, is);      
460       long nbins = PolarScan_getNbins(scan);
461       long nrays = PolarScan_getNrays(scan);
462       double rscale = PolarScan_getRscale(scan);
463       double elangleForThisScan = PolarScan_getElangle(scan);
464       
465       if (elangleForThisScan * RAD2DEG >= self->emin) { /* We only do the calculation for scans with elangle >= the minimum one */
466         RaveDateTime_t* startDTofThisScan = WrwpInternal_getStartDateTimeFromScan(scan);
467         RaveDateTime_t* endDTofThisScan = WrwpInternal_getEndDateTimeFromScan(scan);        
468         RaveAttribute_t* malfuncattr = PolarScan_getAttribute(scan, "how/malfunc");
469         if (malfuncattr != NULL) {
470           RaveAttribute_getString(malfuncattr, &malfuncString); /* Set the malfuncString if attr is not NULL */
471           RAVE_OBJECT_RELEASE(malfuncattr);
472         }
473         if (malfuncString == NULL || strcmp(malfuncString, "False") == 0) { /* Assuming malfuncString = NULL means no malfunc */
474           countAcceptedScans = countAcceptedScans + 1;
475           if (firstInit == 0 && iz==0) { /* We only perform the date time calc at first iz-iteration*/
476             /* Initialize using the first accepted scan and define 2 strings for the combined datetime */
477             firstStartDT = RAVE_OBJECT_COPY(startDTofThisScan);
478             lastEndDT = RAVE_OBJECT_COPY(endDTofThisScan);
479             firstInit = 1;
480           }
481           if (firstInit == 1 && countAcceptedScans > 1 && iz==0) {
482             if (RaveDateTime_compare(startDTofThisScan, firstStartDT) < 0) {
483               /* if start datetime of this scan is before the first saved start datetime, save this one instead */
484               RAVE_OBJECT_RELEASE(firstStartDT);
485               firstStartDT = RAVE_OBJECT_COPY(startDTofThisScan);
486             }
487             if (RaveDateTime_compare(endDTofThisScan, lastEndDT) > 0) {
488               /* If end datetime of this scan is after the last saved end datetime, save this one instead */
489               RAVE_OBJECT_RELEASE(lastEndDT);
490               lastEndDT = RAVE_OBJECT_COPY(endDTofThisScan);
491             }
492           }
493           // radial wind scans
494           if (PolarScan_hasParameter(scan, "VRAD") || PolarScan_hasParameter(scan, "VRADH")) {
495             PolarScanParam_t* vrad = NULL;
496             if (PolarScan_hasParameter(scan, "VRAD")) {
497               vrad = PolarScan_getParameter(scan, "VRAD");
498             } else {
499               vrad = PolarScan_getParameter(scan, "VRADH");
500             } 
501             gain = PolarScanParam_getGain(vrad);
502             offset = PolarScanParam_getOffset(vrad);
503             nodata = PolarScanParam_getNodata(vrad);
504             undetect = PolarScanParam_getUndetect(vrad);
505
506             for (ir = 0; ir < nrays; ir++) {
507               for (ib = 0; ib < nbins; ib++) {
508                 PolarNavigator_reToDh(polnav, (ib+0.5)*rscale, elangleForThisScan, &d, &h);
509                 PolarScanParam_getValue(vrad, ib, ir, &val);
510                 if ((h >= iz) &&
511                     (h < iz + self->dz) &&
512                     (d >= self->dmin) &&
513                     (d <= self->dmax) &&
514                     (val != nodata) &&
515                     (val != undetect) &&
516                     (abs(offset + gain * val) >= self->vmin)) {
517                   *(v+nv) = offset+gain*val;
518                   *(az+nv) = 360./nrays*ir*DEG2RAD;
519                   *(A+nv*NOC) = sin(*(az+nv));
520                   *(A+nv*NOC+1) = cos(*(az+nv));
521                   *(A+nv*NOC+2) = 1;
522                   *(b+nv) = *(v+nv);
523                   nv = nv+1;
524                 }
525               }
526             }
527             RAVE_OBJECT_RELEASE(vrad);          
528           }
529
530           // reflectivity scans
531           if (PolarScan_hasParameter(scan, "DBZH")) {
532             PolarScanParam_t* dbz = PolarScan_getParameter(scan, "DBZH");
533             gain = PolarScanParam_getGain(dbz);
534             offset = PolarScanParam_getOffset(dbz);
535             nodata = PolarScanParam_getNodata(dbz);
536             undetect = PolarScanParam_getUndetect(dbz);
537
538             for (ir = 0; ir < nrays; ir++) {
539               for (ib = 0; ib < nbins; ib++) {
540                 PolarNavigator_reToDh(polnav, (ib+0.5)*rscale, elangleForThisScan, &d, &h);
541                 PolarScanParam_getValue (dbz, ib, ir, &val);
542                 if ((h >= iz) &&
543                     (h < iz+self->dz) &&
544                     (d >= self->dmin) &&
545                     (d <= self->dmax) &&
546                     (val != nodata) &&
547                     (val != undetect)) {
548                  *(z+nz) = dBZ2Z(offset+gain*val);
549                  zsum = zsum + *(z+nz);
550                 nz = nz+1;
551                 }
552               }
553             }
554             RAVE_OBJECT_RELEASE(dbz);         
555           }        
556         }
557         RAVE_OBJECT_RELEASE(startDTofThisScan);
558         RAVE_OBJECT_RELEASE(endDTofThisScan);
559       }
560       RAVE_OBJECT_RELEASE(scan); 
561     }
562
563     if (countAcceptedScans == 0) { /* Emergency exit if no accepted scans were found */
564       result = NULL;               /* If this is the case, we don't bother with checking */
565       RAVE_FREE(A);                /* the same thing for the other atmospheric layers, */
566       RAVE_FREE(b);                /* we skip it and return 0 directly */
567       RAVE_FREE(v);
568       RAVE_FREE(z);
569       RAVE_FREE(az);
570       RAVE_INFO0("Could not find any acceptable scans, dropping out");
571       goto done;
572     } 
573
574     /* Perform radial wind calculations and reflectivity calculations */
575     if (nv>3) {
576       //***************************************************************
577       // fitting: y = gamma+alpha*sin(x+beta)                         *
578       // alpha -> amplitude                                           *
579       // beta -> phase shift                                          *
580       // gamma -> consider an y-shift due to the terminal velocity of *
581       //          falling rain drops                                  *
582       //***************************************************************
583
584       /* QR decomposition */
585       /*info = */
586       LAPACKE_dgels(LAPACK_ROW_MAJOR, 'N', NOR, NOC, nrhs, A, lda, b, ldb);
587
588       /* parameter of the wind model */
589       alpha = sqrt(pow(*(b),2) + pow(*(b+1),2));
590       beta = atan2(*(b+1), *b);
591       gamma = *(b+2);
592
593       /* wind velocity */
594       vvel = alpha;
595
596       /* wind direction */
597       vdir = 0;
598       if (alpha < 0) {
599         vdir = (PI/2-beta) * RAD2DEG;
600       } else if (alpha > 0) {
601         vdir = (3*PI/2-beta) * RAD2DEG;
602       }
603       if (vdir < 0) {
604         vdir = vdir + 360;
605       } else if (vdir > 360) {
606         vdir = vdir - 360;
607       }
608
609       /* standard deviation of the wind velocity */
610       for (i=0; i<nv; i++) {
611         vstd = vstd + pow (*(v+i) - (gamma+alpha*sin(*(az+i)+beta)),2);
612       }
613       vstd = sqrt (vstd/(nv-1));
614       
615       /* Calculate the x-component (East) and y-component (North) of the wind 
616          velocity using the wind direction and the magnitude of the wind velocity */
617       vdir_rad = vdir * DEG2RAD;
618       u_wnd_comp = vvel * cos(vdir_rad);
619       v_wnd_comp = vvel * sin(vdir_rad);
620     }
621
622     // reflectivity calculations
623     if (nz > 1) {
624       /* standard deviation of reflectivity */
625       for (i = 0; i < nz; i++) {
626         zstd = zstd + pow(*(z+i) - (zsum/nz),2);
627       }
628       zmean = Z2dBZ(zsum/nz);
629       zstd = sqrt(zstd/(nz-1));
630       zstd = Z2dBZ(zstd);
631     }
632            
633     if ((nv < NMIN) || (nz < NMIN)) {
634       if (nv_field != NULL) RaveField_setValue(nv_field, 0, yindex, 0);
635       if (hght_field != NULL) RaveField_setValue(hght_field, 0, yindex,centerOfLayer);
636       if (uwnd_field != NULL) RaveField_setValue(uwnd_field, 0, yindex,self->nodata_VP);
637       if (vwnd_field != NULL) RaveField_setValue(vwnd_field, 0, yindex,self->nodata_VP);
638       if (ff_field != NULL) RaveField_setValue(ff_field, 0, yindex, self->nodata_VP);
639       if (ff_dev_field != NULL) RaveField_setValue(ff_dev_field, 0, yindex, self->nodata_VP);
640       if (dd_field != NULL) RaveField_setValue(dd_field, 0, yindex, self->nodata_VP);
641       if (dbzh_field != NULL) RaveField_setValue(dbzh_field, 0, yindex, self->nodata_VP);
642       if (dbzh_dev_field != NULL) RaveField_setValue(dbzh_dev_field, 0, yindex, self->nodata_VP);
643       if (nz_field != NULL) RaveField_setValue(nz_field, 0, yindex, 0);
644     } else {
645       if (nv_field != NULL) RaveField_setValue(nv_field, 0, yindex, nv);
646       if (hght_field != NULL) RaveField_setValue(hght_field, 0, yindex,centerOfLayer);
647       if (uwnd_field != NULL) RaveField_setValue(uwnd_field, 0, yindex,u_wnd_comp);
648       if (vwnd_field != NULL) RaveField_setValue(vwnd_field, 0, yindex,v_wnd_comp);
649       if (ff_field != NULL) RaveField_setValue(ff_field, 0, yindex, vvel);
650       if (ff_dev_field != NULL) RaveField_setValue(ff_dev_field, 0, yindex, vstd);
651       if (dd_field != NULL) RaveField_setValue(dd_field, 0, yindex, vdir);
652       if (dbzh_field != NULL) RaveField_setValue(dbzh_field, 0, yindex, zmean);
653       if (dbzh_dev_field != NULL) RaveField_setValue(dbzh_dev_field, 0, yindex, zstd);
654       if (nz_field != NULL) RaveField_setValue(nz_field, 0, yindex, nz);
655     }
656
657     RAVE_FREE(A);
658     RAVE_FREE(b);
659     RAVE_FREE(v);
660     RAVE_FREE(z);
661     RAVE_FREE(az);
662
663     yindex++;   
664   }
665
666   if (uwnd_field) WrwpInternal_addNodataUndetectGainOffset(uwnd_field, self->nodata_VP, self->undetect_VP, self->gain_VP, self->offset_VP);
667   if (vwnd_field) WrwpInternal_addNodataUndetectGainOffset(vwnd_field, self->nodata_VP, self->undetect_VP, self->gain_VP, self->offset_VP);
668   if (hght_field) WrwpInternal_addNodataUndetectGainOffset(hght_field, self->nodata_VP, self->undetect_VP, self->gain_VP, self->offset_VP);
669   if (nv_field) WrwpInternal_addNodataUndetectGainOffset(nv_field, self->nodata_VP, self->undetect_VP, self->gain_VP, self->offset_VP);
670   if (ff_field) WrwpInternal_addNodataUndetectGainOffset(ff_field, self->nodata_VP, self->undetect_VP, self->gain_VP, self->offset_VP);
671   if (ff_dev_field) WrwpInternal_addNodataUndetectGainOffset(ff_dev_field, self->nodata_VP, self->undetect_VP, self->gain_VP, self->offset_VP);
672   if (dd_field) WrwpInternal_addNodataUndetectGainOffset(dd_field, self->nodata_VP, self->undetect_VP, self->gain_VP, self->offset_VP);
673   if (dbzh_field) WrwpInternal_addNodataUndetectGainOffset(dbzh_field, self->nodata_VP, self->undetect_VP, self->gain_VP, self->offset_VP);
674   if (dbzh_dev_field) WrwpInternal_addNodataUndetectGainOffset(dbzh_dev_field, self->nodata_VP, self->undetect_VP, self->gain_VP, self->offset_VP);
675
676   result = RAVE_OBJECT_NEW(&VerticalProfile_TYPE);
677   if (result != NULL) {
678     VerticalProfile_setLevels(result, ysize);
679     if ((uwnd_field != NULL && !VerticalProfile_setUWND(result, uwnd_field)) ||
680         (vwnd_field != NULL && !VerticalProfile_setVWND(result, vwnd_field)) ||
681         (nv_field != NULL && !VerticalProfile_setNV(result, nv_field)) ||
682         (hght_field != NULL && !VerticalProfile_setHGHT(result, hght_field)) ||
683         (ff_field != NULL && !VerticalProfile_setFF(result, ff_field)) ||
684         (ff_dev_field != NULL && !VerticalProfile_setFFDev(result, ff_dev_field)) ||
685         (dd_field != NULL && !VerticalProfile_setDD(result, dd_field)) ||
686         (dbzh_field != NULL && !VerticalProfile_setDBZ(result, dbzh_field)) ||
687         (dbzh_dev_field != NULL && !VerticalProfile_setDBZDev(result, dbzh_dev_field))) {
688       RAVE_ERROR0("Failed to set vertical profile fields");
689       RAVE_OBJECT_RELEASE(result);
690     }
691   }
692   VerticalProfile_setLongitude(result, PolarVolume_getLongitude(inobj));
693   VerticalProfile_setLatitude(result, PolarVolume_getLatitude(inobj));
694   VerticalProfile_setHeight(result, PolarVolume_getHeight(inobj));
695   VerticalProfile_setSource(result, PolarVolume_getSource(inobj));
696   VerticalProfile_setInterval(result, self->dz);
697   VerticalProfile_setMinheight(result, 0);
698   VerticalProfile_setMaxheight(result, self->hmax);
699   VerticalProfile_setDate(result, PolarVolume_getDate(inobj));
700   VerticalProfile_setTime(result, PolarVolume_getTime(inobj));
701
702   /* Set the times and product, starttime is the starttime for the lowest elev
703      endtime is the endtime for the highest elev. */
704   VerticalProfile_setStartDate(result, RaveDateTime_getDate(firstStartDT));
705   VerticalProfile_setStartTime(result, RaveDateTime_getTime(firstStartDT));
706   VerticalProfile_setEndDate(result, RaveDateTime_getDate(lastEndDT));
707   VerticalProfile_setEndTime(result, RaveDateTime_getTime(lastEndDT));
708   VerticalProfile_setProduct(result, product);
709    
710   /* Need to filter the attribute data with respect to selected min elangle
711      Below attributes satisfy reguirements from E-profile, add other when needed */
712   WrwpInternal_findAndAddAttribute(result, inobj, "how/highprf", self->emin);
713   WrwpInternal_findAndAddAttribute(result, inobj, "how/lowprf", self->emin);
714   WrwpInternal_findAndAddAttribute(result, inobj, "how/pulsewidth", self->emin);
715   WrwpInternal_findAndAddAttribute(result, inobj, "how/wavelength", self->emin);
716   WrwpInternal_findAndAddAttribute(result, inobj, "how/RXbandwidth", self->emin);
717   WrwpInternal_findAndAddAttribute(result, inobj, "how/RXlossH", self->emin);
718   WrwpInternal_findAndAddAttribute(result, inobj, "how/TXlossH", self->emin);
719   WrwpInternal_findAndAddAttribute(result, inobj, "how/antgainH", self->emin);
720   WrwpInternal_findAndAddAttribute(result, inobj, "how/azmethod", self->emin);
721   WrwpInternal_findAndAddAttribute(result, inobj, "how/binmethod", self->emin);
722   WrwpInternal_findAndAddAttribute(result, inobj, "how/malfunc", self->emin);
723   WrwpInternal_findAndAddAttribute(result, inobj, "how/nomTXpower", self->emin);
724   WrwpInternal_findAndAddAttribute(result, inobj, "how/radar_msg", self->emin);
725   WrwpInternal_findAndAddAttribute(result, inobj, "how/radconstH", self->emin);
726   WrwpInternal_findAndAddAttribute(result, inobj, "how/radomelossH", self->emin);
727   WrwpInternal_findAndAddAttribute(result, inobj, "how/rpm", self->emin);
728   WrwpInternal_findAndAddAttribute(result, inobj, "how/software", self->emin);
729   WrwpInternal_findAndAddAttribute(result, inobj, "how/system", self->emin);
730   WrwpInternal_addIntAttribute(result, "how/minrange", Wrwp_getDMIN(self));
731   WrwpInternal_addIntAttribute(result, "how/maxrange", Wrwp_getDMAX(self));
732
733 done:
734   RAVE_OBJECT_RELEASE(polnav);
735   RAVE_OBJECT_RELEASE(ff_field);
736   RAVE_OBJECT_RELEASE(ff_dev_field);
737   RAVE_OBJECT_RELEASE(dd_field);
738   RAVE_OBJECT_RELEASE(dbzh_field);
739   RAVE_OBJECT_RELEASE(dbzh_dev_field);
740   RAVE_OBJECT_RELEASE(nz_field);
741   RAVE_OBJECT_RELEASE(nv_field);
742   RAVE_OBJECT_RELEASE(hght_field);
743   RAVE_OBJECT_RELEASE(uwnd_field);
744   RAVE_OBJECT_RELEASE(vwnd_field);
745   RAVE_OBJECT_RELEASE(firstStartDT);
746   RAVE_OBJECT_RELEASE(lastEndDT);
747   RaveList_freeAndDestroy(&wantedFields);
748
749   return result;
750 }
751
752 /*@} End of Interface functions */
753
754 RaveCoreObjectType Wrwp_TYPE = {
755     "Wrwp",
756     sizeof(Wrwp_t),
757     Wrwp_constructor,
758     Wrwp_destructor
759 };
760