modified gauss fit and fft calculation in projxyzt

This commit is contained in:
Alberto Mariani 2013-05-23 17:09:54 +00:00
parent 4faca7994e
commit 892947b1de

View File

@ -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
@ -6144,70 +6146,67 @@ 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)
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))