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

projxyzt modified:in output correct inverse radii of curvature

Browse Source
This commit is contained in:
Alberto Mariani 2013-05-21 15:41:15 +00:00
parent 9df962f420
commit 185920b12d
3 changed files with 146 additions and 54 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Makefile
View File

@ -3,7 +3,7 @@ EXE=gray
# Objects list # Objects list
OBJ=gray.o grayl.o reflections.o green_func_p.o \ OBJ=gray.o grayl.o reflections.o green_func_p.o \
const_and_precisions.o itm_constants.o itm_types.o const_and_precisions.o itm_constants.o itm_types.o singleton.o
# Alternative search paths # Alternative search paths
vpath %.f90 src vpath %.f90 src
@ -11,7 +11,7 @@ vpath %.f src
# Fortran compiler name and flags # Fortran compiler name and flags
FC=gfortran FC=gfortran
FFLAGS=-O3 FFLAGS=-O3 -fcheck=all
DIRECTIVES = -DREVISION="'$(shell svnversion src)'" DIRECTIVES = -DREVISION="'$(shell svnversion src)'"
@ -22,7 +22,7 @@ $(EXE): $(OBJ)
$(FC) $(FFLAGS) -o $@ $^ $(FC) $(FFLAGS) -o $@ $^
# Dependencies on modules # Dependencies on modules
gray.o: green_func_p.o reflections.o gray.o: green_func_p.o reflections.o singleton.o
green_func_p.o: const_and_precisions.o green_func_p.o: const_and_precisions.o
const_and_precisions.o: itm_types.o itm_constants.o const_and_precisions.o: itm_types.o itm_constants.o
itm_constants.o: itm_types.o itm_constants.o: itm_types.o

195
src/gray.f
View File

