From 8e593fdb1a7ad3a76a5618c25349db3da957d7af Mon Sep 17 00:00:00 2001
From: Lorenzo Figini <lorenzo.figini@istp.cnr.it>
Date: Tue, 14 Oct 2014 14:14:31 +0000
Subject: [PATCH] corrected detection of turning point for rays missing inner
 wall

---
 src/gray.f | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/src/gray.f b/src/gray.f
index ed56d1c..598c1c5 100644
--- a/src/gray.f
+++ b/src/gray.f
@@ -178,7 +178,7 @@ c
       dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6),y(6),dery(6)
       dimension ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3)
       dimension xv(3),anv(3),xvrfl(3),anvrfl(3),xvvac(3),anw(3),anwcl(3)
-      dimension xvjk(3,jmx,kmx),anvjk(3,jmx,kmx)
+      dimension xwcl(3),xvjk(3,jmx,kmx),anvjk(3,jmx,kmx)
       dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
       real*8 rlim(nbb),zlim(nbb)
       logical inside_plasma
@@ -195,7 +195,7 @@ c
 c
       common/iiv/iiv
       common/iov/iop,iow,ihcd,istore
-      common/refln/anwcl,jclosest
+      common/refln/anwcl,xwcl,jclosest
       common/psjki/psjki
       common/psival/psinv
       common/psinv11/psinv11
@@ -309,6 +309,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
               if (j.lt.jclosest) then
                 jclosest=j
                 anwcl=anw
+                xwcl=xvrfl
               end if
               xv=xvrfl
               anv=anvrfl
@@ -343,7 +344,11 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
           kkk=nrayth
           if(j.eq.1) kkk=1
           do k=1,kkk
-            akdotn=dot_product(anvjk(:,j,k),anwcl)
+            anwclr=(anwcl(1)*xwcl(1)+anwcl(2)*xwcl(2))
+     .            /sqrt(xwcl(1)**2+xwcl(2)**2)
+            anvjkr=(anvjk(1,j,k)*xvjk(1,j,k)+anvjk(2,j,k)*xvjk(2,j,k))
+     .            /sqrt(xvjk(1,j,k)**2+xvjk(2,j,k)**2)
+            akdotn=anwclr*anvjkr+anwcl(3)*anvjk(3,j,k)
             if(akdotn.lt.aknmin) aknmin=akdotn
           end do
         end do
@@ -2698,7 +2703,7 @@ c   spline interpolation of H(lambda,rhop) and dH/dlambda
       dimension pdjki(jmx,kmx,nmx),ppabs(jmx,kmx,nmx),tau1v(jmx,kmx)
       dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx)
       dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx)
-      dimension istore(jmx,kmx),anwcl(3)
+      dimension istore(jmx,kmx),anwcl(3),xwcl(3)
       parameter(tmax=5,npts=500)
       real*8 ttv(npts+1),extv(npts+1)
 
@@ -2715,13 +2720,14 @@ c
       common/nray/nrayr,nrayth
       common/nstep/nstep
       common/tau1v/tau1v
-      common/refln/anwcl,jclosest
+      common/refln/anwcl,xwcl,jclosest
 c
       if(nstep.gt.nmx) nstep=nmx
       if(nrayr.gt.jmx) nrayr=jmx
       if(nrayth.gt.kmx) nrayth=kmx
       jclosest=nrayr+1
       anwcl(1:3)=0.0d0
+      xwcl(1:3)=0.0d0
 c
       do i=1,nstep
         do k=1,nrayth