From 33f9dd613047340cc779ad9c064ced8424d77d21 Mon Sep 17 00:00:00 2001
From: Daniele Micheletti <daniele.micheletti@istp.cnr.it>
Date: Tue, 9 Feb 2016 11:18:47 +0000
Subject: [PATCH] fixed bug in read_beams2 to correctly locate the position of
 a point outside the beam grid

---
 src/beams.f90 | 95 ++++++++++++++++++++++++++++++---------------------
 1 file changed, 57 insertions(+), 38 deletions(-)

diff --git a/src/beams.f90 b/src/beams.f90
index 6f360e4..ea14446 100644
--- a/src/beams.f90
+++ b/src/beams.f90
@@ -1,5 +1,5 @@
 module beams
-  use const_and_precisions, only : wp_
+  use const_and_precisions, only : wp_, one
   implicit none
 
 contains
@@ -181,7 +181,7 @@ contains
        txx0, tyx0, txy0, tyy0, txz0, tyz0, txycoord, cycoord, cwaist1,  &
        cwaist2, crci1, crci2, cphi1,cphi2, cx0, cy0, cz0, w, wrk2,      &
        xpolyg, ypolyg, xpolygA, ypolygA, xpolygB, ypolygB, xpolygC,     &
-       ypolygC, xpolygD, ypolygD
+       ypolygC, xpolygD, ypolygD, xoutA, youtA, xoutB, youtB, xoutC, youtC
     real(wp_), DIMENSION(4) :: xvert, yvert
     real(wp_), dimension(1) :: fi
     integer, parameter :: kspl=1
@@ -499,7 +499,9 @@ contains
         if(incheck.eq.0) then
           allocate(xpolygA(nxcoord), ypolygA(nxcoord), xpolygC(nxcoord),   &
           ypolygC(nxcoord), xpolygB(nycoord), ypolygB(nycoord),            &
-          xpolygD(nycoord), ypolygD(nycoord))
+          xpolygD(nycoord), ypolygD(nycoord),                              &
+          xoutA(nxcoord+3), youtA(nxcoord+3), xoutB(nycoord+3),            &
+          youtB(nycoord+3), xoutC(nxcoord+3), youtC(nxcoord+3))
 !     coordinates of vertices v1,v2,v3,v4
           xvert(1) = xpolyg(1)
           xvert(2) = xpolyg(nxcoord)
@@ -520,42 +522,59 @@ contains
           xpolygD(nycoord) = xpolyg(1)
           ypolygD(1:nycoord-1) = ypolyg(2*nxcoord+nycoord-2:npolyg)
           ypolygD(nycoord) = ypolyg(1)
-!    c----------------------------------------------------------------------------------
-!     search for position of xcoord0, ycoord0 with respect to (alpha,beta) data grid
+!     contour of outside regions A (1,2), B(3,4), C(5,6)
+          xoutA = huge(one)
+          xoutA(1:nxcoord) = xpolygA
+          xoutA(nxcoord+3) = xvert(1)
+          youtA = -huge(one)
+          youtA(1:nxcoord) = ypolygA
+          youtA(nxcoord+1) = yvert(2)
+          xoutB = huge(one)
+          xoutB(1:nycoord) = xpolygB
+          xoutB(nycoord+1) = xvert(3)
+          youtB = huge(one)
+          youtB(1:nycoord) = ypolygB
+          youtB(nycoord+3) = yvert(2)
+          xoutC = -huge(one)
+          xoutC(1:nxcoord) = xpolygC
+          xoutC(nxcoord+3) = xvert(3)
+          youtC = huge(one)
+          youtC(1:nxcoord) = ypolygC
+          youtC(nxcoord+1) = yvert(4)
+!     c----------------------------------------------------------------------------------
+!      search for position of xcoord0, ycoord0 with respect to (alpha,beta) data grid
+!      
+!      (6)   |        (5)       |        (4)    
+!      _ _ _v4__________________v3 _ _ _ _ _  
+!             \        C        \
+!              \                 \              (1)->(8) outside regions
+!      (7)    D \                 \ B    (3)    v1->v4 grid vertices
+!                \                 \            A-D grid sides   
+!      _ _ _ _ _ _\_________________\_ _ _ _
+!                v1         A       v2         
+!      (8)        |      (1)        |    (2)     
 !     
-!            |                 |
-!     (6)            (5)            (4)    
-!            |                 |           
-!    _ _ _ v4 _________________v3_ _ _ _   
-!            |        C        |           (1)->(8) outside regions
-!            |                 |           
-!     (7)  D |                 | B  (3)    v1->v4 grid vertices
-!            |                 |           
-!    _ _ _ _ |_________________| _ _ _ _   A-D grid sides
-!          v1         A        v2         
-!            |                 |          
-!     (8)            (1)            (2)    
-!            |                 |           
-!     
-          if(xcoord0.gt.xvert(1).and.xcoord0.lt.xvert(2).and.ycoord0.le.maxval(ypolygA)) then
-            in=1
-          else if(ycoord0.gt.yvert(2).and.ycoord0.lt.yvert(3).and.xcoord0.ge.minval(xpolygB)) then
-            in=3
-          else if(xcoord0.lt.xvert(3).and.xcoord0.gt.xvert(4).and.ycoord0.ge.minval(ypolygC)) then
-            in=5
-          else if(ycoord0.lt.yvert(4).and.ycoord0.gt.yvert(1).and.xcoord0.le.maxval(xpolygD)) then
-            in=7
-          else if(xcoord0.ge.xvert(2).and.ycoord0.le.yvert(2)) then
-            in=2
-          else if(xcoord0.ge.xvert(3).and.ycoord0.ge.yvert(3)) then
-            in=4
-          else if(xcoord0.le.xvert(4).and.ycoord0.ge.yvert(4)) then
-            in=6
-          else if(xcoord0.le.xvert(1).and.ycoord0.le.yvert(1)) then
-            in=8
-          endif
-!    c----------------------------------------------------------------------------------
-!
+          if(inside(xoutA,youtA,nxcoord+3,xcoord0,ycoord0)) then
+            in = 1
+            if(xcoord0.GT.xvert(2)) then
+              in = 2
+            end if
+          else if(inside(xoutB,youtB,nycoord+3,xcoord0,ycoord0)) then
+            in = 3
+            if(ycoord0.GT.yvert(3)) then
+              in = 4
+            end if
+          else if(inside(xoutC,youtC,nxcoord+3,xcoord0,ycoord0)) then
+            in = 5
+            if(xcoord0.LT.xvert(4)) then
+              in = 6
+            end if
+          else
+            in = 7
+            if(ycoord0.LT.yvert(1)) then
+              in = 8
+            end if
+          end if
 !    c----------------------------------------------------------------------------------
 !     (xcoord0,ycoord0) is set to its nearest point on (alpha, beta) grid border
 !     depending on the region