GCC Code Coverage Report

Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 78.5% 893 / 0 / 1137
Functions: 100.0% 87 / 0 / 87
Branches: 62.6% 320 / 0 / 511
xsect.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // xsect.c
3 //
4 // Project: EPA SWMM5
5 // Version: 5.2
6 // Date: 10/17/22 (Build 5.2.2)
7 // Author: L. Rossman
8 // M. Tryby (EPA)
9 //
10 // Cross section geometry functions.
11 //
12 // The primary functions are:
13 // getAofY -- returns area given depth
14 // getWofY -- returns top width given depth
15 // getRofY -- returns hyd. radius given depth
16 // getYofA -- returns flow depth given area
17 // getRofA -- returns hyd. radius given area
18 // getSofA -- returns section factor given area
19 // getAofS -- returns area given section factor
20 // getdSdA -- returns derivative of section factor w.r.t. area
21 // where
22 // Y = flow depth
23 // A = flow area
24 // R = hyd. radius
25 // S = section factor = A*R^(2/3)
26 //
27 // Update History
28 // ==============
29 // Build 5.1.012:
30 // - Height at max. width for Modified Baskethandle shape corrected.
31 // Build 5.1.013:
32 // - Width at full height set to 0 for closed rectangular shape.
33 // Build 5.2.0:
34 // - Support added for Street cross sections.
35 // Build 5.2.2:
36 // - Feasibility check added to Mod. Baskethandle & Rect.-Round shapes.
37 //-----------------------------------------------------------------------------
38 #define _CRT_SECURE_NO_DEPRECATE
39
40 #include <math.h>
41 #include "headers.h"
42 #include "findroot.h"
43
44 #define RECT_ALFMAX 0.97
45 #define RECT_TRIANG_ALFMAX 0.98
46 #define RECT_ROUND_ALFMAX 0.98
47
48 #include "xsect.dat" // File containing geometry tables for rounded shapes
49
50 //-----------------------------------------------------------------------------
51 // Constants
52 //-----------------------------------------------------------------------------
53 // Ratio of area at max. flow to full area
54 // (= 1.0 for open shapes, < 1.0 for closed shapes)
55 double Amax[] = {
56 1.0, // DUMMY
57 0.9756, // CIRCULAR
58 0.9756, // FILLED_CIRCULAR
59 0.97, // RECT_CLOSED
60 1.0, // RECT_OPEN
61 1.0, // TRAPEZOIDAL
62 1.0, // TRIANGULAR
63 1.0, // PARABOLIC
64 1.0, // POWERFUNC
65 0.98, // RECT_TRIANG
66 0.98, // RECT_ROUND
67 0.96, // MOD_BASKET
68 0.96, // HORIZ_ELLIPSE
69 0.96, // VERT_ELLIPSE
70 0.92, // ARCH
71 0.96, // EGGSHAPED
72 0.96, // HORSESHOE
73 0.96, // GOTHIC
74 0.98, // CATENARY
75 0.98, // SEMIELLIPTICAL
76 0.96, // BASKETHANDLE
77 0.96, // SEMICIRCULAR
78 1.0, // IRREGULAR
79 0.96, // CUSTOM
80 0.9756, // FORCE_MAIN
81 1.0}; // STREET_XSECT
82
83 //-----------------------------------------------------------------------------
84 // Shared variables
85 //-----------------------------------------------------------------------------
86 typedef struct
87 {
88 double s; // section factor
89 double qc; // critical flow
90 TXsect* xsect; // pointer to a cross section object
91 } TXsectStar;
92
93 //-----------------------------------------------------------------------------
94 // External functions (declared in funcs.h)
95 //-----------------------------------------------------------------------------
96 // xsect_isOpen
97 // xsect_setParams
98 // xsect_setIrregXsectParams
99 // xsect_setStreetXsectParams
100 // xsect_setCustomXsectParams
101 // xsect_getAmax
102 // xsect_getSofA
103 // xsect_getYofA
104 // xsect_getRofA
105 // xsect_getAofS
106 // xsect_getdSdA
107 // xsect_getAofY
108 // xsect_getRofY
109 // xsect_getWofY
110 // xsect_getYcrit
111
112 //-----------------------------------------------------------------------------
113 // Local functions
114 //-----------------------------------------------------------------------------
115 static void getTransectParams(TXsect *xsect, TTransect *transect);
116
117 static double generic_getAofS(TXsect* xsect, double s);
118 static void evalSofA(double a, double* f, double* df, void* p);
119 static double tabular_getdSdA(TXsect* xsect, double a, double *table, int nItems);
120 static double generic_getdSdA(TXsect* xsect, double a);
121 static double lookup(double x, double *table, int nItems);
122 static double invLookup(double y, double *table, int nItems);
123 static int locate(double y, double *table, int nItems);
124
125 static double rect_closed_getSofA(TXsect* xsect, double a);
126 static double rect_closed_getdSdA(TXsect* xsect, double a);
127 static double rect_closed_getRofA(TXsect* xsect, double a);
128
129 static double rect_open_getSofA(TXsect* xsect, double a);
130 static double rect_open_getdSdA(TXsect* xsect, double a);
131
132 static double rect_triang_getYofA(TXsect* xsect, double a);
133 static double rect_triang_getRofA(TXsect* xsect, double a);
134 static double rect_triang_getSofA(TXsect* xsect, double a);
135 static double rect_triang_getdSdA(TXsect* xsect, double a);
136 static double rect_triang_getAofY(TXsect* xsect, double y);
137 static double rect_triang_getRofY(TXsect* xsect, double y);
138 static double rect_triang_getWofY(TXsect* xsect, double y);
139
140 static double rect_round_getYofA(TXsect* xsect, double a);
141 static double rect_round_getRofA(TXsect* xsect, double a);
142 static double rect_round_getSofA(TXsect* xsect, double a);
143 static double rect_round_getdSdA(TXsect* xsect, double a);
144 static double rect_round_getAofY(TXsect* xsect, double y);
145 static double rect_round_getRofY(TXsect* xsect, double y);
146 static double rect_round_getWofY(TXsect* xsect, double y);
147
148 static double mod_basket_getYofA(TXsect* xsect, double a);
149 static double mod_basket_getRofA(TXsect* xsect, double a);
150 static double mod_basket_getdSdA(TXsect* xsect, double a);
151 static double mod_basket_getAofY(TXsect* xsect, double y);
152 static double mod_basket_getWofY(TXsect* xsect, double y);
153
154 static double trapez_getYofA(TXsect* xsect, double a);
155 static double trapez_getRofA(TXsect* xsect, double a);
156 static double trapez_getdSdA(TXsect* xsect, double a);
157 static double trapez_getAofY(TXsect* xsect, double y);
158 static double trapez_getRofY(TXsect* xsect, double y);
159 static double trapez_getWofY(TXsect* xsect, double y);
160
161 static double triang_getYofA(TXsect* xsect, double a);
162 static double triang_getRofA(TXsect* xsect, double a);
163 static double triang_getdSdA(TXsect* xsect, double a);
164 static double triang_getAofY(TXsect* xsect, double y);
165 static double triang_getRofY(TXsect* xsect, double y);
166 static double triang_getWofY(TXsect* xsect, double y);
167
168 static double parab_getYofA(TXsect* xsect, double a);
169 static double parab_getRofA(TXsect* xsect, double a);
170 static double parab_getPofY(TXsect* xsect, double y);
171 static double parab_getAofY(TXsect* xsect, double y);
172 static double parab_getRofY(TXsect* xsect, double y);
173 static double parab_getWofY(TXsect* xsect, double y);
174
175 static double powerfunc_getYofA(TXsect* xsect, double a);
176 static double powerfunc_getRofA(TXsect* xsect, double a);
177 static double powerfunc_getPofY(TXsect* xsect, double y);
178 static double powerfunc_getAofY(TXsect* xsect, double y);
179 static double powerfunc_getRofY(TXsect* xsect, double y);
180 static double powerfunc_getWofY(TXsect* xsect, double y);
181
182 static double circ_getYofA(TXsect* xsect, double a);
183 static double circ_getSofA(TXsect* xsect, double a);
184 static double circ_getdSdA(TXsect* xsect, double a);
185 static double circ_getAofS(TXsect* xsect, double s);
186 static double circ_getAofY(TXsect* xsect, double y);
187
188 static double filled_circ_getYofA(TXsect* xsect, double a);
189 static double filled_circ_getAofY(TXsect* xsect, double y);
190 static double filled_circ_getRofY(TXsect* xsect, double y);
191
192 static double getYcircular(double alpha);
193 static double getScircular(double alpha);
194 static double getAcircular(double psi);
195 static double getThetaOfAlpha(double alpha);
196 static double getThetaOfPsi(double psi);
197
198 static double getQcritical(double yc, void* p);
199 static double getYcritEnum(TXsect* xsect, double q, double y0);
200 static double getYcritRidder(TXsect* xsect, double q, double y0);
201
202 //=============================================================================
203
204 118367392 int xsect_isOpen(int type)
205 //
206 // Input: type = type of xsection shape
207 // Output: returns 1 if xsection is open, 0 if not
208 // Purpose: determines if a xsection type is open or closed.
209 //
210 {
211 118367392 return ((Amax[type] >= 1.0) ? 1 : 0);
212 }
213
214 //=============================================================================
215
216 10389 int xsect_setParams(TXsect *xsect, int type, double p[], double ucf)
217 //
218 // Input: xsect = ptr. to a cross section data structure
219 // type = xsection shape type
220 // p[] = vector of xsection parameters
221 // ucf = units correction factor
222 // Output: returns TRUE if successful, FALSE if not
223 // Purpose: assigns parameters to a cross section's data structure.
224 //
225 {
226 int index;
227 double aMax, theta;
228
229
3/4
✓ Branch 0 taken 10387 times.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 10387 times.
 10389 if ( type != DUMMY && p[0] <= 0.0 ) return FALSE;
230 10389 xsect->type = type;
231
16/24
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 8670 times.
✓ Branch 2 taken 134 times.
✓ Branch 3 taken 175 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 616 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✓ Branch 11 taken 196 times.
✓ Branch 12 taken 467 times.
✓ Branch 13 taken 14 times.
✓ Branch 14 taken 1 time.
✓ Branch 15 taken 64 times.
✓ Branch 16 taken 9 times.
✓ Branch 17 taken 15 times.
✓ Branch 18 taken 1 time.
✓ Branch 19 taken 1 time.
✓ Branch 20 taken 14 times.
✗ Branch 21 not taken.
✗ Branch 22 not taken.
✓ Branch 23 taken 10 times.
 10389 switch ( xsect->type )
232 {
233 2 case DUMMY:
234 2 xsect->yFull = TINY;
235 2 xsect->wMax = TINY;
236 2 xsect->aFull = TINY;
237 2 xsect->rFull = TINY;
238 2 xsect->sFull = TINY;
239 2 xsect->sMax = TINY;
240 2 break;
241
242 8670 case CIRCULAR:
243 8670 xsect->yFull = p[0]/ucf;
244 8670 xsect->wMax = xsect->yFull;
245 8670 xsect->aFull = PI / 4.0 * xsect->yFull * xsect->yFull;
246 8670 xsect->rFull = 0.2500 * xsect->yFull;
247 8670 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
248 8670 xsect->sMax = 1.08 * xsect->sFull;
249 8670 xsect->ywMax = 0.5 * xsect->yFull;
250 8670 break;
251
252 134 case FORCE_MAIN:
253 134 xsect->yFull = p[0]/ucf;
254 134 xsect->wMax = xsect->yFull;
255 134 xsect->aFull = PI / 4.0 * xsect->yFull * xsect->yFull;
256 134 xsect->rFull = 0.2500 * xsect->yFull;
257 134 xsect->sFull = xsect->aFull * pow(xsect->rFull, 0.63);
258 134 xsect->sMax = 1.06949 * xsect->sFull;
259 134 xsect->ywMax = 0.5 * xsect->yFull;
260
261 // --- save C-factor or roughness in rBot position
262 134 xsect->rBot = p[1];
263 134 break;
264
265 175 case FILLED_CIRCULAR:
266
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 175 times.
 175 if ( p[1] >= p[0] ) return FALSE;
267
268 // --- initially compute full values for unfilled pipe
269 175 xsect->yFull = p[0]/ucf;
270 175 xsect->wMax = xsect->yFull;
271 175 xsect->aFull = PI / 4.0 * xsect->yFull * xsect->yFull;
272 175 xsect->rFull = 0.2500 * xsect->yFull;
273
274 // --- find:
275 // yBot = depth of filled bottom
276 // aBot = area of filled bottom
277 // sBot = width of filled bottom
278 // rBot = wetted perimeter of filled bottom
279 175 xsect->yBot = p[1]/ucf;
280 175 xsect->aBot = circ_getAofY(xsect, xsect->yBot);
281 175 xsect->sBot = xsect_getWofY(xsect, xsect->yBot);
282 350 xsect->rBot = xsect->aBot / (xsect->rFull *
283 175 lookup(xsect->yBot/xsect->yFull, R_Circ, N_R_Circ));
284
285 // --- revise full values for filled bottom
286 175 xsect->aFull -= xsect->aBot;
287 175 xsect->rFull = xsect->aFull /
288 175 (PI*xsect->yFull - xsect->rBot + xsect->sBot);
289 175 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
290 175 xsect->sMax = 1.08 * xsect->sFull;
291 175 xsect->yFull -= xsect->yBot;
292 175 xsect->ywMax = 0.5 * xsect->yFull;
293 175 break;
294
295 ✗ case EGGSHAPED:
296 ✗ xsect->yFull = p[0]/ucf;
297 ✗ xsect->aFull = 0.5105 * xsect->yFull * xsect->yFull;
298 ✗ xsect->rFull = 0.1931 * xsect->yFull;
299 ✗ xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
300 ✗ xsect->sMax = 1.065 * xsect->sFull;
301 ✗ xsect->wMax = 2./3. * xsect->yFull;
302 ✗ xsect->ywMax = 0.64 * xsect->yFull;
303 ✗ break;
304
305 616 case HORSESHOE:
306 616 xsect->yFull = p[0]/ucf;
307 616 xsect->aFull = 0.8293 * xsect->yFull * xsect->yFull;
308 616 xsect->rFull = 0.2538 * xsect->yFull;
309 616 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
310 616 xsect->sMax = 1.077 * xsect->sFull;
311 616 xsect->wMax = 1.0 * xsect->yFull;
312 616 xsect->ywMax = 0.5 * xsect->yFull;
313 616 break;
314
315 ✗ case GOTHIC:
316 ✗ xsect->yFull = p[0]/ucf;
317 ✗ xsect->aFull = 0.6554 * xsect->yFull * xsect->yFull;
318 ✗ xsect->rFull = 0.2269 * xsect->yFull;
319 ✗ xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
320 ✗ xsect->sMax = 1.065 * xsect->sFull;
321 ✗ xsect->wMax = 0.84 * xsect->yFull;
322 ✗ xsect->ywMax = 0.45 * xsect->yFull;
323 ✗ break;
324
325 ✗ case CATENARY:
326 ✗ xsect->yFull = p[0]/ucf;
327 ✗ xsect->aFull = 0.70277 * xsect->yFull * xsect->yFull;
328 ✗ xsect->rFull = 0.23172 * xsect->yFull;
329 ✗ xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
330 ✗ xsect->sMax = 1.05 * xsect->sFull;
331 ✗ xsect->wMax = 0.9 * xsect->yFull;
332 ✗ xsect->ywMax = 0.25 * xsect->yFull;
333 ✗ break;
334
335 ✗ case SEMIELLIPTICAL:
336 ✗ xsect->yFull = p[0]/ucf;
337 ✗ xsect->aFull = 0.785 * xsect->yFull * xsect->yFull;
338 ✗ xsect->rFull = 0.242 * xsect->yFull;
339 ✗ xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
340 ✗ xsect->sMax = 1.045 * xsect->sFull;
341 ✗ xsect->wMax = 1.0 * xsect->yFull;
342 ✗ xsect->ywMax = 0.15 * xsect->yFull;
343 ✗ break;
344
345 ✗ case BASKETHANDLE:
346 ✗ xsect->yFull = p[0]/ucf;
347 ✗ xsect->aFull = 0.7862 * xsect->yFull * xsect->yFull;
348 ✗ xsect->rFull = 0.2464 * xsect->yFull;
349 ✗ xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
350 ✗ xsect->sMax = 1.06078 * xsect->sFull;
351 ✗ xsect->wMax = 0.944 * xsect->yFull;
352 ✗ xsect->ywMax = 0.2 * xsect->yFull;
353 ✗ break;
354
355 ✗ case SEMICIRCULAR:
356 ✗ xsect->yFull = p[0]/ucf;
357 ✗ xsect->aFull = 1.2697 * xsect->yFull * xsect->yFull;
358 ✗ xsect->rFull = 0.2946 * xsect->yFull;
359 ✗ xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
360 ✗ xsect->sMax = 1.06637 * xsect->sFull;
361 ✗ xsect->wMax = 1.64 * xsect->yFull;
362 ✗ xsect->ywMax = 0.15 * xsect->yFull;
363 ✗ break;
364
365 196 case RECT_CLOSED:
366
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 196 times.
 196 if ( p[1] <= 0.0 ) return FALSE;
367 196 xsect->yFull = p[0]/ucf;
368 196 xsect->wMax = p[1]/ucf;
369 196 xsect->aFull = xsect->yFull * xsect->wMax;
370 196 xsect->rFull = xsect->aFull / (2.0 * (xsect->yFull + xsect->wMax));
371 196 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
372 196 aMax = RECT_ALFMAX * xsect->aFull;
373 196 xsect->sMax = aMax * pow(rect_closed_getRofA(xsect, aMax), 2./3.);
374 196 xsect->ywMax = xsect->yFull;
375 196 break;
376
377 467 case RECT_OPEN:
378
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 467 times.
 467 if ( p[1] <= 0.0 ) return FALSE;
379 467 xsect->yFull = p[0]/ucf;
380 467 xsect->wMax = p[1]/ucf;
381
2/4
✓ Branch 0 taken 467 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 467 times.
 467 if (p[2] < 0.0 || p[2] > 2.0) return FALSE; //# sides to ignore
382 467 xsect->sBot = p[2];
383 467 xsect->aFull = xsect->yFull * xsect->wMax;
384 467 xsect->rFull = xsect->aFull / ((2.0 - xsect->sBot) *
385 467 xsect->yFull + xsect->wMax);
386 467 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
387 467 xsect->sMax = xsect->sFull;
388 467 xsect->ywMax = xsect->yFull;
389 467 break;
390
391 14 case RECT_TRIANG:
392
2/4
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 14 times.
 14 if ( p[1] <= 0.0 || p[2] <= 0.0 ) return FALSE;
393 14 xsect->yFull = p[0]/ucf;
394 14 xsect->wMax = p[1]/ucf;
395 14 xsect->yBot = p[2]/ucf;
396 14 xsect->ywMax = xsect->yFull;
397
398 // --- area of bottom triangle
399 14 xsect->aBot = xsect->yBot * xsect->wMax / 2.0;
400
401 // --- slope of bottom side wall
402 14 xsect->sBot = xsect->wMax / xsect->yBot / 2.0;
403
404 // --- length of side wall per unit of depth
405 14 xsect->rBot = sqrt( 1. + xsect->sBot * xsect->sBot );
406
407 14 xsect->aFull = xsect->wMax * (xsect->yFull - xsect->yBot / 2.0);
408 14 xsect->rFull = xsect->aFull / (2.0 * xsect->yBot * xsect->rBot + 2.0 *
409 14 (xsect->yFull - xsect->yBot) + xsect->wMax);
410 14 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
411 14 aMax = RECT_TRIANG_ALFMAX * xsect->aFull;
412 14 xsect->sMax = aMax * pow(rect_triang_getRofA(xsect, aMax), 2./3.);
413 14 break;
414
415 1 case RECT_ROUND:
416
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( p[1] <= 0.0 ) return FALSE;
417
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( p[2] < p[1]/2.0 ) p[2] = p[1]/2.0;
418 1 xsect->yFull = p[0]/ucf;
419 1 xsect->wMax = p[1]/ucf;
420 1 xsect->rBot = p[2]/ucf;
421
422 // --- angle of circular arc
423 1 theta = 2.0 * asin(xsect->wMax / 2.0 / xsect->rBot);
424
425 // --- area of circular bottom
426 1 xsect->aBot = xsect->rBot * xsect->rBot /
427 1 2.0 * (theta - sin(theta));
428
429 // --- section factor for circular bottom
430 1 xsect->sBot = PI * xsect->rBot * xsect->rBot *
431 1 pow(xsect->rBot/2.0, 2./3.);
432
433 // --- depth of circular bottom
434 1 xsect->yBot = xsect->rBot * (1.0 - cos(theta/2.0));
435
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if (xsect->yBot > xsect->yFull) return FALSE;
436 1 xsect->ywMax = xsect->yFull;
437
438 1 xsect->aFull = xsect->wMax * (xsect->yFull - xsect->yBot) + xsect->aBot;
439 1 xsect->rFull = xsect->aFull / (xsect->rBot * theta + 2.0 *
440 1 (xsect->yFull - xsect->yBot) + xsect->wMax);
441 1 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
442 1 aMax = RECT_ROUND_ALFMAX * xsect->aFull;
443 1 xsect->sMax = aMax * pow(rect_round_getRofA(xsect, aMax), 2./3.);
444 1 break;
445
446 64 case MOD_BASKET:
447
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 64 times.
 64 if ( p[1] <= 0.0 ) return FALSE;
448
2/2
✓ Branch 0 taken 63 times.
✓ Branch 1 taken 1 time.
 64 if ( p[2] < p[1]/2.0 ) p[2] = p[1]/2.0;
449 64 xsect->yFull = p[0]/ucf;
450 64 xsect->wMax = p[1]/ucf;
451
452 // --- radius of circular arc
453 64 xsect->rBot = p[2]/ucf;
454
455 // --- angle of circular arc
456 64 theta = 2.0 * asin(xsect->wMax / 2.0 / xsect->rBot);
457 64 xsect->sBot = theta;
458
459 // --- height of circular arc
460 64 xsect->yBot = xsect->rBot * (1.0 - cos(theta/2.0));
461
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 64 times.
 64 if (xsect->yBot > xsect->yFull) return FALSE;
462 64 xsect->ywMax = xsect->yFull - xsect->yBot;
463
464 // --- area of circular arc
465 64 xsect->aBot = xsect->rBot * xsect->rBot /
466 64 2.0 * (theta - sin(theta));
467
468 // --- full area
469 64 xsect->aFull = (xsect->yFull - xsect->yBot) * xsect->wMax +
470 64 xsect->aBot;
471
472 // --- full hydraulic radius & section factor
473 64 xsect->rFull = xsect->aFull / (xsect->rBot * theta + 2.0 *
474 64 (xsect->yFull - xsect->yBot) + xsect->wMax);
475 64 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
476
477 // --- area corresponding to max. section factor
478 64 xsect->sMax = xsect_getSofA(xsect, Amax[MOD_BASKET]*xsect->aFull);
479 64 break;
480
481 9 case TRAPEZOIDAL:
482
3/6
✓ Branch 0 taken 9 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 9 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 9 times.
 9 if ( p[1] < 0.0 || p[2] < 0.0 || p[3] < 0.0 ) return FALSE;
483 9 xsect->yFull = p[0]/ucf;
484 9 xsect->ywMax = xsect->yFull;
485
486 // --- bottom width
487 9 xsect->yBot = p[1]/ucf;
488
489 // --- avg. slope of side walls
490 9 xsect->sBot = ( p[2] + p[3] )/2.0;
491
3/4
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
 9 if ( xsect->yBot == 0.0 && xsect->sBot == 0.0 ) return FALSE;
492
493 // --- length of side walls per unit of depth
494 9 xsect->rBot = sqrt( 1.0 + p[2]*p[2] ) + sqrt( 1.0 + p[3]*p[3] );
495
496 // --- top width
497 9 xsect->wMax = xsect->yBot + xsect->yFull * (p[2] + p[3]);
498
499 9 xsect->aFull = ( xsect->yBot + xsect->sBot * xsect->yFull ) * xsect->yFull;
500 9 xsect->rFull = xsect->aFull / (xsect->yBot + xsect->yFull * xsect->rBot);
501 9 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
502 9 xsect->sMax = xsect->sFull;
503 9 break;
504
505 15 case TRIANGULAR:
506
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
 15 if ( p[1] <= 0.0 ) return FALSE;
507 15 xsect->yFull = p[0]/ucf;
508 15 xsect->wMax = p[1]/ucf;
509 15 xsect->ywMax = xsect->yFull;
510
511 // --- slope of side walls
512 15 xsect->sBot = xsect->wMax / xsect->yFull / 2.;
513
514 // --- length of side wall per unit of depth
515 15 xsect->rBot = sqrt( 1. + xsect->sBot * xsect->sBot );
516
517 15 xsect->aFull = xsect->yFull * xsect->yFull * xsect->sBot;
518 15 xsect->rFull = xsect->aFull / (2.0 * xsect->yFull * xsect->rBot);
519 15 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
520 15 xsect->sMax = xsect->sFull;
521 15 break;
522
523 1 case PARABOLIC:
524
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( p[1] <= 0.0 ) return FALSE;
525 1 xsect->yFull = p[0]/ucf;
526 1 xsect->wMax = p[1]/ucf;
527 1 xsect->ywMax = xsect->yFull;
528
529 // --- rBot :: 1/c^.5, where y = c*x^2 is eqn. of parabolic shape
530 1 xsect->rBot = xsect->wMax / 2.0 / sqrt(xsect->yFull);
531
532 1 xsect->aFull = (2./3.) * xsect->yFull * xsect->wMax;
533 1 xsect->rFull = xsect_getRofY(xsect, xsect->yFull);
534 1 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
535 1 xsect->sMax = xsect->sFull;
536 1 break;
537
538 1 case POWERFUNC:
539
2/4
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
 1 if ( p[1] <= 0.0 || p[2] <= 0.0 ) return FALSE;
540 1 xsect->yFull = p[0]/ucf;
541 1 xsect->wMax = p[1]/ucf;
542 1 xsect->ywMax = xsect->yFull;
543 1 xsect->sBot = 1.0 / p[2];
544 1 xsect->rBot = xsect->wMax / (xsect->sBot + 1) /
545 1 pow(xsect->yFull, xsect->sBot);
546 1 xsect->aFull = xsect->yFull * xsect->wMax / (xsect->sBot+1);
547 1 xsect->rFull = xsect_getRofY(xsect, xsect->yFull);
548 1 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
549 1 xsect->sMax = xsect->sFull;
550 1 break;
551
552 14 case HORIZ_ELLIPSE:
553
1/2
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
 14 if ( p[1] == 0.0 ) p[2] = p[0];
554
1/2
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
 14 if ( p[2] > 0.0 ) // std. ellipse pipe
555 {
556 14 index = (int)floor(p[2]) - 1; // size code
557
1/2
✓ Branch 0 taken 14 times.
✗ Branch 1 not taken.
 14 if ( index < 0 ||
558
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14 times.
 14 index >= NumCodesEllipse ) return FALSE;
559 14 xsect->yFull = MinorAxis_Ellipse[index]/12.;
560 14 xsect->wMax = MajorAxis_Ellipse[index]/12.;
561 14 xsect->aFull = Afull_Ellipse[index];
562 14 xsect->rFull = Rfull_Ellipse[index];
563 }
564 else
565 {
566 // --- length of minor axis
567 ✗ xsect->yFull = p[0]/ucf;
568
569 // --- length of major axis
570 ✗ if ( p[1] < 0.0 ) return FALSE;
571 ✗ xsect->wMax = p[1]/ucf;
572 ✗ xsect->aFull = 1.2692 * xsect->yFull * xsect->yFull;
573 ✗ xsect->rFull = 0.3061 * xsect->yFull;
574 }
575 14 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
576 14 xsect->sMax = xsect->sFull;
577 14 xsect->ywMax = 0.48 * xsect->yFull;
578 14 break;
579
580 ✗ case VERT_ELLIPSE:
581 ✗ if ( p[1] == 0.0 ) p[2] = p[0];
582 ✗ if ( p[2] > 0.0 ) // std. ellipse pipe
583 {
584 ✗ index = (int)floor(p[2]) - 1; // size code
585 ✗ if ( index < 0 ||
586 ✗ index >= NumCodesEllipse ) return FALSE;
587 ✗ xsect->yFull = MajorAxis_Ellipse[index]/12.;
588 ✗ xsect->wMax = MinorAxis_Ellipse[index]/12.;
589 ✗ xsect->aFull = Afull_Ellipse[index];
590 ✗ xsect->rFull = Rfull_Ellipse[index];
591 }
592 else
593 {
594 // --- length of major axis
595 ✗ if ( p[1] < 0.0 ) return FALSE;
596
597 // --- length of minor axis
598 ✗ xsect->yFull = p[0]/ucf;
599 ✗ xsect->wMax = p[1]/ucf;
600 ✗ xsect->aFull = 1.2692 * xsect->wMax * xsect->wMax;
601 ✗ xsect->rFull = 0.3061 * xsect->wMax;
602 }
603 ✗ xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
604 ✗ xsect->sMax = xsect->sFull;
605 ✗ xsect->ywMax = 0.48 * xsect->yFull;
606 ✗ break;
607
608 ✗ case ARCH:
609 ✗ if ( p[1] == 0.0 ) p[2] = p[0];
610 ✗ if ( p[2] > 0.0 ) // std. arch pipe
611 {
612 ✗ index = (int)floor(p[2]) - 1; // size code
613 ✗ if ( index < 0 ||
614 ✗ index >= NumCodesArch ) return FALSE;
615 ✗ xsect->yFull = Yfull_Arch[index]/12.; // Yfull units are inches
616 ✗ xsect->wMax = Wmax_Arch[index]/12.; // Wmax units are inches
617 ✗ xsect->aFull = Afull_Arch[index];
618 ✗ xsect->rFull = Rfull_Arch[index];
619 }
620 else // non-std. arch pipe
621 {
622 ✗ if ( p[1] < 0.0 ) return FALSE;
623 ✗ xsect->yFull = p[0]/ucf;
624 ✗ xsect->wMax = p[1]/ucf;
625 ✗ xsect->aFull = 0.7879 * xsect->yFull * xsect->wMax;
626 ✗ xsect->rFull = 0.2991 * xsect->yFull;
627 }
628 ✗ xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
629 ✗ xsect->sMax = xsect->sFull;
630 ✗ xsect->ywMax = 0.28 * xsect->yFull;
631 ✗ break;
632 }
633 10389 return TRUE;
634 }
635
636 //=============================================================================
637
638 1 void xsect_setIrregXsectParams(TXsect *xsect)
639 //
640 // Input: xsect = ptr. to a cross section data structure
641 // Output: none
642 // Purpose: assigns transect parameters to an irregular shaped cross section.
643 //
644 {
645 1 int index = xsect->transect;
646 1 getTransectParams(xsect, &Transect[index]);
647 1 }
648
649 //=============================================================================
650
651 12 void xsect_setStreetXsectParams(TXsect *xsect)
652 //
653 // Input: xsect = ptr. to a cross section data structure
654 // Output: none
655 // Purpose: assigns transect parameters to a street cross section.
656 //
657 {
658 12 int index = xsect->transect;
659 12 getTransectParams(xsect, &Street[index].transect);
660 12 }
661
662 //=============================================================================
663
664 10 void xsect_setCustomXsectParams(TXsect *xsect)
665 //
666 // Input: xsect = ptr. to a cross section data structure
667 // Output: none
668 // Purpose: assigns parameters to a custom-shaped cross section.
669 //
670 {
671 10 int index = Curve[xsect->transect].refersTo;
672 10 double yFull = xsect->yFull;
673 int i, iMax;
674 double wMax;
675 10 double* wTbl = Shape[index].widthTbl;
676
677 10 xsect->wMax = Shape[index].wMax * yFull;
678 10 xsect->aFull = Shape[index].aFull * yFull * yFull;
679 10 xsect->rFull = Shape[index].rFull * yFull;
680 10 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2./3.);
681 10 xsect->sMax = Shape[index].sMax * yFull * yFull * pow(yFull, 2./3.);
682 10 xsect->aBot = Shape[index].aMax * yFull * yFull;
683
684 // Search shape's width table up to point where width decreases
685 10 iMax = 0;
686 10 wMax = wTbl[0];
687
2/2
✓ Branch 0 taken 382 times.
✓ Branch 1 taken 1 time.
 383 for (i = 1; i < N_SHAPE_TBL; i++)
