-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathread_polygons.f90
202 lines (127 loc) · 6.3 KB
/
read_polygons.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
module read_polygons
! reads the GMT formatted multisegment files that contain the starting and ending polygons.
! required: file called "params.txt" that contains the number of polygons for each step and the file names of a GMT-style multi-segment
! files (where the individual polygons are separated by a "greater than" character > ) corresponding to those steps. See below:
!
! 1 2
! step1.txt
! step2.txt
implicit none
integer :: max_polygons
integer, parameter :: max_points = 1000000, step_unit = 20, param_unit=10
integer, dimension(2) :: number_polygons
integer, dimension(:,:), allocatable :: polygon_points
double precision, dimension(:,:,:), allocatable :: x_coordinates, y_coordinates
double precision, dimension(:,:), allocatable :: poly_x_min, poly_y_min, poly_y_max, poly_x_max
double precision, save :: fining_increment
character (len=256), dimension(2) :: step_file
contains
subroutine read_params()
open(unit=param_unit, file="params.txt", access="sequential", form="formatted", status="old")
read(param_unit,*) number_polygons(1), number_polygons(2)
read(param_unit,*) step_file(1)
read(param_unit,*) step_file(2)
read(param_unit,*) fining_increment
close(unit=param_unit)
end subroutine read_params
subroutine read_polygons_init()
implicit none
integer :: istat, counter, polygon_counter, add_points, add_counter
double precision :: x, y, distance, angle
character (len=1) :: divider ! make sure that there are no spaces before the ">" character or this will mess up
character (len=1), parameter :: divider_character = ">", ignore_character = "#"
call read_params()
max_polygons = max(number_polygons(1),number_polygons(2))
allocate(x_coordinates(2,max_polygons,max_points), y_coordinates(2,max_polygons,max_points), &
polygon_points(2,max_polygons),poly_x_min(2,max_polygons), poly_y_min(2,max_polygons),&
poly_y_max(2,max_polygons), poly_x_max(2,max_polygons), stat=istat)
if(istat /=0) THEN
write(6,*) "problem allocating arrays in read_polygons"
stop
endif
polygon_points = 0
read_files: do counter = 1, 2, 1
open(unit=step_unit, file=step_file(counter), access="sequential", form="formatted", status="old")
polygon_counter = 0
read_polygons: do
read(step_unit,*, iostat=istat) divider
if(istat /=0) THEN
exit read_polygons
endif
if(divider == divider_character) THEN
polygon_counter = polygon_counter + 1
polygon_points(counter,polygon_counter) = 0
cycle read_polygons
elseif (divider == ignore_character) THEN ! anything starting with '#' should be ignored
cycle read_polygons
else
backspace(unit=step_unit)
polygon_points(counter,polygon_counter) = polygon_points(counter,polygon_counter) + 1
call check_array(polygon_points(counter,polygon_counter))
read(step_unit,*) x, y
if (polygon_points(counter,polygon_counter) > 1) THEN ! check if points need to be added
distance = sqrt(&
(x_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter)-1)-x)**2 + &
(y_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter)-1)-y)**2)
if(distance > fining_increment) THEN ! add points
! currently this just does linear interpolation, probably better in the future to do a spline
add_points = int(distance / fining_increment)
angle = atan2(&
y - y_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter)-1), &
x - x_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter)-1))
do add_counter = 1, add_points
x_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter)) = &
x_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter)-add_counter) + &
cos(angle) * dble(add_counter) / dble(add_points+1) * distance
y_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter)) = &
y_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter)-add_counter) + &
sin(angle) * dble(add_counter) / dble(add_points+1) * distance
polygon_points(counter,polygon_counter) = polygon_points(counter,polygon_counter) + 1
call check_array(polygon_points(counter,polygon_counter))
end do
endif
end if
x_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter)) = x
y_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter)) = y
endif
end do read_polygons
close(unit=step_unit)
end do read_files
! check to see if the start point is the same as the end point. The point_in_polygon routine does not work if it is, so this loop removes it
! Also determines the extreme values for the polygon
do counter = 1, 2, 1
do polygon_counter = 1, number_polygons(counter), 1
if (x_coordinates(counter,polygon_counter,1) == &
x_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter)) .and. &
y_coordinates(counter,polygon_counter,1) == &
y_coordinates(counter,polygon_counter,polygon_points(counter,polygon_counter))) THEN
polygon_points(counter,polygon_counter) = polygon_points(counter,polygon_counter) - 1
endif
poly_x_min(counter,polygon_counter) = minval(x_coordinates(counter,polygon_counter,:))
poly_y_min(counter,polygon_counter) = minval(y_coordinates(counter,polygon_counter,:))
poly_y_max(counter,polygon_counter) = maxval(y_coordinates(counter,polygon_counter,:))
poly_x_max(counter,polygon_counter) = maxval(x_coordinates(counter,polygon_counter,:))
end do
end do
end subroutine read_polygons_init
subroutine read_polygons_clear()
! clears the memory from this subroutine
implicit none
integer :: istat
deallocate(x_coordinates, y_coordinates, polygon_points, poly_x_min, poly_y_min, poly_y_max, poly_x_max, stat=istat)
if(istat /=0) THEN
write(6,*) "problem deallocating arrays in read_polygons"
stop
endif
end subroutine read_polygons_clear
subroutine check_array(point_count)
integer, intent(in) :: point_count
if(point_count > max_points) THEN
write(6,*) "Number of points in polygon exceeds internal memory"
write(6,*) "If you want to proceed, you must increase max_points"
write(6,*) "and recompile:", point_count, "> ", max_points
close(unit=step_unit)
stop
end if
end subroutine check_array
end module read_polygons