1 | """Functions for geometries related to culvert flows |
---|
2 | """ |
---|
3 | |
---|
4 | # Import necessary modules |
---|
5 | from math import sqrt |
---|
6 | from Numeric import array, sum |
---|
7 | |
---|
8 | def create_culvert_polygons(end_point0, |
---|
9 | end_point1, |
---|
10 | width, height=None, |
---|
11 | enquiry_gap_factor=1.0, |
---|
12 | enquiry_shape_factor=2.0): |
---|
13 | """Create polygons at the end of a culvert inlet and outlet. |
---|
14 | At either end two polygons will be created; one for the actual flow to pass through and one a little further away |
---|
15 | for enquiring the total energy at both ends of the culvert and transferring flow. |
---|
16 | |
---|
17 | Input (mandatory): |
---|
18 | end_point0 - one end of the culvert (x,y) |
---|
19 | end_point1 - other end of the culvert (x,y) |
---|
20 | width - culvert width |
---|
21 | |
---|
22 | Input (optional): |
---|
23 | height - culvert height, defaults to width making a square culvert |
---|
24 | enquiry_gap_factor - sets the distance to the enquiry polygon |
---|
25 | enquiry_shape_factor - sets the shape of the enquiry polygon |
---|
26 | (large value widens polygon but reduces height |
---|
27 | to preserve same area as exchange polygon) |
---|
28 | |
---|
29 | Output: |
---|
30 | |
---|
31 | Dictionary of four polygons. The dictionary keys are: |
---|
32 | 'exchange_polygon0' - polygon defining the flow area at end_point0 |
---|
33 | 'exchange_polygon1' - polygon defining the flow area at end_point1 |
---|
34 | 'enquiry_polygon0' - polygon defining the enquiry field near end_point0 |
---|
35 | 'enquiry_polygon1' - polygon defining the enquiry field near end_point1 |
---|
36 | """ |
---|
37 | |
---|
38 | |
---|
39 | # Input check |
---|
40 | if height is None: |
---|
41 | height = width |
---|
42 | |
---|
43 | # Dictionary for calculated polygons |
---|
44 | culvert_polygons = {} |
---|
45 | |
---|
46 | |
---|
47 | # Calculate geometry |
---|
48 | x0, y0 = end_point0 |
---|
49 | x1, y1 = end_point1 |
---|
50 | |
---|
51 | dx = x1-x0 |
---|
52 | dy = y1-y0 |
---|
53 | |
---|
54 | vector = array([dx, dy]) |
---|
55 | length = sqrt(sum(vector**2)) |
---|
56 | |
---|
57 | # Unit direction vector and normal |
---|
58 | vector /= length |
---|
59 | normal = array([-dy, dx])/length |
---|
60 | |
---|
61 | |
---|
62 | # Short hands |
---|
63 | w = width/2*normal # Perpendicular vector of 1/2 width |
---|
64 | h = height*vector # Vector of length=height in the |
---|
65 | # direction of the culvert |
---|
66 | |
---|
67 | # Build exchange polygon 0 |
---|
68 | p0 = end_point0 + w |
---|
69 | p1 = end_point0 - w |
---|
70 | p2 = p1 - h |
---|
71 | p3 = p0 - h |
---|
72 | culvert_polygons['exchange_polygon0'] = array([p0,p1,p2,p3]) |
---|
73 | |
---|
74 | |
---|
75 | # Build exchange polygon 1 |
---|
76 | p0 = end_point1 + w |
---|
77 | p1 = end_point1 - w |
---|
78 | p2 = p1 + h |
---|
79 | p3 = p0 + h |
---|
80 | culvert_polygons['exchange_polygon1'] = array([p0,p1,p2,p3]) |
---|
81 | |
---|
82 | |
---|
83 | |
---|
84 | # Redefine shorthands for enquiry polygons |
---|
85 | w = w*enquiry_shape_factor |
---|
86 | h = h/enquiry_shape_factor |
---|
87 | gap = (enquiry_gap_factor + h)*vector |
---|
88 | |
---|
89 | # Build enquiry polygon 0 |
---|
90 | p0 = end_point0 + w - gap |
---|
91 | p1 = end_point0 - w - gap |
---|
92 | p2 = p1 - h |
---|
93 | p3 = p0 - h |
---|
94 | culvert_polygons['enquiry_polygon0'] = array([p0,p1,p2,p3]) |
---|
95 | |
---|
96 | # Build enquiry polygon 1 |
---|
97 | p0 = end_point1 + w + gap |
---|
98 | p1 = end_point1 - w + gap |
---|
99 | p2 = p1 + h |
---|
100 | p3 = p0 + h |
---|
101 | culvert_polygons['enquiry_polygon1'] = array([p0,p1,p2,p3]) |
---|
102 | |
---|
103 | # Return results |
---|
104 | culvert_polygons['vector'] = vector |
---|
105 | culvert_polygons['length'] = length |
---|
106 | culvert_polygons['normal'] = normal |
---|
107 | return culvert_polygons |
---|