istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

corrected detection of turning point for rays missing inner wall

Browse Source
This commit is contained in:
Lorenzo Figini 2014-10-14 14:14:31 +00:00
parent 8f06be1281
commit 8e593fdb1a
1 changed files with 11 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
src/gray.f
View File

@ -178,7 +178,7 @@ c
dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6),y(6),dery(6) dimension tau1v(jmx,kmx),yyrfl(jmx,kmx,6),y(6),dery(6)
dimension ext(jmx,kmx,0:3),eyt(jmx,kmx,0:3) 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 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) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
real*8 rlim(nbb),zlim(nbb) real*8 rlim(nbb),zlim(nbb)
logical inside_plasma logical inside_plasma
@ -195,7 +195,7 @@ c
c c
common/iiv/iiv common/iiv/iiv
common/iov/iop,iow,ihcd,istore common/iov/iop,iow,ihcd,istore
common/refln/anwcl,jclosest common/refln/anwcl,xwcl,jclosest
common/psjki/psjki common/psjki/psjki
common/psival/psinv common/psival/psinv
common/psinv11/psinv11 common/psinv11/psinv11
@ -309,6 +309,7 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
if (j.lt.jclosest) then if (j.lt.jclosest) then
jclosest=j jclosest=j
anwcl=anw anwcl=anw
xwcl=xvrfl
end if end if
xv=xvrfl xv=xvrfl
anv=anvrfl anv=anvrfl
@ -343,7 +344,11 @@ cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
kkk=nrayth kkk=nrayth
if(j.eq.1) kkk=1 if(j.eq.1) kkk=1
do k=1,kkk 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 if(akdotn.lt.aknmin) aknmin=akdotn
end do end do
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 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 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 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) parameter(tmax=5,npts=500)
real*8 ttv(npts+1),extv(npts+1) real*8 ttv(npts+1),extv(npts+1)
@ -2715,13 +2720,14 @@ c
common/nray/nrayr,nrayth common/nray/nrayr,nrayth
common/nstep/nstep common/nstep/nstep
common/tau1v/tau1v common/tau1v/tau1v
common/refln/anwcl,jclosest common/refln/anwcl,xwcl,jclosest
c c
if(nstep.gt.nmx) nstep=nmx if(nstep.gt.nmx) nstep=nmx
if(nrayr.gt.jmx) nrayr=jmx if(nrayr.gt.jmx) nrayr=jmx
if(nrayth.gt.kmx) nrayth=kmx if(nrayth.gt.kmx) nrayth=kmx
jclosest=nrayr+1 jclosest=nrayr+1
anwcl(1:3)=0.0d0 anwcl(1:3)=0.0d0
xwcl(1:3)=0.0d0
c c
do i=1,nstep do i=1,nstep
do k=1,nrayth do k=1,nrayth