local modifications in gaussfit copied in oxb version

src/gray.f

@@ -578,9 +578,11 @@ c
      .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt'
       write(8,*) ' #istep j k xt yt zt rt psin'
       write(9,*) ' #istep j k xt yt zt rt psin'
+      write(82,*) ' #istep j k xt yt zwspl zwparab'
       write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1'
       write(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
-      write(12,*) ' #i sst psi w1 w2'
+      write(12,*) ' #i sst w1 w2 phiw rci1 rci2 phir errw errr '//
+     .            'dk1 dk2 dkpar phik dnpar' 
       write(7,*)'#Icd Pa Jphimx dPdVmx '//
      .'rhotj rhotjava rhotp rhotpav '//
      .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt'
@@ -624,6 +626,7 @@ c
       common/warm/iwarm,ilarm
       common/ieccd/ieccd
       common/idst/idst
+      common/iplane/iplane
 c
       common/filesn/filenmeqq,filenmprf,filenmbm
 c
@@ -665,6 +668,8 @@ c
       common/mode/sox
       common/angles/alpha0,beta0
       common/scal/iscal
+      
+      common/waist/w0csi,w0eta
 c
       open(2,file='gray.data',status= 'unknown')
 c
@@ -742,8 +747,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
@@ -2953,6 +2959,7 @@ c
       hh=h*0.5d0
       h6=h/6.0d0
 c
+      ttest=0
       do j=1,nrayr
         kkk=nrayth
         if(j.eq.1) kkk=1
@@ -2970,17 +2977,17 @@ c
               ddgr(iv,jv)=ggri(iv,jv,j,k)
             end do
           end do
-          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
-          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
-          call fwork(yy,fk4)
+          call fwork(j,k,3,yy,fk4)
 c
           do ieq=1,ndim
             ywrk(ieq,j,k)=y(ieq)
@@ -3034,7 +3041,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)
@@ -3047,14 +3054,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 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
@@ -3072,6 +3083,10 @@ c
       common/anv/anv
       common/xv/xv
       common/idst/idst
+c
+      common/iplane/iplane
+c      
+      common/dery0/dery0,dery0mod
 c
       xx=y(1)
       yy=y(2)
@@ -3179,6 +3194,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
@@ -3195,7 +3218,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))
+        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
@@ -3208,6 +3240,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
@@ -5742,34 +5782,177 @@ c     gg=F(u)/u with F(u) as in Cohen paper
 
 
       subroutine projxyzt(iproj,nfile)
-      implicit real*8 (a-h,o-z)
+      implicit real*8 (a-h,o-z) 
-      parameter(jmx=31,kmx=36)
+      parameter(jmx=31,kmx=36,nmx=8000)
+      parameter(pi=3.14159265358979d0)
       dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx)
+      complex*16 ui,aac,bbc,ccc,ddc,aak,bbk,cck
+      parameter(ui=(0.0d0,1.0d0))
+      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 avn(3,jmx,kmx)
+      dimension aplane(3,jmx,kmx),asip(jmx,kmx)
+      dimension gri(3,jmx,kmx)
+      
+c      parameter(nxmax=2*jmx-1)
+      parameter(nxmax=2*kmx)
+      parameter(kspl=3,nxest=nxmax+4)
+      parameter(km=kspl+1,nu=nxest-km)
+      parameter(ne=nxest,kb1=kspl*nu+km,kb2=kb1+nu-kspl)
+      parameter(lwrk1=nu*nu*(2+kb1+kb2)+
+     .                2*(2*nu+km*(nrmax+ne)+ne-2*kspl)+kb2+1)
+      parameter(lwrk2=nu*nu*(kb2+1)+kb2)
+      parameter(kwrk=nrmax+(nxest-kspl-km)*(nxest-kspl-km))
+      dimension xgridv(nxmax),ygridv(nxmax)
+      dimension zwint(nxmax*nxmax),ccexp(nxmax*nxmax)
+      dimension tx(nxest),ty(nxest),wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk)
+      parameter(lwrkbsp=8*nxmax,liwrkbsp=2*nxmax)
+      dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp)
 c
       common/nray/nrayr,nrayth
       common/wrk/ywrk,ypwrk
-      common/psinv11/psinv11
+      common/derip11/dpsi11dr,dpsi11dz,psinv11,bvx11,bvy11,bvz11
       common/istep/istep