688 {
689
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 373 times.
 382 if ( wTbl[i] < wMax ) break;
690 373 wMax = wTbl[i];
691 373 iMax = i;
692 }
693
694 // Determine height at lowest widest point
695 10 xsect->ywMax = yFull * (double)iMax / (double)(N_SHAPE_TBL-1);
696 10 }
697
698 //=============================================================================
699
700 1833452 double xsect_getAmax(TXsect* xsect)
701 //
702 // Input: xsect = ptr. to a cross section data structure
703 // Output: returns area (ft2)
704 // Purpose: finds xsection area at maximum flow depth.
705 //
706 {
707
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1833452 times.
 1833452 if ( xsect->type == IRREGULAR ) return xsect->aBot;
708
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1833452 times.
 1833452 else if ( xsect->type == CUSTOM ) return xsect->aBot;
709 1833452 else return Amax[xsect->type] * xsect->aFull;
710 }
711
712 //=============================================================================
713
714 12862513 double xsect_getSofA(TXsect *xsect, double a)
715 //
716 // Input: xsect = ptr. to a cross section data structure
717 // a = area (ft2)
718 // Output: returns section factor (ft^(8/3))
719 // Purpose: computes xsection's section factor at a given area.
720 //
721 {
722 12862513 double alpha = a / xsect->aFull;
723 double r;
724
6/13
✓ Branch 0 taken 145424 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 43707 times.
✓ Branch 9 taken 69877 times.
✓ Branch 10 taken 11049 times.
✓ Branch 11 taken 247515 times.
✓ Branch 12 taken 12344941 times.
 12862513 switch ( xsect->type )
725 {
726 145424 case FORCE_MAIN:
727 case CIRCULAR:
728 145424 return circ_getSofA(xsect, a);
729
730 ✗ case EGGSHAPED:
731 ✗ return xsect->sFull * lookup(alpha, S_Egg, N_S_Egg);
732
733 ✗ case HORSESHOE:
734 ✗ return xsect->sFull * lookup(alpha, S_Horseshoe, N_S_Horseshoe);
735
736 ✗ case GOTHIC:
737 ✗ return xsect->sFull * lookup(alpha, S_Gothic, N_S_Gothic);
738
739 ✗ case CATENARY:
740 ✗ return xsect->sFull * lookup(alpha, S_Catenary, N_S_Catenary);
741
742 ✗ case SEMIELLIPTICAL:
743 ✗ return xsect->sFull * lookup(alpha, S_SemiEllip, N_S_SemiEllip);
744
745 ✗ case BASKETHANDLE:
746 ✗ return xsect->sFull * lookup(alpha, S_BasketHandle, N_S_BasketHandle);
747
748 ✗ case SEMICIRCULAR:
749 ✗ return xsect->sFull * lookup(alpha, S_SemiCirc, N_S_SemiCirc);
750
751 43707 case RECT_CLOSED:
752 43707 return rect_closed_getSofA(xsect, a);
753
754 69877 case RECT_OPEN:
755 69877 return rect_open_getSofA(xsect, a);
756
757 11049 case RECT_TRIANG:
758 11049 return rect_triang_getSofA(xsect, a);
759
760 247515 case RECT_ROUND:
761 247515 return rect_round_getSofA(xsect, a);
762
763 12344941 default:
764
2/2
✓ Branch 0 taken 5297 times.
✓ Branch 1 taken 12339644 times.
 12344941 if (a == 0.0) return 0.0;
765 12339644 r = xsect_getRofA(xsect, a);
766
2/2
✓ Branch 0 taken 34 times.
✓ Branch 1 taken 12339610 times.
 12339644 if ( r < TINY ) return 0.0;
767 12339610 return a * pow(r, 2./3.);
768 }
769 }
770
771 //=============================================================================
772
773 32559862 double xsect_getYofA(TXsect *xsect, double a)
774 //
775 // Input: xsect = ptr. to a cross section data structure
776 // a = area (ft2)
777 // Output: returns depth (ft)
778 // Purpose: computes xsection's depth at a given area.
779 //
780 {
781 32559862 double alpha = a / xsect->aFull;
782
12/25
✓ Branch 0 taken 25584334 times.
✓ Branch 1 taken 5617031 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2997 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 12 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✓ Branch 15 taken 6134 times.
✓ Branch 16 taken 1545 times.
✓ Branch 17 taken 14402 times.
✓ Branch 18 taken 14402 times.
✓ Branch 19 taken 14402 times.
✓ Branch 20 taken 356068 times.
✓ Branch 21 taken 919743 times.
✓ Branch 22 taken 14402 times.
✓ Branch 23 taken 14402 times.
✗ Branch 24 not taken.
 32559862 switch ( xsect->type )
783 {
784 25584334 case FORCE_MAIN:
785 25584334 case CIRCULAR: return circ_getYofA(xsect, a);
786
787 5617031 case FILLED_CIRCULAR:
788 5617031 return filled_circ_getYofA(xsect, a);
789
790 ✗ case EGGSHAPED:
791 ✗ return xsect->yFull * lookup(alpha, Y_Egg, N_Y_Egg);
792
793 2997 case HORSESHOE:
794 2997 return xsect->yFull * lookup(alpha, Y_Horseshoe, N_Y_Horseshoe);
795
796 ✗ case GOTHIC:
797 ✗ return xsect->yFull * lookup(alpha, Y_Gothic, N_Y_Gothic);
798
799 ✗ case CATENARY:
800 ✗ return xsect->yFull * lookup(alpha, Y_Catenary, N_Y_Catenary);
801
802 ✗ case SEMIELLIPTICAL:
803 ✗ return xsect->yFull * lookup(alpha, Y_SemiEllip, N_Y_SemiEllip);
804
805 ✗ case BASKETHANDLE:
806 ✗ return xsect->yFull * lookup(alpha, Y_BasketHandle, N_Y_BasketHandle);
807
808 ✗ case SEMICIRCULAR:
809 ✗ return xsect->yFull * lookup(alpha, Y_SemiCirc, N_Y_SemiCirc);
810
811 ✗ case HORIZ_ELLIPSE:
812 ✗ return xsect->yFull * invLookup(alpha, A_HorizEllipse, N_A_HorizEllipse);
813
814 ✗ case VERT_ELLIPSE:
815 ✗ return xsect->yFull * invLookup(alpha, A_VertEllipse, N_A_VertEllipse);
816
817 ✗ case IRREGULAR:
818 ✗ return xsect->yFull * invLookup(alpha,
819 ✗ Transect[xsect->transect].areaTbl, N_TRANSECT_TBL);
820
821 ✗ case CUSTOM:
822 ✗ return xsect->yFull * invLookup(alpha,
823 ✗ Shape[Curve[xsect->transect].refersTo].areaTbl, N_SHAPE_TBL);
824
825 ✗ case STREET_XSECT:
826 ✗ return xsect->yFull * invLookup(alpha,
827 ✗ Street[xsect->transect].transect.areaTbl,
828 ✗ Street[xsect->transect].transect.nTbl);
829
830 ✗ case ARCH:
831 ✗ return xsect->yFull * invLookup(alpha, A_Arch, N_A_Arch);
832
833 6134 case RECT_CLOSED: return a / xsect->wMax;
834
835 1545 case RECT_TRIANG: return rect_triang_getYofA(xsect, a);
836
837 14402 case RECT_ROUND: return rect_round_getYofA(xsect, a);
838
839 14402 case RECT_OPEN: return a / xsect->wMax;
840
841 14402 case MOD_BASKET: return mod_basket_getYofA(xsect, a);
842
843 356068 case TRAPEZOIDAL: return trapez_getYofA(xsect, a);
844
845 919743 case TRIANGULAR: return triang_getYofA(xsect, a);
846
847 14402 case PARABOLIC: return parab_getYofA(xsect, a);
848
849 14402 case POWERFUNC: return powerfunc_getYofA(xsect, a);
850
851 ✗ default: return 0.0;
852 }
853 }
854
855 //=============================================================================
856
857 243024017 double xsect_getAofY(TXsect *xsect, double y)
858 //
859 // Input: xsect = ptr. to a cross section data structure
860 // y = depth (ft)
861 // Output: returns area (ft2)
862 // Purpose: computes xsection's area at a given depth.
863 //
864 {
865 243024017 double yNorm = y / xsect->yFull;
866
2/2
✓ Branch 0 taken 1912503 times.
✓ Branch 1 taken 241111514 times.
 243024017 if ( y <= 0.0 ) return 0.0;
867
16/25
✓ Branch 0 taken 225640498 times.
✓ Branch 1 taken 3320060 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 339212 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 7573 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 13385 times.
✓ Branch 13 taken 87362 times.
✓ Branch 14 taken 41187 times.
✓ Branch 15 taken 137300 times.
✓ Branch 16 taken 10991 times.
✓ Branch 17 taken 95561 times.
✓ Branch 18 taken 93662 times.
✓ Branch 19 taken 205048 times.
✓ Branch 20 taken 6027742 times.
✓ Branch 21 taken 4962217 times.
✓ Branch 22 taken 64858 times.
✓ Branch 23 taken 64858 times.
✗ Branch 24 not taken.
 241111514 switch ( xsect->type )
868 {
869 225640498 case FORCE_MAIN:
870 case CIRCULAR:
871 225640498 return xsect->aFull * lookup(yNorm, A_Circ, N_A_Circ);
872
873 3320060 case FILLED_CIRCULAR:
874 3320060 return filled_circ_getAofY(xsect, y);
875
876 ✗ case EGGSHAPED:
877 ✗ return xsect->aFull * lookup(yNorm, A_Egg, N_A_Egg);
878
879 339212 case HORSESHOE:
880 339212 return xsect->aFull * lookup(yNorm, A_Horseshoe, N_A_Horseshoe);
881
882 ✗ case GOTHIC:
883 ✗ return xsect->aFull * invLookup(yNorm, Y_Gothic, N_Y_Gothic);
884
885 ✗ case CATENARY:
886 ✗ return xsect->aFull * invLookup(yNorm, Y_Catenary, N_Y_Catenary);
887
888 ✗ case SEMIELLIPTICAL:
889 ✗ return xsect->aFull * invLookup(yNorm, Y_SemiEllip, N_Y_SemiEllip);
890
891 ✗ case BASKETHANDLE:
892 ✗ return xsect->aFull * lookup(yNorm, A_Baskethandle, N_A_Baskethandle);
893
894 ✗ case SEMICIRCULAR:
895 ✗ return xsect->aFull * invLookup(yNorm, Y_SemiCirc, N_Y_SemiCirc);
896
897 7573 case HORIZ_ELLIPSE:
898 7573 return xsect->aFull * lookup(yNorm, A_HorizEllipse, N_A_HorizEllipse);
899
900 ✗ case VERT_ELLIPSE:
901 ✗ return xsect->aFull * lookup(yNorm, A_VertEllipse, N_A_VertEllipse);
902
903 ✗ case ARCH:
904 ✗ return xsect->aFull * lookup(yNorm, A_Arch, N_A_Arch);
905
906 13385 case IRREGULAR:
907 13385 return xsect->aFull * lookup(yNorm,
908 13385 Transect[xsect->transect].areaTbl, N_TRANSECT_TBL);
909
910 87362 case CUSTOM:
911 87362 return xsect->aFull * lookup(yNorm,
912 87362 Shape[Curve[xsect->transect].refersTo].areaTbl, N_SHAPE_TBL);
913
914 41187 case STREET_XSECT:
915 41187 return xsect->aFull * lookup(yNorm,
916 41187 Street[xsect->transect].transect.areaTbl,
917 41187 Street[xsect->transect].transect.nTbl);
918
919 137300 case RECT_CLOSED: return y * xsect->wMax;
920
921 10991 case RECT_TRIANG: return rect_triang_getAofY(xsect, y);
922
923 95561 case RECT_ROUND: return rect_round_getAofY(xsect, y);
924
925 93662 case RECT_OPEN: return y * xsect->wMax;
926
927 205048 case MOD_BASKET: return mod_basket_getAofY(xsect, y);
928
929 6027742 case TRAPEZOIDAL: return trapez_getAofY(xsect, y);
930
931 4962217 case TRIANGULAR: return triang_getAofY(xsect, y);
932
933 64858 case PARABOLIC: return parab_getAofY(xsect, y);
934
935 64858 case POWERFUNC: return powerfunc_getAofY(xsect, y);
936
937 ✗ default: return 0.0;
938 }
939 }
940
941 //=============================================================================
942
943 219835684 double xsect_getWofY(TXsect *xsect, double y)
944 //
945 // Input: xsect = ptr. to a cross section data structure
946 // y = depth ft)
947 // Output: returns top width (ft)
948 // Purpose: computes xsection's top width at a given depth.
949 //
950 {
951 219835684 double yNorm = y / xsect->yFull;
952
16/25
✓ Branch 0 taken 203757100 times.
✓ Branch 1 taken 2938068 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 273670 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 6551 times.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✓ Branch 12 taken 13034 times.
✓ Branch 13 taken 105837 times.
✓ Branch 14 taken 59681 times.
✓ Branch 15 taken 81521 times.
✓ Branch 16 taken 9276 times.
✓ Branch 17 taken 88596 times.
✓ Branch 18 taken 190675 times.
✓ Branch 19 taken 151099 times.
✓ Branch 20 taken 6016695 times.
✓ Branch 21 taken 5999861 times.
✓ Branch 22 taken 72010 times.
✓ Branch 23 taken 72010 times.
✗ Branch 24 not taken.
 219835684 switch ( xsect->type )
953 {
954 203757100 case FORCE_MAIN:
955 case CIRCULAR:
956 203757100 return xsect->wMax * lookup(yNorm, W_Circ, N_W_Circ);
957
958 2938068 case FILLED_CIRCULAR:
959 2938068 yNorm = (y + xsect->yBot) / (xsect->yFull + xsect->yBot);
960 2938068 return xsect->wMax * lookup(yNorm, W_Circ, N_W_Circ);
961
962 ✗ case EGGSHAPED:
963 ✗ return xsect->wMax * lookup(yNorm, W_Egg, N_W_Egg);
964
965 273670 case HORSESHOE:
966 273670 return xsect->wMax * lookup(yNorm, W_Horseshoe, N_W_Horseshoe);
967
968 ✗ case GOTHIC:
969 ✗ return xsect->wMax * lookup(yNorm, W_Gothic, N_W_Gothic);
970
971 ✗ case CATENARY:
972 ✗ return xsect->wMax * lookup(yNorm, W_Catenary, N_W_Catenary);
973
974 ✗ case SEMIELLIPTICAL:
975 ✗ return xsect->wMax * lookup(yNorm, W_SemiEllip, N_W_SemiEllip);
976
977 ✗ case BASKETHANDLE:
978 ✗ return xsect->wMax * lookup(yNorm, W_BasketHandle, N_W_BasketHandle);
979
980 ✗ case SEMICIRCULAR:
981 ✗ return xsect->wMax * lookup(yNorm, W_SemiCirc, N_W_SemiCirc);
982
983 6551 case HORIZ_ELLIPSE:
984 6551 return xsect->wMax * lookup(yNorm, W_HorizEllipse, N_W_HorizEllipse);
985
986 ✗ case VERT_ELLIPSE:
987 ✗ return xsect->wMax * lookup(yNorm, W_VertEllipse, N_W_VertEllipse);
988
989 ✗ case ARCH:
990 ✗ return xsect->wMax * lookup(yNorm, W_Arch, N_W_Arch);
991
992 13034 case IRREGULAR:
993 13034 return xsect->wMax * lookup(yNorm,
994 13034 Transect[xsect->transect].widthTbl, N_TRANSECT_TBL);
995
996 105837 case CUSTOM:
997 105837 return xsect->wMax * lookup(yNorm,
998 105837 Shape[Curve[xsect->transect].refersTo].widthTbl, N_SHAPE_TBL);
999
1000 59681 case STREET_XSECT:
1001 59681 return xsect->wMax * lookup(yNorm,
1002 59681 Street[xsect->transect].transect.widthTbl,
1003 59681 Street[xsect->transect].transect.nTbl);
1004
1005 81521 case RECT_CLOSED:
1006
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 81521 times.
 81521 if (yNorm == 1.0) return 0.0;
1007 81521 return xsect->wMax;
1008
1009 9276 case RECT_TRIANG: return rect_triang_getWofY(xsect, y);
1010
1011 88596 case RECT_ROUND: return rect_round_getWofY(xsect, y);
1012
1013 190675 case RECT_OPEN: return xsect->wMax;
1014
1015 151099 case MOD_BASKET: return mod_basket_getWofY(xsect, y);
1016
1017 6016695 case TRAPEZOIDAL: return trapez_getWofY(xsect, y);
1018
1019 5999861 case TRIANGULAR: return triang_getWofY(xsect, y);
1020
1021 72010 case PARABOLIC: return parab_getWofY(xsect, y);
1022
1023 72010 case POWERFUNC: return powerfunc_getWofY(xsect, y);
1024
1025 ✗ default: return 0.0;
1026 }
1027 }
1028
1029 //=============================================================================
1030
1031 104432339 double xsect_getRofY(TXsect *xsect, double y)
1032 //
1033 // Input: xsect = ptr. to a cross section data structure
1034 // y = depth (ft)
1035 // Output: returns hydraulic radius (ft)
1036 // Purpose: computes xsection's hydraulic radius at a given depth.
1037 //
1038 {
1039 104432339 double yNorm = y / xsect->yFull;
1040
14/18
✓ Branch 0 taken 93970107 times.
✓ Branch 1 taken 6187077 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 150014 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 3380 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 4678 times.
✓ Branch 9 taken 48974 times.
✓ Branch 10 taken 13390 times.
✓ Branch 11 taken 4782 times.
✓ Branch 12 taken 28804 times.
✓ Branch 13 taken 1438399 times.
✓ Branch 14 taken 2406740 times.
✓ Branch 15 taken 28805 times.
✓ Branch 16 taken 28805 times.
✓ Branch 17 taken 118384 times.
 104432339 switch ( xsect->type )
1041 {
1042 93970107 case FORCE_MAIN:
1043 case CIRCULAR:
1044 93970107 return xsect->rFull * lookup(yNorm, R_Circ, N_R_Circ);
1045
1046 6187077 case FILLED_CIRCULAR:
1047
2/2
✓ Branch 0 taken 48484 times.
✓ Branch 1 taken 6138593 times.
 6187077 if ( xsect->yBot == 0.0 )
1048 48484 return xsect->rFull * lookup(yNorm, R_Circ, N_R_Circ);
1049 6138593 return filled_circ_getRofY(xsect, y);
1050
1051 ✗ case EGGSHAPED:
1052 ✗ return xsect->rFull * lookup(yNorm, R_Egg, N_R_Egg);
1053
1054 150014 case HORSESHOE:
1055 150014 return xsect->rFull * lookup(yNorm, R_Horseshoe, N_R_Horseshoe);
1056
1057 ✗ case BASKETHANDLE:
1058 ✗ return xsect->rFull * lookup(yNorm, R_Baskethandle, N_R_Baskethandle);
1059
1060 3380 case HORIZ_ELLIPSE:
1061 3380 return xsect->rFull * lookup(yNorm, R_HorizEllipse, N_R_HorizEllipse);
1062
1063 ✗ case VERT_ELLIPSE:
1064 ✗ return xsect->rFull * lookup(yNorm, R_VertEllipse, N_R_VertEllipse);
1065
1066 ✗ case ARCH:
1067 ✗ return xsect->rFull * lookup(yNorm, R_Arch, N_R_Arch);
1068
1069 4678 case IRREGULAR:
1070 4678 return xsect->rFull * lookup(yNorm,
1071 4678 Transect[xsect->transect].hradTbl, N_TRANSECT_TBL);
1072
1073 48974 case CUSTOM:
1074 48974 return xsect->rFull * lookup(yNorm,
1075 48974 Shape[Curve[xsect->transect].refersTo].hradTbl, N_SHAPE_TBL);
1076
1077 13390 case STREET_XSECT:
1078 13390 return xsect->rFull * lookup(yNorm,
1079 13390 Street[xsect->transect].transect.hradTbl,
1080 13390 Street[xsect->transect].transect.nTbl);
1081
1082 4782 case RECT_TRIANG: return rect_triang_getRofY(xsect, y);
1083
1084 28804 case RECT_ROUND: return rect_round_getRofY(xsect, y);
1085
1086 1438399 case TRAPEZOIDAL: return trapez_getRofY(xsect, y);
1087
1088 2406740 case TRIANGULAR: return triang_getRofY(xsect, y);
1089
1090 28805 case PARABOLIC: return parab_getRofY(xsect, y);
1091
1092 28805 case POWERFUNC: return powerfunc_getRofY(xsect, y);
1093
1094 118384 default: return xsect_getRofA( xsect, xsect_getAofY(xsect, y) );
1095 }
1096 }
1097
1098 //=============================================================================
1099
1100 12615319 double xsect_getRofA(TXsect *xsect, double a)
1101 //
1102 // Input: xsect = ptr. to a cross section data structure
1103 // a = area (ft2)
1104 // Output: returns hydraulic radius (ft)
1105 // Purpose: computes xsection's hydraulic radius at a given area.
1106 //
1107 {
1108 double cathy;
1109
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12615319 times.
 12615319 if ( a <= 0.0 ) return 0.0;
1110
8/11
✓ Branch 0 taken 5191285 times.
✓ Branch 1 taken 132990 times.
✓ Branch 2 taken 98681 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 140739 times.
✓ Branch 6 taken 1995203 times.
✓ Branch 7 taken 4668773 times.
✓ Branch 8 taken 193854 times.
✓ Branch 9 taken 193794 times.
✗ Branch 10 not taken.
 12615319 switch ( xsect->type )
1111 {
1112 5191285 case HORIZ_ELLIPSE:
1113 case VERT_ELLIPSE:
1114 case ARCH:
1115 case IRREGULAR:
1116 case FILLED_CIRCULAR:
1117 case CUSTOM:
1118 case STREET_XSECT:
1119 5191285 return xsect_getRofY( xsect, xsect_getYofA(xsect, a) );
1120
1121 132990 case RECT_CLOSED: return rect_closed_getRofA(xsect, a);
1122
1123 98681 case RECT_OPEN: return a / (xsect->wMax +
1124 98681 (2. - xsect->sBot) * a / xsect->wMax);
1125
1126 ✗ case RECT_TRIANG: return rect_triang_getRofA(xsect, a);
1127
1128 ✗ case RECT_ROUND: return rect_round_getRofA(xsect, a);
1129
1130 140739 case MOD_BASKET: return mod_basket_getRofA(xsect, a);
1131
1132 1995203 case TRAPEZOIDAL: return trapez_getRofA(xsect, a);
1133
1134 4668773 case TRIANGULAR: return triang_getRofA(xsect, a);
1135
1136 193854 case PARABOLIC: return parab_getRofA(xsect, a);
1137
1138 193794 case POWERFUNC: return powerfunc_getRofA(xsect, a);
1139
1140 ✗ default:
1141 ✗ cathy = xsect_getSofA(xsect, a);
1142 ✗ if ( cathy < TINY || a < TINY ) return 0.0;
1143 ✗ return pow(cathy/a, 3./2.);
1144 }
1145 }
1146
1147 //=============================================================================
1148
1149 15376361 double xsect_getAofS(TXsect* xsect, double s)
1150 //
1151 // Input: xsect = ptr. to a cross section data structure
1152 // s = section factor (ft^(8/3))
1153 // Output: returns area (ft2)
1154 // Purpose: computes xsection's area at a given section factor.
1155 //
1156 {
1157 15376361 double psi = s / xsect->sFull;
1158
2/2
✓ Branch 0 taken 11315709 times.
✓ Branch 1 taken 4060652 times.
 15376361 if ( s <= 0.0 ) return 0.0;
1159
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4060652 times.
 4060652 if ( s > xsect->sMax ) s = xsect->sMax;
1160
3/10
✗ Branch 0 not taken.
✓ Branch 1 taken 2276409 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2997 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✓ Branch 9 taken 1781246 times.
 4060652 switch ( xsect->type )
1161 {
1162 ✗ case DUMMY: return 0.0;
1163
1164 2276409 case FORCE_MAIN:
1165 2276409 case CIRCULAR: return circ_getAofS(xsect, s);
1166
1167 ✗ case EGGSHAPED:
1168 ✗ return xsect->aFull * invLookup(psi, S_Egg, N_S_Egg);
1169
1170 2997 case HORSESHOE:
1171 2997 return xsect->aFull * invLookup(psi, S_Horseshoe, N_S_Horseshoe);
1172
1173 ✗ case GOTHIC:
1174 ✗ return xsect->aFull * invLookup(psi, S_Gothic, N_S_Gothic);
1175
1176 ✗ case CATENARY:
1177 ✗ return xsect->aFull * invLookup(psi, S_Catenary, N_S_Catenary);
1178
1179 ✗ case SEMIELLIPTICAL:
1180 ✗ return xsect->aFull * invLookup(psi, S_SemiEllip, N_S_SemiEllip);
1181
1182 ✗ case BASKETHANDLE:
1183 ✗ return xsect->aFull * invLookup(psi, S_BasketHandle, N_S_BasketHandle);
1184
1185 ✗ case SEMICIRCULAR:
1186 ✗ return xsect->aFull * invLookup(psi, S_SemiCirc, N_S_SemiCirc);
1187
1188 1781246 default: return generic_getAofS(xsect, s);
1189 }
1190 }
1191
1192 //=============================================================================
1193
1194 8876256 double xsect_getdSdA(TXsect* xsect, double a)
1195 //
1196 // Input: xsect = ptr. to a cross section data structure
1197 // a = area (ft2)
1198 // Output: returns derivative of section factor w.r.t. area (ft^2/3)
1199 // Purpose: computes xsection's derivative of its section factor with
1200 // respect to area at a given area.
1201 //
1202 {
1203
9/16
✓ Branch 0 taken 75865 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✓ Branch 8 taken 43707 times.
✓ Branch 9 taken 69877 times.
✓ Branch 10 taken 11049 times.
✓ Branch 11 taken 82505 times.
✓ Branch 12 taken 67867 times.
✓ Branch 13 taken 1995203 times.
✓ Branch 14 taken 4668773 times.
✓ Branch 15 taken 1861410 times.
 8876256 switch ( xsect->type )
1204 {
1205 75865 case FORCE_MAIN:
1206 case CIRCULAR:
1207 75865 return circ_getdSdA(xsect, a);
1208
1209 ✗ case EGGSHAPED:
1210 ✗ return tabular_getdSdA(xsect, a, S_Egg, N_S_Egg);
1211
1212 ✗ case HORSESHOE:
1213 ✗ return tabular_getdSdA(xsect, a, S_Horseshoe, N_S_Horseshoe);
1214
1215 ✗ case GOTHIC:
1216 ✗ return tabular_getdSdA(xsect, a, S_Gothic, N_S_Gothic);
1217
1218 ✗ case CATENARY:
1219 ✗ return tabular_getdSdA(xsect, a, S_Catenary, N_S_Catenary);
1220
1221 ✗ case SEMIELLIPTICAL:
1222 ✗ return tabular_getdSdA(xsect, a, S_SemiEllip, N_S_SemiEllip);
1223
1224 ✗ case BASKETHANDLE:
1225 ✗ return tabular_getdSdA(xsect, a, S_BasketHandle, N_S_BasketHandle);
1226
1227 ✗ case SEMICIRCULAR:
1228 ✗ return tabular_getdSdA(xsect, a, S_SemiCirc, N_S_SemiCirc);
1229
1230 43707 case RECT_CLOSED:
1231 43707 return rect_closed_getdSdA(xsect, a);
1232
1233 69877 case RECT_OPEN:
1234 69877 return rect_open_getdSdA(xsect, a);
1235
1236 11049 case RECT_TRIANG:
1237 11049 return rect_triang_getdSdA(xsect, a);
1238
1239 82505 case RECT_ROUND:
1240 82505 return rect_round_getdSdA(xsect, a);
1241
1242 67867 case MOD_BASKET:
1243 67867 return mod_basket_getdSdA(xsect, a);
1244
1245 1995203 case TRAPEZOIDAL:
1246 1995203 return trapez_getdSdA(xsect, a);
1247
1248 4668773 case TRIANGULAR:
1249 4668773 return triang_getdSdA(xsect, a);
1250
1251 1861410 default: return generic_getdSdA(xsect, a);
1252 }
1253 }
1254
1255 //=============================================================================
1256
1257 3763349 double xsect_getYcrit(TXsect* xsect, double q)
1258 //
1259 // Input: xsect = ptr. to a cross section data structure
1260 // q = flow rate (cfs)
1261 // Output: returns critical depth (ft)
1262 // Purpose: computes critical depth at a specific flow rate.
1263 //
1264 {
1265 3763349 double q2g = SQR(q) / GRAVITY;
1266 double y, r;
1267
1268
2/2
✓ Branch 0 taken 374578 times.
✓ Branch 1 taken 3388771 times.
 3763349 if ( q2g == 0.0 ) return 0.0;
1269
5/6
✗ Branch 0 not taken.
✓ Branch 1 taken 20536 times.
✓ Branch 2 taken 919743 times.
✓ Branch 3 taken 14402 times.
✓ Branch 4 taken 14402 times.
✓ Branch 5 taken 2419688 times.
 3388771 switch ( xsect->type )
1270 {
1271 ✗ case DUMMY:
1272 ✗ return 0.0;
1273
1274 20536 case RECT_OPEN:
1275 case RECT_CLOSED:
1276 // --- analytical expression for yCritical is
1277 // y = (q2g / w^2)^(1/3) where w = width
1278 20536 y = pow(q2g / SQR(xsect->wMax), 1./3.);
1279 20536 break;
1280
1281 919743 case TRIANGULAR:
1282 // --- analytical expression for yCritical is
1283 // y = (2 * q2g / s^2)^(1/5) where s = side slope
1284 919743 y = pow(2.0 * q2g / SQR(xsect->sBot), 1./5.);
1285 919743 break;
1286
1287 14402 case PARABOLIC:
1288 // --- analytical expression for yCritical is
1289 // y = (27/32 * q2g * c)^(1/4) where y = c*x^2
1290 // is eqn. for parabola and 1/sqrt(c) = rBot
1291 14402 y = pow(27./32. * q2g / SQR(xsect->rBot), 1./4.);
1292 14402 break;
1293
1294 14402 case POWERFUNC:
1295 14402 y = 1. / (2.0 * xsect->sBot + 3.0);
1296 14402 y = pow( q2g * (xsect->sBot + 1.0) / SQR(xsect->rBot), y);
1297 14402 break;
1298
1299 2419688 default:
1300 // --- first estimate yCritical for an equivalent circular conduit
1301 // using 1.01 * (q2g / yFull)^(1/4)
1302 2419688 y = 1.01 * pow(q2g / xsect->yFull, 1./4.);
1303
2/2
✓ Branch 0 taken 206713 times.
✓ Branch 1 taken 2212975 times.
 2419688 if (y >= xsect->yFull) y = 0.97 * xsect->yFull;
1304
1305 // --- then find ratio of conduit area to equiv. circular area
1306 2419688 r = xsect->aFull / (PI / 4.0 * SQR(xsect->yFull));
1307
1308 // --- use interval enumeration method to find yCritical if
1309 // area ratio not too far from 1.0
1310
3/4
✓ Branch 0 taken 2419688 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2063620 times.
✓ Branch 3 taken 356068 times.
 2419688 if ( r >= 0.5 && r <= 2.0 )
1311 2063620 y = getYcritEnum(xsect, q, y);
1312
1313 // --- otherwise use Ridder's root finding method
1314 356068 else y = getYcritRidder(xsect, q, y);
1315 }
1316
1317 // --- do not allow yCritical to be > yFull
1318
1/2
✓ Branch 0 taken 3388771 times.
✗ Branch 1 not taken.
 3388771 return MIN(y, xsect->yFull);
1319 }
1320
1321 //=============================================================================
1322
1323 13 void getTransectParams(TXsect *xsect, TTransect *transect)
1324 //
1325 // Input: xsect = ptr. to a cross section data structure
1326 // transect = ptr. to a transect data structure
1327 // Output: none
1328 // Purpose: gets a cross section's properties from its transect's properties.
1329 //
1330 {
1331 int i, iMax;
1332 double wMax;
1333 13 double* wTbl = transect->widthTbl;
1334
1335 13 xsect->yFull = transect->yFull;
1336 13 xsect->wMax = transect->wMax;
1337 13 xsect->aFull = transect->aFull;
1338 13 xsect->rFull = transect->rFull;
1339 13 xsect->sFull = xsect->aFull * pow(xsect->rFull, 2. / 3.);
1340 13 xsect->sMax = transect->sMax;
1341 13 xsect->aBot = transect->aMax;
1342
1343 // Search transect's width table up to point where width decreases
1344 13 iMax = 0;
1345 13 wMax = wTbl[0];
1346
2/2
✓ Branch 0 taken 650 times.
✓ Branch 1 taken 13 times.
 663 for (i = 1; i < transect->nTbl; i++)
1347 {
1348
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 650 times.
 650 if (wTbl[i] < wMax) break;
1349 650 wMax = wTbl[i];
1350 650 iMax = i;
1351 }
1352
1353 // Determine height at lowest widest point
1354 13 xsect->ywMax = xsect->yFull * (double)iMax / ((double)(transect->nTbl) - 1);
1355 13 }
1356
1357 //=============================================================================
1358
1359 1781246 double generic_getAofS(TXsect* xsect, double s)
1360 //
1361 // Input: xsect = ptr. to a cross section data structure
1362 // s = section factor (ft^8/3)
1363 // Output: returns area (ft2)
1364 // Purpose: finds area given section factor by
1365 // solving S = A*(A/P(A))^(2/3) using Newton-Raphson iterations.
1366 //
1367 {
1368 double a, a1, a2, tol;
1369 TXsectStar xsectStar;
1370
1371
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1781246 times.
 1781246 if (s <= 0.0) return 0.0;
1372
1373 // --- if S is between sMax and sFull then
1374 // bracket A between aFull and aMax
1375
3/4
✓ Branch 0 taken 1781246 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 90031 times.
✓ Branch 3 taken 1691215 times.
 1781246 if ( (s <= xsect->sMax && s >= xsect->sFull)
1376
2/2
✓ Branch 0 taken 27173 times.
✓ Branch 1 taken 62858 times.
 90031 && xsect->sMax != xsect->sFull )
1377 {
1378 27173 a1 = xsect->aFull; // do this because sFull < sMax
1379 27173 a2 = xsect_getAmax(xsect);
1380 }
1381
1382 // --- otherwise bracket A between 0 and aMax
1383 else
1384 {
1385 1754073 a1 = 0.0;
1386 1754073 a2 = xsect_getAmax(xsect);
1387 }
1388
1389 // --- place S & xsect in xsectStar for access by evalSofA function
1390 1781246 xsectStar.xsect = xsect;
1391 1781246 xsectStar.s = s;
1392
1393 // --- compute starting guess for A
1394 1781246 a = 0.5 * (a1 + a2);
1395
1396 // use the Newton-Raphson root finder function to find A
1397 1781246 tol = 0.0001 * xsect->aFull;
1398 1781246 findroot_Newton(a1, a2, &a, tol, evalSofA, &xsectStar);
1399 1781246 return a;
1400 }
1401
1402 //=============================================================================
1403
1404 8800391 void evalSofA(double a, double* f, double* df, void* p)
1405 //
1406 // Input: a = area
1407 // Output: f = root finding function
1408 // df = derivative of root finding function
1409 // Purpose: function used in conjunction with getAofS() that evaluates
1410 // f = S(a) - s and df = dS(a)/dA.
1411 //
1412 {
1413 TXsectStar* xsectStar;
1414 double s;
1415
1416 8800391 xsectStar = (TXsectStar *)p;
1417 8800391 s = xsect_getSofA(xsectStar->xsect, a);
1418 8800391 *f = s - xsectStar->s;
1419 8800391 *df = xsect_getdSdA(xsectStar->xsect, a);
1420 8800391 }
1421
1422 //=============================================================================
1423
1424 60432 double tabular_getdSdA(TXsect* xsect, double a, double *table, int nItems)
1425 //
1426 // Input: xsect = ptr. to cross section data structure
1427 // a = area (ft2)
1428 // table = ptr. to table of section factor v. normalized area
1429 // nItems = number of equally spaced items in table
1430 // Output: returns derivative of section factor w.r.t. area (ft^2/3)
1431 // Purpose: computes derivative of section factor w.r.t area
1432 // using geometry tables.
1433 //
1434 {
1435 int i;
1436 60432 double alpha = a / xsect->aFull;
1437 60432 double delta = 1.0 / ((double)nItems-1);
1438 double dSdA;
1439
1440 // --- find which segment of table contains alpha
1441 60432 i = (int)(alpha / delta);
1442
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 60431 times.
 60432 if ( i >= nItems - 1 ) i = nItems - 2;
1443
1444 // --- compute slope from this interval of table
1445 60432 dSdA = (table[i+1] - table[i]) / delta;
1446
1447 // --- convert slope to un-normalized value
1448 60432 return dSdA * xsect->sFull / xsect->aFull;
1449 }
1450
1451 //=============================================================================
1452
1453 1958317 double generic_getdSdA(TXsect* xsect, double a)
1454 //
1455 // Input: xsect = ptr. to cross section data structure
1456 // a = area (ft2)
1457 // Output: returns derivative of section factor w.r.t. area (ft^2/3)
1458 // Purpose: computes derivative of section factor w.r.t area
1459 // using central difference approximation.
1460 //
1461 {
1462 double a1, a2;
1463 1958317 double alpha = a / xsect->aFull;
1464 1958317 double alpha1 = alpha - 0.001;
1465 1958317 double alpha2 = alpha + 0.001;
1466
2/2
✓ Branch 0 taken 5336 times.
✓ Branch 1 taken 1952981 times.
 1958317 if ( alpha1 < 0.0 ) alpha1 = 0.0;
1467 1958317 a1 = alpha1 * xsect->aFull;
1468 1958317 a2 = alpha2 * xsect->aFull;
1469 1958317 return (xsect_getSofA(xsect, a2) - xsect_getSofA(xsect, a1)) / (a2 - a1);
1470 }
1471
1472 //=============================================================================
1473
1474 549957237 double lookup(double x, double *table, int nItems)
1475 //
1476 // Input: x = value of independent variable in a geometry table
1477 // table = ptr. to geometry table
1478 // nItems = number of equally spaced items in table
1479 // Output: returns value of dependent table variable
1480 // Purpose: looks up a value in a geometry table (i.e., finds y given x).
1481 //
1482 {
1483 double delta, x0, x1, y, y2;
1484 int i;
1485
1486 // --- find which segment of table contains x
1487 549957237 delta = 1.0 / ((double)nItems-1);
1488 549957237 i = (int)(x / delta);
1489
2/2
✓ Branch 0 taken 454606 times.
✓ Branch 1 taken 549502631 times.
 549957237 if ( i >= nItems - 1 ) return table[nItems-1];
1490
1491 // --- compute x at start and end of segment
1492 549502631 x0 = i * delta;
1493 549502631 x1 = ((double)i+1) * delta;
1494
1495 // --- linearly interpolate a y-value
1496 549502631 y = table[i] + (x - x0) * (table[i+1] - table[i]) / delta;
1497
1498 // --- use quadratic interpolation for low x value
1499
2/2
✓ Branch 0 taken 69241543 times.
✓ Branch 1 taken 480261088 times.
 549502631 if ( i < 2 )
1500 {
1501 69241543 y2 = y + (x - x0) * (x - x1) / (delta*delta) *
1502 69241543 (table[i]/2.0 - table[i+1] + table[i+2]/2.0) ;
1503
2/2
✓ Branch 0 taken 68700118 times.
✓ Branch 1 taken 541425 times.
 69241543 if ( y2 > 0.0 ) y = y2;
1504 }
1505
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 549502631 times.
 549502631 if ( y < 0.0 ) y = 0.0;
1506 549502631 return y;
1507 }
1508
1509 //=============================================================================
1510
1511 832898 double invLookup(double y, double *table, int nItems)
1512 //
1513 // Input: y = value of dependent variable in a geometry table
1514 // table = ptr. to geometry table
1515 // nItems = number of equally spaced items in table
1516 // Output: returns value of independent table variable
1517 // Purpose: performs inverse lookup in a geometry table (i.e., finds
1518 // x given y).
1519 //
1520 // Notes: This function assumes that the geometry table has either strictly
1521 // increasing entries or that the maximum entry is always third
1522 // from the last (which is true for all section factor tables). In
1523 // the latter case, the location of a large y can be ambiguous
1524 // -- it can be both below and above the location of the maximum.
1525 // In such cases this routine searches only the interval above
1526 // the maximum (i.e., the last 2 segments of the table).
1527 //
1528 // nItems-1 is the highest subscript for the table's data.
1529 //
1530 // The x value's in a geometry table lie between 0 and 1.
1531 //
1532 {
1533 double dx; // x-increment of table
1534 double x, x0, dy; // interpolation variables
1535 int n; // # items in increasing portion of table
1536 int i; // lower table index that brackets y
1537
1538 // --- compute table's uniform x-increment
1539 832898 dx = 1.0 / (double)((double)nItems-1);
1540
1541 // --- truncate item count if last 2 table entries are decreasing
1542 832898 n = nItems;
1543
1/2
✓ Branch 0 taken 832898 times.
✗ Branch 1 not taken.
 832898 if ( table[n-3] > table[n-1] ) n = n - 2;
1544
1545 // --- check if y falls in decreasing portion of table
1546
2/4
✓ Branch 0 taken 832898 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 832898 times.
 832898 if ( n < nItems && y > table[nItems-1])
1547 {
1548 ✗ if ( y >= table[nItems-3] ) return ((double)n-1) * dx;
1549 ✗ if ( y <= table[nItems-2] ) i = nItems - 2;
1550 ✗ else i = nItems - 3;
1551 }
1552
1553 // --- otherwise locate the interval where y falls in the table
1554 832898 else i = locate(y, table, n-1);
1555
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 832898 times.
 832898 if ( i >= n - 1 ) return ((double)n-1) * dx;
1556
1557 // --- compute x at start and end of segment
1558 832898 x0 = i * dx;
1559
1560 // --- linearly interpolate an x value
1561 832898 dy = table[i+1] - table[i];
1562
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 832898 times.
 832898 if ( dy == 0.0 ) x = x0;
1563 832898 else x = x0 + (y - table[i]) * dx / dy;
1564
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 832898 times.
 832898 if ( x < 0.0 ) x = 0.0;
1565
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 832898 times.
 832898 if ( x > 1.0 ) x = 1.0;
1566 832898 return x;
1567 }
1568
1569 //=============================================================================
1570
1571 832898 int locate(double y, double *table, int jLast)
1572 //
1573 // Input: y = value being located in table
1574 // table = ptr. to table with monotonically increasing entries
1575 // jLast = highest table entry index to search over
1576 // Output: returns index j of table such that table[j] <= y <= table[j+1]
1577 // Purpose: uses bisection method to locate the highest table index whose
1578 // table entry does not exceed a given value.
1579 //
1580 // Notes: This function is only used in conjunction with invLookup().
1581 //
1582 {
1583 int j;
1584 832898 int j1 = 0;
1585 832898 int j2 = jLast;
1586
1587 // Check if value <= first table entry
1588
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 832898 times.
 832898 if ( y <= table[0] ) return 0;
1589
1590 // Check if value >= the last entry
1591
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 832898 times.
 832898 if ( y >= table[jLast] ) return jLast;
1592
1593 // While a portion of the table still remains
1594
2/2
✓ Branch 0 taken 4733685 times.
✓ Branch 1 taken 832898 times.
 5566583 while ( j2 - j1 > 1)
1595 {
1596 // Find midpoint of remaining portion of table
1597 4733685 j = (j1 + j2) >> 1;
1598
1599 // Value is greater or equal to midpoint: search from midpoint to j2
1600
2/2
✓ Branch 0 taken 1791301 times.
✓ Branch 1 taken 2942384 times.
 4733685 if ( y >= table[j] ) j1 = j;
1601
1602 // Value is less than midpoint: search from j1 to midpoint
1603 2942384 else j2 = j;
1604 }
1605
1606 // Return the lower index of the remaining interval,
1607 832898 return j1;
1608 }
1609
1610 //=============================================================================
1611
1612 7345275 double getQcritical(double yc, void* p)
1613 //
1614 // Input: yc = critical depth (ft)
1615 // p = pointer to a TXsectStar object
1616 // Output: returns flow difference value (cfs)
1617 // Purpose: finds difference between critical flow at depth yc and
1618 // some target value.
1619 //
1620 {
1621 double a, w, qc;
1622 TXsectStar* xsectStar;
1623
1624 7345275 xsectStar = (TXsectStar *)p;
1625 7345275 a = xsect_getAofY(xsectStar->xsect, yc);
1626 7345275 w = xsect_getWofY(xsectStar->xsect, yc);
1627 7345275 qc = -xsectStar->qc;
1628
2/2
✓ Branch 0 taken 6799086 times.
✓ Branch 1 taken 546189 times.
 7345275 if ( w > 0.0 ) qc = a * sqrt(GRAVITY * a / w) - xsectStar->qc;
1629 7345275 return qc;
1630 }
1631
1632 //=============================================================================
1633
1634 2063620 double getYcritEnum(TXsect* xsect, double q, double y0)
1635 //
1636 // Input: xsect = ptr. to cross section data structure
1637 // q = critical flow rate (cfs)
1638 // y0 = estimate of critical depth (ft)
1639 // Output: returns true critical depth (ft)
1640 // Purpose: solves a * sqrt(a(y)*g / w(y)) - q for y using interval
1641 // enumeration with starting guess of y0.
1642 //
1643 {
1644 double q0, dy, qc, yc;
1645 int i1, i;
1646 TXsectStar xsectStar;
1647
1648 // --- divide cross section depth into 25 increments and
1649 // locate increment corresponding to initial guess y0
1650 2063620 dy = xsect->yFull / 25.;
1651 2063620 i1 = (int)(y0 / dy);
1652
1653 // --- evaluate critical flow at this increment
1654 2063620 xsectStar.xsect = xsect;
1655 2063620 xsectStar.qc = 0.0;
1656 2063620 q0 = getQcritical(i1*dy, &xsectStar);
1657
1658 // --- initial flow lies below target flow
1659
2/2
✓ Branch 0 taken 1534906 times.
✓ Branch 1 taken 528714 times.
 2063620 if ( q0 < q )
1660 {
1661 // --- search each successive higher depth increment
1662 1534906 yc = xsect->yFull;
1663
1/2
✓ Branch 0 taken 1542397 times.
✗ Branch 1 not taken.
 1542397 for ( i = i1+1; i <= 25; i++)
1664 {
1665 // --- if critical flow at current depth is above target
1666 // then use linear interpolation to compute critical depth
1667 1542397 qc = getQcritical(i*dy, &xsectStar);
1668
2/2
✓ Branch 0 taken 1534906 times.
✓ Branch 1 taken 7491 times.
 1542397 if ( qc >= q )
1669 {
1670 1534906 yc = ( (q-q0) / (qc - q0) + ((double)i-1) ) * dy;
1671 1534906 break;
1672 }
1673 7491 q0 = qc;
1674 }
1675 }
1676
1677 // --- initial flow lies above target flow
1678 else
1679 {
1680 // --- search each successively lower depth increment
1681 528714 yc = 0.0;
1682
1/2
✓ Branch 0 taken 831017 times.
✗ Branch 1 not taken.
 831017 for ( i = i1-1; i >= 0; i--)
1683 {
1684 // --- if critical flow at current depth is below target
1685 // then use linear interpolation to compute critical depth
1686 831017 qc = getQcritical(i*dy, &xsectStar);
1687
2/2
✓ Branch 0 taken 528714 times.
✓ Branch 1 taken 302303 times.
 831017 if ( qc < q )
1688 {
1689 528714 yc = ( (q-qc) / (q0-qc) + (double)i ) * dy;
1690 528714 break;
1691 }
1692 302303 q0 = qc;
1693 }
1694 }
1695 2063620 return yc;
1696 }
1697
1698 //=============================================================================
1699
1700 356068 double getYcritRidder(TXsect* xsect, double q, double y0)
1701 //
1702 // Input: xsect = ptr. to cross section data structure
1703 // q = critical flow rate (cfs)
1704 // y0 = estimate of critical depth (ft)
1705 // Output: returns true critical depth (ft)
1706 // Purpose: solves a * sqrt(a(y)*g / w(y)) - q for y using Ridder's
1707 // root finding method with starting guess of y0.
1708 //
1709 {
1710 356068 double y1 = 0.0;
1711 356068 double y2 = 0.99 * xsect->yFull;
1712 double yc;
1713 double q0, q1, q2;
1714 TXsectStar xsectStar;
1715
1716 // --- store reference to cross section in global pointer
1717 356068 xsectStar.xsect = xsect;
1718 356068 xsectStar.qc = 0.0;
1719
1720 // --- check if critical flow at (nearly) full depth < target flow
1721 356068 q2 = getQcritical(y2, &xsectStar);
1722
2/2
✓ Branch 0 taken 61119 times.
✓ Branch 1 taken 294949 times.
 356068 if (q2 < q ) return xsect->yFull;
1723
1724 // --- evaluate critical flow at initial depth guess y0
1725 // and at 1/2 of full depth
1726 294949 q0 = getQcritical(y0, &xsectStar);
1727 294949 q1 = getQcritical(0.5*xsect->yFull, &xsectStar);
1728
1729 // --- adjust search interval on depth so it contains flow q
1730
2/2
✓ Branch 0 taken 289203 times.
✓ Branch 1 taken 5746 times.
 294949 if ( q0 > q )
1731 {
1732 289203 y2 = y0;
1733
2/2
✓ Branch 0 taken 109466 times.
✓ Branch 1 taken 179737 times.
 289203 if ( q1 < q ) y1 = 0.5*xsect->yFull;
1734 }
1735 else
1736 {
1737 5746 y1 = y0;
1738
2/2
✓ Branch 0 taken 2860 times.
✓ Branch 1 taken 2886 times.
 5746 if ( q1 > q ) y2 = 0.5*xsect->yFull;
1739 }
1740
1741 // --- save value of target critical flow in global variable
1742 294949 xsectStar.qc = q;
1743
1744 // --- call Ridder root finding procedure with error tolerance
1745 // of 0.001 ft. to find critical depth yc
1746 294949 yc = findroot_Ridder(y1, y2, 0.001, getQcritical, &xsectStar);
1747 294949 return yc;
1748 }
1749
1750
1751 //=============================================================================
1752 // RECT_CLOSED fuctions
1753 //=============================================================================
1754
1755 43707 double rect_closed_getSofA(TXsect* xsect, double a)
1756 {
1757 // --- if a > area corresponding to Smax then
1758 // interpolate between sMax and Sfull
1759 43707 double alfMax = RECT_ALFMAX;
1760
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 43707 times.
 43707 if ( a / xsect->aFull > alfMax )
1761 {
1762 ✗ return xsect->sMax + (xsect->sFull - xsect->sMax) *
1763 ✗ (a/xsect->aFull - alfMax) / (1.0 - alfMax);
1764 }
1765
1766 // --- otherwise use regular formula
1767 43707 return a * pow(xsect_getRofA(xsect, a), 2./3.);
1768 }
1769
1770 43707 double rect_closed_getdSdA(TXsect* xsect, double a)
1771 {
1772 double alpha, alfMax, r;
1773
1774 // --- if above level corresponding to sMax, then
1775 // use slope between sFull & sMax
1776 43707 alfMax = RECT_ALFMAX;
1777 43707 alpha = a / xsect->aFull;
1778
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 43707 times.
 43707 if ( alpha > alfMax )
1779 {
1780 ✗ return (xsect->sFull - xsect->sMax) /
1781 ✗ ((1.0 - alfMax) * xsect->aFull);
1782 }
1783
1784 // --- for small a/aFull use generic central difference formula
1785
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 43707 times.
 43707 if ( alpha <= 1.0e-30 ) return generic_getdSdA(xsect, a);
1786
1787 // --- otherwise evaluate dSdA = [5/3 - (2/3)(dP/dA)R]R^(2/3)
1788 // (where P = wetted perimeter & dPdA = 2/width)
1789 43707 r = xsect_getRofA(xsect, a);
1790 43707 return (5./3. - (2./3.) * (2.0/xsect->wMax) * r) * pow(r, 2./3.);
1791 }
1792
1793 133186 double rect_closed_getRofA(TXsect* xsect, double a)
1794 {
1795 double p;
1796
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 133186 times.
 133186 if ( a <= 0.0 ) return 0.0;
1797 133186 p = xsect->wMax + 2.*a/xsect->wMax; // Wetted Perim = width + 2*area/width
1798
2/2
✓ Branch 0 taken 70 times.
✓ Branch 1 taken 133116 times.
 133186 if ( a/xsect->aFull > RECT_ALFMAX )
1799 {
1800 70 p += (a/xsect->aFull - RECT_ALFMAX) / (1.0 - RECT_ALFMAX) * xsect->wMax;
1801 }
1802 133186 return a / p;
1803 }
1804
1805
1806 //=============================================================================
1807 // RECT_OPEN fuctions
1808 //=============================================================================
1809
1810 69877 double rect_open_getSofA(TXsect* xsect, double a)
1811 {
1812 69877 double y = a / xsect->wMax;
1813 69877 double r = a / ((2.0-xsect->sBot)*y + xsect->wMax);
1814 69877 return a * pow(r, 2./3.);
1815 }
1816
1817
1818 69877 double rect_open_getdSdA(TXsect* xsect, double a)
1819 {
1820 double r, dPdA;
1821
1822 // --- for small a/aFull use generic central difference formula
1823
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 69877 times.
 69877 if ( a / xsect->aFull <= 1.0e-30 ) return generic_getdSdA(xsect, a);
1824
1825 // --- otherwise evaluate dSdA = [5/3 - (2/3)(dP/dA)R]R^(2/3)
1826 // (where P = wetted perimeter)
1827 69877 r = xsect_getRofA(xsect, a);
1828 69877 dPdA = (2.0 - xsect->sBot) / xsect->wMax; // since P = geom2 + 2a/geom2
1829 69877 return (5./3. - (2./3.) * dPdA * r) * pow(r, 2./3.);
1830 }
1831
1832
1833 //=============================================================================
1834 // RECT_TRIANG fuctions
1835 //=============================================================================
1836
1837 23657 double rect_triang_getYofA(TXsect* xsect, double a)
1838 {
1839 // below upper section
1840
2/2
✓ Branch 0 taken 14373 times.
✓ Branch 1 taken 9284 times.
 23657 if ( a <= xsect->aBot ) return sqrt(a / xsect->sBot);
1841
1842 // above bottom section
1843 9284 else return xsect->yBot + (a - xsect->aBot) / xsect->wMax;
1844 }
1845
1846 22112 double rect_triang_getRofA(TXsect* xsect, double a)
1847 {
1848 double y;
1849 double p, alf;
1850
1851
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 22112 times.
 22112 if ( a <= 0.0 ) return 0.0;
1852 22112 y = rect_triang_getYofA(xsect, a);
1853
1854 // below upper section
1855
2/2
✓ Branch 0 taken 12828 times.
✓ Branch 1 taken 9284 times.
 22112 if ( y <= xsect->yBot ) return a / (2. * y * xsect->rBot);
1856
1857 // wetted perimeter without contribution of top surface
1858 9284 p = 2. * xsect->yBot * xsect->rBot + 2. * (y - xsect->yBot);
1859
1860 // top-surface contribution
1861 9284 alf = (a / xsect->aFull) - RECT_TRIANG_ALFMAX;
1862
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9284 times.
 9284 if ( alf > 0.0 ) p += alf / (1.0 - RECT_TRIANG_ALFMAX) * xsect->wMax;
1863 9284 return a / p;
1864 }
1865
1866 11049 double rect_triang_getSofA(TXsect* xsect, double a)
1867 {
1868 // --- if a > area corresponding to sMax, then
1869 // interpolate between sMax and Sfull
1870 11049 double alfMax = RECT_TRIANG_ALFMAX;
1871
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11049 times.
 11049 if ( a / xsect->aFull > alfMax )
1872 ✗ return xsect->sMax + (xsect->sFull - xsect->sMax) *
1873 ✗ (a/xsect->aFull - alfMax) / (1.0 - alfMax);
1874
1875 // --- otherwise use regular formula
1876 11049 else return a * pow(rect_triang_getRofA(xsect, a), 2./3.);
1877 }
1878
1879 11049 double rect_triang_getdSdA(TXsect* xsect, double a)
1880 {
1881 double alpha, alfMax, dPdA, r;
1882
1883 // --- if a > area corresponding to sMax, then
1884 // use slope between sFull & sMax
1885 11049 alfMax = RECT_TRIANG_ALFMAX;
1886 11049 alpha = a / xsect->aFull;
1887
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11049 times.
 11049 if ( alpha > alfMax )
1888 ✗ return (xsect->sFull - xsect->sMax) / ((1.0 - alfMax) * xsect->aFull);
1889
1890 // --- use generic central difference method for very small a
1891
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11049 times.
 11049 if ( alpha <= 1.0e-30 ) return generic_getdSdA(xsect, a);
1892
1893 // --- find deriv. of wetted perimeter
1894
2/2
✓ Branch 0 taken 4635 times.
✓ Branch 1 taken 6414 times.
 11049 if ( a > xsect->aBot ) dPdA = 2.0 / xsect->wMax; // for upper rectangle
1895 6414 else dPdA = xsect->rBot / sqrt(a * xsect->sBot); // for triang. bottom
1896
1897 // --- get hyd. radius & evaluate section factor derivative formula
1898 11049 r = rect_triang_getRofA(xsect, a);
1899 11049 return (5./3. - (2./3.) * dPdA * r) * pow(r, 2./3.);
1900 }
1901
1902 10991 double rect_triang_getAofY(TXsect* xsect, double y)
1903 {
1904
2/2
✓ Branch 0 taken 5631 times.
✓ Branch 1 taken 5360 times.
 10991 if ( y <= xsect->yBot ) return y * y * xsect->sBot; // below upper section
1905 5360 else return xsect->aBot + (y - xsect->yBot) * xsect->wMax; // above bottom section
1906 }
1907
1908 4782 double rect_triang_getRofY(TXsect* xsect, double y)
1909 {
1910 double p, a, alf;
1911
1912 // y is below upper rectangular section
1913
2/2
✓ Branch 0 taken 3235 times.
✓ Branch 1 taken 1547 times.
 4782 if ( y <= xsect->yBot ) return y * xsect->sBot / (2. * xsect->rBot);
1914
1915 // area
1916 1547 a = xsect->aBot + (y - xsect->yBot) * xsect->wMax;
1917
1918 // wetted perimeter without contribution of top surface
1919 1547 p = 2. * xsect->yBot * xsect->rBot + 2. * (y - xsect->yBot);
1920
1921 // top-surface contribution
1922 1547 alf = (a / xsect->aFull) - RECT_TRIANG_ALFMAX;
1923
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1547 times.
 1547 if ( alf > 0.0 ) p += alf / (1.0 - RECT_TRIANG_ALFMAX) * xsect->wMax;
1924 1547 return a / p;
1925 }
1926
1927 9276 double rect_triang_getWofY(TXsect* xsect, double y)
1928 {
1929
2/2
✓ Branch 0 taken 4637 times.
✓ Branch 1 taken 4639 times.
 9276 if ( y <= xsect->yBot ) return 2.0 * xsect->sBot * y; // below upper section
1930 4639 else return xsect->wMax; // above bottom section
1931 }
1932
1933
1934 //=============================================================================
1935 // RECT_ROUND fuctions
1936 //=============================================================================
1937
1938 14402 double rect_round_getYofA(TXsect* xsect, double a)
1939 {
1940 double alpha;
1941
1942 // --- if above circular bottom:
1943
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14402 times.
 14402 if ( a > xsect->aBot )
1944 ✗ return xsect->yBot + (a - xsect->aBot) / xsect->wMax;
1945
1946 // --- otherwise use circular xsection method to find height
1947 14402 alpha = a / (PI * xsect->rBot * xsect->rBot);
1948
2/2
✓ Branch 0 taken 6075 times.
✓ Branch 1 taken 8327 times.
 14402 if ( alpha < 0.04 ) return (2.0 * xsect->rBot) * getYcircular(alpha);
1949 8327 return (2.0 * xsect->rBot) * lookup(alpha, Y_Circ, N_Y_Circ);
1950 }
1951
1952 1 double rect_round_getRofA(TXsect* xsect, double a)
1953 {
1954 double y1, theta1, p, arg;
1955
1956 // --- if above circular invert ...
1957
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( a <= 0.0 ) return 0.0;
1958
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
 1 if ( a > xsect->aBot )
1959 {
1960 // wetted perimeter without contribution of top surface
1961 1 y1 = (a - xsect->aBot) / xsect->wMax;
1962 1 theta1 = 2.0 * asin(xsect->wMax/2.0/xsect->rBot);
1963 1 p = xsect->rBot*theta1 + 2.0*y1;
1964
1965 // top-surface contribution
1966 1 arg = (a / xsect->aFull) - RECT_ROUND_ALFMAX;
1967
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( arg > 0.0 ) p += arg / (1.0 - RECT_ROUND_ALFMAX) * xsect->wMax;
1968 1 return a / p;
1969 }
1970
1971 // --- if within circular invert ...
1972 ✗ y1 = rect_round_getYofA(xsect, a);
1973 ✗ theta1 = 2.0*acos(1.0 - y1/xsect->rBot);
1974 ✗ p = xsect->rBot * theta1;
1975 ✗ return a / p;
1976 }
1977
1978 247515 double rect_round_getSofA(TXsect* xsect, double a)
1979 {
1980 double alpha, aFull, sFull;
1981
1982 // --- if a > area corresponding to sMax,
1983 // interpolate between sMax and sFull
1984 247515 double alfMax = RECT_ROUND_ALFMAX;
1985
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 247515 times.
 247515 if ( a / xsect->aFull > alfMax )
1986 {
1987 ✗ return xsect->sMax + (xsect->sFull - xsect->sMax) *
1988 ✗ (a / xsect->aFull - alfMax) / (1.0 - alfMax);
1989 }
1990
1991 // --- if above circular invert, use generic function
1992
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 247515 times.
 247515 else if ( a > xsect->aBot )
1993 {
1994 ✗ return a * pow(xsect_getRofA(xsect, a), 2./3.);
1995 }
1996
1997 // --- otherwise use circular xsection function applied
1998 // to full circular shape of bottom section
1999 else
2000 {
2001 247515 aFull = PI * xsect->rBot * xsect->rBot;
2002 247515 alpha = a / aFull;
2003 247515 sFull = xsect->sBot;
2004
2005 // --- use special function for small a/aFull
2006
2/2
✓ Branch 0 taken 57848 times.
✓ Branch 1 taken 189667 times.
 247515 if ( alpha < 0.04 ) return sFull * getScircular(alpha);
2007
2008 // --- otherwise use table
2009 189667 else return sFull * lookup(alpha, S_Circ, N_S_Circ);
2010 }
2011 }
2012
2013 82505 double rect_round_getdSdA(TXsect* xsect, double a)
2014 {
2015 double alfMax, r, dPdA;
2016
2017 // --- if a > area corresponding to sMax, then
2018 // use slope between sFull & sMax
2019 82505 alfMax = RECT_ROUND_ALFMAX;
2020
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 82505 times.
 82505 if ( a / xsect->aFull > alfMax )
2021 {
2022 ✗ return (xsect->sFull - xsect->sMax) /
2023 ✗ ((1.0 - alfMax) * xsect->aFull);
2024 }
2025
2026 // --- if above circular invert, use analytical function for dS/dA
2027
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 82505 times.
 82505 else if ( a > xsect->aBot )
2028 {
2029 ✗ r = rect_round_getRofA(xsect, a);
2030 ✗ dPdA = 2.0 / xsect->wMax; // d(wet perim)/dA for rect.
2031 ✗ return (5./3. - (2./3.) * dPdA * r) * pow(r, 2./3.);
2032 }
2033
2034 // --- otherwise use generic finite difference function
2035 82505 else return generic_getdSdA(xsect, a);
2036 }
2037
2038 95561 double rect_round_getAofY(TXsect* xsect, double y)
2039 {
2040 double theta1;
2041
2042 // --- if above circular invert...
2043
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 95561 times.
 95561 if ( y > xsect->yBot )
2044 ✗ return xsect->aBot + (y - xsect->yBot) * xsect->wMax;
2045
2046 // --- find area of circular section
2047 95561 theta1 = 2.0*acos(1.0 - y/xsect->rBot);
2048 95561 return 0.5 * xsect->rBot * xsect->rBot * (theta1 - sin(theta1));
2049 }
2050
2051 28804 double rect_round_getRofY(TXsect* xsect, double y)
2052 {
2053 double theta1;
2054
2055 // --- if above top of circular bottom, use RofA formula
2056
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 28804 times.
 28804 if ( y <= 0.0 ) return 0.0;
2057
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 28804 times.
 28804 if ( y > xsect->yBot )
2058 ✗ return rect_round_getRofA( xsect, rect_round_getAofY(xsect, y) );
2059
2060 // --- find hyd. radius of circular section
2061 28804 theta1 = 2.0*acos(1.0 - y/xsect->rBot);
2062 28804 return 0.5 * xsect->rBot * (1.0 - sin(theta1)) / theta1;
2063 }
2064
2065 88596 double rect_round_getWofY(TXsect* xsect, double y)
2066 {
2067 // --- return width if depth above circular bottom section
2068
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 88596 times.
 88596 if ( y > xsect->yBot ) return xsect->wMax;
2069
2070 // --- find width of circular section
2071 88596 return 2.0 * sqrt( y * (2.0*xsect->rBot - y) );
2072 }
2073
2074
2075 //=============================================================================
2076 // MOD_BASKETHANDLE fuctions
2077 //=============================================================================
2078
2079 // Note: the variables rBot, yBot, and aBot refer to properties of the
2080 // circular top portion of the cross-section (not the bottom)
2081
2082 57672 double mod_basket_getYofA(TXsect* xsect, double a)
2083 {
2084 double alpha, y1;
2085
2086 // --- water level below top of rectangular bottom
2087
2/2
✓ Branch 0 taken 14402 times.
✓ Branch 1 taken 43270 times.
 57672 if ( a <= xsect->aFull - xsect->aBot ) return a / xsect->wMax;
2088
2089 // --- find unfilled top area / area of full circular top
2090 43270 alpha = (xsect->aFull - a) / (PI * xsect->rBot * xsect->rBot);
2091
2092 // --- find unfilled height
2093
2/2
✓ Branch 0 taken 57 times.
✓ Branch 1 taken 43213 times.
 43270 if ( alpha < 0.04 ) y1 = getYcircular(alpha);
2094 43213 else y1 = lookup(alpha, Y_Circ, N_Y_Circ);
2095 43270 y1 = 2.0 * xsect->rBot * y1;
2096
2097 // --- return difference between full height & unfilled height
2098 43270 return xsect->yFull - y1;
2099 }
2100
2101 140739 double mod_basket_getRofA(TXsect* xsect, double a)
2102 {
2103 double y1, p, theta1;
2104
2105 // --- water level is below top of rectangular bottom;
2106 // return hyd. radius of rectangle
2107
2/2
✓ Branch 0 taken 97469 times.
✓ Branch 1 taken 43270 times.
 140739 if ( a <= xsect->aFull - xsect->aBot )
2108 97469 return a / (xsect->wMax + 2.0 * a / xsect->wMax);
2109
2110 // --- find height of empty area
2111 43270 y1 = xsect->yFull - mod_basket_getYofA(xsect, a);
2112
2113 // --- find angle of circular arc corresponding to this height
2114 43270 theta1 = 2.0 * acos(1.0 - y1 / xsect->rBot);
2115
2116 // --- find perimeter of wetted portion of circular arc
2117 // (angle of full circular opening was stored in sBot)
2118 43270 p = (xsect->sBot - theta1) * xsect->rBot;
2119
2120 // --- add on wetted perimeter of bottom rectangular area
2121 43270 y1 = xsect->yFull - xsect->yBot;
2122 43270 p = p + 2.0*y1 + xsect->wMax;
2123
2124 // --- return area / wetted perimeter
2125 43270 return a / p;
2126 }
2127
2128 67867 double mod_basket_getdSdA(TXsect* xsect, double a)
2129 {
2130 double r, dPdA;
2131
2132 // --- if water level below top of rectangular bottom but not
2133 // empty then use same code as for rectangular xsection
2134
3/4
✓ Branch 0 taken 53465 times.
✓ Branch 1 taken 14402 times.
✓ Branch 2 taken 53465 times.
✗ Branch 3 not taken.
 67867 if ( a <= xsect->aFull - xsect->aBot && a/xsect->aFull > 1.0e-30 )
2135 {
2136 53465 r = a / (xsect->wMax + 2.0 * a / xsect->wMax);
2137 53465 dPdA = 2.0 / xsect->wMax;
2138 53465 return (5./3. - (2./3.) * dPdA * r) * pow(r, 2./3.);
2139 }
2140
2141 // --- otherwise use generic function
2142 14402 else return generic_getdSdA(xsect, a);
2143 }
2144
2145 205048 double mod_basket_getAofY(TXsect* xsect, double y)
2146 {
2147 double a1, theta1, y1;
2148
2149 // --- if water level is below top of rectangular bottom
2150 // return depth * width
2151
2/2
✓ Branch 0 taken 173442 times.
✓ Branch 1 taken 31606 times.
 205048 if ( y <= xsect->yFull - xsect->yBot ) return y * xsect->wMax;
2152
2153 // --- find empty top circular area
2154 31606 y1 = xsect->yFull - y;
2155 31606 theta1 = 2.0*acos(1.0 - y1/xsect->rBot);
2156 31606 a1 = 0.5 * xsect->rBot * xsect->rBot * (theta1 - sin(theta1));
2157
2158 // --- return difference between full and empty areas
2159 31606 return xsect->aFull - a1;
2160 }
2161
2162 151099 double mod_basket_getWofY(TXsect* xsect, double y)
2163 {
2164 double y1;
2165
2166 // --- if water level below top of rectangular bottom then return width
2167
2/2
✓ Branch 0 taken 1862 times.
✓ Branch 1 taken 149237 times.
 151099 if ( y <= 0.0 ) return 0.0;
2168
2/2
✓ Branch 0 taken 117631 times.
✓ Branch 1 taken 31606 times.
 149237 if ( y <= xsect->yFull - xsect->yBot ) return xsect->wMax;
2169
2170 // --- find width of empty top circular section
2171 31606 y1 = xsect->yFull - y;
2172 31606 return 2.0 * sqrt( y1 * (2.0*xsect->rBot - y1) );
2173 }
2174
2175
2176 //=============================================================================
2177 // TRAPEZOIDAL fuctions
2178 //
2179 // Note: yBot = width of bottom
2180 // sBot = avg. of side slopes
2181 // rBot = length of sides per unit of depth
2182 //=============================================================================
2183
2184 4346474 double trapez_getYofA(TXsect* xsect, double a)
2185 {
2186
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4346474 times.
 4346474 if ( xsect->sBot == 0.0 ) return a / xsect->yBot;
2187 4346474 return ( sqrt( xsect->yBot*xsect->yBot + 4.*xsect->sBot*a )
2188 4346474 - xsect->yBot )/(2. * xsect->sBot);
2189 }
2190
2191 3990406 double trapez_getRofA(TXsect* xsect, double a)
2192 {
2193 3990406 return a / (xsect->yBot + trapez_getYofA(xsect, a) * xsect->rBot);
2194 }
2195
2196 1995203 double trapez_getdSdA(TXsect* xsect, double a)
2197 {
2198 double r, dPdA;
2199 // --- use generic central difference method for very small a
2200
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1995203 times.
 1995203 if ( a/xsect->aFull <= 1.0e-30 ) return generic_getdSdA(xsect, a);
2201
2202 // --- otherwise use analytical formula:
2203 // dSdA = [5/3 - (2/3)(dP/dA)R]R^(2/3)
2204 1995203 r = trapez_getRofA(xsect, a);
2205 1995203 dPdA = xsect->rBot /
2206 1995203 sqrt( xsect->yBot * xsect->yBot + 4. * xsect->sBot * a );
2207 1995203 return (5./3. - (2./3.) * dPdA * r) * pow(r, 2./3.);
2208 }
2209
2210 7466141 double trapez_getAofY(TXsect* xsect, double y)
2211 {
2212 7466141 return ( xsect->yBot + xsect->sBot * y ) * y;
2213 }
2214
2215 1438399 double trapez_getRofY(TXsect* xsect, double y)
2216 {
2217
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1438399 times.
 1438399 if ( y == 0.0 ) return 0.0;
2218 1438399 return trapez_getAofY(xsect, y) / (xsect->yBot + y * xsect->rBot);
2219 }
2220
2221 6016695 double trapez_getWofY(TXsect* xsect, double y)
2222 {
2223 6016695 return xsect->yBot + 2.0 * y * xsect->sBot;
2224 }
2225
2226
2227 //=============================================================================
2228 // TRIANGULAR fuctions
2229 //=============================================================================
2230
2231 10257289 double triang_getYofA(TXsect* xsect, double a)
2232 {
2233 10257289 return sqrt(a / xsect->sBot);
2234 }
2235
2236 9337546 double triang_getRofA(TXsect* xsect, double a)
2237 {
2238 9337546 return a / (2. * triang_getYofA(xsect, a) * xsect->rBot);
2239 }
2240
2241 4668773 double triang_getdSdA(TXsect* xsect, double a)
2242 {
2243 double r, dPdA;
2244 // --- use generic finite difference method for very small 'a'
2245
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4668773 times.
 4668773 if ( a/xsect->aFull <= 1.0e-30 ) return generic_getdSdA(xsect, a);
2246
2247 // --- evaluate dSdA = [5/3 - (2/3)(dP/dA)R]R^(2/3)
2248 4668773 r = triang_getRofA(xsect, a);
2249 4668773 dPdA = xsect->rBot / sqrt(a * xsect->sBot);
2250 4668773 return (5./3. - (2./3.) * dPdA * r) * pow(r, 2./3.);
2251 }
2252
2253 4962217 double triang_getAofY(TXsect* xsect, double y)
2254 {
2255 4962217 return y * y * xsect->sBot;
2256 }
2257
2258 2406740 double triang_getRofY(TXsect* xsect, double y)
2259 {
2260 2406740 return (y * xsect->sBot) / (2. * xsect->rBot);
2261 }
2262
2263 5999861 double triang_getWofY(TXsect* xsect, double y)
2264 {
2265 5999861 return 2.0 * xsect->sBot * y;
2266 }
2267
2268
2269 //=============================================================================
2270 // PARABOLIC fuctions
2271 //=============================================================================
2272
2273 208256 double parab_getYofA(TXsect* xsect, double a)
2274 {
2275 208256 return pow( (3./4.) * a / xsect->rBot, 2./3. );
2276 }
2277
2278 193854 double parab_getRofA(TXsect* xsect, double a)
2279 {
2280
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 193854 times.
 193854 if ( a <= 0.0 ) return 0.0;
2281 193854 return a / parab_getPofY( xsect, parab_getYofA(xsect, a) );
2282 }
2283
2284 222659 double parab_getPofY(TXsect* xsect, double y)
2285 {
2286 222659 double x = 2. * sqrt(y) / xsect->rBot;
2287 222659 double t = sqrt(1.0 + x * x);
2288 222659 return 0.5 * xsect->rBot * xsect->rBot * ( x * t + log(x + t) );
2289 }
2290
2291 93663 double parab_getAofY(TXsect* xsect, double y)
2292 {
2293 93663 return (4./3. * xsect->rBot * y * sqrt(y));
2294 }
2295
2296 28805 double parab_getRofY(TXsect* xsect, double y)
2297 {
2298
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 28805 times.
 28805 if ( y <= 0.0 ) return 0.0;
2299 28805 return parab_getAofY(xsect, y) / parab_getPofY(xsect, y);
2300 }
2301
2302 72010 double parab_getWofY(TXsect* xsect, double y)
2303 {
2304 72010 return 2.0 * xsect->rBot * sqrt(y);
2305 }
2306
2307
2308 //=============================================================================
2309 // POWERFUNC fuctions
2310 //=============================================================================
2311
2312 208196 double powerfunc_getYofA(TXsect* xsect, double a)
2313 {
2314 208196 return pow(a / xsect->rBot, 1.0 / (xsect->sBot + 1.0));
2315 }
2316
2317 193794 double powerfunc_getRofA(TXsect* xsect, double a)
2318 {
2319
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 193794 times.
 193794 if ( a <= 0.0 ) return 0.0;
2320 193794 return a / powerfunc_getPofY(xsect, powerfunc_getYofA(xsect, a));
2321 }
2322
2323 222599 double powerfunc_getPofY(TXsect* xsect, double y)
2324 {
2325 222599 double dy1 = 0.02 * xsect->yFull;
2326 222599 double h = (xsect->sBot + 1.0) * xsect->rBot / 2.0;
2327 222599 double m = xsect->sBot;
2328 222599 double p = 0.0;
2329 222599 double y1 = 0.0;
2330 222599 double x1 = 0.0;
2331 double x2, y2, dx, dy;
2332 do
2333 {
2334 3523489 y2 = y1 + dy1;
2335
2/2
✓ Branch 0 taken 222599 times.
✓ Branch 1 taken 3300890 times.
 3523489 if ( y2 > y ) y2 = y;
2336 3523489 x2 = h * pow(y2, m);
2337 3523489 dx = x2 - x1;
2338 3523489 dy = y2 - y1;
2339 3523489 p += sqrt(dx*dx + dy*dy);
2340 3523489 x1 = x2;
2341 3523489 y1 = y2;
2342
2/2
✓ Branch 0 taken 3300890 times.
✓ Branch 1 taken 222599 times.
 3523489 } while ( y2 < y );
2343 222599 return 2.0 * p;
2344 }
2345
2346 93663 double powerfunc_getAofY(TXsect* xsect, double y)
2347 {
2348 93663 return xsect->rBot * pow(y, xsect->sBot + 1.0);
2349 }
2350
2351 28805 double powerfunc_getRofY(TXsect* xsect, double y)
2352 {
2353
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 28805 times.
 28805 if ( y <= 0.0 ) return 0.0;
2354 28805 return powerfunc_getAofY(xsect, y) / powerfunc_getPofY(xsect, y);
2355 }
2356
2357 72010 double powerfunc_getWofY(TXsect* xsect, double y)
2358 {
2359 72010 return (xsect->sBot + 1.0) * xsect->rBot * pow(y, xsect->sBot);
2360 }
2361
2362
2363 //=============================================================================
2364 // CIRCULAR functions
2365 //=============================================================================
2366
2367 31201365 double circ_getYofA(TXsect* xsect, double a)
2368 {
2369 31201365 double alpha = a / xsect->aFull;
2370
2371 // --- use special function for small a/aFull
2372
2/2
✓ Branch 0 taken 24720495 times.
✓ Branch 1 taken 6480870 times.
 31201365 if ( alpha < 0.04 ) return xsect->yFull * getYcircular(alpha);
2373
2374 // --- otherwise use table
2375 6480870 else return xsect->yFull * lookup(alpha, Y_Circ, N_Y_Circ);
2376 }
2377
2378 2276409 double circ_getAofS(TXsect* xsect, double s)
2379 {
2380 2276409 double psi = s / xsect->sFull;
2381
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2276409 times.
 2276409 if (psi == 0.0) return 0.0;
2382
2/2
✓ Branch 0 taken 94346 times.
✓ Branch 1 taken 2182063 times.
 2276409 if (psi >= 1.0) return xsect->aFull;
2383
2384 // --- use special function for small s/sFull
2385
2/2
✓ Branch 0 taken 1352162 times.
✓ Branch 1 taken 829901 times.
 2182063 if (psi <= 0.015) return xsect->aFull * getAcircular(psi);
2386
2387 // --- otherwise use table
2388 829901 else return xsect->aFull * invLookup(psi, S_Circ, N_S_Circ);
2389 }
2390
2391 145424 double circ_getSofA(TXsect* xsect, double a)
2392 {
2393 145424 double alpha = a / xsect->aFull;
2394
2395 // --- use special function for small a/aFull
2396
2/2
✓ Branch 0 taken 33042 times.
✓ Branch 1 taken 112382 times.
 145424 if ( alpha < 0.04 ) return xsect->sFull * getScircular(alpha);
2397
2398 // --- otherwise use table
2399 else
2400 112382 return xsect->sFull * lookup(alpha, S_Circ, N_S_Circ);
2401 }
2402
2403 75865 double circ_getdSdA(TXsect* xsect, double a)
2404 {
2405 double alpha, theta, p, r, dPdA;
2406
2407 // --- for near-zero area, use generic central difference formula
2408 75865 alpha = a / xsect->aFull;
2409
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 75825 times.
 75865 if ( alpha <= 1.0e-30 ) return 1.0e-30; //generic_getdSdA(xsect, a);
2410
2411 // --- for small a/aFull use analytical derivative
2412
2/2
✓ Branch 0 taken 15393 times.
✓ Branch 1 taken 60432 times.
 75825 else if ( alpha < 0.04 )
2413 {
2414 15393 theta = getThetaOfAlpha(alpha);
2415 15393 p = theta * xsect->yFull / 2.0;
2416 15393 r = a / p;
2417 15393 dPdA = 4.0 / xsect->yFull / (1. - cos(theta));
2418 15393 return (5./3. - (2./3.) * dPdA * r) * pow(r, 2./3.);
2419 }
2420
2421 // --- otherwise use generic tabular getdSdA
2422 60432 else return tabular_getdSdA(xsect, a, S_Circ, N_S_Circ);
2423 }
2424
2425 ////////////////////////////////////////////////
2426 // This is an alternate method used in SWMM 4.4.
2427 ////////////////////////////////////////////////
2428 /*
2429 double circ_getdSdA(TXsect* xsect, double a)
2430 {
2431 double alpha, a1, a2, da, s1, s2, ds;
2432 alpha = a / xsect->aFull;
2433 if ( alpha <= 1.0e-30 ) return 1.0e-30;
2434 da = 0.002;
2435 a1 = alpha - 0.001;
2436 a2 = alpha + 0.001;
2437 if ( a1 < 0.0 )
2438 {
2439 a1 = 0.0;
2440 da = alpha + 0.001;
2441 }
2442 s1 = getScircular(a1);
2443 s2 = getScircular(a2);
2444 ds = (s2 - s1) / da;
2445 if ( ds <= 1.0e-30 ) ds = 1.0e-30;
2446 return xsect->sFull * ds / xsect->aFull;
2447 }
2448 */
2449
2450 9458828 double circ_getAofY(TXsect* xsect, double y)
2451 {
2452 double yNorm;
2453 9458828 yNorm = y / xsect->yFull;
2454 9458828 return xsect->aFull * lookup(yNorm, A_Circ, N_A_Circ);
2455 }
2456
2457
2458 //=============================================================================
2459 // FILLED_CIRCULAR functions
2460 //=============================================================================
2461
2462 5617031 double filled_circ_getYofA(TXsect* xsect, double a)
2463 {
2464 double y;
2465
2466 // --- temporarily remove filled portion of circle
2467 5617031 xsect->yFull += xsect->yBot;
2468 5617031 xsect->aFull += xsect->aBot;
2469 5617031 a += xsect->aBot;
2470
2471 // --- find depth in unfilled circle
2472 5617031 y = circ_getYofA(xsect, a);
2473
2474 // --- restore original values
2475 5617031 y -= xsect->yBot;
2476 5617031 xsect->yFull -= xsect->yBot;
2477 5617031 xsect->aFull -= xsect->aBot;
2478 5617031 return y;
2479 }
2480
2481 3320060 double filled_circ_getAofY(TXsect* xsect, double y)
2482 {
2483 double a;
2484
2485 // --- temporarily remove filled portion of circle
2486 3320060 xsect->yFull += xsect->yBot;
2487 3320060 xsect->aFull += xsect->aBot;
2488 3320060 y += xsect->yBot;
2489
2490 // --- find area of unfilled circle
2491 3320060 a = circ_getAofY(xsect, y);
2492
2493 // --- restore original values
2494 3320060 a -= xsect->aBot;
2495 3320060 xsect->yFull -= xsect->yBot;
2496 3320060 xsect->aFull -= xsect->aBot;
2497 3320060 return a;
2498 }
2499
2500 6138593 double filled_circ_getRofY(TXsect* xsect, double y)
2501 {
2502 double a, r, p;
2503
2504 // --- temporarily remove filled portion of circle
2505 6138593 xsect->yFull += xsect->yBot;
2506 6138593 xsect->aFull += xsect->aBot;
2507 6138593 y += xsect->yBot;
2508
2509 // --- get area, hyd. radius & wetted perimeter of unfilled circle
2510 6138593 a = circ_getAofY(xsect, y);
2511 6138593 r = 0.25 * xsect->yFull * lookup(y/xsect->yFull, R_Circ, N_R_Circ);
2512 6138593 p = (a/r);
2513
2514 // --- reduce area and wetted perimeter by amount of filled circle
2515 // (rBot = filled perimeter, sBot = filled width)
2516 6138593 a = a - xsect->aBot;
2517 6138593 p = p - xsect->rBot + xsect->sBot;
2518
2519 // --- compute actual hyd. radius & restore xsect parameters
2520 6138593 r = a / p;
2521 6138593 xsect->yFull -= xsect->yBot;
2522 6138593 xsect->aFull -= xsect->aBot;
2523 6138593 return r;
2524 }
2525
2526
2527 //=============================================================================
2528 // Special functions for circular cross sections
2529 //=============================================================================
2530
2531 24726627 double getYcircular(double alpha)
2532 {
2533 double theta;
2534
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 24726627 times.
 24726627 if ( alpha >= 1.0 ) return 1.0;
2535
2/2
✓ Branch 0 taken 22629875 times.
✓ Branch 1 taken 2096752 times.
 24726627 if ( alpha <= 0.0 ) return 0.0;
2536
2/2
✓ Branch 0 taken 215 times.
✓ Branch 1 taken 2096537 times.
 2096752 if ( alpha <= 1.0e-5 )
2537 {
2538 215 theta = pow(37.6911*alpha, 1./3.);
2539 215 return theta * theta / 16.0;
2540 }
2541 2096537 theta = getThetaOfAlpha(alpha);
2542 2096537 return (1.0 - cos(theta/2.)) / 2.0;
2543 }
2544
2545 90890 double getScircular(double alpha)
2546 {
2547 double theta;
2548
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 90890 times.
 90890 if ( alpha >= 1.0 ) return 1.0;
2549
2/2
✓ Branch 0 taken 147 times.
✓ Branch 1 taken 90743 times.
 90890 if ( alpha <= 0.0 ) return 0.0;
2550
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 90742 times.
 90743 if ( alpha <= 1.0e-5 )
2551 {
2552 1 theta = pow(37.6911*alpha, 1./3.);
2553 1 return pow(theta, 13./3.) / 124.4797;
2554 }
2555 90742 theta = getThetaOfAlpha(alpha);
2556 90742 return pow((theta - sin(theta)), 5./3.) / (2.0 * PI) / pow(theta, 2./3.);
2557 }
2558
2559 1352162 double getAcircular(double psi)
2560 {
2561 double theta;
2562
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1352162 times.
 1352162 if ( psi >= 1.0 ) return 1.0;
2563
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1352162 times.
 1352162 if ( psi <= 0.0 ) return 0.0;
2564
2/2
✓ Branch 0 taken 4093 times.
✓ Branch 1 taken 1348069 times.
 1352162 if ( psi <= 1.0e-6 )
2565 {
2566 4093 theta = pow(124.4797*psi, 3./13.);
2567 4093 return theta*theta*theta / 37.6911;
2568 }
2569 1348069 theta = getThetaOfPsi(psi);
2570 1348069 return (theta - sin(theta)) / (2.0 * PI);
2571 }
2572
2573 2202672 double getThetaOfAlpha(double alpha)
2574 {
2575 int k;
2576 double theta, theta1, ap, d;
2577
2578
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2202672 times.
 2202672 if ( alpha > 0.04 ) theta = 1.2 + 5.08 * (alpha - 0.04) / 0.96;
2579 2202672 else theta = 0.031715 - 12.79384 * alpha + 8.28479 * sqrt(alpha);
2580 2202672 theta1 = theta;
2581 2202672 ap = (2.0*PI) * alpha;
2582
1/2
✓ Branch 0 taken 6852469 times.
✗ Branch 1 not taken.
 6852469 for (k = 1; k <= 40; k++ )
2583 {
2584 6852469 d = - (ap - theta + sin(theta)) / (1.0 - cos(theta));
2585 // --- modification to improve convergence for large theta
2586
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 6852469 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 6852469 if ( d > 1.0 ) d = SIGN( 1.0, d );
2587 6852469 theta = theta - d;
2588
2/2
✓ Branch 0 taken 2202672 times.
✓ Branch 1 taken 4649797 times.
 6852469 if ( fabs(d) <= 0.0001 ) return theta;
2589 }
2590 ✗ return theta1;
2591 }
2592
2593 1348069 double getThetaOfPsi(double psi)
2594 {
2595 int k;
2596 double theta, theta1, ap, tt, tt23, t3, d;
2597
2598
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1348069 times.
 1348069 if (psi > 0.90) theta = 4.17 + 1.12 * (psi - 0.90) / 0.176;
2599
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1348069 times.
 1348069 else if (psi > 0.5) theta = 3.14 + 1.03 * (psi - 0.5) / 0.4;
2600
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1348069 times.
 1348069 else if (psi > 0.015) theta = 1.2 + 1.94 * (psi - 0.015) / 0.485;
2601 1348069 else theta = 0.12103 - 55.5075 * psi +
2602 1348069 15.62254 * sqrt(psi);
2603 1348069 theta1 = theta;
2604 1348069 ap = (2.0*PI) * psi;
2605
2606
1/2
✓ Branch 0 taken 5309152 times.
✗ Branch 1 not taken.
 5309152 for (k = 1; k <= 40; k++)
2607 {
2608 5309152 theta = fabs(theta);
2609 5309152 tt = theta - sin(theta);
2610 5309152 tt23 = pow(tt, 2./3.);
2611 5309152 t3 = pow(theta, 1./3.);
2612 5309152 d = ap * theta / t3 - tt * tt23;
2613 5309152 d = d / ( ap*(2./3.)/t3 - (5./3.)*tt23*(1.0-cos(theta)) );
2614 5309152 theta = theta - d;
2615
2/2
✓ Branch 0 taken 1348069 times.
✓ Branch 1 taken 3961083 times.
 5309152 if ( fabs(d) <= 0.0001 ) return theta;
2616 }
2617 ✗ return theta1;
2618 }
2619