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

Changed arguments in fwork. Added projection on plane when integrating in Sr (iplane=2,3) - fix needed in projxyzt (ywrk)

Browse Source
This commit is contained in:
Lorenzo Figini 2013-05-07 14:38:59 +00:00
parent e341062dcb
commit 5c641a55ec
1 changed files with 151 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

190
src/gray.f
View File

@ -2861,7 +2861,6 @@ c
dimension ywrk(ndim,jmx,kmx),ypwrk(ndim,jmx,kmx)
dimension grad2(jmx,kmx),dgrad2v(3,jmx,kmx)
dimension gri(3,jmx,kmx),ggri(3,3,jmx,kmx)
dimension anv11(3,3)
c
common/nray/nrayr,nrayth
common/dsds/dst
@ -2874,9 +2873,6 @@ c
common/gr/gr2
common/dgr/dgr2,dgr,ddgr
common/igrad/igrad
c
common/anv11/anv11
common/iint/iint
c
h=dst
hh=h*0.5d0
@ -2900,35 +2896,17 @@ c
ddgr(iv,jv)=ggri(iv,jv,j,k)
end do
end do
iint=1
if(j.eq.1) then
anv11(1,iint)=yy(4)
anv11(2,iint)=yy(5)
anv11(3,iint)=yy(6)
endif
call fwork(yy,fk2)
call fwork(j,k,1,yy,fk2)
c
do ieq=1,ndim
yy(ieq)=y(ieq)+fk2(ieq)*hh
end do
iint=2
if(j.eq.1) then
anv11(1,iint)=yy(4)
anv11(2,iint)=yy(5)
anv11(3,iint)=yy(6)
endif
call fwork(yy,fk3)
call fwork(j,k,2,yy,fk3)
c
do ieq=1,ndim
yy(ieq)=y(ieq)+fk3(ieq)*h
end do
iint=3
if(j.eq.1) then
anv11(1,iint)=yy(4)
anv11(2,iint)=yy(5)
anv11(3,iint)=yy(6)
endif
call fwork(yy,fk4)
call fwork(j,k,3,yy,fk4)
c
do ieq=1,ndim
ywrk(ieq,j,k)=y(ieq)
@ -2982,7 +2960,7 @@ c
end do
end if
c
call fwork(yy,yyp)
call fwork(j,k,3,yy,yyp)
c
do ieq=1,ndim
ypwrk(ieq,j,k)=yyp(ieq)
@ -2995,15 +2973,18 @@ c
c
c
c
subroutine fwork(y,dery)
subroutine fwork(j,k,isubst,y,dery)
implicit real*8 (a-h,o-z)
parameter(ndim=6)
dimension y(ndim),dery(ndim)
dimension xv(3),anv(3),vgv(3),bv(3),derbv(3,3),derxg(3),deryg(3)
dimension derdxv(3),danpldxv(3),derdnv(3)
dimension dgr2(3),dgr(3),ddgr(3,3),dbgr(3)
dimension anv11(3,3)
dimension yy11(6,3),derdxv11(3,3)
dimension dery0(3)
c
save yy11,derdxv11
common/gr/gr2
common/dgr/dgr2,dgr,ddgr
common/ddd/dd,an2s,an2,fdia,bdotgr,ddi,ddr11
@ -3023,8 +3004,8 @@ c
common/idst/idst
c
common/iplane/iplane
common/anv11/anv11
common/iint/iint
c
common/dery0/dery0,dery0mod
c
xx=y(1)
yy=y(2)
@ -3132,6 +3113,14 @@ c
c
derdnm=derdnm+derdnv(iv)**2
end do
if(j.eq.1.and.k.eq.1) then
do ii=1,3
yy11(ii,isubst)=y(ii)
yy11(ii+3,isubst)=y(ii+3)
derdxv11(ii,isubst)=derdxv(ii)
enddo
endif
c
derdnm=sqrt(derdnm)
c
@ -3148,12 +3137,16 @@ c integration variable: c*t
denom=derdom
else
c integration variable: Sr
denom=-(anv(1)*derdnv(1)+anv(2)*derdnv(2)+anv(3)*derdnv(3))
end if
c
if(iplane.eq.1.and.idst.eq.2) then
denom=-(anv11(1,iint)*derdnv(1)+anv11(2,iint)*derdnv(2)
. +anv11(3,iint)*derdnv(3))
if (iplane.eq.1.and.j.gt.1) then
c advance outer rays to the plane through x_11 and perp to N_11
denom=-((yy11(4,isubst)*derdnv(1)+yy11(5,isubst)*derdnv(2)
. +yy11(6,isubst)*derdnv(3))
. -(derdxv11(1,isubst)*(y(1)-yy11(1,isubst))+
. derdxv11(2,isubst)*(y(2)-yy11(2,isubst))+
. derdxv11(3,isubst)*(y(3)-yy11(3,isubst))))
else
denom=-(anv(1)*derdnv(1)+anv(2)*derdnv(2)+anv(3)*derdnv(3))
end if
end if
c
c coefficient for integration in s
@ -3166,6 +3159,14 @@ c
dery(4) = derdxv(1)/denom
dery(5) = derdxv(2)/denom
dery(6) = derdxv(3)/denom
if(j.eq.1) then
do ll=1,3
dery0(ll)=dery(ll)
enddo
dery0mod=sqrt(dery0(1)*dery0(1)+
. dery0(2)*dery0(2)+dery0(3)*dery0(3))
endif
c
c vgv : ~ group velocity
c
@ -5709,6 +5710,14 @@ c gg=F(u)/u with F(u) as in Cohen paper
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
parameter(nrmax=(jmx-1)*kmx+1)
dimension xtiv(nrmax),ytiv(nrmax),zwjv(nrmax),w(nrmax)
dimension pvett(3),pvettn(3),dery0(3),dery0n(3)
dimension deltax(3,jmx,kmx),avn(3,jmx,kmx),avmod(jmx,kmx)
dimension deltaxdotn0(jmx,kmx),avndotn0(jmx,kmx)
dimension aalpha(jmx,kmx),aplane(3,jmx,kmx)
dimension gri(3,jmx,kmx)
dimension aincr(3,jmx,kmx),asip(jmx,kmx)
c parameter(nxmax=2*jmx-1)
parameter(nxmax=2*kmx)
parameter(kspl=3,nxest=nxmax+4)
@ -5733,6 +5742,10 @@ c
common/dnpar/dnpara
common/atjki/tauv,alphav
common/waist/w0csi,w0eta
common/dery0/dery0,dery0mod
common/iplane/iplane
common/gradjk/gri
c
x4m=0.0d0
x3ym=0.0d0
@ -5750,12 +5763,87 @@ c
c initialize grid dimension for spline interpolation
xmaxgrid=2*max(w0csi,w0eta)
iray=0
c
if(iplane.le.1) then
do j=1,nrayr
kktx=nrayth
if(j.eq.1) kktx=1
do k=1,kktx
asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2
enddo
enddo
else
if(iplane.eq.2) then
c prjection parallel to vg on the plane perpendicular to n0 passing through x11
do j=2,nrayr
do k=1,nrayth
avmod(j,k)=sqrt(ypwrk(1,j,k)*ypwrk(1,j,k)+
. ypwrk(2,j,k)*ypwrk(2,j,k)+ypwrk(3,j,k)*ypwrk(3,j,k))
do ii=1,3
deltax(ii,j,k)=ywrk(ii,1,1)-ywrk(ii,j,k)
avn(ii,j,k)=ypwrk(ii,j,k)/avmod(j,k)
enddo
deltaxdotn0(j,k)=deltax(1,j,k)*ywrk(4,1,1)+
. deltax(2,j,k)*ywrk(5,1,1)+
. deltax(3,j,k)*ywrk(6,1,1)
avndotn0(j,k)=avn(1,j,k)*ywrk(4,1,1)+
. avn(2,j,k)*ywrk(5,1,1)+
. avn(3,j,k)*ywrk(6,1,1)
aalpha(j,k)=deltaxdotn0(j,k)/avndotn0(j,k)
do ll=1,3
aplane(ll,j,k)=ywrk(ll,j,k)+aalpha(j,k)*avn(ll,j,k)
enddo
enddo
enddo
c ortogonal projection on the plane perpendicular to n0 passing through x11
else
do j=2,nrayr
do k=1,nrayth
do ii=1,3
deltax(ii,j,k)=ywrk(ii,1,1)-ywrk(ii,j,k)
enddo
deltaxdotn0(j,k)=deltax(1,j,k)*ywrk(4,1,1)+
. deltax(2,j,k)*ywrk(5,1,1)+
. deltax(3,j,k)*ywrk(6,1,1)
an02=ywrk(4,1,1)*ywrk(4,1,1)+ywrk(5,1,1)*ywrk(5,1,1)+
. ywrk(6,1,1)*ywrk(6,1,1)
aalpha(j,k)=deltaxdotn0(j,k)/an02
do ll=1,3
aplane(ll,j,k)=ywrk(ll,j,k)+aalpha(j,k)*ywrk(ll+3,1,1)
enddo
enddo
enddo
end if
do ll=1,3
aplane(ll,1,1)=ywrk(ll,1,1)
enddo
do j=1,nrayr
kktx=nrayth
if(j.eq.1) kktx=1
do k=1,kktx
do ii=1,3
aincr(ii,j,k)=aplane(ii,j,k)-ywrk(ii,j,k)
ywrk(ii,j,k)=aplane(ii,j,k)
enddo
enddo
enddo
c Si evaluation on the projection plane(Taylor, first order)
do j=1,nrayr
kktx=nrayth
if(j.eq.1) kktx=1
do k=1,kktx
asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2+
. gri(1,j,k)*aincr(1,j,k)+
. gri(2,j,k)*aincr(2,j,k)+
. gri(3,j,k)*aincr(3,j,k)
enddo
enddo
endif
c
do j=1,nrayr
kktx=nrayth
if(j.eq.1) kktx=1
do k=1,kktx
zwj=(dble(j-1)*rwmax/dble(nrayr-1))**2+
zwj=asip(j,k)+
. 0.5*(tauv(j,k,istep)-tauv(1,1,istep))
c
dx=ywrk(1,j,k)-ywrk(1,1,1)
@ -5766,6 +5854,30 @@ c
diry=ywrk(5,j,k)
dirz=ywrk(6,j,k)
dir=sqrt(dirx*dirx+diry*diry+dirz*dirz)
if (j>1) then
k2=mod(k+kktx/4-1,kktx)+1
dx2=ywrk(1,j,k2)-ywrk(1,1,1)
dy2=ywrk(2,j,k2)-ywrk(2,1,1)
dz2=ywrk(3,j,k2)-ywrk(3,1,1)
pvett(1)=dy*dz2-dy2*dz
pvett(2)=dz*dx2-dz2*dx
pvett(3)=dx*dy2-dx2*dy
pvettmod=sqrt(pvett(1)*pvett(1)+
. pvett(2)*pvett(2)+pvett(3)*pvett(3))
do ll=1,3
pvettn(ll)=pvett(ll)/pvettmod
dery0n(ll)=dery0(ll)/dery0mod
enddo
c write(*,*) 'dotn0',j,k,(pvettn(1)*dirx+
c . pvettn(2)*diry+pvettn(3)*dirz)/dir
c write(*,*) 'dotvg0',j,k,(pvettn(1)*dery0(1)+
c . pvettn(2)*dery0(2)+pvettn(3)*dery0(3))/dery0mod
endif
c dirx=dery0(1)
c diry=dery0(2)
c dirz=dery0(3)
c dir=dery0mod
c
if(j.eq.1.and.k.eq.1) then
csth1=dirz/dir
@ -6753,7 +6865,7 @@ c ivac=-1 second interface vacuum-plasma WITHOUT wall reflection
y(4)=anv(1)
y(5)=anv(2)
y(6)=anv(3)
call fwork(y,dery)
call fwork(1,1,3,y,dery)
if (st.ge.smax.or.(psinv.gt.0.0d0.and.psinv.lt.psdbnd)) exit
i=i+1
end do