source: branches/publications/ORCHIDEE_gmd_2018_MICT-LEAK/src_global/polygones.f90 @ 6892

Last change on this file since 6892 was 3953, checked in by albert.jornet, 8 years ago

Merge: from trunk revisions [3917:3946/trunk/ORCHIDEE]
File size: 22.9 KB
Line 
1
2MODULE polygones
3  !
4  ! This module is there to handle polygones. It is an essential tool in interpolation methods.
5  !
6  USE defprec
7
8  IMPLICIT NONE
9
10  PRIVATE
11  PUBLIC :: polygones_pointinside, polygones_extend, polygones_intersection, polygones_cleanup, &
12       &    polygones_area, polygones_crossing, polygones_convexhull
13  !
14  REAL(r_std),PARAMETER                     :: zero=0.0
15  !
16CONTAINS
17!
18!=======================================================================================================
19!
20  SUBROUTINE polygones_pointinside(nvert_in, poly, point_x, point_y, inside)
21    !
22    ! Routine to determine of point (point_x, point_y) is inside of polygone poly
23    ! Mathode based on docume http://www.ecse.rpi.edu/Homepages/wrf/Research/Short_Notes/pnpoly.html
24    !
25    ! Arguments
26    INTEGER(i_std), INTENT(in)              :: nvert_in
27    REAL(r_std), DIMENSION(:,:), INTENT(in) :: poly
28    REAL(r_std), INTENT(in)                 :: point_x, point_y
29    LOGICAL, INTENT(out)                    :: inside
30    !
31    ! Local
32    INTEGER(i_std)                          :: i, j, nvert
33    REAL(r_std)                             :: xonline
34    !
35    !
36    inside = .FALSE.
37    !
38    nvert=SIZE(poly,DIM=1)
39    IF ( nvert < nvert_in) THEN
40       WRITE(*,*) "Input polygone is smaller than the size proposed in the interface."
41       STOP "polygones_pointinside"
42    ENDIF
43    IF (SIZE(poly,DIM=2) .NE. 2) THEN
44       WRITE(*,*) "This cannot be a polygone :", SIZE(poly)
45       STOP "polygones_pointinside"
46    ENDIF
47    !
48    j = nvert_in
49    DO i=1,nvert_in
50       IF ( (poly(i,2) > point_y) .NEQV. (poly(j,2) > point_y) ) THEN
51          xonline = (poly(j,1)-poly(i,1))*(point_y-poly(i,2))/(poly(j,2)-poly(i,2))+poly(i,1)
52          IF ( point_x < xonline ) THEN
53             inside = (.NOT. inside)
54          ENDIF
55       ENDIF
56       j = i
57    ENDDO
58    !
59  END SUBROUTINE polygones_pointinside
60!
61!=======================================================================================================
62!
63  SUBROUTINE polygones_lineintersect(la, lb, intersection, point_x, point_y)
64    !
65    ! Simple alogorith tho determine the intersections of 2 lines.
66    !
67    ! ARGUMENTS
68    REAL(r_std), DIMENSION(2,2), INTENT(in)    :: la, lb
69    LOGICAL, INTENT(out)                       :: intersection
70    REAL(r_std), INTENT(out)                   :: point_x, point_y
71    !
72    ! Local
73    REAL(r_std)                                :: den, ua, ub
74    REAL(r_std)                                :: x1, x2, x3, x4, y1, y2, y3, y4
75    !
76    intersection = .FALSE.
77    !
78     x1 = la(1,1)
79     y1 = la(1,2)
80     x2 = la(2,1)
81     y2 = la(2,2)
82
83     x3 = lb(1,1)
84     y3 = lb(1,2)
85     x4 = lb(2,1)
86     y4 = lb(2,2)
87
88     den = (y4-y3)*(x2-x1)-(x4-x3)*(y2-y1)
89
90     IF ( ABS(den) > EPSILON(den) ) THEN
91        ua = ((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3))/den
92        ub = ((x2-x1)*(y1-y3)-(y2-y1)*(x1-x3))/den
93        IF ( ua >= 0 .AND. ua <= 1 .AND. ub >= 0 .AND. ub <= 1) THEN
94           intersection = .TRUE.
95           point_x = x1 + ua*(x2-x1)
96           point_y = y1 + ua*(y2-y1)
97        ENDIF
98     ENDIF
99
100  END SUBROUTINE polygones_lineintersect
101!
102!=======================================================================================================
103!
104  SUBROUTINE polygones_extend(nvert_in, poly_in, nbdots, nvert_out, poly_out)
105    !
106    ! This function simply add nbdots onto each side of the polygone.
107    !
108    ! Arguments
109    INTEGER(i_std), INTENT(in)               :: nvert_in
110    REAL(r_std), DIMENSION(:,:), INTENT(in)  :: poly_in
111    INTEGER(i_std), INTENT(in)               :: nbdots
112    INTEGER(i_std), INTENT(out)              :: nvert_out
113    REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out
114    !
115    ! Local
116    !
117    INTEGER(i_std)                        :: nvert, nvert_tmp, i, j, id, ipos
118    REAL(r_std)                           :: xs, xe, ys, ye
119    !
120    nvert=SIZE(poly_in,DIM=1)
121    IF ( nvert < nvert_in ) THEN
122       WRITE(*,*) "Input polygone is smaller than the size proposed in the interface."
123       STOP "polygones_extend"
124    ENDIF
125    IF (SIZE(poly_in,DIM=2) .NE. 2) THEN
126       WRITE(*,*) "Poly_in cannot be a polygone :", SIZE(poly_in)
127       STOP "polygones_extend"
128    ENDIF
129    nvert_tmp=SIZE(poly_out,DIM=1)
130    IF (SIZE(poly_out,DIM=2) .NE. 2) THEN
131       WRITE(*,*) "Poly_out cannot be a polygone :", SIZE(poly_out)
132       STOP "polygones_extend"
133    ENDIF
134    !
135    IF ( nvert_tmp < nvert_in*nbdots) THEN
136       WRITE(*,*) "Poly_out is too small for extending it :", SIZE(poly_out), " we would need :", nvert_in*nbdots
137       STOP "polygones_extend"
138    ENDIF
139    nvert_out = nvert_in*nbdots
140    !
141    j = nvert_in
142    ipos = 0
143    DO i=1,nvert_in
144       xs = poly_in(j,1)
145       xe = poly_in(i,1)
146       ys = poly_in(j,2)
147       ye = poly_in(i,2)
148       IF (xs == xe ) THEN
149          DO id=1,nbdots 
150             ipos = ipos + 1
151             poly_out(ipos,1) = xs
152             poly_out(ipos,2) = ys + (id-1)*(ye-ys)/nbdots
153          ENDDO
154       ELSE IF (ys == ye ) THEN
155          DO id=1,nbdots
156             ipos = ipos + 1
157             poly_out(ipos,1) = xs + (id-1)*(xe-xs)/nbdots
158             poly_out(ipos,2) = ys
159          ENDDO
160       ELSE
161          DO id=1,nbdots
162             ipos = ipos + 1
163             poly_out(ipos,2) = ys + (id-1)*(ye-ys)/nbdots
164             poly_out(ipos,1) = (xs-xe)*(poly_out(ipos,2)-ye)/(ys-ye)+xe
165          ENDDO
166       ENDIF
167       j = i
168    ENDDO
169    !
170  END SUBROUTINE polygones_extend
171!
172!=======================================================================================================
173!
174  SUBROUTINE polygones_intersection(nvert_a, poly_a, nvert_b, poly_b, nbpts, poly_out)
175    !
176    ! Find the polygon resulting from the intersection of poly_a and poly_b
177    !
178    ! Arguments
179    INTEGER(i_std), INTENT(in)               :: nvert_a
180    REAL(r_std), DIMENSION(:,:), INTENT(in)  :: poly_a
181    INTEGER(i_std), INTENT(in)               :: nvert_b
182    REAL(r_std), DIMENSION(:,:), INTENT(in)  :: poly_b
183    INTEGER(i_std), INTENT(out)              :: nbpts
184    REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out
185    !
186    ! Local
187    INTEGER(i_std)                            :: nvert_tmpa, nvert_tmpb, nvert_out, i, j
188    LOGICAL                                   :: inside
189    INTEGER(i_std)                            :: in_a, in_b
190    INTEGER(i_std), ALLOCATABLE, DIMENSION(:) :: ptin_a, ptin_b
191    !
192    nvert_tmpa=SIZE(poly_a,DIM=1)
193    IF ( nvert_tmpa < nvert_a) THEN
194       WRITE(*,*) "Input polygone (poly_a) is smaller than the size proposed in the interface."
195       STOP "polygones_intersection"
196    ENDIF
197    IF (SIZE(poly_a,DIM=2) .NE. 2) THEN
198       WRITE(*,*) "Poly_a cannot be a polygone :", SIZE(poly_a)
199       STOP "polygones_intersection"
200    ENDIF
201    ALLOCATE(ptin_a(nvert_a))
202    nvert_tmpb=SIZE(poly_b,DIM=1)
203    IF ( nvert_tmpb < nvert_a) THEN
204       WRITE(*,*) "Input polygone (poly_b) is smaller than the size proposed in the interface."
205       STOP "polygones_intersection"
206    ENDIF
207    IF (SIZE(poly_b,DIM=2) .NE. 2) THEN
208       WRITE(*,*) "Poly_b cannot be a polygone :", SIZE(poly_b)
209       STOP "polygones_intersection"
210    ENDIF
211    ALLOCATE(ptin_b(nvert_b))
212    !
213    !
214    in_a = 0
215    DO i=1,nvert_a
216       CALL polygones_pointinside(nvert_b, poly_b, poly_a(i,1), poly_a(i,2), inside)
217       IF ( inside ) THEN
218          in_a = in_a + 1
219          ptin_a(in_a) = i
220       ENDIF
221    ENDDO
222    !
223    in_b = 0
224    DO i=1,nvert_b
225       CALL polygones_pointinside(nvert_a, poly_a, poly_b(i,1), poly_b(i,2), inside)
226       IF ( inside ) THEN
227          in_b = in_b + 1
228          ptin_b(in_b) = i
229       ENDIF
230    ENDDO
231    !
232    nbpts = in_a + in_b
233    !
234    nvert_out=SIZE(poly_out,DIM=1)
235    IF (SIZE(poly_out,DIM=2) .NE. 2) THEN
236       WRITE(*,*) "Poly_out cannot be a polygone :", SIZE(poly_out)
237       STOP "polygones_intersection"
238    ENDIF
239    IF ( nbpts <= nvert_out ) THEN
240       i = 0
241       DO j=1,in_a
242          i = i+1
243          poly_out(i,:) = poly_a(ptin_a(j),:)
244       ENDDO
245       DO j=1,in_b
246          i = i+1
247          poly_out(i,:) = poly_b(ptin_b(j),:)
248       ENDDO
249    ELSE
250       WRITE(*,*) "The intersection polygone has", nbpts, "points and that is more than available in poly_out", SIZE(poly_out)
251       STOP "polygones_intersection"
252    ENDIF
253    !
254    DEALLOCATE(ptin_a, ptin_b)
255    !
256  END SUBROUTINE polygones_intersection
257!
258!=======================================================================================================
259!
260  SUBROUTINE polygones_cleanup(nvert_in, poly_in, nvert_out, poly_out)
261    !
262    ! Clean_up the polygone by deleting points which are redundant and ordering based
263    ! on the proximity of the points. Beware this does not give a convex hull to the surface
264    ! inside the polygone.
265    !
266    ! Arguments
267    INTEGER(i_std), INTENT(in)               :: nvert_in
268    REAL(r_std), DIMENSION(:,:), INTENT(in)  :: poly_in
269    INTEGER(i_std), INTENT(out)              :: nvert_out
270    REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out
271    !
272    ! Local
273    !
274    INTEGER(i_std)                           :: nvert, nvert_tmpin, nvert_tmp
275    INTEGER(i_std)                           :: i, j
276    INTEGER(i_std), DIMENSION(1)             :: ismin
277    REAL(r_std), ALLOCATABLE, DIMENSION(:)   :: dist
278    REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: poly_tmp
279    !
280    nvert_tmpin=SIZE(poly_in,DIM=1)
281    IF ( nvert_tmpin < nvert_in ) THEN
282       WRITE(*,*) "Input polygone is smaller than the size proposed in the interface."
283       STOP "polygones_cleanup"
284    ENDIF
285    IF (SIZE(poly_in,DIM=2) .NE. 2) THEN
286       WRITE(*,*) "Poly_in cannot be a polygone :", SIZE(poly_in)
287       STOP "polygones_cleanup"
288    ENDIF
289    nvert_tmp=SIZE(poly_out,DIM=1)
290    IF (SIZE(poly_out,DIM=2) .NE. 2) THEN
291       WRITE(*,*) "Poly_out cannot be a polygone :", SIZE(poly_out)
292       STOP "polygones_cleanup"
293    ENDIF
294    !
295    ALLOCATE(poly_tmp(nvert_in,2))
296    ALLOCATE(dist(nvert_in))
297    !
298    nvert = nvert_in-1
299    poly_tmp(1:nvert,:) = poly_in(2:nvert_in,:)
300    poly_out(1,:) = poly_in(1,:)
301    nvert_out = 1
302    !
303    DO WHILE ( nvert > 0 ) 
304       !
305       DO j=1,nvert
306          dist(j) = SQRT((poly_out(nvert_out,1)-poly_tmp(j,1))**2+(poly_out(nvert_out,2)-poly_tmp(j,2))**2)
307       ENDDO
308       !
309       ismin = MINLOC(dist(1:nvert))
310       !
311       IF ( dist(ismin(1)) > EPSILON(dist(1)) ) THEN
312          !
313          ! The smallest distance between 2 vertices is larger than zero :
314          ! Add the vertex to poly_out and delete it from poly_tmp
315          !
316          nvert_out = nvert_out + 1
317          IF ( nvert_out > nvert_tmp) THEN
318             WRITE(*,*) "Output polygone too small"
319             STOP "polygones_cleanup"
320          ENDIF
321          poly_out(nvert_out,:) = poly_tmp(ismin(1),:)
322          poly_tmp(ismin(1):nvert-1,:) = poly_tmp(ismin(1)+1:nvert,:)
323       ELSE
324          !
325          ! Else the vertex already exists in poly_out and thus we just need
326          ! to delete it from poly_tmp
327          !
328          poly_tmp(ismin(1):nvert-1,:) = poly_tmp(ismin(1)+1:nvert,:)
329       ENDIF
330       !
331       ! One less vertices exists in poly_tmp
332       !
333       nvert = nvert-1
334       !
335    ENDDO
336
337    DEALLOCATE(poly_tmp, dist)
338
339  END SUBROUTINE polygones_cleanup
340!
341!=======================================================================================================
342!
343  SUBROUTINE polygones_area(nvert_in, poly_in, dx, dy, area)
344    !
345    ! Find the polygon resulting from the intersection of poly_a and poly_b
346    !
347    ! Arguments
348    INTEGER(i_std), INTENT(in)               :: nvert_in
349    REAL(r_std), DIMENSION(:,:), INTENT(in)  :: poly_in
350    REAL(r_std), INTENT(in)                  :: dx, dy
351    REAL(r_std), INTENT(out)                 :: area
352    !
353    ! Local
354    !
355    INTEGER(i_std)     :: nvert, i, j
356    !
357    nvert = SIZE(poly_in,DIM=1)
358    IF ( nvert < nvert_in) THEN
359       WRITE(*,*) "Input polygone is smaller than the size proposed in the interface."
360       STOP "polygones_area"
361    ENDIF
362    IF (SIZE(poly_in,DIM=2) .NE. 2) THEN
363       WRITE(*,*) "Poly_in cannot be a polygone :", SIZE(poly_in)
364       STOP "polygones_area"
365    ENDIF
366    !
367    area = zero
368    j = nvert_in
369    DO i=1,nvert_in
370!!       area = area + (dy*(poly_in(j,2)+poly_in(i,2))/2.0)*(dx*(poly_in(j,1)-poly_in(i,1)))
371       area = area + dy*dx/2.0*(poly_in(j,2)+poly_in(i,2))*(poly_in(j,1)-poly_in(i,1))
372       j = i
373    ENDDO
374    !
375    area = ABS(area)
376    !
377  END SUBROUTINE polygones_area
378!
379!=======================================================================================================
380!
381  SUBROUTINE polygones_crossing(nvert_a, poly_a, nvert_b, poly_b, nbpts, crossings)
382    !
383    ! Find the polygon resulting from the intersection of poly_a and poly_b
384    !
385    ! Arguments
386    INTEGER(i_std), INTENT(in)               :: nvert_a
387    REAL(r_std), DIMENSION(:,:), INTENT(in)  :: poly_a
388    INTEGER(i_std), INTENT(in)               :: nvert_b
389    REAL(r_std), DIMENSION(:,:), INTENT(in)  :: poly_b
390    INTEGER(i_std), INTENT(out)              :: nbpts
391    REAL(r_std), DIMENSION(:,:), INTENT(out) :: crossings
392    !
393    ! Local
394    !
395    INTEGER(i_std)                            :: nvert, ia, ja, ib, jb
396    REAL(r_std), DIMENSION(2,2)               :: la, lb
397    REAL(r_std)                               :: x, y
398    LOGICAL                                   :: intersect
399    !
400    nvert=SIZE(poly_a,DIM=1)
401    IF ( nvert < nvert_a) THEN
402       WRITE(*,*) "Input polygone (poly_a) is smaller than the size proposed in the interface."
403       STOP "polygones_crossing"
404    ENDIF
405    IF (SIZE(poly_a,DIM=2) .NE. 2) THEN
406       WRITE(*,*) "Poly_a cannot be a polygone :", SIZE(poly_a)
407       STOP "polygones_crossing"
408    ENDIF
409    nvert=SIZE(poly_b,DIM=1)
410    IF ( nvert < nvert_b) THEN
411       WRITE(*,*) "Input polygone (poly_b) is smaller than the size proposed in the interface."
412       STOP "polygones_crossing"
413    ENDIF
414    IF (SIZE(poly_b,DIM=2) .NE. 2) THEN
415       WRITE(*,*) "Poly_b cannot be a polygone :", SIZE(poly_b)
416       STOP "polygones_crossing"
417    ENDIF
418    nvert=SIZE(crossings,DIM=1)
419    IF (SIZE(crossings,DIM=2) .NE. 2) THEN
420       WRITE(*,*) "Poly_a cannot be a polygone :", SIZE(poly_a)
421       STOP "polygones_crossing"
422    ENDIF
423    !
424    nbpts = 0
425    !
426    ja = nvert_a
427    DO ia=1,nvert_a
428       !
429       ! Put one side of the polygone a into la
430       !
431       la(1,1) = poly_a(ja,1)
432       la(1,2) = poly_a(ja,2)
433       la(2,1) = poly_a(ia,1)
434       la(2,2) = poly_a(ia,2)
435       !
436       jb = nvert_b
437       DO ib=1,nvert_b
438          !
439          ! Put one side of the polygone b into lb
440          !
441          lb(1,1) = poly_b(jb,1)
442          lb(1,2) = poly_b(jb,2)
443          lb(2,1) = poly_b(ib,1)
444          lb(2,2) = poly_b(ib,2)
445          !
446          CALL polygones_lineintersect(la, lb, intersect, x, y)
447          !
448          IF ( intersect ) THEN
449             nbpts = nbpts + 1
450             IF ( nbpts <= nvert ) THEN
451                crossings(nbpts, 1) = x
452                crossings(nbpts, 2) = y
453             ELSE
454                WRITE(*,*) "Polygone to write the intersection points is too small", nvert, nbpts
455                STOP "polygones_crossing"
456             ENDIF
457          ENDIF
458          !
459          jb = ib
460       ENDDO
461       ja = ia
462    ENDDO
463    !
464
465  END SUBROUTINE polygones_crossing
466!
467!=======================================================================================================
468!
469  SUBROUTINE polygones_convexhull(nvert_in, poly_in, nvert_out, poly_out)
470    !
471    ! Routine orders points in poly_in so that they form a convex shape. This is based on the Graham scan
472    ! algorith as implemented by Alan Miller : http://jblevins.org/mirror/amiller/
473    !
474    ! The output polygone can be smaller than the input one as we are computing the convex envelope.
475    !
476    ! Arguments
477    INTEGER(i_std), INTENT(in)               :: nvert_in
478    REAL(r_std), DIMENSION(:,:), INTENT(in)  :: poly_in
479    INTEGER(i_std), INTENT(out)              :: nvert_out
480    REAL(r_std), DIMENSION(:,:), INTENT(out) :: poly_out
481    !
482    ! Local
483    !
484    INTEGER(i_std)                           :: nvert, nvert_tmp
485    INTEGER(i_std), ALLOCATABLE, DIMENSION(:):: iwk, next, vertex
486    REAL(r_std)                              :: xmin, temp, dy, dx2, dy1, dy2
487    REAL(r_std)                              :: x1, x2, y1, y2
488    REAL(r_std)                              :: xmax, ymin, ymax
489    REAL(r_std)                              :: dx, dx1
490    REAL(r_std)                              :: dmax, dmax1, dmax2
491    REAL(r_std)                              :: dist, dmin
492    INTEGER(i_std)                           :: i, i1, i2, j, jp1, jp2, i2save, i3, i2next
493    LOGICAL                                  :: points_todo
494    INTEGER(i_std), PARAMETER                :: nextinc=20
495    !
496    nvert_tmp=SIZE(poly_in,DIM=1)
497    IF ( nvert_tmp < nvert_in) THEN
498       WRITE(*,*) "Input polygone is smaller than the size proposed in the interface."
499       STOP "polygones_convexhull"
500    ENDIF
501    IF ( nvert_tmp <= 2 ) THEN
502       WRITE(*,*) "The input polygone is too small to compute a convex hull."
503       STOP "polygones_convexhull"
504    ENDIF
505    IF (SIZE(poly_in,DIM=2) .NE. 2) THEN
506       WRITE(*,*) "Poly_a cannot be a polygone :", SIZE(poly_in)
507       STOP "polygones_convexhull"
508    ENDIF
509    IF (SIZE(poly_out,DIM=2) .NE. 2) THEN
510       WRITE(*,*) "Poly_b cannot be a polygone :", SIZE(poly_out)
511       STOP "polygones_convexhull"
512    ENDIF
513    !
514    ALLOCATE(vertex(nvert_in))
515    ALLOCATE(iwk(nvert_in))
516    ALLOCATE(next(nvert_in*nextinc))
517    !
518    !
519    !  Choose the points with smallest & largest x- values as the
520    !  first two vertices of the polygon.
521    !
522    IF (poly_in(1,1) > poly_in(nvert_in,1)) THEN
523       vertex(1) = nvert_in
524       vertex(2) = 1
525       xmin = poly_in(nvert_in,1)
526       xmax = poly_in(1,1)
527    ELSE
528       vertex(1) = 1
529       vertex(2) = nvert_in
530       xmin = poly_in(1,1)
531       xmax = poly_in(nvert_in,1)
532    END IF
533   
534    DO i = 2, nvert_in-1
535       temp = poly_in(i,1)
536       IF (temp < xmin) THEN
537          vertex(1) = i
538          xmin = temp
539       ELSE IF (temp > xmax) THEN
540          vertex(2) = i
541          xmax = temp
542       END IF
543    END DO
544    !
545    !       Special case, xmax = xmin.
546    !
547    IF (xmax == xmin) THEN
548       IF (poly_in(1,2) > poly_in(nvert_in,2)) THEN
549          vertex(1) = nvert_in
550          vertex(2) = 1
551          ymin = poly_in(nvert_in,2)
552          ymax = poly_in(1,2)
553       ELSE
554          vertex(1) = 1
555          vertex(2) = nvert_in
556          ymin = poly_in(1,2)
557          ymax = poly_in(nvert_in,2)
558       END IF
559       
560       DO i = 2, nvert_in-1
561          temp = poly_in(i,2)
562          IF (temp < ymin) THEN
563             vertex(1) = i
564             ymin = temp
565          ELSE IF (temp > ymax) THEN
566             vertex(2) = i
567             ymax = temp
568          END IF
569       END DO
570       
571       nvert = 2
572       IF (ymax == ymin) nvert = 1
573       RETURN
574    END IF
575    !
576    !  Set up two initial lists of points; those points above & those below the
577    !  line joining the first two vertices.    next(i) will hold the pointer to the
578    !  point furthest from the line joining vertex(i) to vertex(i+1) on the left
579    !  hand side.
580    !
581    i1 = vertex(1)
582    i2 = vertex(2)
583    iwk(i1) = -1
584    iwk(i2) = -1
585    dx = xmax - xmin
586    y1 = poly_in(i1,2)
587    dy = poly_in(i2,2) - y1
588    dmax = zero
589    dmin = zero
590    next(1) = -1
591    next(2) = -1
592    !
593    DO i = 1, nvert_in
594       IF (i == vertex(1) .OR. i == vertex(2)) CYCLE
595       dist = (poly_in(i,2) - y1)*dx - (poly_in(i,1) - xmin)*dy
596       IF (dist > zero) THEN
597          iwk(i1) = i
598          i1 = i
599          IF (dist > dmax) THEN
600             next(1) = i
601             dmax = dist
602          END IF
603       ELSE IF (dist < zero) THEN
604          iwk(i2) = i
605          i2 = i
606          IF (dist < dmin) THEN
607             next(2) = i
608             dmin = dist
609          END IF
610       END IF
611    END DO
612    !
613    !  Ends of lists are indicated by pointers to -ve positions.
614    !
615    iwk(i1) = -1
616    iwk(i2) = -1
617    nvert = 2
618    !
619    j = 1
620    !
621    !  Start of main process.
622    !
623    !  Introduce new vertex between vertices j & j+1, if one has been found.
624    !  Otherwise increase j.   Exit if no more vertices.
625    !
626    !
627    points_todo = .TRUE.
628    DO WHILE ( points_todo )
629       !
630       DO WHILE ( points_todo .AND. next(j) < 0 )
631          IF (j == nvert) points_todo = .FALSE.
632          j = j + 1
633       ENDDO
634          !
635       IF ( points_todo ) THEN
636          !
637          jp1 = j + 1
638          IF ( jp1 >= nvert_in*nextinc) THEN
639             STOP "polygones_convexhull : please increase nextinc"
640          ENDIF
641          !
642          DO i = nvert, jp1, -1
643             vertex(i+1) = vertex(i)
644             next(i+1) = next(i)
645          END DO
646          jp2 = jp1 + 1
647          nvert = nvert + 1
648          IF (jp2 > nvert) jp2 = 1
649          i1 = vertex(j)
650          i2 = next(j)
651          i3 = vertex(jp2)
652          vertex(jp1) = i2
653          !
654          !  Process the list of points associated with vertex j.   New list at vertex j
655          !  consists of those points to the left of the line joining it to the new
656          !  vertex (j+1).   Similarly for the list at the new vertex.
657          !  Points on or to the right of these lines are dropped.
658          !
659          x1 = poly_in(i1,1)
660          x2 = poly_in(i2,1)
661          y1 = poly_in(i1,2)
662          y2 = poly_in(i2,2)
663          !
664          dx1 = x2 - x1
665          dx2 = poly_in(i3,1) - x2
666          dy1 = y2 - y1
667          dy2 = poly_in(i3,2) - y2
668          dmax1 = zero
669          dmax2 = zero
670          next(j) = -1
671          next(jp1) = -1
672          i2save = i2
673          i2next = iwk(i2)
674          i = iwk(i1)
675          iwk(i1) = -1
676          iwk(i2) = -1
677          !
678          DO WHILE ( i > 0 )
679             IF (i /= i2save) THEN
680                dist = (poly_in(i,2) - y1)*dx1 - (poly_in(i,1) - x1)*dy1
681                IF (dist > zero) THEN
682                   iwk(i1) = i
683                   i1 = i
684                   IF (dist > DMAX1) THEN
685                      next(j) = i
686                      dmax1 = dist
687                   END IF
688                ELSE
689                   dist = (poly_in(i,2) - y2)*dx2 - (poly_in(i,1) - x2)*dy2
690                   IF (dist > zero) THEN
691                      iwk(i2) = i
692                      i2 = i
693                      IF (dist > dmax2) THEN
694                         next(jp1) = i
695                         dmax2 = dist
696                      END IF
697                   END IF
698                END IF
699                i = iwk(i)
700             ELSE
701                i = i2next
702             END IF
703             !
704             !  Get next point from old list at vertex j.
705             !
706          ENDDO
707          !
708          !  End lists with -ve values.
709          !
710          iwk(i1) = -1
711          iwk(i2) = -1
712          !
713       ENDIF
714    ENDDO
715    !
716    ! Copy the polygone in he right order to the output variable
717    !
718    nvert_out = nvert
719    nvert_tmp=SIZE(poly_out,DIM=1)
720    IF ( nvert_tmp < nvert_out) THEN
721       WRITE(*,*) "Ouptput polygone is smaller than the size proposed in the interface."
722       STOP "polygones_convexhull"
723    ENDIF
724    !
725    DO i=1,nvert
726       poly_out(i,:) = poly_in(vertex(i),:)
727    ENDDO
728    !
729    DEALLOCATE(vertex, iwk, next)
730    !
731  END SUBROUTINE polygones_convexhull
732
733END MODULE polygones
Note: See TracBrowser for help on using the repository browser.