transect.c
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | //----------------------------------------------------------------------------- | ||
| 2 | // transect.c | ||
| 3 | // | ||
| 4 | // Project: EPA SWMM5 | ||
| 5 | // Version: 5.2 | ||
| 6 | // Date: 07/13/23 (Build 5.2.4) | ||
| 7 | // Author: L. Rossman | ||
| 8 | // | ||
| 9 | // Geometry processing for irregular cross-section transects. | ||
| 10 | // | ||
| 11 | // Update History | ||
| 12 | // ============== | ||
| 13 | // Build 5.2.0: | ||
| 14 | // - Function added to create a transect for a Street cross-section. | ||
| 15 | // Build 5.2.4: | ||
| 16 | // - Corrected street transect points in transect_createStreetTransect. | ||
| 17 | //----------------------------------------------------------------------------- | ||
| 18 | #define _CRT_SECURE_NO_DEPRECATE | ||
| 19 | |||
| 20 | #include <stdlib.h> | ||
| 21 | #include <string.h> | ||
| 22 | #include <math.h> | ||
| 23 | #include "headers.h" | ||
| 24 | |||
| 25 | //----------------------------------------------------------------------------- | ||
| 26 | // Constants | ||
| 27 | //----------------------------------------------------------------------------- | ||
| 28 | #define MAXSTATION 1500 // max. number of stations in a transect | ||
| 29 | |||
| 30 | //----------------------------------------------------------------------------- | ||
| 31 | // Shared variables | ||
| 32 | //----------------------------------------------------------------------------- | ||
| 33 | static int Ntransects; // total number of transects | ||
| 34 | static int Nstations; // number of stations in current transect | ||
| 35 | static double Station[MAXSTATION+1]; // x-coordinate of each station | ||
| 36 | static double Elev[MAXSTATION+1]; // elevation of each station | ||
| 37 | static double Nleft; // Manning's n for left overbank | ||
| 38 | static double Nright; // Manning's n for right overbank | ||
| 39 | static double Nchannel; // Manning's n for main channel | ||
| 40 | static double Xleftbank; // station where left overbank ends | ||
| 41 | static double Xrightbank; // station where right overbank begins | ||
| 42 | static double Xfactor; // multiplier for station spacing | ||
| 43 | static double Yfactor; // factor added to station elevations | ||
| 44 | static double Lfactor; // main channel/flood plain length | ||
| 45 | |||
| 46 | //----------------------------------------------------------------------------- | ||
| 47 | // External functions (declared in funcs.h) | ||
| 48 | //----------------------------------------------------------------------------- | ||
| 49 | // transect_create (called by createObjects in project.c) | ||
| 50 | // transect_delete (called by deleteObjects in project.c) | ||
| 51 | // transect_readParams (called by parseLine in input.c) | ||
| 52 | // transect_validate (called by input_readData) | ||
| 53 | // transect_createStreetTransect (called by street_readparams) | ||
| 54 | |||
| 55 | //----------------------------------------------------------------------------- | ||
| 56 | // Local functions | ||
| 57 | //----------------------------------------------------------------------------- | ||
| 58 | static int setParams(int transect, char* id, double x[]); | ||
| 59 | static int setManning(double n[]); | ||
| 60 | static int addStation(double x, double y); | ||
| 61 | static double getFlow(int k, double a, double wp, int findFlow); | ||
| 62 | static void createTables(TTransect *transect, double ymin, double ymax); | ||
| 63 | static void getGeometry(TTransect *transect, int i, double y); | ||
| 64 | static void getSliceGeom(int k, double y, double yu, double yd, double *w, | ||
| 65 | double *a, double *wp); | ||
| 66 | static void setMaxSectionFactor(TTransect *transect); | ||
| 67 | |||
| 68 | //============================================================================= | ||
| 69 | |||
| 70 | 58 | int transect_create(int n) | |
| 71 | // | ||
| 72 | // Input: n = number of transect objects to create | ||
| 73 | // Output: returns an error code | ||
| 74 | // Purpose: creates an array of cross-section transects. | ||
| 75 | // | ||
| 76 | { | ||
| 77 | 58 | Ntransects = n; | |
| 78 |
2/2✓ Branch 0 taken 57 times.
✓ Branch 1 taken 1 time.
|
58 | if ( n == 0 ) return 0; |
| 79 | 1 | Transect = (TTransect *) calloc(Ntransects, sizeof(TTransect)); | |
| 80 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( Transect == NULL ) return ERR_MEMORY; |
| 81 | 1 | Nchannel = 0.0; | |
| 82 | 1 | Nleft = 0.0; | |
| 83 | 1 | Nright = 0.0; | |
| 84 | 1 | Nstations = 0; | |
| 85 | 1 | return 0; | |
| 86 | } | ||
| 87 | |||
| 88 | //============================================================================= | ||
| 89 | |||
| 90 | 58 | void transect_delete(void) | |
| 91 | // | ||
| 92 | // Input: none | ||
| 93 | // Output: none | ||
| 94 | // Purpose: deletes memory allocated for all transects. | ||
| 95 | // | ||
| 96 | { | ||
| 97 |
2/2✓ Branch 0 taken 57 times.
✓ Branch 1 taken 1 time.
|
58 | if ( Ntransects == 0 ) return; |
| 98 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | FREE(Transect); |
| 99 | 1 | Ntransects = 0; | |
| 100 | } | ||
| 101 | |||
| 102 | //============================================================================= | ||
| 103 | |||
| 104 | 3 | int transect_readParams(int* count, char* tok[], int ntoks) | |
| 105 | // | ||
| 106 | // Input: count = transect index | ||
| 107 | // tok[] = array of string tokens | ||
| 108 | // ntoks = number of tokens | ||
| 109 | // Output: updated value of count, | ||
| 110 | // returns an error code | ||
| 111 | // Purpose: read parameters of a transect from a tokenized line of input data. | ||
| 112 | // | ||
| 113 | // Format of transect data follows that used for HEC-2 program: | ||
| 114 | // NC nLeft nRight nChannel | ||
| 115 | // X1 name nSta xLeftBank xRightBank 0 0 0 xFactor yFactor | ||
| 116 | // GR Elevation Station ... | ||
| 117 | // | ||
| 118 | { | ||
| 119 | int i, k; | ||
| 120 | 3 | int index = *count; // transect index | |
| 121 | int errcode; // error code | ||
| 122 | double x[10]; // parameter values | ||
| 123 | char* id; // transect ID name | ||
| 124 | |||
| 125 | // --- match first token to a transect keyword | ||
| 126 | 3 | k = findmatch(tok[0], TransectKeyWords); | |
| 127 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if ( k < 0 ) return error_setInpError(ERR_KEYWORD, tok[0]); |
| 128 | |||
| 129 | // --- read parameters associated with keyword | ||
| 130 |
3/4✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
|
3 | switch ( k ) |
| 131 | { | ||
| 132 | // --- NC line: Manning n values | ||
| 133 | 1 | case 0: | |
| 134 | |||
| 135 | // --- finish processing the previous transect | ||
| 136 | 1 | transect_validate(index - 1); | |
| 137 | |||
| 138 | // --- read Manning's n values | ||
| 139 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( ntoks < 4 ) return error_setInpError(ERR_ITEMS, ""); |
| 140 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
|
4 | for (i = 1; i <= 3; i++) |
| 141 | { | ||
| 142 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
|
3 | if ( ! getDouble(tok[i], &x[i]) ) |
| 143 | ✗ | return error_setInpError(ERR_NUMBER, tok[i]); | |
| 144 | } | ||
| 145 | 1 | return setManning(x); | |
| 146 | |||
| 147 | // --- X1 line: identifies start of next transect | ||
| 148 | 1 | case 1: | |
| 149 | |||
| 150 | // --- check that transect was already added to project | ||
| 151 | // (by input_countObjects) | ||
| 152 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( ntoks < 10 ) return error_setInpError(ERR_ITEMS, ""); |
| 153 | 1 | id = project_findID(TRANSECT, tok[1]); | |
| 154 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( id == NULL ) return error_setInpError(ERR_NAME, tok[1]); |
| 155 | |||
| 156 | // --- read in rest of numerical values on data line | ||
| 157 |
2/2✓ Branch 0 taken 8 times.
✓ Branch 1 taken 1 time.
|
9 | for ( i = 2; i < 10; i++ ) |
| 158 | { | ||
| 159 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 8 times.
|
8 | if ( ! getDouble(tok[i], &x[i]) ) |
| 160 | ✗ | return error_setInpError(ERR_NUMBER, tok[i]); | |
| 161 | } | ||
| 162 | |||
| 163 | // --- update total transect count | ||
| 164 | 1 | *count = index + 1; | |
| 165 | |||
| 166 | // --- transfer parameter values to transect's properties | ||
| 167 | 1 | return setParams(index, id, x); | |
| 168 | |||
| 169 | // --- GR line: station elevation & location data | ||
| 170 | 1 | case 2: | |
| 171 | |||
| 172 | // --- check that line contains pairs of data values | ||
| 173 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( (ntoks - 1) % 2 > 0 ) return error_setInpError(ERR_ITEMS, ""); |
| 174 | |||
| 175 | // --- parse each pair of Elevation-Station values | ||
| 176 | 1 | i = 1; | |
| 177 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
|
4 | while ( i < ntoks ) |
| 178 | { | ||
| 179 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
|
3 | if ( ! getDouble(tok[i], &x[1]) ) |
| 180 | ✗ | return error_setInpError(ERR_NUMBER, tok[i]); | |
| 181 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
|
3 | if ( ! getDouble(tok[i+1], &x[2]) ) |
| 182 | ✗ | return error_setInpError(ERR_NUMBER, tok[i+1]); | |
| 183 | 3 | errcode = addStation(x[1], x[2]); | |
| 184 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if ( errcode ) return errcode; |
| 185 | 3 | i += 2; | |
| 186 | } | ||
| 187 | 1 | return 0; | |
| 188 | } | ||
| 189 | ✗ | return 0; | |
| 190 | } | ||
| 191 | |||
| 192 | //============================================================================= | ||
| 193 | |||
| 194 | 9 | void transect_validate(int j) | |
| 195 | // | ||
| 196 | // Input: j = transect index | ||
| 197 | // Output: none | ||
| 198 | // Purpose: validates transect data and creates its geometry tables. | ||
| 199 | // | ||
| 200 | { | ||
| 201 | int i; | ||
| 202 | double ymin, ymax; | ||
| 203 | 9 | double oldNchannel = Nchannel; | |
| 204 | |||
| 205 | // --- check for valid transect data | ||
| 206 |
3/4✓ Branch 0 taken 1 time.
✓ Branch 1 taken 8 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
|
9 | if ( j < 0 || j >= Ntransects ) return; |
| 207 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( Nstations < 2 ) |
| 208 | { | ||
| 209 | ✗ | report_writeErrorMsg(ERR_TRANSECT_TOO_FEW, Transect[j].ID); | |
| 210 | ✗ | return; | |
| 211 | } | ||
| 212 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( Nstations >= MAXSTATION ) |
| 213 | { | ||
| 214 | ✗ | report_writeErrorMsg(ERR_TRANSECT_TOO_MANY, Transect[j].ID); | |
| 215 | ✗ | return; | |
| 216 | } | ||
| 217 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( Nchannel <= 0.0 ) |
| 218 | { | ||
| 219 | ✗ | report_writeErrorMsg(ERR_TRANSECT_MANNING, Transect[j].ID); | |
| 220 | ✗ | return; | |
| 221 | } | ||
| 222 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( Xleftbank > Xrightbank ) |
| 223 | { | ||
| 224 | ✗ | report_writeErrorMsg(ERR_TRANSECT_OVERBANK, Transect[j].ID); | |
| 225 | ✗ | return; | |
| 226 | } | ||
| 227 | |||
| 228 | // --- adjust main channel's Mannings n to make its equivalent | ||
| 229 | // length equal to that of entire flood plain | ||
| 230 | 1 | Nchannel = Nchannel * sqrt(Lfactor); | |
| 231 | 1 | Transect[j].lengthFactor = Lfactor; | |
| 232 | |||
| 233 | // --- find max. depth across transect | ||
| 234 | 1 | ymax = Elev[1]; | |
| 235 | 1 | ymin = Elev[1]; | |
| 236 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | for (i = 2; i <= Nstations; i++) |
| 237 | { | ||
| 238 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | ymax = MAX(Elev[i], ymax); |
| 239 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 1 time.
|
2 | ymin = MIN(Elev[i], ymin); |
| 240 | } | ||
| 241 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( ymin >= ymax ) |
| 242 | { | ||
| 243 | ✗ | report_writeErrorMsg(ERR_TRANSECT_NO_DEPTH, Transect[j].ID); | |
| 244 | ✗ | return; | |
| 245 | } | ||
| 246 | |||
| 247 | // --- add vertical sides to transect to reach full ht. on both ends | ||
| 248 | 1 | Station[0] = Station[1]; | |
| 249 | 1 | Elev[0] = ymax; | |
| 250 | 1 | Nstations++; | |
| 251 | 1 | Station[Nstations] = Station[Nstations-1]; | |
| 252 | 1 | Elev[Nstations] = Elev[0]; | |
| 253 | |||
| 254 | // --- create geometry tables | ||
| 255 | 1 | Transect[j].nTbl = N_TRANSECT_TBL; | |
| 256 | 1 | createTables(&Transect[j], ymin, ymax); | |
| 257 | |||
| 258 | // --- save unadjusted main channel roughness | ||
| 259 | 1 | Transect[j].roughness = oldNchannel; | |
| 260 | } | ||
| 261 | |||
| 262 | //============================================================================= | ||
| 263 | |||
| 264 | 4 | void createTables(TTransect *transect, double ymin, double ymax) | |
| 265 | { | ||
| 266 | int i, nLast; | ||
| 267 | double dy, y; | ||
| 268 | |||
| 269 | 4 | transect->yFull = ymax - ymin; | |
| 270 | 4 | transect->wMax = 0.0; | |
| 271 | |||
| 272 | // --- set 1st table entries to zero | ||
| 273 | 4 | transect->areaTbl[0] = 0.0; | |
| 274 | 4 | transect->hradTbl[0] = 0.0; | |
| 275 | 4 | transect->widthTbl[0] = 0.0; | |
| 276 | |||
| 277 | // --- compute geometry for each depth increment | ||
| 278 | 4 | dy = (ymax - ymin) / ((double)(transect->nTbl) - 1); | |
| 279 | 4 | y = ymin; | |
| 280 |
2/2✓ Branch 0 taken 200 times.
✓ Branch 1 taken 4 times.
|
204 | for (i = 1; i < transect->nTbl; i++) |
| 281 | { | ||
| 282 | 200 | y += dy; | |
| 283 | 200 | transect->areaTbl[i] = 0.0; | |
| 284 | 200 | transect->hradTbl[i] = 0.0; | |
| 285 | 200 | transect->widthTbl[i] = 0.0; | |
| 286 | 200 | getGeometry(transect, i, y); | |
| 287 | } | ||
| 288 | |||
| 289 | // --- determine max. section factor | ||
| 290 | 4 | setMaxSectionFactor(transect); | |
| 291 | |||
| 292 | // --- normalize geometry table entries | ||
| 293 | // (full cross-section values are last table entries) | ||
| 294 | 4 | nLast = transect->nTbl - 1; | |
| 295 | 4 | transect->aFull = transect->areaTbl[nLast]; | |
| 296 | 4 | transect->rFull = transect->hradTbl[nLast]; | |
| 297 | 4 | transect->wMax = transect->widthTbl[nLast]; | |
| 298 | |||
| 299 |
2/2✓ Branch 0 taken 200 times.
✓ Branch 1 taken 4 times.
|
204 | for (i = 1; i <= nLast; i++) |
| 300 | { | ||
| 301 | 200 | transect->areaTbl[i] /= transect->aFull; | |
| 302 | 200 | transect->hradTbl[i] /= transect->rFull; | |
| 303 | 200 | transect->widthTbl[i] /= transect->wMax; | |
| 304 | } | ||
| 305 | |||
| 306 | // --- set width at 0 height equal to width at 4% of max. height | ||
| 307 | 4 | transect->widthTbl[0] = transect->widthTbl[1]; | |
| 308 | 4 | } | |
| 309 | |||
| 310 | //============================================================================= | ||
| 311 | |||
| 312 | 1 | int setManning(double n[]) | |
| 313 | // | ||
| 314 | // Input: n[] = array of Manning's n values | ||
| 315 | // Output: returns an error code | ||
| 316 | // Purpose: sets Manning's n for overbanks and main channel of a transect. | ||
| 317 | // | ||
| 318 | { | ||
| 319 | int i; | ||
| 320 |
2/2✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
|
4 | for (i=1; i<=3; i++) |
| 321 | { | ||
| 322 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if ( n[i] < 0.0 ) return ERR_NUMBER; |
| 323 | } | ||
| 324 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if ( n[1] > 0.0 ) Nleft = n[1]; |
| 325 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if ( n[2] > 0.0 ) Nright = n[2]; |
| 326 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if ( n[3] > 0.0 ) Nchannel = n[3]; |
| 327 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( Nleft == 0.0 ) Nleft = Nchannel; |
| 328 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | if ( Nright == 0.0 ) Nright = Nchannel; |
| 329 | 1 | return 0; | |
| 330 | } | ||
| 331 | |||
| 332 | //============================================================================= | ||
| 333 | |||
| 334 | 1 | int setParams(int j, char* id, double x[]) | |
| 335 | // | ||
| 336 | // Input: j = transect index | ||
| 337 | // id = transect ID name | ||
| 338 | // x[] = array of parameter values | ||
| 339 | // Output: returns an error code | ||
| 340 | // Purpose: assigns parameter values to current transect being processed. | ||
| 341 | // | ||
| 342 | { | ||
| 343 |
2/4✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
|
1 | if ( j < 0 || j >= Ntransects ) return ERR_NUMBER; |
| 344 | 1 | Transect[j].ID = id; // ID name | |
| 345 | 1 | Xleftbank = x[3] / UCF(LENGTH); // left overbank location | |
| 346 | 1 | Xrightbank = x[4] / UCF(LENGTH); // right overbank location | |
| 347 | 1 | Lfactor = x[7]; // channel/bank length | |
| 348 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if ( Lfactor == 0.0 ) Lfactor = 1.0; |
| 349 | 1 | Xfactor = x[8]; // station location multiplier | |
| 350 |
1/2✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
|
1 | if ( Xfactor == 0.0 ) Xfactor = 1.0; |
| 351 | 1 | Xleftbank *= Xfactor; // adjusted left bank | |
| 352 | 1 | Xrightbank *= Xfactor; // adjusted right bank | |
| 353 | 1 | Yfactor = x[9] / UCF(LENGTH); // elevation offset | |
| 354 | 1 | Nstations = 0; | |
| 355 | 1 | return 0; | |
| 356 | } | ||
| 357 | |||
| 358 | //============================================================================= | ||
| 359 | |||
| 360 | 3 | int addStation(double y, double x) | |
| 361 | // | ||
| 362 | // Input: y = station elevation value | ||
| 363 | // x = station distance value | ||
| 364 | // Output: returns an error code | ||
| 365 | // Purpose: adds a new station to the transect currently being processed. | ||
| 366 | // | ||
| 367 | { | ||
| 368 | // --- check for valid number of stations | ||
| 369 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if ( Nstations < 0 ) return ERR_TRANSECT_UNKNOWN; |
| 370 | 3 | Nstations++; | |
| 371 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if ( Nstations >= MAXSTATION ) return 0; |
| 372 | |||
| 373 | // --- add station distance, modified by distance multiplier | ||
| 374 | 3 | Station[Nstations] = x * Xfactor / UCF(LENGTH); | |
| 375 | |||
| 376 | // --- add station elevation, modified by offset elevation | ||
| 377 | 3 | Elev[Nstations] = (y + Yfactor) / UCF(LENGTH); | |
| 378 | |||
| 379 | // --- check if station distances are non-increasing | ||
| 380 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 1 time.
|
3 | if ( Nstations > 1 |
| 381 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
|
2 | && Station[Nstations] < Station[Nstations-1] ) |
| 382 | ✗ | return ERR_TRANSECT_SEQUENCE; | |
| 383 | 3 | return 0; | |
| 384 | } | ||
| 385 | |||
| 386 | //============================================================================= | ||
| 387 | |||
| 388 | 200 | void getGeometry(TTransect *transect, int i, double y) | |
| 389 | // | ||
| 390 | // Input: transect = transect being analyzed | ||
| 391 | // i = index of current entry in geometry tables | ||
| 392 | // y = depth of current entry in geometry tables | ||
| 393 | // Output: none | ||
| 394 | // Purpose: computes entries in a transect's geometry tables at a given depth. | ||
| 395 | // | ||
| 396 | { | ||
| 397 | int k; // station index | ||
| 398 | double ylo, // lower elev. of transect slice | ||
| 399 | yhi, // higher elev. of transect slice | ||
| 400 | w, // top width of transect slice | ||
| 401 | wp, // wetted perimeter of transect slice | ||
| 402 | wpSum, // total wetted perimeter across transect | ||
| 403 | a, // area of transect slice | ||
| 404 | aSum, // total area across transect | ||
| 405 | q, // flow across transect slices with same roughness | ||
| 406 | qSum; // total flow across transect | ||
| 407 | int findFlow; // true if flow thru area slice needs updating | ||
| 408 | |||
| 409 | // --- initialize | ||
| 410 | 200 | wpSum = 0.0; | |
| 411 | 200 | aSum = 0.0; | |
| 412 | 200 | qSum = 0.0; | |
| 413 | |||
| 414 | // --- examine each horizontal station from left to right | ||
| 415 |
2/2✓ Branch 0 taken 1250 times.
✓ Branch 1 taken 200 times.
|
1450 | for (k = 1; k <= Nstations; k++) |
| 416 | { | ||
| 417 | // --- determine low & high elevations for transect sub-section | ||
| 418 |
2/2✓ Branch 0 taken 900 times.
✓ Branch 1 taken 350 times.
|
1250 | if ( Elev[k-1] >= Elev[k] ) |
| 419 | { | ||
| 420 | 900 | yhi = Elev[k-1]; | |
| 421 | 900 | ylo = Elev[k]; | |
| 422 | } | ||
| 423 | else | ||
| 424 | { | ||
| 425 | 350 | yhi = Elev[k]; | |
| 426 | 350 | ylo = Elev[k-1]; | |
| 427 | } | ||
| 428 | |||
| 429 | // --- skip station if its totally dry | ||
| 430 |
2/2✓ Branch 0 taken 343 times.
✓ Branch 1 taken 907 times.
|
1250 | if ( ylo >= y ) continue; |
| 431 | |||
| 432 | // --- get top width, area & wetted perimeter values for transect | ||
| 433 | // slice between station k and k-1 | ||
| 434 | 907 | getSliceGeom(k, y, ylo, yhi, &w, &a, &wp); | |
| 435 | |||
| 436 | // --- update total transect values | ||
| 437 | 907 | wpSum += wp; | |
| 438 | 907 | aSum += a; | |
| 439 | 907 | transect->areaTbl[i] += a; | |
| 440 | 907 | transect->widthTbl[i] += w; | |
| 441 | |||
| 442 | // --- must update flow if station elevation is above water level | ||
| 443 |
2/2✓ Branch 0 taken 196 times.
✓ Branch 1 taken 711 times.
|
907 | if ( Elev[k] >= y ) findFlow = TRUE; |
| 444 | 711 | else findFlow = FALSE; | |
| 445 | |||
| 446 | // --- update flow across transect if called for | ||
| 447 | 907 | q = getFlow(k, aSum, wpSum, findFlow); | |
| 448 |
2/2✓ Branch 0 taken 200 times.
✓ Branch 1 taken 707 times.
|
907 | if ( q > 0.0 ) |
| 449 | { | ||
| 450 | 200 | qSum += q; | |
| 451 | 200 | aSum = 0.0; | |
| 452 | 200 | wpSum = 0.0; | |
| 453 | } | ||
| 454 | |||
| 455 | } // next station k | ||
| 456 | |||
| 457 | // --- find hyd. radius table entry solving Manning eq. with | ||
| 458 | // total flow, total area, and main channel n | ||
| 459 | 200 | aSum = transect->areaTbl[i]; | |
| 460 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 200 times.
|
200 | if ( aSum == 0.0 ) |
| 461 | ✗ | transect->hradTbl[i] = transect->hradTbl[i-1]; | |
| 462 | else | ||
| 463 | 200 | transect->hradTbl[i] = pow(qSum * Nchannel / 1.49 / aSum, 1.5); | |
| 464 | 200 | } | |
| 465 | |||
| 466 | //============================================================================= | ||
| 467 | |||
| 468 | 907 | void getSliceGeom(int k, double y, double ylo, double yhi, double *w, | |
| 469 | double *a, double *wp) | ||
| 470 | // | ||
| 471 | // Input: k = station index | ||
| 472 | // y = water elevation | ||
| 473 | // ylo = transect elevation on low side of slice | ||
| 474 | // yhi = transect elevation on high side of slice | ||
| 475 | // Output w = width of transect slice | ||
| 476 | // a = area of transect slice | ||
| 477 | // wp = wetted perimeter of transect slice | ||
| 478 | // Purpose: finds area, width & wetted perim. for slice of transect that | ||
| 479 | // is covered by given water depth. | ||
| 480 | // | ||
| 481 | // yhi | | ||
| 482 | // | | ||
| 483 | // y |********** | ||
| 484 | // |********** --> slice of transect being analyzed | ||
| 485 | // ylo |**********| | ||
| 486 | // |**********| | ||
| 487 | // |**********| | ||
| 488 | // Station Station | ||
| 489 | // k-1 k | ||
| 490 | // | ||
| 491 | { | ||
| 492 | double width, ratio; | ||
| 493 | |||
| 494 | // --- compute width & wetted perimeter of transect slice | ||
| 495 | 907 | width = fabs(Station[k] - Station[k-1]); | |
| 496 | 907 | (*w) = width; | |
| 497 | 907 | (*wp) = sqrt(width * width + (yhi - ylo) * (yhi - ylo)); | |
| 498 | 907 | (*a) = 0.0; | |
| 499 | |||
| 500 | // --- find area for completely submerged slice | ||
| 501 |
2/2✓ Branch 0 taken 515 times.
✓ Branch 1 taken 392 times.
|
907 | if ( y > yhi ) |
| 502 | { | ||
| 503 | 515 | (*a) = width * ( (y - yhi) + (y - ylo) ) / 2.0; | |
| 504 | } | ||
| 505 | |||
| 506 | // --- otherwise find area and adjust width & wetted perim. for | ||
| 507 | // partly submerged slice | ||
| 508 |
1/2✓ Branch 0 taken 392 times.
✗ Branch 1 not taken.
|
392 | else if ( yhi > ylo ) |
| 509 | { | ||
| 510 | 392 | ratio = (y - ylo) / (yhi - ylo); | |
| 511 | 392 | (*a) = width * (yhi - ylo) / 2.0 * ratio * ratio; | |
| 512 | 392 | (*w) *= ratio; | |
| 513 | 392 | (*wp) *= ratio; | |
| 514 | } | ||
| 515 | 907 | } | |
| 516 | |||
| 517 | //============================================================================= | ||
| 518 | |||
| 519 | 907 | double getFlow(int k, double a, double wp, int findFlow) | |
| 520 | // | ||
| 521 | // Input: k = index of station at end of transect sub-section | ||
| 522 | // a = flow area of sub-section | ||
| 523 | // wp = wetted perimeter of flow area of sub-section | ||
| 524 | // findFlow = TRUE if flow needs updating | ||
| 525 | // Output: returns normal flow (per unit of slope) | ||
| 526 | // Purpose: finds flow through a sub-section of a transect. | ||
| 527 | // | ||
| 528 | { | ||
| 529 | double n; // Manning's n | ||
| 530 | |||
| 531 |
2/2✓ Branch 0 taken 711 times.
✓ Branch 1 taken 196 times.
|
907 | if ( findFlow == FALSE) |
| 532 | { | ||
| 533 | // --- flow needs updating if we are at last station | ||
| 534 |
2/2✓ Branch 0 taken 53 times.
✓ Branch 1 taken 658 times.
|
711 | if ( k == Nstations - 1 ) findFlow = TRUE; |
| 535 | |||
| 536 | // --- flow needs updating if we are at end of left overbank and | ||
| 537 | // there is a change in Manning's n and section not vertical | ||
| 538 |
2/2✓ Branch 0 taken 303 times.
✓ Branch 1 taken 355 times.
|
658 | else if ( Station[k] == Xleftbank ) |
| 539 | { | ||
| 540 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 303 times.
|
303 | if ( Nleft != Nchannel && |
| 541 | ✗ | Station[k] != Station[k-1] ) findFlow = TRUE; | |
| 542 | } | ||
| 543 | |||
| 544 | // --- flow needs updating if we are at start of right overbank and | ||
| 545 | // there is a change in Manning's n and section not vertical | ||
| 546 |
2/2✓ Branch 0 taken 203 times.
✓ Branch 1 taken 152 times.
|
355 | else if ( Station[k] == Xrightbank ) |
| 547 | { | ||
| 548 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 203 times.
|
203 | if ( Nright != Nchannel && |
| 549 | ✗ | Station[k] != Station[k+1] ) findFlow = TRUE; | |
| 550 | } | ||
| 551 | } | ||
| 552 | |||
| 553 | // --- if flow needs updating | ||
| 554 |
2/2✓ Branch 0 taken 249 times.
✓ Branch 1 taken 658 times.
|
907 | if ( findFlow ) |
| 555 | { | ||
| 556 | // --- find value of Manning's n to use | ||
| 557 | 249 | n = Nchannel; | |
| 558 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 249 times.
|
249 | if ( Station[k-1] < Xleftbank ) n = Nleft; |
| 559 |
2/2✓ Branch 0 taken 50 times.
✓ Branch 1 taken 199 times.
|
249 | if ( Station[k] > Xrightbank ) n = Nright; |
| 560 | |||
| 561 | // --- compute flow through flow area | ||
| 562 | // (PHI is the Manning Eqn. constant defined in consts.h) | ||
| 563 | 249 | return PHI / n * a * pow(a/wp, 2./3.); | |
| 564 | } | ||
| 565 | 658 | return 0.0; | |
| 566 | } | ||
| 567 | |||
| 568 | //============================================================================= | ||
| 569 | |||
| 570 | 4 | void setMaxSectionFactor(TTransect *transect) | |
| 571 | // | ||
| 572 | // Input: transect = transect being analyzed | ||
| 573 | // Output: none | ||
| 574 | // Purpose: determines the maximum section factor for a transect and the | ||
| 575 | // area where this maxumum occurs. | ||
| 576 | // | ||
| 577 | { | ||
| 578 | int i; | ||
| 579 | double sf; | ||
| 580 | |||
| 581 | 4 | transect->aMax = 0.0; | |
| 582 | 4 | transect->sMax = 0.0; | |
| 583 |
2/2✓ Branch 0 taken 200 times.
✓ Branch 1 taken 4 times.
|
204 | for (i = 1; i < transect->nTbl; i++) |
| 584 | { | ||
| 585 | 200 | sf = transect->areaTbl[i] * pow(transect->hradTbl[i], 2. / 3.); | |
| 586 |
1/2✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
|
200 | if (sf > transect->sMax) |
| 587 | { | ||
| 588 | 200 | transect->sMax = sf; | |
| 589 | 200 | transect->aMax = transect->areaTbl[i]; | |
| 590 | } | ||
| 591 | } | ||
| 592 | 4 | } | |
| 593 | |||
| 594 | //============================================================================= | ||
| 595 | |||
| 596 | 3 | void transect_createStreetTransect(TStreet* street) | |
| 597 | // | ||
| 598 | { | ||
| 599 | double ymin, ymax, y1, y3, y4; | ||
| 600 | double w1, w2, w3, w4; | ||
| 601 | |||
| 602 | // Point 0 = top of backing | ||
| 603 | // Point 1 = top of curb | ||
| 604 | // Point 2 = bottom of curb | ||
| 605 | // Point 3 = bottom of depressed gutter | ||
| 606 | // Point 4 = top of depressed gutter | ||
| 607 | // Point 5 = street crown | ||
| 608 | |||
| 609 | // --- assign height (y) and width (w) to road & gutter sections | ||
| 610 | 3 | ymin = 0.0; | |
| 611 | 3 | w1 = street->backWidth; | |
| 612 | 3 | w2 = street->gutterWidth; | |
| 613 | 3 | w3 = street->width; | |
| 614 | 3 | w4 = w3 - w2; | |
| 615 | 3 | y3 = street->gutterDepression + street->slope * w2; | |
| 616 | 3 | y1 = street->curbHeight + street->gutterDepression; | |
| 617 | 3 | ymax = street->backSlope * street->backWidth + y1; | |
| 618 | 3 | y4 = y3 + street->slope * w4; | |
| 619 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | ymax = MAX(ymax, y4); |
| 620 | |||
| 621 | // --- assign Station,Elevation points to the street's sections | ||
| 622 | 3 | Station[0] = 0.0; | |
| 623 | 3 | Elev[0] = ymax; | |
| 624 | 3 | Station[1] = w1; | |
| 625 | 3 | Elev[1] = y1; | |
| 626 | 3 | Station[2] = w1; | |
| 627 | 3 | Elev[2] = 0.0; | |
| 628 | 3 | Station[3] = w1 + w2; | |
| 629 | 3 | Elev[3] = y3; | |
| 630 | 3 | Station[4] = w1 + w3; | |
| 631 | 3 | Elev[4] = y4; | |
| 632 | |||
| 633 | // --- a half street ends here | ||
| 634 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 2 times.
|
3 | if (street->sides == 1) |
| 635 | { | ||
| 636 | 1 | Station[5] = Station[4]; | |
| 637 | 1 | Elev[5] = ymax; | |
| 638 | 1 | Nstations = 5; | |
| 639 | 1 | street->transect.nTbl = N_TRANSECT_TBL; | |
| 640 | } | ||
| 641 | |||
| 642 | // --- the right side of a full street mirrors the left side | ||
| 643 | else | ||
| 644 | { | ||
| 645 | 2 | Station[5] = Station[4] + w4; | |
| 646 | 2 | Elev[5] = y3; | |
| 647 | 2 | Station[6] = Station[5] + w2; | |
| 648 | 2 | Elev[6] = 0.0; | |
| 649 | 2 | Station[7] = Station[6]; | |
| 650 | 2 | Elev[7] = y1; | |
| 651 | 2 | Station[8] = Station[7] + w1; | |
| 652 | 2 | Elev[8] = ymax; | |
| 653 | 2 | Nstations = 8; | |
| 654 | 2 | street->transect.nTbl = N_TRANSECT_TBL; | |
| 655 | } | ||
| 656 | |||
| 657 | // --- assign Manning's N to street | ||
| 658 | 3 | Nchannel = street->roughness; | |
| 659 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | if (street->backWidth == 0.0) |
| 660 | { | ||
| 661 | 3 | Nleft = Nchannel; | |
| 662 | 3 | Nright = Nchannel; | |
| 663 | 3 | Xleftbank = Station[0]; | |
| 664 | 3 | Xrightbank = Station[Nstations]; | |
| 665 | } | ||
| 666 | else | ||
| 667 | { | ||
| 668 | ✗ | Nleft = street->backRoughness; | |
| 669 | ✗ | Nright = Nleft; | |
| 670 | ✗ | Xleftbank = Station[1]; | |
| 671 | ✗ | if (street->sides == 2) | |
| 672 | ✗ | Xrightbank = Station[Nstations - 1]; | |
| 673 | else | ||
| 674 | ✗ | Xrightbank = Station[Nstations]; | |
| 675 | } | ||
| 676 | |||
| 677 | // --- create the street's geometry tables | ||
| 678 | 3 | createTables(&(street->transect), ymin, ymax); | |
| 679 | 3 | street->transect.roughness = street->roughness; | |
| 680 | 3 | } | |
| 681 |