istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

projxyzt fixed. now ywrk is left unchanged

Browse Source
This commit is contained in:
Lorenzo Figini 2013-05-09 09:05:36 +00:00
parent 5c641a55ec
commit 9df962f420
1 changed files with 56 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

90
src/gray.f
View File

@ -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
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))
avmod=sqrt(ypwrk(1,j,k)**2+ypwrk(2,j,k)**2+ypwrk(3,j,k)**2)
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)
avn(ii,j,k)=ypwrk(ii,j,k)/avmod
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)+
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(j,k)=deltaxdotn0(j,k)/avndotn0(j,k)
aalpha=deltaxdotn0/avndotn0
do ll=1,3
aplane(ll,j,k)=ywrk(ll,j,k)+aalpha(j,k)*avn(ll,j,k)
aplane(ll,j,k)=ywrk(ll,j,k)+aalpha*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 ll=1,3
aplane(ll,1,1)=ywrk(ll,1,1)
enddo
@ -5820,11 +5800,26 @@ 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)
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
@ -5832,9 +5827,9 @@ c Si evaluation on the projection plane(Taylor, first order)
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