f944c9612e69e80990a03d3c0d130f962a7cefb0
[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  * @file
20  * @author Gunther Haase, SMHI
21  * @date 2013-02-06
22  */
23
24 #include "wrwp.h"
25
26 PolarVolume_t* wrwp(PolarVolume_t* inobj) {
27         PolarVolume_t* result = NULL;
28         PolarScan_t* scan = NULL;
29         PolarScanParam_t* vrad = NULL;
30         PolarScanParam_t* dbz = NULL;
31         PolarNavigator_t* polnav = NULL;
32         int nrhs = NRHS, lda = LDA, ldb = LDB, info;
33         int nscans = 0, nbins = 0, nrays, nv, nz, i, iz, is, ib, ir;
34         double rscale, elangle, gain, offset, nodata, undetect, val;
35         double d, h;
36         double alpha, beta, gamma, vvel, vdir, vstd, zsum, zmean, zstd;
37
38         result = RAVE_OBJECT_COPY (inobj);
39         polnav = RAVE_OBJECT_NEW (&PolarNavigator_TYPE);
40         nscans = PolarVolume_getNumberOfScans (inobj);
41
42         // loop over atmospheric layers
43         for (iz=0; iz<HMAX; iz+=DZ) {
44
45                 /* allocate memory and initialize with zeros */
46                 double *A = RAVE_CALLOC ((size_t)(NOR*NOC), sizeof (double));
47                 double *b = RAVE_CALLOC ((size_t)(NOR), sizeof (double));
48                 double *v = RAVE_CALLOC ((size_t)(NOR), sizeof (double));
49                 double *z = RAVE_CALLOC ((size_t)(NOR), sizeof (double));
50                 double *az = RAVE_CALLOC ((size_t)(NOR), sizeof (double));
51
52                 vdir = NAN;
53                 vvel = NAN;
54                 vstd = 0.0;
55                 zsum = 0.0;
56                 zmean = NAN;
57                 zstd = 0.0;
58                 nv = 0;
59                 nz = 0;
60
61                 // loop over scans
62                 for (is=0; is<nscans; is++) {
63                         scan = PolarVolume_getScan (inobj, is);
64                         nbins = PolarScan_getNbins (scan);
65                         nrays = PolarScan_getNrays (scan);
66                         rscale = PolarScan_getRscale (scan);
67                         elangle = PolarScan_getElangle (scan);
68
69                         // radial wind scans
70                         if (PolarScan_hasParameter (scan, "VRAD")) {
71                                 vrad = PolarScan_getParameter (scan, "VRAD");
72                                 gain = PolarScanParam_getGain (vrad);
73                                 offset = PolarScanParam_getOffset (vrad);
74                                 nodata = PolarScanParam_getNodata (vrad);
75                                 undetect = PolarScanParam_getUndetect (vrad);
76
77                                 for (ir=0; ir<nrays; ir++) {
78                                         for (ib=0; ib<nbins; ib++) {
79                                                 PolarNavigator_reToDh (polnav, (ib+0.5)*rscale, elangle, &d, &h);
80                                                 PolarScanParam_getValue (vrad, ib, ir, &val);
81                                                 if ((h>=iz) && (h<iz+DZ) && (d>=DMIN) && (d<=DMAX) && (elangle*RAD2DEG>=EMIN) &&
82                                                         (val!=nodata) && (val!=undetect) && (abs(offset+gain*val)>=VMIN)) {
83                                                         *(v+nv) = offset+gain*val;
84                                                         *(az+nv) = 360./nrays*ir*DEG2RAD;
85                                                         *(A+nv*NOC) = sin(*(az+nv));
86                                                         *(A+nv*NOC+1) = cos(*(az+nv));
87                                                         *(A+nv*NOC+2) = 1;
88                                                         *(b+nv) = *(v+nv);
89                                                         nv = nv+1;
90                                                 }
91                                         }
92                                 }
93                         }
94
95                         // reflectivity scans
96                         if (PolarScan_hasParameter (scan, "DBZH")) {
97                                 dbz = PolarScan_getParameter (scan, "DBZH");
98                                 gain = PolarScanParam_getGain (dbz);
99                                 offset = PolarScanParam_getOffset (dbz);
100                                 nodata = PolarScanParam_getNodata (dbz);
101                                 undetect = PolarScanParam_getUndetect (dbz);
102
103                                 for (ir=0; ir<nrays; ir++) {
104                                         for (ib=0; ib<nbins; ib++) {
105                                                 PolarNavigator_reToDh (polnav, (ib+0.5)*rscale, elangle, &d, &h);
106                                                 PolarScanParam_getValue (dbz, ib, ir, &val);
107                                                 if ((h>=iz) && (h<iz+DZ) && (d>=DMIN) && (d<=DMAX) && (elangle*RAD2DEG>=EMIN) &&
108                                                         (val!=nodata) && (val!=undetect)) {
109                                                         *(z+nz) = dBZ2Z (offset+gain*val);
110                                                         zsum = zsum + *(z+nz);
111                                                         nz = nz+1;
112                                                 }
113                                         }
114                                 }
115                         }
116                 }
117
118                 // radial wind calculations
119                 if (nv>3) {
120
121                 //***************************************************************
122                 // fitting: y = gamma+alpha*sin(x+beta)                         *
123                 // alpha -> amplitude                                           *
124                 // beta -> phase shift                                          *
125                 // gamma -> consider an y-shift due to the terminal velocity of *
126                 //          falling rain drops                                  *
127                 //***************************************************************
128
129                         /* QR decomposition */
130                         info = LAPACKE_dgels (LAPACK_ROW_MAJOR, 'N', NOR, NOC, nrhs, A, lda, b, ldb);
131
132                         /* parameter of the wind model */
133                         alpha = sqrt (pow (*(b),2) + pow (*(b+1),2));
134                         beta = atan2 (*(b+1), *b);
135                         gamma = *(b+2);
136
137                         /* wind velocity */
138                         vvel = alpha;
139
140                 /* wind direction */
141                 vdir = 0;
142                 if (alpha < 0) vdir = (PI/2-beta) * RAD2DEG;
143                 else if (alpha > 0) vdir = (3*PI/2-beta) * RAD2DEG;
144                 if (vdir < 0) vdir = vdir + 360;
145                 else if (vdir > 360) vdir = vdir - 360;
146
147                 /* standard deviation of the wind velocity */
148                 for (i=0; i<nv; i++) {
149                         vstd = vstd + pow (*(v+i) - (gamma+alpha*sin (*(az+i)+beta)),2);
150                 }
151                 vstd = sqrt (vstd/(nv-1));
152                 }
153
154                 // reflectivity calculations
155                 if (nz>1) {
156
157                         /* standard deviation of reflectivity */
158                         for (i=0; i<nz; i++) {
159                                 zstd = zstd + pow (*(z+i) - (zsum/nz),2);
160                         }
161                         zmean = Z2dBZ (zsum/nz);
162                         zstd = sqrt (zstd/(nz-1));
163                         zstd = Z2dBZ (zstd);
164                 }
165
166                 if ((nv<NMIN) || (nz<NMIN))
167                         printf("%6d %6d %6.2f %6.2f %6.2f %6.2f %6.2f %6.2f \n", iz+DZ/2, 0, -9999., -9999., -9999., -9999., -9999., -9999.);
168                 else
169                         printf("%6d %6d %8.2f %8.2f %8.2f %8.2f %8.2f %8.2f \n", iz+DZ/2, nv, vvel, vstd, vdir, -9999., zmean, zstd);
170
171                 RAVE_FREE (A);
172                 RAVE_FREE (b);
173                 RAVE_FREE (v);
174                 RAVE_FREE (z);
175                 RAVE_FREE (az);
176
177         }
178         return result;
179 }