+      common/rwmax/rwmax
       common/ss/st
+      common/dnpar/dnpara
+      common/atjki/tauv,alphav
+      common/waist/w0csi,w0eta
+      
+      common/dery0/dery0,dery0mod
+      common/iplane/iplane
+      common/gradjk/gri
 c
-      rtimn=1.d+30
+      x4m=0.0d0
-      rtimx=-1.d-30
+      x3ym=0.0d0
+      x2y2m=0.0d0
+      xy3m=0.0d0
+      y4m=0.0d0
+      x2zrm=0.0d0
+      xyzrm=0.0d0
+      y2zrm=0.0d0
+      x2zwm=0.0d0
+      xyzwm=0.0d0
+      y2zwm=0.0d0
+      z2wm=0.0d0
+      z2rm=0.0d0
+c initialize grid dimension for spline interpolation
+      xmaxgrid=2*max(w0csi,w0eta)
+      iray=0
+      if(iplane.le.1) then
+        do j=1,nrayr     
+          kktx=nrayth
+          if(j.eq.1) kktx=1
+          do k=1,kktx
+            do ll=1,3
+              aplane(ll,j,k)=ywrk(ll,j,k)
+            enddo
+            asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2
+          enddo
+        enddo  
+      endif  
+      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=sqrt(ypwrk(1,j,k)**2+ypwrk(2,j,k)**2+ypwrk(3,j,k)**2)
+            do ii=1,3
+              avn(ii,j,k)=ypwrk(ii,j,k)/avmod
+            enddo 
+            deltaxdotn0=(ywrk(1,1,1)-ywrk(1,j,k))*ywrk(4,1,1)+
+     .                  (ywrk(2,1,1)-ywrk(2,j,k))*ywrk(5,1,1)+
+     .                  (ywrk(3,1,1)-ywrk(3,j,k))*ywrk(6,1,1)
+            avndotn0=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=deltaxdotn0/avndotn0
+            do ll=1,3
+              aplane(ll,j,k)=ywrk(ll,j,k)+aalpha*avn(ll,j,k)
+            enddo  
+          enddo    
+        enddo
+        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
+            asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2
+          enddo
+        enddo  
+      endif
+      if(iplane.eq.3) then  
+c     ortogonal projection on the plane perpendicular to n0 passing through x11
+        do j=2,nrayr
+          do k=1,nrayth
+            deltaxdotn0=(ywrk(1,1,1)-ywrk(1,j,k))*ywrk(4,1,1)+
+     .                  (ywrk(2,1,1)-ywrk(2,j,k))*ywrk(5,1,1)+
+     .                  (ywrk(3,1,1)-ywrk(3,j,k))*ywrk(6,1,1)
+            an02=ywrk(4,1,1)**2+ywrk(5,1,1)**2+ywrk(6,1,1)**2     
+            aalpha=deltaxdotn0/an02
+            do ll=1,3
+              aplane(ll,j,k)=ywrk(ll,j,k)+aalpha*ywrk(ll+3,1,1)
+            enddo  
+          enddo    
+        enddo
+        do ll=1,3
+          aplane(ll,1,1)=ywrk(ll,1,1)
+        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)*(aplane(1,j,k)-ywrk(1,j,k))+
+     .                gri(2,j,k)*(aplane(2,j,k)-ywrk(2,j,k))+
+     .                gri(3,j,k)*(aplane(3,j,k)-ywrk(3,j,k))
+          enddo
+        enddo  
+      endif 
+c      
+      do j=1,nrayr     
+        kktx=nrayth
+        if(j.eq.1) kktx=1
+        do k=1,kktx
+          zwj=asip(j,k)+
+     .    0.5*(tauv(j,k,istep)-tauv(1,1,istep))
 c
-      jd=1
+          dx=aplane(1,j,k)-aplane(1,1,1)
-      if(iproj.eq.0) jd=nrayr-1
+          dy=aplane(2,j,k)-aplane(2,1,1)
-      do j=1,nrayr,jd
+          dz=aplane(3,j,k)-aplane(3,1,1)
-        kkk=nrayth
-        if(j.eq.1) kkk=1
-        do k=1,kkk
-c
-          dx=ywrk(1,j,k)-ywrk(1,1,1)
-          dy=ywrk(2,j,k)-ywrk(2,1,1)
-          dz=ywrk(3,j,k)-ywrk(3,1,1)
 c
           dirx=ywrk(4,j,k)
           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=aplane(1,j,k2)-aplane(1,1,1)
