GCC Code Coverage Report


Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 82.4% 285 / 0 / 346
Functions: 90.5% 19 / 0 / 21
Branches: 57.5% 130 / 0 / 226

table.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // table.c
3 //
4 // Project: EPA SWMM5
5 // Version: 5.2
6 // Date: 10/17/22 (Build 5.2.2)
7 // Author: L. Rossman
8 //
9 // Table (curve and time series) functions.
10 //
11 // Curve and Time Series objects in SWMM 5 are both modeled with
12 // TTable data structures.
13 //
14 // The table_getFirstEntry and table_getNextEntry functions, as well as the
15 // Time Series functions that use them, are not thread safe.
16 //
17 // Update History
18 // ==============
19 // Build 5.1.008:
20 // - The lookup functions used for Curve tables (table_lookup, table_lookupEx,
21 // table_intervalLookup, table_inverseLookup, table_getSlope, table_getMaxY,
22 // table_getArea, and table_getInverseArea) were made thread-safe (thanks to
23 // suggestions by CHI).
24 // Build 5.2.0:
25 // - First line of Curve's input data can contain just the curve name and type.
26 // - The table_getArea function was renamed table_getStorageVolume and was
27 // - refactored.
28 // - The table_getInverseArea function was renamed table_getStorageDepth and
29 // was refactored.
30 // - Support added for relative file names.
31 // Build 5.2.2:
32 // - Prevent re-reading a time series file from start once end is reached.
33 //-----------------------------------------------------------------------------
34 #define _CRT_SECURE_NO_DEPRECATE
35
36 #include <stdlib.h>
37 #include <math.h>
38 #include <string.h>
39 #include "headers.h"
40
41 //-----------------------------------------------------------------------------
42 // Local functions
43 //-----------------------------------------------------------------------------
44 int table_getNextFileEntry(TTable* table, double* x, double* y);
45 int table_parseFileLine(char* line, TTable* table, double* x, double* y);
46 double table_interpolate(double x, double x1, double y1, double x2, double y2);
47
48
49 //=============================================================================
50
51 688123 double table_interpolate(double x, double x1, double y1, double x2, double y2)
52 //
53 // Input: x = x value being interpolated
54 // x1, x2 = x values on either side of x
55 // y1, y2 = y values corrresponding to x1 and x2, respectively
56 // Output: returns the y value corresponding to x
57 // Purpose: interpolates a y value for a given x value.
58 //
59 {
60 688123 double dx = x2 - x1;
61
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 688123 times.
688123 if ( fabs(dx) < 1.0e-20 ) return (y1 + y2) / 2.;
62 688123 return y1 + (x - x1) * (y2 - y1) / dx;
63 }
64
65 //=============================================================================
66
67 1554 int table_readCurve(char* tok[], int ntoks)
68 //
69 // Input: tok[] = array of string tokens
70 // ntoks = number of tokens
71 // Output: returns an error code
72 // Purpose: reads a tokenized line of data for a curve table.
73 //
74 {
75 1554 int j, m, k, k1 = 1;
76 double x, y;
77
78 // --- check for minimum number of tokens
79
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1554 times.
1554 if ( ntoks < 2 ) return error_setInpError(ERR_ITEMS, "");
80
81 // --- check that curve exists in database
82 1554 j = project_findObject(CURVE, tok[0]);
83
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1554 times.
1554 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
84
85 // --- check if this is first line of curve's data
86 // (curve's ID will not have been assigned yet)
87
2/2
✓ Branch 0 taken 221 times.
✓ Branch 1 taken 1333 times.
1554 if ( Curve[j].ID == NULL )
88 {
89 // --- assign ID pointer & curve type
90 221 Curve[j].ID = project_findID(CURVE, tok[0]);
91 221 m = findmatch(tok[1], CurveTypeWords);
92
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 221 times.
221 if ( m < 0 ) return error_setInpError(ERR_KEYWORD, tok[1]);
93 221 Curve[j].curveType = m;
94
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 221 times.
221 if (ntoks == 2) return 0;
95 221 k1 = 2;
96 }
97
98 // --- start reading pairs of X-Y value tokens
99
2/2
✓ Branch 0 taken 1554 times.
✓ Branch 1 taken 1554 times.
3108 for ( k = k1; k < ntoks; k = k+2)
100 {
101
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1554 times.
1554 if ( k+1 >= ntoks ) return error_setInpError(ERR_ITEMS, "");
102
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1554 times.
1554 if ( ! getDouble(tok[k], &x) )
103 return error_setInpError(ERR_NUMBER, tok[k]);
104
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1554 times.
1554 if ( ! getDouble(tok[k+1], &y) )
105 return error_setInpError(ERR_NUMBER, tok[k+1]);
106 1554 table_addEntry(&Curve[j], x, y);
107 }
108 1554 return 0;
109 }
110
111 //=============================================================================
112
113 40970 int table_readTimeseries(char* tok[], int ntoks)
114 //
115 // Input: tok[] = array of string tokens
116 // ntoks = number of tokens
117 // Output: returns an error code
118 // Purpose: reads a tokenized line of data for a time series table.
119 //
120 {
121 int j; // time series index
122 int k; // token index
123 int state; // 1: next token should be a date
124 // 2: next token should be a time
125 // 3: next token should be a value
126 double x, y; // time & value table entries
127 DateTime d; // day portion of date/time value
128 DateTime t; // time portion of date/time value
129 char fname[MAXFNAME + 1];
130
131 // --- check for minimum number of tokens
132
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40970 times.
40970 if ( ntoks < 3 ) return error_setInpError(ERR_ITEMS, "");
133
134 // --- check that time series exists in database
135 40970 j = project_findObject(TSERIES, tok[0]);
136
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40970 times.
40970 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
137
138 // --- if first line of data, assign ID pointer
139
2/2
✓ Branch 0 taken 102 times.
✓ Branch 1 taken 40868 times.
40970 if ( Tseries[j].ID == NULL )
140 102 Tseries[j].ID = project_findID(TSERIES, tok[0]);
141
142 // --- check if time series data is in an external file
143
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 40969 times.
40970 if ( strcomp(tok[1], w_FILE ) )
144 {
145 1 sstrncpy(fname, tok[2], MAXFNAME);
146 1 sstrncpy(Tseries[j].file.name, addAbsolutePath(fname), MAXFNAME);
147 1 Tseries[j].file.mode = USE_FILE;
148 1 return 0;
149 }
150
151 // --- parse each token of input line
152 40969 x = 0.0;
153 40969 k = 1;
154 40969 state = 1; // start off looking for a date
155
2/2
✓ Branch 0 taken 122907 times.
✓ Branch 1 taken 40969 times.
204845 while ( k < ntoks )
156 {
157
3/4
✗ Branch 0 not taken.
✓ Branch 1 taken 40969 times.
✓ Branch 2 taken 40969 times.
✓ Branch 3 taken 40969 times.
122907 switch(state)
158 {
159 40969 case 1: // look for a date entry
160
2/2
✓ Branch 1 taken 40137 times.
✓ Branch 2 taken 832 times.
40969 if ( datetime_strToDate(tok[k], &d) )
161 {
162 40137 Tseries[j].lastDate = d;
163 40137 k++;
164 }
165
166 // --- next token must be a time
167 40969 state = 2;
168 40969 break;
169
170 40969 case 2: // look for a time entry
171
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40969 times.
40969 if ( k >= ntoks ) return error_setInpError(ERR_ITEMS, "");
172
173 // --- first check for decimal hours format
174
2/2
✓ Branch 1 taken 143 times.
✓ Branch 2 taken 40826 times.
40969 if ( getDouble(tok[k], &t) ) t /= 24.0;
175
176 // --- then for an hrs:min format
177
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 40826 times.
40826 else if ( !datetime_strToTime(tok[k], &t) )
178 return error_setInpError(ERR_NUMBER, tok[k]);
179
180 // --- save date + time in x
181 40969 x = Tseries[j].lastDate + t;
182
183 // --- next token must be a numeric value
184 40969 k++;
185 40969 state = 3;
186 40969 break;
187
188 40969 case 3:
189 // --- extract a numeric value from token
190
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40969 times.
40969 if ( k >= ntoks ) return error_setInpError(ERR_ITEMS, "");
191
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 40969 times.
40969 if ( ! getDouble(tok[k], &y) )
192 return error_setInpError(ERR_NUMBER, tok[k]);
193
194 // --- add date/time & value to time series
195 40969 table_addEntry(&Tseries[j], x, y);
196
197 // --- start over looking first for a date
198 40969 k++;
199 40969 state = 1;
200 40969 break;
201 }
202 }
203 40969 return 0;
204 }
205
206 //=============================================================================
207
208 42523 int table_addEntry(TTable* table, double x, double y)
209 //
210 // Input: table = pointer to a TTable structure
211 // x = x value
212 // y = y value
213 // Output: returns TRUE if successful, FALSE if not
214 // Purpose: adds a new x/y entry to a table.
215 //
216 {
217 TTableEntry *entry;
218 42523 entry = (TTableEntry *) malloc(sizeof(TTableEntry));
219
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 42523 times.
42523 if ( !entry ) return FALSE;
220 42523 entry->x = x;
221 42523 entry->y = y;
222 42523 entry->next = NULL;
223
2/2
✓ Branch 0 taken 322 times.
✓ Branch 1 taken 42201 times.
42523 if ( table->firstEntry == NULL ) table->firstEntry = entry;
224 42201 else table->lastEntry->next = entry;
225 42523 table->lastEntry = entry;
226 42523 return TRUE;
227 }
228
229 //=============================================================================
230
231 323 void table_deleteEntries(TTable *table)
232 //
233 // Input: table = pointer to a TTable structure
234 // Output: none
235 // Purpose: deletes all x/y entries in a table.
236 //
237 {
238 TTableEntry *entry;
239 TTableEntry *nextEntry;
240 323 entry = table->firstEntry;
241
2/2
✓ Branch 0 taken 42523 times.
✓ Branch 1 taken 323 times.
42846 while (entry)
242 {
243 42523 nextEntry = entry->next;
244 42523 free(entry);
245 42523 entry = nextEntry;
246 }
247 323 table->firstEntry = NULL;
248 323 table->lastEntry = NULL;
249 323 table->thisEntry = NULL;
250
251
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 322 times.
323 if (table->file.file)
252 {
253 1 fclose(table->file.file);
254 1 table->file.file = NULL;
255 }
256 323 }
257
258 //=============================================================================
259
260 323 void table_init(TTable *table)
261 //
262 // Input: table = pointer to a TTable structure
263 // Output: none
264 // Purpose: initializes properties when table is first created.
265 //
266 {
267 323 table->ID = NULL;
268 323 table->refersTo = -1;
269 323 table->firstEntry = NULL;
270 323 table->lastEntry = NULL;
271 323 table->thisEntry = table->firstEntry;
272 323 table->lastDate = 0.0;
273 323 table->x1 = 0.0;
274 323 table->x2 = 0.0;
275 323 table->y1 = 0.0;
276 323 table->y2 = 0.0;
277 323 table->dxMin = 0.0;
278 323 table->file.mode = NO_FILE;
279 323 table->file.file = NULL;
280 323 table->curveType = -1;
281 323 }
282
283 //=============================================================================
284
285 323 int table_validate(TTable *table)
286 //
287 // Input: table = pointer to a TTable structure
288 // Output: returns error code
289 // Purpose: checks that table's x-values are in ascending order.
290 //
291 {
292 int result;
293 double x1, x2, y1, y2;
294 323 double dx, dxMin = BIG;
295
296 // --- open external file if used as the table's data source
297
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 322 times.
323 if ( table->file.mode == USE_FILE )
298 {
299 1 table->file.file = fopen(table->file.name, "rt");
300
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
1 if ( table->file.file == NULL ) return ERR_TABLE_FILE_OPEN;
301 }
302
303 // --- retrieve the first data entry in the table
304 323 result = table_getFirstEntry(table, &x1, &y1);
305
306 // --- return error condition if external file has no valid data
307
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 323 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
323 if ( !result && table->file.mode == USE_FILE ) return ERR_TABLE_FILE_READ;
308
309 // --- retrieve successive table entries and check for non-increasing x-values
310
2/2
✓ Branch 1 taken 152474 times.
✓ Branch 2 taken 323 times.
152797 while ( table_getNextEntry(table, &x2, &y2) )
311 {
312 152474 dx = x2 - x1;
313
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 152474 times.
152474 if ( dx <= 0.0 )
314 {
315 table->x2 = x2;
316 return ERR_CURVE_SEQUENCE;
317 }
318
2/2
✓ Branch 0 taken 152023 times.
✓ Branch 1 taken 451 times.
152474 dxMin = MIN(dxMin, dx);
319 152474 x1 = x2;
320 }
321 323 table->dxMin = dxMin;
322
323 // --- return error if external file could not be read completely
324
3/4
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 322 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
323 if ( table->file.mode == USE_FILE && !feof(table->file.file) )
325 return ERR_TABLE_FILE_READ;
326 323 return 0;
327 }
328
329 //=============================================================================
330
331 641895 int table_getFirstEntry(TTable *table, double *x, double *y)
332 //
333 // Input: table = pointer to a TTable structure
334 // Output: x = x-value of first table entry
335 // y = y-value of first table entry
336 // returns TRUE if successful, FALSE if not
337 // Purpose: retrieves the first x/y entry in a table.
338 //
339 // NOTE: also moves the current position pointer (thisEntry) to the 1st entry.
340 //
341 {
342 TTableEntry *entry;
343 641895 *x = 0;
344 641895 *y = 0.0;
345
346
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 641893 times.
641895 if ( table->file.mode == USE_FILE )
347 {
348
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
2 if ( table->file.file == NULL ) return FALSE;
349 2 rewind(table->file.file);
350 2 return table_getNextFileEntry(table, x, y);
351 }
352
353 641893 entry = table->firstEntry;
354
1/2
✓ Branch 0 taken 641893 times.
✗ Branch 1 not taken.
641893 if ( entry )
355 {
356 641893 *x = entry->x;
357 641893 *y = entry->y;
358 641893 table->thisEntry = entry;
359 641893 return TRUE;
360 }
361 else return FALSE;
362 }
363
364 //=============================================================================
365
366 3210951 int table_getNextEntry(TTable *table, double *x, double *y)
367 //
368 // Input: table = pointer to a TTable structure
369 // Output: x = x-value of next table entry
370 // y = y-value of next table entry
371 // returns TRUE if successful, FALSE if not
372 // Purpose: retrieves the next x/y entry in a table.
373 //
374 // NOTE: also updates the current position pointer (thisEntry).
375 //
376 {
377 TTableEntry *entry;
378
379
2/2
✓ Branch 0 taken 127598 times.
✓ Branch 1 taken 3083353 times.
3210951 if ( table->file.mode == USE_FILE )
380 127598 return table_getNextFileEntry(table, x, y);
381
382 3083353 entry = table->thisEntry->next;
383
2/2
✓ Branch 0 taken 2444495 times.
✓ Branch 1 taken 638858 times.
3083353 if ( entry )
384 {
385 2444495 *x = entry->x;
386 2444495 *y = entry->y;
387 2444495 table->thisEntry = entry;
388 2444495 return TRUE;
389 }
390 638858 else return FALSE;
391 }
392
393 //=============================================================================
394
395 40992 double table_lookup(TTable *table, double x)
396 //
397 // Input: table = pointer to a TTable structure
398 // x = an x-value
399 // Output: returns a y-value
400 // Purpose: retrieves the y-value corresponding to an x-value in a table,
401 // using interploation if necessary.
402 //
403 // NOTE: if x is below the first table entry, then the first y-value is
404 // returned; if x is above the last entry, then the last y-value is
405 // returned.
406 //
407 {
408 double x1,y1,x2,y2;
409 TTableEntry* entry;
410
411 40992 entry = table->firstEntry;
412
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40992 times.
40992 if ( entry == NULL ) return 0.0;
413 40992 x1 = entry->x;
414 40992 y1 = entry->y;
415
2/2
✓ Branch 0 taken 27136 times.
✓ Branch 1 taken 13856 times.
40992 if ( x <= x1 ) return y1;
416
1/2
✓ Branch 0 taken 38890 times.
✗ Branch 1 not taken.
38890 while ( entry->next )
417 {
418 38890 entry = entry->next;
419 38890 x2 = entry->x;
420 38890 y2 = entry->y;
421
2/2
✓ Branch 0 taken 13856 times.
✓ Branch 1 taken 25034 times.
38890 if ( x <= x2 ) return table_interpolate(x, x1, y1, x2, y2);
422 25034 x1 = x2;
423 25034 y1 = y2;
424 }
425 return y1;
426 }
427
428 //=============================================================================
429
430 27122 double table_getSlope(TTable *table, double x)
431 //
432 // Input: table = pointer to a TTable structure
433 // x = an x-value
434 // Output: returns the slope of the curve at x
435 // Purpose: retrieves the slope of the curve at the line segment containing x.
436 //
437 {
438 double x1,y1,x2,y2;
439 double dx;
440 TTableEntry* entry;
441
442 27122 entry = table->firstEntry;
443
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 27122 times.
27122 if ( entry == NULL ) return 0.0;
444 27122 x1 = entry->x;
445 27122 y1 = entry->y;
446 27122 x2 = x1;
447 27122 y2 = y1;
448
1/2
✓ Branch 0 taken 27122 times.
✗ Branch 1 not taken.
27122 while ( entry->next )
449 {
450 27122 entry = entry->next;
451 27122 x2 = entry->x;
452 27122 y2 = entry->y;
453
1/2
✓ Branch 0 taken 27122 times.
✗ Branch 1 not taken.
27122 if ( x <= x2 ) break;
454 x1 = x2;
455 y1 = y2;
456 }
457 27122 dx = x2 - x1;
458
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 27122 times.
27122 if ( dx == 0.0 ) return 0.0;
459 27122 return (y2 - y1) / dx;
460 }
461
462 //=============================================================================
463
464 67179 double table_lookupEx(TTable *table, double x)
465 //
466 // Input: table = pointer to a TTable structure
467 // x = an x-value
468 // Output: returns a y-value
469 // Purpose: retrieves the y-value corresponding to an x-value in a table,
470 // using interploation if necessary within the table and linear
471 // extrapolation outside of the table.
472 //
473 {
474 double x1,y1,x2,y2;
475 67179 double s = 0.0;
476 TTableEntry* entry;
477
478 67179 entry = table->firstEntry;
479
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 67179 times.
67179 if (entry == NULL ) return 0.0;
480 67179 x1 = entry->x;
481 67179 y1 = entry->y;
482
2/2
✓ Branch 0 taken 2463 times.
✓ Branch 1 taken 64716 times.
67179 if ( x <= x1 )
483 {
484
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2463 times.
2463 if (x1 > 0.0 ) return x/x1*y1;
485 2463 else return y1;
486 }
487
1/2
✓ Branch 0 taken 98783 times.
✗ Branch 1 not taken.
98783 while ( entry->next )
488 {
489 98783 entry = entry->next;
490 98783 x2 = entry->x;
491 98783 y2 = entry->y;
492
1/2
✓ Branch 0 taken 98783 times.
✗ Branch 1 not taken.
98783 if ( x2 != x1 ) s = (y2 - y1) / (x2 - x1);
493
2/2
✓ Branch 0 taken 64716 times.
✓ Branch 1 taken 34067 times.
98783 if ( x <= x2 ) return table_interpolate(x, x1, y1, x2, y2);
494 34067 x1 = x2;
495 34067 y1 = y2;
496 }
497 if ( s < 0.0 ) s = 0.0;
498 return y1 + s*(x - x1);
499 }
500
501 //=============================================================================
502
503 4320 double table_intervalLookup(TTable *table, double x)
504 //
505 // Input: table = pointer to a TTable structure
506 // x = an x-value
507 // Output: returns a y-value
508 // Purpose: retrieves the y-value corresponding to the first table entry
509 // whose x-value is > x.
510 //
511 {
512 TTableEntry* entry;
513
514 4320 entry = table->firstEntry;
515
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4320 times.
4320 if (entry == NULL ) return 0.0;
516
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4320 times.
4320 if ( x < entry->x ) return entry->y;
517
1/2
✓ Branch 0 taken 9895 times.
✗ Branch 1 not taken.
9895 while ( entry->next )
518 {
519 9895 entry = entry->next;
520
2/2
✓ Branch 0 taken 4320 times.
✓ Branch 1 taken 5575 times.
9895 if ( x < entry->x ) return entry->y;
521 }
522 return entry->y;
523 }
524
525 //=============================================================================
526
527 double table_inverseLookup(TTable *table, double y)
528 //
529 // Input: table = pointer to a TTable structure
530 // y = a y-value
531 // Output: returns an x-value
532 // Purpose: retrieves the x-value corresponding to an y-value in a table,
533 // using interploation if necessary.
534 //
535 // NOTE: if y is below the first table entry, then the first x-value is
536 // returned; if y is above the last entry, then the last x-value is
537 // returned.
538 //
539 {
540 double x1,y1,x2,y2;
541 TTableEntry* entry;
542
543 entry = table->firstEntry;
544 if (entry == NULL ) return 0.0;
545 x1 = entry->x;
546 y1 = entry->y;
547 if ( y <= y1 ) return x1;
548 while ( entry->next )
549 {
550 entry = entry->next;
551 x2 = entry->x;
552 y2 = entry->y;
553 if ( y <= y2 ) return table_interpolate(y, y1, x1, y2, x2);
554 x1 = x2;
555 y1 = y2;
556 }
557 return x1;
558 }
559
560 //=============================================================================
561
562 double table_getMaxY(TTable *table, double x)
563 //
564 // Input: table = pointer to a TTable structure
565 // x = an x-value
566 // Output: returns the maximum y-value for x-values below x.
567 // Purpose: finds the largest y value in the initial non-decreasing
568 // portion of a table that appear before value x.
569 //
570 {
571 double ymax;
572 TTableEntry* entry;
573
574 entry = table->firstEntry;
575 ymax = entry->y;
576 while ( x > entry->x && entry->next )
577 {
578 entry = entry->next;
579 if ( entry->y < ymax ) return ymax;
580 ymax = entry->y;
581 }
582 return 0.0;
583 }
584
585 //=============================================================================
586
587 44927 double table_getStorageVolume(TTable *table, double x)
588 //
589 // Input: table = pointer to a TTable structure
590 // x = a depth value
591 // Output: returns a storage volume
592 // Purpose: finds volume for a given depth in a Storage Curve table.
593 //
594 {
595 44927 double a, a1, x1, v, dx = 0.0, dy = 0.0, s;
596 TTableEntry* entry;
597
598 // --- get first entry in table
599 44927 v = 0.0;
600 44927 entry = table->firstEntry;
601
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 44927 times.
44927 if (entry == NULL) return 0.0;
602 44927 x1 = entry->x;
603 44927 a1 = entry->y;
604
605 // --- target depth is below first tabulated depth
606
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 44927 times.
44927 if (x <= x1)
607 {
608 if (x1 < 1.e-6) return 0.0;
609 return (a1/x1) * x * x / 2.0;
610 }
611
612 // --- otherwise traverse table entries until target depth is bracketed
613
2/2
✓ Branch 0 taken 68059 times.
✓ Branch 1 taken 44 times.
68103 while (entry->next)
614 {
615 68059 entry = entry->next;
616 // --- target is bracketed - apply end area method to interpolated area
617
2/2
✓ Branch 0 taken 44883 times.
✓ Branch 1 taken 23176 times.
68059 if (entry->x >= x)
618 {
619 44883 a = table_interpolate(x, x1, a1, entry->x, entry->y);
620 44883 return v + (a1 + a) / 2.0 * (x - x1);
621 }
622 // --- target not yet bracketed so update volume using end area method
623 else
624 {
625 23176 dx = entry->x - x1;
626 23176 dy = entry->y - a1;
627 23176 v = v + (a1 + entry->y) / 2.0 * dx;
628 23176 x1 = entry->x;
629 23176 a1 = entry->y;
630 }
631 }
632
633 // --- extrapolate area if table limit exceeded
634
1/2
✓ Branch 0 taken 44 times.
✗ Branch 1 not taken.
44 if (dx > 1.0e-6)
635 {
636 44 s = dy / dx;
637 44 a = a1 + s * (x - x1);
638 // --- don't extrapolate below 0 in case s is negative
639
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 44 times.
44 if (a < 0.0)
640 {
641 v = v - a1 * a1 / s / 2.0;
642 }
643 // --- apply end area method to extrapolated area
644 44 else v = v + (a1 + a) / 2.0 * (x - x1);
645 }
646 44 return v;
647 }
648
649 //=============================================================================
650
651 5669 double table_getStorageDepth(TTable *table, double v)
652 //
653 // Input: table = pointer to a TTable structure
654 // v = a storage volume
655 // Output: returns a storage depth
656 // Purpose: finds depth for a given volume in a Storage Curve table.
657 //
658 {
659 5669 double a1, a2, d1, d2, dd = 0.0, da = 0.0, v1, v2, s;
660 TTableEntry* entry;
661
662 // --- see if target volume is below that of 1st table entry
663
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5669 times.
5669 if (v == 0.0) return 0.0;
664 5669 entry = table->firstEntry;
665
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5669 times.
5669 if (entry == NULL) return 0.0;
666 5669 d1 = entry->x;
667 5669 a1 = entry->y;
668 5669 v1 = a1 * d1 / 2.0;
669
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5669 times.
5669 if (v <= v1)
670 {
671 if (a1 > 0.0) return sqrt(2.0 * v * d1 / a1);
672 else return 0.0;
673 }
674
675 // --- add next table entry to volume until target volume is bracketed
676
1/2
✓ Branch 0 taken 10477 times.
✗ Branch 1 not taken.
10477 while (entry->next)
677 {
678 10477 entry = entry->next;
679 10477 d2 = entry->x;
680 10477 a2 = entry->y;
681 10477 dd = d2 - d1;
682 10477 da = a2 - a1;
683 10477 v2 = v1 + (a1 + a2) / 2.0 * dd;
684
685 // target volume is bracketed
686
2/2
✓ Branch 0 taken 5669 times.
✓ Branch 1 taken 4808 times.
10477 if (v <= v2)
687 {
688 // --- target coincides with point on curve
689
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5669 times.
5669 if (dd <= 0.0) return d1;
690
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5669 times.
5669 if (da == 0.0)
691 {
692 if (fabs(v2 - v1) < 1.e-6) return d1;
693 else return d1 + dd * (v - v1) / (v2 - v1);
694 }
695 // --- if area decreases with depth then replace point 1 with point 2
696
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5669 times.
5669 if (da < 0.0)
697 {
698 d1 = d2;
699 a1 = a2;
700 v1 = v2;
701 }
702 // --- interpolate between volumes derived from curve
703 5669 s = da / dd;
704 5669 return d1 + (sqrt(a1*a1 + 2.0*s*(v-v1)) - a1) / s;
705 }
706
707 // --- replace point 1 with point 2
708 4808 d1 = d2;
709 4808 a1 = a2;
710 4808 v1 = v2;
711 }
712
713 // --- extrapolate volume if table limit exceeded
714 if (dd == 0.0 || da == 0.0)
715 {
716 if (a1 > 0.0) dd = (v - v1) / a1;
717 else dd = 0.0;
718 }
719 else
720 {
721 s = da / dd;
722 dd = (sqrt(a1*a1 + 2.0*s*(v - v1)) - a1) / s;
723 if (dd < 0.0) dd = 0.0;
724 }
725 return d1 + dd;
726 }
727
728 //=============================================================================
729
730 111 void table_tseriesInit(TTable *table)
731 //
732 // Input: table = pointer to a TTable structure
733 // Output: none
734 // Purpose: initializes the time bracket within a time series table.
735 //
736 {
737 111 table_getFirstEntry(table, &(table->x1), &(table->y1));
738 111 table->x2 = table->x1;
739 111 table->y2 = table->y1;
740 111 table_getNextEntry(table, &(table->x2), &(table->y2));
741 111 }
742
743 //=============================================================================
744
745 1203092 double table_tseriesLookup(TTable *table, double x, char extend)
746 //
747 // Input: table = pointer to a TTable structure
748 // x = a date/time value
749 // extend = TRUE if time series extended on either end
750 // Output: returns a y-value
751 // Purpose: retrieves the y-value corresponding to a time series date,
752 // using interploation if necessary.
753 //
754 // NOTE: if extend is FALSE and date x is outside the range of the table
755 // then 0 is returned; if TRUE then the first or last value is
756 // returned.
757 //
758 {
759 // --- x lies within current time bracket
760
1/2
✓ Branch 0 taken 1203092 times.
✗ Branch 1 not taken.
1203092 if ( table->x1 <= x
761
2/2
✓ Branch 0 taken 564509 times.
✓ Branch 1 taken 638583 times.
1203092 && table->x2 >= x
762
1/2
✓ Branch 0 taken 564509 times.
✗ Branch 1 not taken.
564509 && table->x1 != table->x2 )
763 564509 return table_interpolate(x, table->x1, table->y1, table->x2, table->y2);
764
765 // --- end of external time series file has been reached
766
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 638581 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
638583 if ( table->file.mode == USE_FILE && feof(table->file.file) )
767 {
768 if (extend == TRUE) return table->y1;
769 else return 0;
770 }
771
772 // --- x lies before current time bracket:
773 // move to start of time series
774
3/4
✓ Branch 0 taken 178 times.
✓ Branch 1 taken 638405 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 178 times.
638583 if ( table->x1 == table->x2 || x < table->x1 )
775 {
776 638405 table_getFirstEntry(table, &(table->x1), &(table->y1));
777
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 638405 times.
638405 if ( x < table->x1 )
778 {
779 if ( extend == TRUE ) return table->y1;
780 else return 0;
781 }
782 }
783
784 // --- x lies beyond current time bracket:
785 // update start of next time bracket
786 638583 table->x1 = table->x2;
787 638583 table->y1 = table->y2;
788
789 // --- get end of next time bracket
790
2/2
✓ Branch 1 taken 2413117 times.
✓ Branch 2 taken 638424 times.
3051541 while ( table_getNextEntry(table, &(table->x2), &(table->y2)) )
791 {
792 // --- x lies within the bracket
793
2/2
✓ Branch 0 taken 159 times.
✓ Branch 1 taken 2412958 times.
2413117 if ( x <= table->x2 )
794 159 return table_interpolate(x, table->x1, table->y1,
795 table->x2, table->y2);
796 // --- otherwise move to next time bracket
797 2412958 table->x1 = table->x2;
798 2412958 table->y1 = table->y2;
799 }
800
801 // --- return last value or 0 if beyond last data value
802
2/2
✓ Branch 0 taken 24 times.
✓ Branch 1 taken 638400 times.
638424 if ( extend == TRUE ) return table->y1;
803 638400 else return 0.0;
804 }
805
806 //=============================================================================
807
808 127600 int table_getNextFileEntry(TTable* table, double* x, double* y)
809 //
810 // Input: table = pointer to a TTable structure
811 // x = pointer to a date (as decimal days)
812 // y = pointer to a time series value
813 // Output: updates values of x and y;
814 // returns TRUE if successful, FALSE if not
815 // Purpose: retrieves the next date and value for a time series
816 // table stored in an external file.
817 //
818 {
819 char line[MAXLINE+1];
820 int code;
821
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 127600 times.
127600 if ( table->file.file == NULL ) return FALSE;
822
3/4
✓ Branch 1 taken 127602 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 127601 times.
✓ Branch 5 taken 1 time.
127602 while ( !feof(table->file.file) && fgets(line, MAXLINE, table->file.file) != NULL )
823 {
824 127601 code = table_parseFileLine(line, table, x, y);
825
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 127599 times.
127601 if ( code < 0 ) continue; //skip blank & comment lines
826 127599 return code;
827 }
828 1 return FALSE;
829 }
830
831 //=============================================================================
832
833 127601 int table_parseFileLine(char* line, TTable* table, double* x, double* y)
834 //
835 // Input: table = pointer to a TTable structure
836 // x = pointer to a date (as decimal days)
837 // y = pointer to a time series value
838 // Output: updates values of x and y;
839 // returns -1 if line was a comment,
840 // TRUE if line successfully parsed,
841 // FALSE if line could not be parsed
842 // Purpose: parses a line of time series data from an external file.
843 //
844 {
845 int n;
846 char s1[50],
847 s2[50],
848 s3[50];
849 char* tStr; // time as string
850 char* yStr; // value as string
851 double yy; // value as double
852 DateTime d; // day portion of date/time value
853 DateTime t; // time portion of date/time value
854
855 // --- get 3 string tokens from line and check if its a comment
856 127601 n = sscanf(line, "%s %s %s", s1, s2, s3);
857
858 // --- return if line is blank or is a comment
859 127601 tStr = strtok(line, SEPSTR);
860
3/4
✓ Branch 0 taken 127601 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 127599 times.
127601 if ( tStr == NULL || *tStr == ';' ) return -1;
861
862 // --- line only has a time and a value
863
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 127599 times.
127599 if ( n == 2 )
864 {
865 // --- calendar date is same as last recorded date
866 d = table->lastDate;
867 tStr = s1;
868 yStr = s2;
869 }
870
871 // --- line has date, time and a value
872
1/2
✓ Branch 0 taken 127599 times.
✗ Branch 1 not taken.
127599 else if ( n == 3 )
873 {
874 // --- convert date string to numeric value
875
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 127599 times.
127599 if ( !datetime_strToDate(s1, &d) ) return FALSE;
876
877 // --- update last recorded calendar date
878 127599 table->lastDate = d;
879 127599 tStr = s2;
880 127599 yStr = s3;
881 }
882 else return FALSE;
883
884 // --- convert time string to numeric value
885
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 127599 times.
127599 if ( getDouble(tStr, &t) ) t /= 24.0;
886
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 127599 times.
127599 else if ( !datetime_strToTime(tStr, &t) ) return FALSE;
887
888 // --- convert value string to numeric value
889
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 127599 times.
127599 if ( !getDouble(yStr, &yy) ) return FALSE;
890
891 // --- assign values to current date and value
892 127599 *x = d + t;
893 127599 *y = yy;
894 127599 return TRUE;
895 }
896