source: branches/publications/ORCHIDEE_IPSLCM5A2.1.r5307/ORCHIDEE/src_global/polygones.f90 @ 7683

Last change on this file since 7683 was 3447, checked in by josefine.ghattas, 8 years ago

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