From 9df962f4207e7c3411e0568f554f7fa7081e6c73 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Thu, 9 May 2013 09:05:36 +0000 Subject: [PATCH] projxyzt fixed. now ywrk is left unchanged --- src/gray.f | 118 +++++++++++++++++++++++++---------------------------- 1 file changed, 56 insertions(+), 62 deletions(-) diff --git a/src/gray.f b/src/gray.f index bdeae1b..98a3004 100644 --- a/src/gray.f +++ b/src/gray.f @@ -5712,11 +5712,9 @@ c gg=F(u)/u with F(u) as in Cohen paper dimension xtiv(nrmax),ytiv(nrmax),zwjv(nrmax),w(nrmax) dimension pvett(3),pvettn(3),dery0(3),dery0n(3) - dimension deltax(3,jmx,kmx),avn(3,jmx,kmx),avmod(jmx,kmx) - dimension deltaxdotn0(jmx,kmx),avndotn0(jmx,kmx) - dimension aalpha(jmx,kmx),aplane(3,jmx,kmx) + dimension avn(3,jmx,kmx) + dimension aplane(3,jmx,kmx),asip(jmx,kmx) dimension gri(3,jmx,kmx) - dimension aincr(3,jmx,kmx),asip(jmx,kmx) c parameter(nxmax=2*jmx-1) parameter(nxmax=2*kmx) @@ -5768,51 +5766,33 @@ c initialize grid dimension for spline interpolation 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 - else - if(iplane.eq.2) then + 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(j,k)=sqrt(ypwrk(1,j,k)*ypwrk(1,j,k)+ - . ypwrk(2,j,k)*ypwrk(2,j,k)+ypwrk(3,j,k)*ypwrk(3,j,k)) - do ii=1,3 - deltax(ii,j,k)=ywrk(ii,1,1)-ywrk(ii,j,k) - avn(ii,j,k)=ypwrk(ii,j,k)/avmod(j,k) - enddo - deltaxdotn0(j,k)=deltax(1,j,k)*ywrk(4,1,1)+ - . deltax(2,j,k)*ywrk(5,1,1)+ - . deltax(3,j,k)*ywrk(6,1,1) - avndotn0(j,k)=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(j,k)=deltaxdotn0(j,k)/avndotn0(j,k) - do ll=1,3 - aplane(ll,j,k)=ywrk(ll,j,k)+aalpha(j,k)*avn(ll,j,k) - enddo - enddo - enddo -c ortogonal projection on the plane perpendicular to n0 passing through x11 - else - do j=2,nrayr - do k=1,nrayth - do ii=1,3 - deltax(ii,j,k)=ywrk(ii,1,1)-ywrk(ii,j,k) - enddo - deltaxdotn0(j,k)=deltax(1,j,k)*ywrk(4,1,1)+ - . deltax(2,j,k)*ywrk(5,1,1)+ - . deltax(3,j,k)*ywrk(6,1,1) - an02=ywrk(4,1,1)*ywrk(4,1,1)+ywrk(5,1,1)*ywrk(5,1,1)+ - . ywrk(6,1,1)*ywrk(6,1,1) - aalpha(j,k)=deltaxdotn0(j,k)/an02 - do ll=1,3 - aplane(ll,j,k)=ywrk(ll,j,k)+aalpha(j,k)*ywrk(ll+3,1,1) - enddo - enddo - enddo - end if + 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 @@ -5820,21 +5800,36 @@ c ortogonal projection on the plane perpendicular to n0 passing through x11 kktx=nrayth if(j.eq.1) kktx=1 do k=1,kktx - do ii=1,3 - aincr(ii,j,k)=aplane(ii,j,k)-ywrk(ii,j,k) - ywrk(ii,j,k)=aplane(ii,j,k) - enddo + asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2 enddo - enddo -c Si evaluation on the projection plane(Taylor, first order) + 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)*aincr(1,j,k)+ - . gri(2,j,k)*aincr(2,j,k)+ - . gri(3,j,k)*aincr(3,j,k) + . 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 @@ -5846,9 +5841,9 @@ c zwj=asip(j,k)+ . 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) - dz=ywrk(3,j,k)-ywrk(3,1,1) + dx=aplane(1,j,k)-aplane(1,1,1) + dy=aplane(2,j,k)-aplane(2,1,1) + dz=aplane(3,j,k)-aplane(3,1,1) c dirx=ywrk(4,j,k) diry=ywrk(5,j,k) @@ -5857,14 +5852,13 @@ c if (j>1) then k2=mod(k+kktx/4-1,kktx)+1 - dx2=ywrk(1,j,k2)-ywrk(1,1,1) - dy2=ywrk(2,j,k2)-ywrk(2,1,1) - dz2=ywrk(3,j,k2)-ywrk(3,1,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)*pvett(1)+ - . pvett(2)*pvett(2)+pvett(3)*pvett(3)) + 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