@ -505,7 +505,7 @@ c
.'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt'
write(8,*) ' #istep j k xt yt zt rt psin' write(8,*) ' #istep j k xt yt zt rt psin'
write(9,*) ' #istep j k xt yt zt rt psin' write(9,*) ' #istep j k xt yt zt rt psin'
write(82,*) ' #istep j k xt yt zwspl zwparab' write(82,*) ' #istep j k xt yt zwspl zwparab srspl,etre,etim'
write(17,*) ' #sst Dr_11 Dr_Nr1 Di_Nr1' 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(33,*) ' #i j k sst x y R z psi tauv Npl alpha index_rt'
write(12,*) ' #i sst w1 w2 phiw rci1 rci2 phir errw errr '// write(12,*) ' #i sst w1 w2 phiw rci1 rci2 phir errw errr '//
@ -541,7 +541,7 @@ c
character*24 filenmeqq,filenmprf,filenmbm character*24 filenmeqq,filenmprf,filenmbm
parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10) parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10)
parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0) parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0)
parameter(nmx=8000,nbb=1000) parameter(nmx=8000,nbb=5000)
real*8 rlim(nbb),zlim(nbb) real*8 rlim(nbb),zlim(nbb)
c c
common/xgcn/xgcn common/xgcn/xgcn
@ -1068,7 +1068,7 @@ c
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(nnw=501,nnh=501) parameter(nnw=501,nnh=501)
parameter(pi=3.14159265358979d0) parameter(pi=3.14159265358979d0)
parameter(nbb=1000) parameter(nbb=5000)
character*48 stringa character*48 stringa
dimension fpol(nnw),pres(nnw),qpsi(nnw) dimension fpol(nnw),pres(nnw),qpsi(nnw)
dimension ffprim(nnw),pprim(nnw) dimension ffprim(nnw),pprim(nnw)
@ -5701,6 +5701,7 @@ c gg=F(u)/u with F(u) as in Cohen paper
subroutine projxyzt(iproj,nfile) subroutine projxyzt(iproj,nfile)
use singleton, only:fft
implicit real*8 (a-h,o-z) implicit real*8 (a-h,o-z)
parameter(jmx=31,kmx=36,nmx=8000) parameter(jmx=31,kmx=36,nmx=8000)
parameter(pi=3.14159265358979d0) parameter(pi=3.14159265358979d0)
@ -5709,11 +5710,12 @@ c gg=F(u)/u with F(u) as in Cohen paper
parameter(ui=(0.0d0,1.0d0)) parameter(ui=(0.0d0,1.0d0))
dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx)
parameter(nrmax=(jmx-1)*kmx+1) parameter(nrmax=(jmx-1)*kmx+1)
dimension xtiv(nrmax),ytiv(nrmax),zwjv(nrmax),w(nrmax) dimension xtiv(nrmax),ytiv(nrmax),zwjv(nrmax),srv(nrmax),w(nrmax)
dimension pvett(3),pvettn(3),dery0(3),dery0n(3) dimension pvett(3),pvettn(3),dery0(3),dery0n(3)
dimension avn(3,jmx,kmx) dimension avn(3,jmx,kmx)
dimension aplane(3,jmx,kmx),asip(jmx,kmx) dimension aplane(3,jmx,kmx),asip(jmx,kmx),asrp(jmx,kmx),
. taup(jmx,kmx)
dimension gri(3,jmx,kmx) dimension gri(3,jmx,kmx)
c parameter(nxmax=2*jmx-1) c parameter(nxmax=2*jmx-1)
@ -5727,7 +5729,16 @@ c parameter(nxmax=2*jmx-1)
parameter(kwrk=nrmax+(nxest-kspl-km)*(nxest-kspl-km)) parameter(kwrk=nrmax+(nxest-kspl-km)*(nxest-kspl-km))
dimension xgridv(nxmax),ygridv(nxmax) dimension xgridv(nxmax),ygridv(nxmax)
dimension zwint(nxmax*nxmax),ccexp(nxmax*nxmax) dimension zwint(nxmax*nxmax),ccexp(nxmax*nxmax)
complex*16 aecompl(nxmax,nxmax),aemodel(nxmax,nxmax)
complex*16 adiff(nxmax,nxmax)
complex*16 trick(2*nxmax-1,2*nxmax-1)
complex*16 trickdiff(2*nxmax-1,2*nxmax-1)
complex*16 aetcompl(2*nxmax-1,2*nxmax-1)
complex*16 adifft(2*nxmax-1,2*nxmax-1)
dimension akk(nxmax)
dimension tx(nxest),ty(nxest),wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk) dimension tx(nxest),ty(nxest),wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk)
dimension srint(nxmax*nxmax),ccexps(nxmax*nxmax)
dimension txs(nxest),tys(nxest)
parameter(lwrkbsp=8*nxmax,liwrkbsp=2*nxmax) parameter(lwrkbsp=8*nxmax,liwrkbsp=2*nxmax)
dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp) dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp)
c c
@ -5770,6 +5781,7 @@ c initialize grid dimension for spline interpolation
aplane(ll,j,k)=ywrk(ll,j,k) aplane(ll,j,k)=ywrk(ll,j,k)
enddo enddo
asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2 asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2
asrp(j,k)=0
enddo enddo
enddo enddo
endif endif
@ -5801,6 +5813,15 @@ c prjection parallel to vg on the plane perpendicular to n0 passing through
if(j.eq.1) kktx=1 if(j.eq.1) kktx=1
do k=1,kktx do k=1,kktx
asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2 asip(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))
if(istep.eq.0) then
taup(j,k)=0
else
taup(j,k)=(tauv(j,k,istep)-tauv(1,1,istep))+
. alphav(j,k,istep)*adsp
endif
enddo enddo
enddo enddo
endif endif
@ -5826,10 +5847,22 @@ c Si evaluation on the projection plane(Taylor, first order)
kktx=nrayth kktx=nrayth
if(j.eq.1) kktx=1 if(j.eq.1) kktx=1
do k=1,kktx do k=1,kktx
asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2+ asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2
. gri(1,j,k)*(aplane(1,j,k)-ywrk(1,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(2,j,k)*(aplane(2,j,k)-ywrk(2,j,k))
. gri(3,j,k)*(aplane(3,j,k)-ywrk(3,j,k)) . +gri(3,j,k)*(aplane(3,j,k)-ywrk(3,j,k))
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))
adsp=Sqrt((aplane(1,j,k)-ywrk(1,j,k))**2+
. (aplane(2,j,k)-ywrk(2,j,k))**2+
. (aplane(3,j,k)-ywrk(3,j,k))**2)
if(istep.eq.0) then
taup(j,k)=0
else
taup(j,k)=(tauv(j,k,istep)-tauv(1,1,istep))+
. alphav(j,k,istep)*adsp
endif
enddo enddo
enddo enddo
endif endif
@ -5838,40 +5871,18 @@ c
kktx=nrayth kktx=nrayth
if(j.eq.1) kktx=1 if(j.eq.1) kktx=1
do k=1,kktx do k=1,kktx
zwj=asip(j,k)+ zwj=(dble(j-1)*rwmax/dble(nrayr-1))**2
. 0.5*(tauv(j,k,istep)-tauv(1,1,istep)) zwjsp=asip(j,k)
c . +0.5*taup(j,k)
c c
dx=aplane(1,j,k)-aplane(1,1,1) dx=ywrk(1,j,k)-ywrk(1,1,1)
dy=aplane(2,j,k)-aplane(2,1,1) dy=ywrk(2,j,k)-ywrk(2,1,1)
dz=aplane(3,j,k)-aplane(3,1,1) dz=ywrk(3,j,k)-ywrk(3,1,1)
c c
dirx=ywrk(4,j,k) dirx=ywrk(4,j,k)
diry=ywrk(5,j,k) diry=ywrk(5,j,k)
dirz=ywrk(6,j,k) dirz=ywrk(6,j,k)
dir=sqrt(dirx*dirx+diry*diry+dirz*dirz) dir=sqrt(dirx*dirx+diry*diry+dirz*dirz)
if (j>1) then
k2=mod(k+kktx/4-1,kktx)+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)**2+pvett(2)**2+pvett(3)**2)
do ll=1,3
pvettn(ll)=pvett(ll)/pvettmod
dery0n(ll)=dery0(ll)/dery0mod
enddo
c write(*,*) 'dotn0',j,k,(pvettn(1)*dirx+
c . pvettn(2)*diry+pvettn(3)*dirz)/dir
c write(*,*) 'dotvg0',j,k,(pvettn(1)*dery0(1)+
c . pvettn(2)*dery0(2)+pvettn(3)*dery0(3))/dery0mod
endif
c dirx=dery0(1)
c diry=dery0(2)
c dirz=dery0(3)
c dir=dery0mod
c c
if(j.eq.1.and.k.eq.1) then if(j.eq.1.and.k.eq.1) then
csth1=dirz/dir csth1=dirz/dir
@ -5891,7 +5902,8 @@ c store x,y,z values for spline interpolation
iray=iray+1 iray=iray+1
xtiv(iray)=xti xtiv(iray)=xti
ytiv(iray)=yti ytiv(iray)=yti
zwjv(iray)=zwj zwjv(iray)=zwjsp
srv(iray)=asrp(j,k)
c initialize grid dimension for spline interpolation c initialize grid dimension for spline interpolation
xmaxgrid=max(xmaxgrid,rti) xmaxgrid=max(xmaxgrid,rti)
c c
@ -6054,8 +6066,10 @@ c
c in common dnpara=dnpar/2 to match westerhof Delta function c in common dnpara=dnpar/2 to match westerhof Delta function
c Delta -> exp[-(...)^2/(2*DeltaQ)] c Delta -> exp[-(...)^2/(2*DeltaQ)]
c c
write(12,99) istep,st,
. wcsi,weta,phiwdeg,rcicsi,rcieta,phirdeg,errw,errr,
. dk1,dk2,dkpar,phik*180.0d0/pi,dnpar
c
c Spline fit c Spline fit
if(nrayr.gt.1) then if(nrayr.gt.1) then
@ -6083,7 +6097,7 @@ c call interp spline
. nkntx,tx,nknty,ty,ccexp,resid,wrk1,lwrk1,wrk2,lwrk2, . nkntx,tx,nknty,ty,ccexp,resid,wrk1,lwrk1,wrk2,lwrk2,
. iwrk,kwrk,ierr) . iwrk,kwrk,ierr)
if (ierr.gt.0) then if (ierr.gt.0) then
print*, 'surfit:',istep,nkntx,nknty,ierr,resid print*, 'surfit zw:',istep,nkntx,nknty,ierr,resid
do j=1,nxgrid do j=1,nxgrid
do i=1,nxgrid do i=1,nxgrid
zwint(nxgrid*(i-1)+j)=0.0d0 zwint(nxgrid*(i-1)+j)=0.0d0
@ -6095,24 +6109,101 @@ c call eval spline
. ygridv,nxgrid,zwint,wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ierr) . ygridv,nxgrid,zwint,wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ierr)
if (ierr.ne.0) print*, 'bispev:',istep,ierr if (ierr.ne.0) print*, 'bispev:',istep,ierr
end if end if
call surfit(iopt,npoints,xtiv,ytiv,srv,w,xmin,xmax,ymin,ymax,
. kspl,kspl,sspl,nxest,nxest,nxest,eps,
. nkntxs,txs,nkntys,tys,ccexps,resid,wrk1,lwrk1,
. wrk2,lwrk2,iwrk,kwrk,ierr)
if (ierr.gt.0) then
print*, 'surfit sr:',istep,nkntxs,nkntys,ierr,resid
do j=1,nxgrid
do i=1,nxgrid
srint(nxgrid*(i-1)+j)=0.0d0
end do
end do
else
c call eval spline
call bispev(txs,nkntxs,tys,nkntys,ccexps,kspl,kspl,xgridv,
. nxgrid,ygridv,nxgrid,srint,wrkbsp,lwrkbsp,iwrkbsp,
. liwrkbsp,ierr)
if (ierr.ne.0) print*, 'bispev:',istep,ierr
end if
end if
c call fast Fourier transform 2D
deltak=2*pi/(xgridv(2)-xgridv(1))/(2*nxgrid-1)
do i=1,nxgrid
akk(i)=-(nxgrid-1)*deltak/2+(i-1)*deltak
enddo
do j=1,nxgrid do j=1,nxgrid
do i=1,nxgrid do i=1,nxgrid
ij=nxgrid*(i-1)+j ij=nxgrid*(i-1)+j
write(82,111) istep,i,j,xgridv(i),ygridv(j),zwint(ij), aecompl(i,j)=exp(ui*(ak0*srint(ij)+ui*zwint(ij)))
. aaw*xgridv(i)**2+ccw*xgridv(i)*ygridv(j)+bbw*ygridv(j)**2 c . *0.25*(1-cos((2*pi*(i-1))/(nxgrid-1)))*
c . (1-cos((2*pi*(j-1))/(nxgrid-1)))
aemodel(i,j)=exp(ui*(aac*xgridv(i)**2+bbc*ygridv(j)**2+
. ccc*xgridv(i)*ygridv(j)))
adiff(i,j)=aemodel(i,j)-aecompl(i,j)
end do end do
write(82,*) ''
end do end do
write(82,*) '' do j=1,2*nxgrid-1
do i=1,2*nxgrid-1
trick(i,j)=0
trickdiff(i,j)=0
end do
end do
do j=1,nxgrid
do i=1,nxgrid
trick(i+(nxgrid-1)/2,j+(nxgrid-1)/2)=aecompl(i,j)
trickdiff(i+(nxgrid-1)/2,j+(nxgrid-1)/2)=adiff(i,j)
end do
end do
aetcompl(1:2*nxgrid-1,1:2*nxgrid-1)=
. fft(trick(1:2*nxgrid-1,1:2*nxgrid-1))
adifft(1:2*nxgrid-1,1:2*nxgrid-1)=
. fft(trickdiff(1:2*nxgrid-1,1:2*nxgrid-1))
areaetc=real(aetcompl(1,1))
aimaetc=aimag(aetcompl(1,1))
acentral=sqrt(areaetc**2+aimaetc**2)
nnp=5
nindex=0
adev=0
do j=1,nxgrid
do i=1,nxgrid
ij=nxgrid*(i-1)+j
c areaetc=real(aetcompl(mod(i-1+3*(nxgrid-1)/2,2*nxgrid-1)+1
c . ,mod(j-1+3*(nxgrid-1)/2,2*nxgrid-1)+1))
c aimaetc=aimag(aetcompl(mod(i-1+3*(nxgrid-1)/2,2*nxgrid-1)+1,
c . mod(j-1+3*(nxgrid-1)/2,2*nxgrid-1)+1))
areaetc=real(aetcompl(mod(i+3*(nxgrid-1)/2,2*nxgrid-1)+1
. ,mod(j+3*(nxgrid-1)/2,2*nxgrid-1)+1))
aimaetc=aimag(aetcompl(mod(i+3*(nxgrid-1)/2,2*nxgrid-1)+1,
. mod(j+3*(nxgrid-1)/2,2*nxgrid-1)+1))
adifftr=real(adifft(mod(i+3*(nxgrid-1)/2,2*nxgrid-1)+1
. ,mod(j+3*(nxgrid-1)/2,2*nxgrid-1)+1))
adiffti=aimag(adifft(mod(i+3*(nxgrid-1)/2,2*nxgrid-1)+1,
. mod(j+3*(nxgrid-1)/2,2*nxgrid-1)+1))
write(82,111) istep,i,j,xgridv(i),ygridv(j),akk(i),akk(j),
. zwint(ij),
. aaw*xgridv(i)**2+ccw*xgridv(i)*ygridv(j)+bbw*ygridv(j)**2,
. srint(ij),areaetc,aimaetc,sqrt(areaetc**2+aimaetc**2)
. /acentral,exp(-(akw*akk(i)**2+bkw*akk(j)**2
. +ckw*akk(i)*akk(j))),sqrt(adifftr**2+adiffti**2),
. aecompl(i,j),aemodel(i,j),adiff(i,j)
if(abs(i-(nxgrid+1)/2).le.nnp.and.
. abs(j-(nxgrid+1)/2).le.nnp) then
nindex=nindex+1
adev=adev+(log(sqrt(areaetc**2+aimaetc**2)/acentral)
. +(akw*akk(i)**2+bkw*akk(j)**2+ckw*akk(i)*akk(j)))**2
endif endif
end do
write(12,99) istep,st, write(82,*) ''
. wcsi,weta,phiwdeg,rcicsi,rcieta,phirdeg,errw,errr, end do
. dk1,dk2,dkpar,phik*180.0d0/pi,dnpar write(82,*) ''
c adev=sqrt(adev/nindex)
write(83,*) istep,adev,akw,bkw,ckw,dk1,dk2
return return
99 format(i5,22(1x,e16.8e3)) 99 format(i5,22(1x,e16.8e3))
111 format(3i5,12(1x,e16.8e3)) 111 format(3i5,30(1x,e16.8e3))
end end
c c
@ -6808,7 +6899,7 @@ c wave vector and electric field after reflection in lab frame
implicit none implicit none
integer*4 ivac integer*4 ivac
integer nbb,nlim,i,imax integer nbb,nlim,i,imax
parameter(nbb=1000) parameter(nbb=5000)
real*8 st,rrm,zzm,psinv,dst,psdbnd,dstvac,smax real*8 st,rrm,zzm,psinv,dst,psdbnd,dstvac,smax
real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6) real*8 anv(3),xvstart(3),xvend(3),walln(3),y(6),dery(6)
real*8 xv0(3) real*8 xv0(3)

1
src/grayl.f
View File

@ -10902,3 +10902,4 @@ c
errest = 2.0d0*errest errest = 2.0d0*errest
go to 82 go to 82
end end