diff --git a/src/gray.f b/src/gray.f index 34785f0..acb5fd7 100644 --- a/src/gray.f +++ b/src/gray.f @@ -5731,10 +5731,10 @@ c parameter(nxmax=2*jmx-1) 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) + complex*16 atrasfmod,atrasfmod0,adifft0 dimension akk(nxmax) dimension tx(nxest),ty(nxest),wrk1(lwrk1),wrk2(lwrk2),iwrk(kwrk) dimension srint(nxmax*nxmax),ccexps(nxmax*nxmax) @@ -5772,7 +5772,7 @@ c z2wm=0.0d0 z2rm=0.0d0 c initialize grid dimension for spline interpolation - xmaxgrid=2*max(w0csi,w0eta) + xmaxgrid=0 !2*max(w0csi,w0eta) iray=0 if(iplane.le.1) then do j=1,nrayr @@ -5907,7 +5907,7 @@ c store x,y,z values for spline interpolation zwjv(iray)=zwjsp srv(iray)=asrp(j,k) c initialize grid dimension for spline interpolation - xmaxgrid=max(xmaxgrid,rti) + xmaxgrid=max(xmaxgrid,1.1*rti) c dirxt= (dirx*csps1-diry*snps1)/dir diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir @@ -6136,6 +6136,8 @@ c call fast Fourier transform 2D do i=1,nxgrid akk(i)=-(nxgrid-1)*deltak/2+(i-1)*deltak enddo + rhomax=xmaxgrid + rho0=0.5*rhomax do j=1,nxgrid do i=1,nxgrid ij=nxgrid*(i-1)+j @@ -6143,71 +6145,68 @@ c call fast Fourier transform 2D 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) + . ccc*xgridv(i)*ygridv(j))) + argcos=pi*(max(rho0,min(sqrt(xgridv(i)**2+ygridv(j)**2) + . ,rhomax))-rho0)/(rhomax-rho0) + adiff(i,j)=(aecompl(i,j)-aemodel(i,j))*0.5*(1+cos(argcos)) +c . *0.25*(1-cos((2*pi*(i-1))/(nxgrid-1))) +c . *(1-cos((2*pi*(j-1))/(nxgrid-1))) 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 + nnp=10 nindex=0 adev=0 adev2=0 + atrasfmod0=1.0d0/(-ui*Sqrt(-ddc)) +c controllare offset griglia + adifft0=adifft(1,1) 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) + adifftij=adifft(mod(i+3*(nxgrid-1)/2,2*nxgrid-1)+1 + . ,mod(j+3*(nxgrid-1)/2,2*nxgrid-1)+1) + atrasfmod=exp(-ui*(bbc*akk(i)**2+aac*akk(j)**2 + . -ccc*akk(i)*akk(j))/(-ddc))/(-ui*Sqrt(-ddc)) + write(81,99) istep,xgridv(i),ygridv(j), + . aemodel(i,j),aecompl(i,j),adiff(i,j) + write(82,99) istep,akk(i),akk(j), + . sqrt(real(atrasfmod/atrasfmod0)**2+ + . aimag(atrasfmod/atrasfmod0)**2), + . sqrt(real((atrasfmod+adifftij)/(atrasfmod0+adifft0))**2+ + . aimag((atrasfmod+adifftij)/(atrasfmod0+adifft0))**2) if(abs(i-(nxgrid+1)/2).le.nnp.and. . abs(j-(nxgrid+1)/2).le.nnp) then nindex=nindex+1 - adev2=adev2+(log(sqrt(areaetc**2+aimaetc**2)/acentral) - . +(akw*akk(i)**2+bkw*akk(j)**2+ckw*akk(i)*akk(j)))**2 - adev=adev+log(sqrt(areaetc**2+aimaetc**2)/acentral) - . +(akw*akk(i)**2+bkw*akk(j)**2+ckw*akk(i)*akk(j)) + adev2=adev2+(sqrt(real(atrasfmod/atrasfmod0)**2 + . +aimag(atrasfmod/atrasfmod0)**2) + . -sqrt(real((atrasfmod+adifftij)/(atrasfmod0+adifft0))**2 + . +aimag((atrasfmod+adifftij)/(atrasfmod0+adifft0))**2))**2 + adev=adev+sqrt(real(atrasfmod/atrasfmod0)**2 + . +aimag(atrasfmod/atrasfmod0)**2) + . -sqrt(real((atrasfmod+adifftij)/(atrasfmod0+adifft0))**2 + . +aimag((atrasfmod+adifftij)/(atrasfmod0+adifft0))**2) endif end do + write(81,*) '' write(82,*) '' end do + write(81,*) '' write(82,*) '' adev2=sqrt(adev2/nindex) adev=adev/nindex write(83,*) istep,adev2,adev,aar,bbr,ccr,aaw,bbw,ccw,akw,bkw,ckw, - . dk1,dk2 + . dk1,dk2,wcsi,weta,rcicsi,rcieta return 99 format(i5,22(1x,e16.8e3)) 111 format(3i5,30(1x,e16.8e3))