added option for ray integration on planes perp to central k

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