+            dy2=aplane(2,j,k2)-aplane(2,1,1)
+            dz2=aplane(3,j,k2)-aplane(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)**2+pvett(2)**2+pvett(3)**2)
+            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
@@ -5784,39 +5967,235 @@ c
           xti= dx*csps1-dy*snps1
           yti=(dx*snps1+dy*csps1)*csth1-dz*snth1
           zti=(dx*snps1+dy*csps1)*snth1+dz*csth1
-          rti=sqrt(xti**2+yti**2)
+          rti=sqrt(xti**2+yti**2)     
+c store x,y,z values for spline interpolation
+          iray=iray+1
+          xtiv(iray)=xti
+          ytiv(iray)=yti 
+          zwjv(iray)=zwj
+c initialize grid dimension for spline interpolation
+          xmaxgrid=max(xmaxgrid,rti)
 c
           dirxt= (dirx*csps1-diry*snps1)/dir
           diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir
           dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir
+c
+          bxti= bvx11*csps1-bvy11*snps1
+          byti=(bvx11*snps1+bvy11*csps1)*csth1-bvz11*snth1
+c          bzti=(bvx11*snps1+bvy11*csps1)*snth1+bvz11*csth1
+c
+          rr11=sqrt(ywrk(1,1,1)**2+ywrk(2,1,1)**2)
+          dpsx=dpsi11dr*ywrk(1,1,1)/rr11
+          dpsy=dpsi11dr*ywrk(2,1,1)/rr11
+          dpsz=dpsi11dz
+          xdpsi=dpsx*csps1-dpsy*snps1
+          ydpsi=(dpsx*snps1+dpsy*csps1)*csth1-dpsz*snth1
+          zdpsi=(dpsx*snps1+dpsy*csps1)*snth1+dpsz*csth1
+          sq=sqrt(xdpsi**2+ydpsi**2+zdpsi**2)
+          if(sq.gt.0.0d0) then
+            xdpsi=dx*xdpsi/sq
+            ydpsi=dx*ydpsi/sq
+            zdpsi=dx*zdpsi/sq
+          end if  
 c
           if(k.eq.1) then
             xti1=xti
             yti1=yti
             zti1=zti
             rti1=rti
-          end if
+            xdpsi1=xdpsi
+            ydpsi1=ydpsi
+            zdpsi1=zdpsi
+          end if  
+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(.not.(iproj.eq.0.and.j.eq.1))
+        if((iproj.eq.1).or.(iproj.eq.0.and.j.eq.nrayr))
-     .    write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
+     .    write(nfile,111) istep,j,k,xti,yti,zti,rti,
-c
+     .    xdpsi,ydpsi,zdpsi
-          if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti
-          if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti
 c
+          x4m=x4m+xti**4
+          x3ym=x3ym+xti**3*yti
+          x2y2m=x2y2m+xti**2*yti**2
+          xy3m=xy3m+xti*yti**3
+          y4m=y4m+yti**4
+          x2zrm=x2zrm+xti**2*zti
+          xyzrm=xyzrm+xti*yti*zti
+          y2zrm=y2zrm+yti**2*zti
+          x2zwm=x2zwm+xti**2*zwj
+          xyzwm=xyzwm+xti*yti*zwj
+          y2zwm=y2zwm+yti**2*zwj
+          z2wm=z2wm+zwj*zwj
+          z2rm=z2rm+zti*zti
         end do
-c
+c        
-c        if(.not.(iproj.eq.0.and.j.eq.1))
+        if((iproj.eq.1.and.j.gt.1).or.(iproj.eq.0.and.j.eq.nrayr))
-c     .    write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11
+     .   write(nfile,111) istep,j,k,xti1,yti1,zti1,rti1,
+     .  xdpsi1,ydpsi1,zdpsi1
         if(iproj.eq.1) write(nfile,*) ' '
       end do
-c
+c      
       write(nfile,*) ' '
 c
