diff --git a/Makefile b/Makefile index 173e322..66c7d4e 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ EXE=gray # Objects list 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 vpath %.f90 src @@ -11,7 +11,7 @@ vpath %.f src # Fortran compiler name and flags FC=gfortran -FFLAGS=-O3 +FFLAGS=-O3 -fcheck=all DIRECTIVES = -DREVISION="'$(shell svnversion src)'" @@ -22,7 +22,7 @@ $(EXE): $(OBJ) $(FC) $(FFLAGS) -o $@ $^ # 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 const_and_precisions.o: itm_types.o itm_constants.o itm_constants.o: itm_types.o diff --git a/src/gray.f b/src/gray.f index 98a3004..48360c7 100644 --- a/src/gray.f +++ b/src/gray.f @@ -505,7 +505,7 @@ c .'N Npl ki alpha tau Pt Jcd dIds nh iohkw index_rt' write(8,*) ' #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(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 '// @@ -541,7 +541,7 @@ c character*24 filenmeqq,filenmprf,filenmbm parameter(qe=4.8032d-10,me=9.1095d-28,vc=2.9979d+10) parameter(pi=3.14159265358979d0,cvdr=pi/180.0d0) - parameter(nmx=8000,nbb=1000) + parameter(nmx=8000,nbb=5000) real*8 rlim(nbb),zlim(nbb) c common/xgcn/xgcn @@ -1068,7 +1068,7 @@ c implicit real*8 (a-h,o-z) parameter(nnw=501,nnh=501) parameter(pi=3.14159265358979d0) - parameter(nbb=1000) + parameter(nbb=5000) character*48 stringa dimension fpol(nnw),pres(nnw),qpsi(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) + use singleton, only:fft implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36,nmx=8000) 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)) dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) 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 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) 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)) dimension xgridv(nxmax),ygridv(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 srint(nxmax*nxmax),ccexps(nxmax*nxmax) + dimension txs(nxest),tys(nxest) parameter(lwrkbsp=8*nxmax,liwrkbsp=2*nxmax) dimension wrkbsp(lwrkbsp),iwrkbsp(liwrkbsp) c @@ -5770,6 +5781,7 @@ c initialize grid dimension for spline interpolation aplane(ll,j,k)=ywrk(ll,j,k) enddo asip(j,k)=(dble(j-1)*rwmax/dble(nrayr-1))**2 + asrp(j,k)=0 enddo enddo endif @@ -5801,6 +5813,15 @@ c prjection parallel to vg on the plane perpendicular to n0 passing through if(j.eq.1) kktx=1 do k=1,kktx 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 endif @@ -5826,10 +5847,22 @@ 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+ - . 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)) + asip(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)) + 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 endif @@ -5838,40 +5871,18 @@ c kktx=nrayth if(j.eq.1) kktx=1 do k=1,kktx - zwj=asip(j,k)+ - . 0.5*(tauv(j,k,istep)-tauv(1,1,istep)) + zwj=(dble(j-1)*rwmax/dble(nrayr-1))**2 + zwjsp=asip(j,k) +c . +0.5*taup(j,k) c - 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) + 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) c dirx=ywrk(4,j,k) diry=ywrk(5,j,k) dirz=ywrk(6,j,k) 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 if(j.eq.1.and.k.eq.1) then csth1=dirz/dir @@ -5891,7 +5902,8 @@ c store x,y,z values for spline interpolation iray=iray+1 xtiv(iray)=xti ytiv(iray)=yti - zwjv(iray)=zwj + zwjv(iray)=zwjsp + srv(iray)=asrp(j,k) c initialize grid dimension for spline interpolation xmaxgrid=max(xmaxgrid,rti) c @@ -6054,8 +6066,10 @@ c c in common dnpara=dnpar/2 to match westerhof Delta function c Delta -> exp[-(...)^2/(2*DeltaQ)] 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 if(nrayr.gt.1) then @@ -6083,7 +6097,7 @@ c call interp spline . nkntx,tx,nknty,ty,ccexp,resid,wrk1,lwrk1,wrk2,lwrk2, . iwrk,kwrk,ierr) 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 i=1,nxgrid zwint(nxgrid*(i-1)+j)=0.0d0 @@ -6095,24 +6109,101 @@ c call eval spline . ygridv,nxgrid,zwint,wrkbsp,lwrkbsp,iwrkbsp,liwrkbsp,ierr) if (ierr.ne.0) print*, 'bispev:',istep,ierr 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 i=1,nxgrid ij=nxgrid*(i-1)+j - write(82,111) istep,i,j,xgridv(i),ygridv(j),zwint(ij), - . aaw*xgridv(i)**2+ccw*xgridv(i)*ygridv(j)+bbw*ygridv(j)**2 + aecompl(i,j)=exp(ui*(ak0*srint(ij)+ui*zwint(ij))) +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 + 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 end do write(82,*) '' end do write(82,*) '' - end if - - write(12,99) istep,st, - . wcsi,weta,phiwdeg,rcicsi,rcieta,phirdeg,errw,errr, - . dk1,dk2,dkpar,phik*180.0d0/pi,dnpar -c + adev=sqrt(adev/nindex) + write(83,*) istep,adev,akw,bkw,ckw,dk1,dk2 return 99 format(i5,22(1x,e16.8e3)) -111 format(3i5,12(1x,e16.8e3)) +111 format(3i5,30(1x,e16.8e3)) end c @@ -6808,7 +6899,7 @@ c wave vector and electric field after reflection in lab frame implicit none integer*4 ivac integer nbb,nlim,i,imax - parameter(nbb=1000) + parameter(nbb=5000) 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 xv0(3) diff --git a/src/grayl.f b/src/grayl.f index d1aa7b9..a1ba10d 100644 --- a/src/grayl.f +++ b/src/grayl.f @@ -10902,3 +10902,4 @@ c errest = 2.0d0*errest go to 82 end +