VMEC 8.52
3D Equilibrium solver with nested flux surfaces.
Loading...
Searching...
No Matches
line_segment.f
Go to the documentation of this file.
1
3
4!*******************************************************************************
5! File line_segment.f
6!
7! Module is part of the libstell.a. This module containes code to create a
8! profile constructed of line sigments. These line segments are assumed to be
9! specified such that xx(i) < xx(i + 1)
10!
11!*******************************************************************************
12
16 USE stel_kinds
18
20 PRIVATE get_indices, y_value, y_value_int, slope, offset, check
21
22 CONTAINS
23
24!*******************************************************************************
25! Public Functions and Subroutines.
26!*******************************************************************************
27!*******************************************************************************
28! Find the linearly interpolated value of y at position x from line sigments
29! specified by yy(xx).
30!*******************************************************************************
31 SUBROUTINE line_seg(x, y, xx, yy, n)
32 IMPLICIT none
33
34!-------------------------------------------------------------------------------
35! Variable declarations.
36!-------------------------------------------------------------------------------
37 INTEGER, INTENT(in) :: n
38
39 REAL(rprec), INTENT(out) :: y
40 REAL(rprec), INTENT(in) :: x
41 REAL(rprec), DIMENSION(n), INTENT(in) :: xx, yy
42!-------------------------------------------------------------------------------
43! Local Variable declarations.
44!-------------------------------------------------------------------------------
45 INTEGER :: ilow, ihigh
46
47!-------------------------------------------------------------------------------
48! Start of executable code
49!-------------------------------------------------------------------------------
50! Check for valid line segments.
51 IF (n .le. 1) THEN
52 stop "Line sigments require at least two points"
53 ELSE IF (xx(1) .ge. xx(2)) THEN
54 stop "Line sigments must be specified in increasing order of x"
55 END IF
56
57! If x is outside of array, linearly interprolate from the last two points.
58 IF (x .lt. xx(1)) THEN
59 ilow = 1; ihigh = 2
60 ELSE IF (x .gt. xx(n)) THEN
61 ilow = n - 1; ihigh = n
62 ELSE
63 CALL get_indices(x, xx, 1, n, ilow, ihigh)
64 END IF
65
66 y = y_value(x, yy, xx, ilow, ihigh)
67
68 END SUBROUTINE
69
70!*******************************************************************************
71! Find the linearly integrated value of y at position x from line sigments
72! specified by yy(xx).
73!*******************************************************************************
74 SUBROUTINE line_seg_int(x, y, xx, yy, n)
75 IMPLICIT none
76
77!-------------------------------------------------------------------------------
78! Variable declarations.
79!-------------------------------------------------------------------------------
80 INTEGER, INTENT(in) :: n
81
82 REAL(rprec), INTENT(out) :: y
83 REAL(rprec), INTENT(in) :: x
84 REAL(rprec), DIMENSION(n), INTENT(in) :: xx, yy
85!-------------------------------------------------------------------------------
86! Local Variable declarations.
87!-------------------------------------------------------------------------------
88 INTEGER :: ilow, ihigh, i0low, i0high, i
89
90!-------------------------------------------------------------------------------
91! Start of executable code
92!-------------------------------------------------------------------------------
93! Check for valid line segments.
94 IF (n .le. 1) THEN
95 stop "Line sigments require at least two points"
96 ELSE IF (xx(1) .ge. xx(2)) THEN
97 stop "Line sigments must be specified in increasing order of x"
98 END IF
99
100! If x is outside of array, linearly interprolate from the last two points.
101 IF (x .lt. xx(1)) THEN
102 ilow = 1; ihigh = 2
103 ELSE IF (x .gt. xx(n)) THEN
104 ilow = n - 1; ihigh = n
105 ELSE
106 CALL get_indices(x, xx, 1, n, ilow, ihigh)
107 END IF
108
109! Determine where zero is.
110 IF (zero .lt. xx(1)) THEN
111 i0low = 1; i0high = 2
112 ELSE IF (zero .gt. xx(n)) THEN
113 i0low = n - 1; i0high = n
114 ELSE
115 CALL get_indices(zero, xx, 1, n, i0low, i0high)
116 END IF
117
118! If 0 and x have the same low and high indices, then the integrate directly
119! from 0 to x. Other wise Integrate from:
120! 0 to xx(i0high); xx(i0high) to xx(ilow); xx(ilow) to x
121 IF (ilow .eq. i0low) THEN
122 y = y_value_int(zero, x, yy, xx, ilow, ihigh)
123 ELSE
124 y = y_value_int(zero, xx(i0high), yy, xx, i0low, i0high)
125 DO i = i0high, ilow - 1
126 y = y + y_value_int(xx(i), xx(i+1), yy, xx, i, i+1)
127 END DO
128 y = y + y_value_int(xx(ilow), x, yy, xx, ilow, ihigh)
129 END IF
130
131 END SUBROUTINE
132
133!*******************************************************************************
134! Private Functions and Subroutines.
135!*******************************************************************************
136!*******************************************************************************
137! Find bounding indices of the xx array
138!*******************************************************************************
139 PURE RECURSIVE &
140 &SUBROUTINE get_indices(x, xx, lBound, uBound, ilow, ihigh)
141 IMPLICIT none
142
143!-------------------------------------------------------------------------------
144! Variable declarations.
145!-------------------------------------------------------------------------------
146 REAL(rprec), INTENT(in) :: x
147 REAL(rprec), DIMENSION(:), INTENT(in) :: xx
148 INTEGER, INTENT(in) :: lbound, ubound ! Lower and upper bounds of the xx array
149 INTEGER, INTENT(out) :: ilow, ihigh
150!-------------------------------------------------------------------------------
151! Local Variable declarations.
152!-------------------------------------------------------------------------------
153 INTEGER :: hbound ! Half way index of the xx array.
154
155!-------------------------------------------------------------------------------
156! Start of executable code
157!-------------------------------------------------------------------------------
158! Perform a tree search to find the interval that x falls between. This
159! algorthim works by dividing the array into two sub arrays then determining if
160! x is in the lower or upper array. Then the algortrim recurses on the sub array
161! until the size of the array equals 2. This assumes that the array xx is
162! correctly sorted.
163 IF (ubound - lbound .eq. 1) THEN
164 ilow = lbound; ihigh = ubound
165 RETURN
166 END IF
167
168 hbound = (ubound + lbound)/2
169 IF (x .ge. xx(lbound) .and. x .lt. xx(hbound)) THEN
170 CALL get_indices(x, xx, lbound, hbound, ilow, ihigh)
171 ELSE
172 CALL get_indices(x, xx, hbound, ubound, ilow, ihigh)
173 END IF
174 END SUBROUTINE
175
176!*******************************************************************************
177! Find m of y = m*x + b
178!*******************************************************************************
179 PURE FUNCTION slope(yy, xx, ilow, ihigh)
180!-------------------------------------------------------------------------------
181! Variable declarations.
182!-------------------------------------------------------------------------------
183 REAL(rprec), DIMENSION(:), INTENT(in) :: xx, yy
184 INTEGER, INTENT(in) :: ilow, ihigh
185 REAL(rprec) :: slope
186!-------------------------------------------------------------------------------
187! Start of executable code
188!-------------------------------------------------------------------------------
189 slope = (yy(ihigh) - yy(ilow))/(xx(ihigh) - xx(ilow))
190
191 END FUNCTION
192
193!*******************************************************************************
194! Find b of y = m*x + b
195!*******************************************************************************
196 PURE FUNCTION offset(yy, xx, ilow, ihigh)
197!-------------------------------------------------------------------------------
198! Variable declarations.
199!-------------------------------------------------------------------------------
200 REAL(rprec), DIMENSION(:), INTENT(in) :: xx, yy
201 INTEGER, INTENT(in) :: ilow, ihigh
202 REAL(rprec) :: offset
203!-------------------------------------------------------------------------------
204! Start of executable code
205!-------------------------------------------------------------------------------
206 offset = (xx(ihigh)*yy(ilow) - xx(ilow)*yy(ihigh)) &
207 & / (xx(ihigh) - xx(ilow))
208
209 END FUNCTION
210
211!*******************************************************************************
212! Find y(x)
213!*******************************************************************************
214 PURE FUNCTION y_value(x, yy, xx, ilow, ihigh)
215!-------------------------------------------------------------------------------
216! Variable declarations.
217!-------------------------------------------------------------------------------
218 REAL(rprec), DIMENSION(:), INTENT(in) :: xx, yy
219 REAL(rprec), INTENT(in) :: x
220 INTEGER, INTENT(in) :: ilow, ihigh
221 REAL(rprec) :: y_value
222!-------------------------------------------------------------------------------
223! Start of executable code
224!-------------------------------------------------------------------------------
225 y_value = slope(yy, xx, ilow, ihigh)*x &
226 & + offset(yy, xx, ilow, ihigh)
227
228 END FUNCTION
229
230!*******************************************************************************
231! Find Int[y(x),x]
232!*******************************************************************************
233 PURE FUNCTION y_value_int(x0, x1, yy, xx, ilow, ihigh)
234!-------------------------------------------------------------------------------
235! Variable declarations.
236!-------------------------------------------------------------------------------
237 REAL(rprec), DIMENSION(:), INTENT(in) :: xx, yy
238 REAL(rprec), INTENT(in) :: x0, x1
239 INTEGER, INTENT(in) :: ilow, ihigh
240 REAL(rprec) :: y_value_int
241!-------------------------------------------------------------------------------
242! Start of executable code
243!-------------------------------------------------------------------------------
244 y_value_int = slope(yy, xx, ilow, ihigh)/2.0*(x1**2.0 - x0**2.0) &
245 & + offset(yy, xx, ilow, ihigh)*(x1 - x0)
246
247 END FUNCTION
248
249!*******************************************************************************
250! Unit Test Functions
251!*******************************************************************************
252!*******************************************************************************
253! Main test function
254!*******************************************************************************
255 FUNCTION line_seg_test()
256!-------------------------------------------------------------------------------
257! Variable declarations.
258!-------------------------------------------------------------------------------
259 LOGICAL :: line_seg_test
260 REAL(rprec) :: result
261
262!-------------------------------------------------------------------------------
263! Start of executable code
264!-------------------------------------------------------------------------------
265! Test slope function. yy(0.0,1.0) and xx(0.0,1.0) should have a slope of 1.0
266 result = slope((/ 0.0d+0, 1.0d+0 /), &
267 & (/ 0.0d+0, 1.0d+0 /), &
268 & 1, 2)
269 line_seg_test = check(1.0d+0, result, 1, "slope()")
270 IF (.not.line_seg_test) RETURN
271! Test slope function. yy(1.0,0.0) and xx(0.0,1.0) should have a slope of -1.0
272 result = slope((/ 1.0d+0, 0.0d+0 /), &
273 & (/ 0.0d+0, 1.0d+0 /), &
274 & 1, 2)
275 line_seg_test = check(-1.0d+0, result, 2, "slope()")
276 IF (.not.line_seg_test) RETURN
277! Test slope function. yy(1.0,1.0) and xx(0.0,1.0) should have a slope of 0.0
278 result = slope((/ 1.0d+0, 1.0d+0 /), &
279 & (/ 0.0d+0, 1.0d+0 /), &
280 & 1, 2)
281 line_seg_test = check(0.0d+0, result, 3, "slope()")
282 IF (.not.line_seg_test) RETURN
283
284! Test slope function. yy(0.0,1.0) and xx(1.0,2.0) should have a slope of 1.0
285 result = slope((/ 0.0d+0, 1.0d+0 /), &
286 & (/ 1.0d+0, 2.0d+0 /), &
287 & 1, 2)
288 line_seg_test = check(1.0d+0, result, 4, "slope()")
289
290! Test offset function. yy(0.0,1.0) and xx(0.0,1.0) should have an offset of 0.0
291 result = offset((/ 0.0d+0, 1.0d+0 /), &
292 & (/ 0.0d+0, 1.0d+0 /), &
293 & 1, 2)
294 line_seg_test = check(0.0d+0, result, 1, "offset()")
295 IF (.not.line_seg_test) RETURN
296! Test offset function. yy(1.0,0.0) and xx(0.0,1.0) should have an offset of 1.0
297 result = offset((/ 1.0d+0, 0.0d+0 /), &
298 & (/ 0.0d+0, 1.0d+0 /), &
299 & 1, 2)
300 line_seg_test = check(1.0d+0, result, 2, "offset()")
301 IF (.not.line_seg_test) RETURN
302! Test offset function. yy(2.0,1.0) and xx(1.0,2.0) should have an offset of 3.0
303 result = offset((/ 2.0d+0, 1.0d+0 /), &
304 & (/ 1.0d+0, 2.0d+0 /), &
305 & 1, 2)
306 line_seg_test = check(3.0d+0, result, 3, "offset()")
307 IF (.not.line_seg_test) RETURN
308
309! Test y_value function. yy(0.0,2.0) and xx(0.0,2.0) should have a value of 1.0
310! at x = 1.0
311 result = y_value(1.0d+0, &
312 & (/ 0.0d+0, 2.0d+0 /), &
313 & (/ 0.0d+0, 2.0d+0 /), &
314 & 1, 2)
315 line_seg_test = check(1.0d+0, result, 1, "y_value()")
316 IF (.not.line_seg_test) RETURN
317! Test y_value function. yy(1.0,2.0) and xx(1.0,2.0) should have a value of 0.0
318! at x = 0.0
319 result = y_value(0.0d+0, &
320 & (/ 1.0d+0, 2.0d+0 /), &
321 & (/ 1.0d+0, 2.0d+0 /), &
322 & 1, 2)
323 line_seg_test = check(0.0d+0, result, 2, "y_value()")
324 IF (.not.line_seg_test) RETURN
325! Test y_value function. yy(0.0,1.0) and xx(0.0,1.0) should have a value of 3.0
326! at x = 3.0
327 result = y_value(3.0d+0, &
328 & (/ 0.0d+0, 1.0d+0 /), &
329 & (/ 0.0d+0, 1.0d+0 /), &
330 & 1, 2)
331 line_seg_test = check(3.0d+0, result, 3, "y_value()")
332 IF (.not.line_seg_test) RETURN
333
334! Test y_value_int function. yy(1.0,1.0) and xx(0.0,1.0) should have a value of
335! 1.0 at x1 = 1.0 to x0 = 0.0
336 result = y_value_int(0.0d+0, 1.0d+0, &
337 & (/ 1.0d+0, 1.0d+0 /), &
338 & (/ 0.0d+0, 1.0d+0 /), &
339 & 1, 2)
340 line_seg_test = check(1.0d+0, result, 1, "y_value_int()")
341 IF (.not.line_seg_test) RETURN
342
343! Test y_value_int function. yy(1.0,1.0) and xx(1.0,2.0) should have a value of
344! 1.0 at x1 = 1.0 to x0 = 0.0
345 result = y_value_int(0.0d+0, 1.0d+0, &
346 & (/ 1.0d+0, 1.0d+0 /), &
347 & (/ 1.0d+0, 2.0d+0 /), &
348 & 1, 2)
349 line_seg_test = check(1.0d+0, result, 2, "y_value_int()")
350 IF (.not.line_seg_test) RETURN
351! Test y_value_int function. yy(1.0,1.0) and xx(0.0,1.0) should have a value of
352! 1.0 at x1 = 2.0 to x0 = 1.0
353 result = y_value_int(1.0d+0, 2.0d+0, &
354 & (/ 1.0d+0, 1.0d+0 /), &
355 & (/ 0.0d+0, 1.0d+0 /), &
356 & 1, 2)
357 line_seg_test = check(1.0d+0, result, 3, "y_value_int()")
358 IF (.not.line_seg_test) RETURN
359! Test y_value_int function. yy(0.0,1.0) and xx(0.0,1.0) should have a value of
360! 0.5 at x1 = 1.0 to x0 = 0.0
361 result = y_value_int(0.0d+0, 1.0d+0, &
362 & (/ 0.0d+0, 1.0d+0 /), &
363 & (/ 0.0d+0, 1.0d+0 /), &
364 & 1, 2)
365 line_seg_test = check(0.5d+0, result, 4, "y_value_int()")
366 IF (.not.line_seg_test) RETURN
367! Test y_value_int function. yy(1.0,2.0) and xx(0.0,1.0) should have a value of
368! 1.5 at x1 = 1.0 to x0 = 0.0
369 result = y_value_int(0.0d+0, 1.0d+0, &
370 & (/ 1.0d+0, 2.0d+0 /), &
371 & (/ 0.0d+0, 1.0d+0 /), &
372 & 1, 2)
373 line_seg_test = check(1.5d+0, result, 5, "y_value_int()")
374 IF (.not.line_seg_test) RETURN
375! Test y_value_int function. yy(1.0,2.0) and xx(1.0,2.0) should have a value of
376! 0.5 at x1 = 1.0 to x0 = 0.0
377 result = y_value_int(0.0d+0, 1.0d+0, &
378 & (/ 1.0d+0, 2.0d+0 /), &
379 & (/ 1.0d+0, 2.0d+0 /), &
380 & 1, 2)
381 line_seg_test = check(0.5d+0, result, 6, "y_value_int()")
382 IF (.not.line_seg_test) RETURN
383! Test y_value_int function. yy(0.0,1.0) and xx(0.0,1.0) should have a value of
384! 1.5 at x1 = 2.0 to x0 = 1.0
385 result = y_value_int(1.0d+0, 2.0d+0, &
386 & (/ 1.0d+0, 2.0d+0 /), &
387 & (/ 1.0d+0, 2.0d+0 /), &
388 & 1, 2)
389 line_seg_test = check(1.5d+0, result, 7, "y_value_int()")
390 IF (.not.line_seg_test) RETURN
391
392! Test line_seg function. yy(1.0,1.0,0.0,0.0) and xx(0.0,1.0,2.0,3.0) should
393! have a value of 1.0 at x = 0.5
394 call line_seg(0.5d+0, result, &
395 & (/ 0.0d+0, 1.0d+0, 2.0d+0, 3.0d+0 /), &
396 & (/ 1.0d+0, 1.0d+0, 0.0d+0, 0.0d+0 /), &
397 & 4)
398 line_seg_test = check(1.0d+0, result, 1, "line_seg()")
399 IF (.not.line_seg_test) RETURN
400! Test line_seg function. yy(1.0,1.0,0.0,0.0) and xx(0.0,1.0,2.0,3.0) should
401! have a value of 0.0 at x = 2.5
402 call line_seg(2.5d+0, result, &
403 & (/ 0.0d+0, 1.0d+0, 2.0d+0, 3.0d+0 /), &
404 & (/ 1.0d+0, 1.0d+0, 0.0d+0, 0.0d+0 /), &
405 & 4)
406 line_seg_test = check(0.0d+0, result, 2, "line_seg()")
407 IF (.not.line_seg_test) RETURN
408! Test line_seg function. yy(1.0,1.0,0.0,0.0) and xx(0.0,1.0,2.0,3.0) should
409! have a value of 0.5 at x = 1.5
410 call line_seg(1.5d+0, result, &
411 & (/ 0.0d+0, 1.0d+0, 2.0d+0, 3.0d+0 /), &
412 & (/ 1.0d+0, 1.0d+0, 0.0d+0, 0.0d+0 /), &
413 & 4)
414 line_seg_test = check(0.5d+0, result, 3, "line_seg()")
415 IF (.not.line_seg_test) RETURN
416! Test line_seg function. yy(1.0,1.0,0.0,0.0) and xx(1.0,2.0,3.0,4.0) should
417! have a value of 1.0 at x = 0.0
418 call line_seg(0.0d+0, result, &
419 & (/ 1.0d+0, 2.0d+0, 3.0d+0, 4.0d+0 /), &
420 & (/ 1.0d+0, 1.0d+0, 0.0d+0, 0.0d+0 /), &
421 & 4)
422 line_seg_test = check(1.0d+0, result, 4, "line_seg()")
423 IF (.not.line_seg_test) RETURN
424! Test line_seg function. yy(1.0,1.0,0.0,0.0) and xx(0.0,1.0,2.0,3.0) should
425! have a value of 0.0 at x = 4.0
426 call line_seg(4.0d+0, result, &
427 & (/ 0.0d+0, 1.0d+0, 2.0d+0, 3.0d+0 /), &
428 & (/ 1.0d+0, 1.0d+0, 0.0d+0, 0.0d+0 /), &
429 & 4)
430 line_seg_test = check(0.0d+0, result, 5, "line_seg()")
431 IF (.not.line_seg_test) RETURN
432
433! Test line_seg_int function. yy(1.0,1.0,0.0,0.0) and xx(0.0,1.0,2.0,3.0) should
434! have a value of 0.5 at x = 0.5
435 call line_seg_int(0.5d+0, result, &
436 & (/ 0.0d+0, 1.0d+0, 2.0d+0, 3.0d+0 /), &
437 & (/ 1.0d+0, 1.0d+0, 0.0d+0, 0.0d+0 /), &
438 & 4)
439 line_seg_test = check(0.5d+0, result, 1, "line_seg_int()")
440 IF (.not.line_seg_test) RETURN
441! Test line_seg_int function. yy(1.0,1.0,0.0,0.0) and xx(0.0,1.0,2.0,3.0) should
442! have a value of 1.375 at x = 1.5
443 call line_seg_int(1.5d+0, result, &
444 & (/ 0.0d+0, 1.0d+0, 2.0d+0, 3.0d+0 /), &
445 & (/ 1.0d+0, 1.0d+0, 0.0d+0, 0.0d+0 /), &
446 & 4)
447 line_seg_test = check(1.375d+0, result, 2, "line_seg_int()")
448 IF (.not.line_seg_test) RETURN
449! Test line_seg_int function. yy(1.0,1.0,0.0,0.0) and xx(0.0,1.0,2.0,3.0) should
450! have a value of 1.5 at x = 2.5
451 call line_seg_int(2.5d+0, result, &
452 & (/ 0.0d+0, 1.0d+0, 2.0d+0, 3.0d+0 /), &
453 & (/ 1.0d+0, 1.0d+0, 0.0d+0, 0.0d+0 /), &
454 & 4)
455 line_seg_test = check(1.5d+0, result, 3, "line_seg_int()")
456 IF (.not.line_seg_test) RETURN
457! Test line_seg_int function. yy(1.0,1.0,0.0,0.0) and xx(1.0,2.0,3.0,4.0) should
458! have a value of 0.5 at x = 0.5
459 call line_seg_int(0.5d+0, result, &
460 & (/ 1.0d+0, 2.0d+0, 3.0d+0, 4.0d+0 /), &
461 & (/ 1.0d+0, 1.0d+0, 0.0d+0, 0.0d+0 /), &
462 & 4)
463 line_seg_test = check(0.5d+0, result, 4, "line_seg_int()")
464 IF (.not.line_seg_test) RETURN
465! Test line_seg_int function. yy(1.0,0.0,0.0,1.0) and xx(0.0,1.0,2.0,3.0) should
466! have a value of 2.5 at x = 4.0
467 call line_seg_int(4.0d+0, result, &
468 & (/ 0.0d+0, 1.0d+0, 2.0d+0, 3.0d+0 /), &
469 & (/ 1.0d+0, 0.0d+0, 0.0d+0, 1.0d+0 /), &
470 & 4)
471 line_seg_test = check(2.5d+0, result, 5, "line_seg_int()")
472 IF (.not.line_seg_test) RETURN
473! Test line_seg_int function. yy(1.0,0.0,0.0,1.0) and xx(1.0,2.0,3.0,4.0) should
474! have a value of 4.0 at x = 5.0
475 call line_seg_int(5.0d+0, result, &
476 & (/ 1.0d+0, 2.0d+0, 3.0d+0, 4.0d+0 /), &
477 & (/ 1.0d+0, 0.0d+0, 0.0d+0, 1.0d+0 /), &
478 & 4)
479 line_seg_test = check(4.0d+0, result, 6, "line_seg_int()")
480 IF (.not.line_seg_test) RETURN
481
482 END FUNCTION
483
484!*******************************************************************************
485! Check Test result
486!*******************************************************************************
487 FUNCTION check(expected, received, testNum, name)
488!-------------------------------------------------------------------------------
489! Variable declarations.
490!-------------------------------------------------------------------------------
491 LOGICAL :: check
492 REAL(rprec), INTENT(in) :: expected, received
493 INTEGER, INTENT(in) :: testnum
494 CHARACTER (LEN=*), INTENT(in) :: name
495!-------------------------------------------------------------------------------
496! Start of executable code
497!-------------------------------------------------------------------------------
498 check = expected .eq. received
499 IF (.not.check) THEN
500 write(*,*) "line_segment_mod.f: ", name, " test", testnum, &
501 & "failed."
502 write(*,*) "Expected", expected, "Received", received
503 END IF
504
505 END FUNCTION
506
507 END MODULE line_segment
This module containes code to create a profile constructed of line segments. These line segments are ...
subroutine, public line_seg_int(x, y, xx, yy, n)
subroutine, public line_seg(x, y, xx, yy, n)
logical function, public line_seg_test()
real(dp), parameter zero
integer, parameter rprec