From e341062dcba6cdf8818917c86b53bbaeebb79830 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Tue, 30 Apr 2013 16:35:51 +0000 Subject: [PATCH] added option for ray integration on planes perp to central k --- src/gray.f | 45 ++++++++++++++++++++++++++++++++++++++++++--- 1 file changed, 42 insertions(+), 3 deletions(-) diff --git a/src/gray.f b/src/gray.f index fa226e7..ef8d60f 100644 --- a/src/gray.f +++ b/src/gray.f @@ -553,6 +553,7 @@ c common/warm/iwarm,ilarm common/ieccd/ieccd common/idst/idst + common/iplane/iplane c common/filesn/filenmeqq,filenmprf,filenmbm c @@ -673,8 +674,9 @@ c from center of mirror and with angular spread c ipass=1/2 1 or 2 passes into plasma c iox=1/2 OM/XM c idst=0/1/2 0 integration in s, 1 integr. in ct, 2 integr. in Sr +c iplane=0/1 1 plane integration wavefronts c - read(2,*) dst,nstep,istpr0,istpl0,idst + read(2,*) dst,nstep,istpr0,istpl0,idst,iplane read(2,*) igrad,ipass,rwallm read(2,*) iox, psipol0,chipol0 c @@ -2859,6 +2861,7 @@ 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 @@ -2871,11 +2874,15 @@ 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 h6=h/6.0d0 c + ttest=0 do j=1,nrayr kkk=nrayth if(j.eq.1) kkk=1 @@ -2893,16 +2900,34 @@ 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) 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) 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) c do ieq=1,ndim @@ -2977,6 +3002,7 @@ c 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) c common/gr/gr2 common/dgr/dgr2,dgr,ddgr @@ -2995,6 +3021,10 @@ c common/anv/anv common/xv/xv common/idst/idst +c + common/iplane/iplane + common/anv11/anv11 + common/iint/iint c xx=y(1) yy=y(2) @@ -3120,6 +3150,11 @@ c integration variable: c*t 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)) + end if c c coefficient for integration in s c ds/dst, where st is the integration variable @@ -5788,6 +5823,10 @@ c c if(istep.eq.0) . write(10,111) istep,j,k,xti,yti,zti,dirxt,diryt,dirzt,dir +c dr=sqrt(dx**2+dy**2+dz**2) +c write(11,111) istep,j,k,dx/dr,dy/dr,dz/dr, +c . snth1*snps1,snth1*csps1,csth1, +c . snth1*snps1*dx/dr+snth1*csps1*dy/dr+csth1*dz/dr if((iproj.eq.1).or.(iproj.eq.0.and.j.eq.nrayr)) . write(nfile,111) istep,j,k,xti,yti,zti,rti, @@ -5920,9 +5959,9 @@ c Spline fit ymin=-xmaxgrid ymax=xmaxgrid nxgrid=2*nrayr-1 - dx=(xmax-xmin)/(nxgrid-1) + dxgrid=(xmax-xmin)/(nxgrid-1) do i=1,nxgrid - xgridv(i)=xmin+(i-1)*dx + xgridv(i)=xmin+(i-1)*dxgrid ygridv(i)=xgridv(i) end do