diff --git a/src/gray.f b/src/gray.f index ef8d60f..bdeae1b 100644 --- a/src/gray.f +++ b/src/gray.f @@ -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