modified projxyzt: added window function and corrected the case iplane=2

src/gray.f

@@ -5714,9 +5714,11 @@ c     gg=F(u)/u with F(u) as in Cohen paper
       
       dimension pvett(3),pvettn(3),dery0(3),dery0n(3)
       dimension avn(3,jmx,kmx)
-      dimension aplane(3,jmx,kmx),asip(jmx,kmx),asrp(jmx,kmx),
+      dimension aplane(3,jmx,kmx),ak0sip(jmx,kmx),asrp(jmx,kmx),
      .          taup(jmx,kmx)
       dimension gri(3,jmx,kmx)
+      dimension dxspv(jmx,kmx),dyspv(jmx,kmx),dzspv(jmx,kmx)
+      dimension zwjvv(jmx,kmx),zwjspvv(jmx,kmx)
       
 c      parameter(nxmax=2*jmx-1)
       parameter(nxmax=2*kmx)
@@ -5782,7 +5784,7 @@ c initialize grid dimension for spline interpolation
             do ll=1,3
               aplane(ll,j,k)=ywrk(ll,j,k)
             enddo
-            asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2
+            ak0sip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2
             asrp(j,k)=0
           enddo
         enddo  
@@ -5814,7 +5816,7 @@ c     prjection parallel to vg on the plane perpendicular to n0 passing through
           kktx=nrayth
           if(j.eq.1) kktx=1
           do k=1,kktx
-            asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2
+            ak0sip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2
             asrp(j,k)=ywrk(4,j,k)*(aplane(1,j,k)-ywrk(1,j,k))+
      .                ywrk(5,j,k)*(aplane(2,j,k)-ywrk(2,j,k))+
      .                ywrk(6,j,k)*(aplane(3,j,k)-ywrk(3,j,k))
@@ -5849,7 +5851,7 @@ c       Si evaluation on the projection plane(Taylor, first order)
           kktx=nrayth
           if(j.eq.1) kktx=1
           do k=1,kktx
-            asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2
+            ak0sip(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))
@@ -5873,10 +5875,51 @@ c
         kktx=nrayth
         if(j.eq.1) kktx=1
         do k=1,kktx
-          zwj=(dble(j-1)*rwmax/dble(nrayr-1))**2
+          zwjspvv(j,k)=ak0sip(j,k)
-          zwjsp=asip(j,k)
 c     .    +0.5*taup(j,k)
-c
+          dxspv(j,k)=aplane(1,j,k)-aplane(1,1,1)
+          dyspv(j,k)=aplane(2,j,k)-aplane(2,1,1)
+          dzspv(j,k)=aplane(3,j,k)-aplane(3,1,1)
+        enddo
+      enddo          
+c 
+      if(iplane.eq.2) then          
+        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
+        zwjvv(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))
+      else
+        do j=1,nrayr     
+          kktx=nrayth
+          if(j.eq.1) kktx=1
+          do k=1,kktx
+            zwjvv(j,k)=zwjspvv(j,k)
+          enddo
+        enddo    
+      endif    
+c      
+      do j=1,nrayr     
+        kktx=nrayth
+        if(j.eq.1) kktx=1
+        do k=1,kktx
+          zwj=zwjvv(j,k)
+          zwjsp=zwjspvv(j,k)       
+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)
@@ -5899,11 +5942,19 @@ 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
+          dxsp=dxspv(j,k)
+          dysp=dyspv(j,k)
+          dzsp=dzspv(j,k)
+          xspti=dxsp*csps1-dysp*snps1
+          yspti=(dxsp*snps1+dysp*csps1)*csth1-dzsp*snth1
+          zspti=(dxsp*snps1+dysp*csps1)*snth1+dzsp*csth1
+          rspti=sqrt(xspti**2+yspti**2)               
 c store x,y,z values for spline interpolation
           iray=iray+1
-          xtiv(iray)=xti
+          xtiv(iray)=xspti
-          ytiv(iray)=yti 
+          ytiv(iray)=yspti 
           zwjv(iray)=zwjsp
           srv(iray)=asrp(j,k)
 c initialize grid dimension for spline interpolation