From d7ced56519da8b20034a942f7ba6800aa4f5abde Mon Sep 17 00:00:00 2001 From: Alberto Mariani Date: Wed, 29 May 2013 13:46:25 +0000 Subject: [PATCH] modified projxyzt: added window function and corrected the case iplane=2 --- src/gray.f | 71 ++++++++++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 61 insertions(+), 10 deletions(-) diff --git a/src/gray.f b/src/gray.f index acb5fd7..dee2966 100644 --- a/src/gray.f +++ b/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 - zwjsp=asip(j,k) + zwjspvv(j,k)=ak0sip(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 - ytiv(iray)=yti + xtiv(iray)=xspti + ytiv(iray)=yspti zwjv(iray)=zwjsp srv(iray)=asrp(j,k) c initialize grid dimension for spline interpolation