-      write(12,99) istep,st,psinv11,rtimn,rtimx
+c     computation of the SI paraboloid
+c
+      denw= x2y2m**3 + x4m*xy3m**2 + x3ym**2*y4m - 
+     .  x2y2m*(2*x3ym*xy3m + x4m*y4m)
+      aaw=  -(x2y2m*xy3m*xyzwm) + x2y2m**2*y2zwm - 
+     .  x3ym*xy3m*y2zwm + x3ym*xyzwm*y4m + 
+     .  x2zwm*(xy3m**2 - x2y2m*y4m)
+      ccw=  x2y2m**2*xyzwm + x4m*xy3m*y2zwm - 
+     .  x2y2m*(x2zwm*xy3m + x3ym*y2zwm) + x2zwm*x3ym*y4m - 
+     .  x4m*xyzwm*y4m
+      bbw=   x2y2m**2*x2zwm - x2zwm*x3ym*xy3m + x4m*xy3m*xyzwm + 
+     .  x3ym**2*y2zwm - x2y2m*(x3ym*xyzwm + x4m*y2zwm)
+      aaw=aaw/denw
+      bbw=bbw/denw
+      ccw=ccw/denw
+      phiw = 0.5d0*atan(ccw/(aaw-bbw+1.0d-32))
+      phiwdeg=phiw*180.d0/pi
+      aaaw = 0.5d0*(aaw+bbw + (aaw-bbw)/cos(2.0d0*phiw))
+      bbbw = aaw+bbw - aaaw
+      errw= 2.0d0*aaw*bbw*x2y2m + ccw**2*x2y2m - 2.0d0*aaw*x2zwm + 
+     -  2.0d0*aaw*ccw*x3ym + aaw**2*x4m + 2.0d0*bbw*ccw*xy3m - 
+     -  2.0d0*ccw*xyzwm - 2.0d0*bbw*y2zwm + bbw**2*y4m + z2wm
+      errw=sqrt(abs(errw)/dble((nrayr-1)*nrayth))
+      wcsi = 1.0d0/sqrt(aaaw)
+      weta = 1.0d0/sqrt(bbbw)
+c
+c     computation of paraboloid Sr=z-(ax^2+cxy+by^2)=const
+c
+      denr= x2y2m**3 + x4m*xy3m**2 + x3ym**2*y4m - 
+     .  x2y2m*(2*x3ym*xy3m + x4m*y4m)
+      aar=  -(x2y2m*xy3m*xyzrm) + x2y2m**2*y2zrm - 
+     .  x3ym*xy3m*y2zrm + x3ym*xyzrm*y4m + 
+     .  x2zrm*(xy3m**2 - x2y2m*y4m)
+      ccr=  x2y2m**2*xyzrm + x4m*xy3m*y2zrm - 
+     .  x2y2m*(x2zrm*xy3m + x3ym*y2zrm) + x2zrm*x3ym*y4m - 
+     .  x4m*xyzrm*y4m
+      bbr=   x2y2m**2*x2zrm - x2zrm*x3ym*xy3m + x4m*xy3m*xyzrm + 
+     .  x3ym**2*y2zrm - x2y2m*(x3ym*xyzrm + x4m*y2zrm)
+      aar=aar/denr
+      bbr=bbr/denr
+      ccr=ccr/denr
+      phir = 0.5d0*atan(ccr/(aar-bbr+1.0d-32))
+      phirdeg=phir*180.d0/pi
+      aaar = 0.5d0*(aar+bbr + (aar-bbr)/cos(2.0d0*phir))
+      bbbr = aar+bbr - aaar
+      errr= 2.0d0*aar*bbr*x2y2m + ccr**2*x2y2m - 2.0d0*aar*x2zrm + 
+     -  2.0d0*aar*ccr*x3ym + aar**2*x4m + 2.0d0*bbr*ccr*xy3m - 
+     -  2.0d0*ccr*xyzrm - 2.0d0*bbr*y2zrm + bbr**2*y4m + z2rm
+      errr=sqrt(abs(errr)/dble((nrayr-1)*nrayth))
+      rcicsi=-2.0d0*aaar
+      rcieta=-2.0d0*bbbr
+c
+c     computation of Fourier Transform of exp[i k0 (Sr + i Si)]
+c     k spectrum 
+c
+      aar=-ak0*aar
+      bbr=-ak0*bbr
+      ccr=-ak0*ccr
+      aac=aar+ui*aaw
+      bbc=bbr+ui*bbw
+      ccc=ccr+ui*ccw
+      ddc=ccc*ccc-4.0d0*aac*bbc
+      aak=bbc/ddc
+      bbk=aac/ddc
+      cck=-ccc/ddc
+      akw=dimag(aak)
+      bkw=dimag(bbk)
+      ckw=dimag(cck)
+      dkpar2=4.0d0*(bxti**2*bkw+byti**2*akw-bxti*byti*ckw)/
+     .        (4.0d0*akw*bkw-ckw*ckw)
+      phik = 0.5d0*atan(ckw/(akw-bkw+1.0d-32))
+      aakw = 0.5d0*(akw+bkw + (akw-bkw)/cos(2.0d0*phik))
+      bbkw = akw+bkw - aakw
+      dk1 = 1.0d0/sqrt(aakw)
+      dk2 = 1.0d0/sqrt(bbkw)
+c      co = cos(phik)
+c      si = sin(phik)
+c      bxtir= co*bxti+si*byti
+c      bytir=-si*bxti+co*byti
+c      dkpar2=bxtir**2/aakw+bytir**2/bbkw
+      dkpar=sqrt(dkpar2)
+      dnpar=dkpar/ak0
+      dnpara=0.5d0*dnpar
+c
+c     spectral distribution of electric field E(k): 
+c     exp[-(N//-N//0)^2/(DNPAR^2)]
+c     dnpar^2 = 2 (Delta N//)^2 ,  Maj !!!
+c
+c     in vacuum dkpar^2 = 4/w0^2
+c
+c     in common dnpara=dnpar/2 to match  westerhof Delta function
+c     Delta -> exp[-(...)^2/(2*DeltaQ)]
+c
+
+
+c     Spline fit
+      if(nrayr.gt.1) then
+
+      npoints=iray
+      xmin=-xmaxgrid
+      xmax=xmaxgrid
+      ymin=-xmaxgrid
+      ymax=xmaxgrid
+      nxgrid=2*nrayr-1
+      dxgrid=(xmax-xmin)/(nxgrid-1)
+      do i=1,nxgrid
+        xgridv(i)=xmin+(i-1)*dxgrid
+        ygridv(i)=xgridv(i)
+      end do
+
+c     call interp spline
+      iopt=0
+      sspl=1.0d-3
+      eps=1.0d-7
+      do iray=1,npoints
+        w(iray) = 1.0d0
+      end do
+      call surfit(iopt,npoints,xtiv,ytiv,zwjv,w,xmin,xmax,ymin,ymax,
+     .            kspl,kspl,sspl,nxest,nxest,nxest,eps,
+     .            nkntx,tx,nknty,ty,ccexp,resid,wrk1,lwrk1,wrk2,lwrk2,
+     .            iwrk,kwrk,ierr)
+      if (ierr.gt.0) then
+        print*, 'surfit:',istep,nkntx,nknty,ierr,resid
+        do j=1,nxgrid
+          do i=1,nxgrid
+            zwint(nxgrid*(i-1)+j)=0.0d0
+          end do
+        end do
+      else
+c     call eval spline
+        call bispev(tx,nkntx,ty,nknty,ccexp,kspl,kspl,xgridv,nxgrid,
+     .     ygridv,nxgrid,zwint,wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ierr)
+        if (ierr.ne.0) print*, 'bispev:',istep,ierr
+      end if
+      do j=1,nxgrid
+        do i=1,nxgrid
+          ij=nxgrid*(i-1)+j
+          write(82,111) istep,i,j,xgridv(i),ygridv(j),zwint(ij),
+     .      aaw*xgridv(i)**2+ccw*xgridv(i)*ygridv(j)+bbw*ygridv(j)**2
+        end do
+        write(82,*) ''
+      end do
+      write(82,*) ''
+      end if
+
+      write(12,99) istep,st,
+     . wcsi,weta,phiwdeg,rcicsi,rcieta,phirdeg,errw,errr,
+     . dk1,dk2,dkpar,phik*180.0d0/pi,dnpar
+c      
       return
- 99   format(i5,12(1x,e16.8e3))
+ 99   format(i5,22(1x,e16.8e3))
 111   format(3i5,12(1x,e16.8e3))
       end
+
 c
 c
 c
@@ -6561,7 +6940,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