projxyzt with gaussian fit

src/gray.f

@@ -507,7 +507,8 @@ c
       write(9,*) ' #istep j k xt yt zt rt psin'
       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'
@@ -5662,24 +5663,43 @@ c     gg=F(u)/u with F(u) as in Cohen paper
 
       subroutine projxyzt(iproj,nfile)
       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)
 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
 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
-      jd=1
+      do j=1,nrayr     
-      if(iproj.eq.0) jd=nrayr-1
+        kktx=nrayth
-      do j=1,nrayr,jd
+        if(j.eq.1) kktx=1
-        kkk=nrayth
+        do k=1,kktx
-        if(j.eq.1) kkk=1
+          zwj=(dble(j-1)*rwmax/dble(nrayr-1))**2+
-        do k=1,kkk
+     .    0.5*(tauv(j,k,istep)-tauv(1,1,istep))
+  
 c
           dx=ywrk(1,j,k)-ywrk(1,1,1)
           dy=ywrk(2,j,k)-ywrk(2,1,1)
@@ -5708,34 +5728,167 @@ 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
+            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
 
-          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        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      
       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
+      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