From 9edc3079ba0612d3dfa1682bfc46da971f6bb3c2 Mon Sep 17 00:00:00 2001 From: Alberto Mariani Date: Thu, 18 Apr 2013 15:42:10 +0000 Subject: [PATCH] projxyzt with gaussian fit --- src/gray.f | 203 ++++++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 178 insertions(+), 25 deletions(-) diff --git a/src/gray.f b/src/gray.f index 46505e2..f2c08e2 100644 --- a/src/gray.f +++ b/src/gray.f @@ -507,7 +507,8 @@ c write(9,*) ' #istep j k xt yt zt rt psin' 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 psi w1 w2' + write(12,*) ' #i sst w1 w2 phiw rci1 rci2 phir errw errr '// + . 'dk1 dk2 dkpar phik dnpar' write(7,*)'#Icd Pa Jphimx dPdVmx '// .'rhotj rhotjava rhotp rhotpav '// .'drhotjava drhotpav ratjbmx stmx psipol chipol index_rt' @@ -5661,25 +5662,44 @@ c gg=F(u)/u with F(u) as in Cohen paper subroutine projxyzt(iproj,nfile) - implicit real*8 (a-h,o-z) - parameter(jmx=31,kmx=36) + implicit real*8 (a-h,o-z) + parameter(jmx=31,kmx=36,nmx=8000) + parameter(pi=3.14159265358979d0) dimension ywrk(6,jmx,kmx),ypwrk(6,jmx,kmx) + complex*16 ui,aac,bbc,ccc,ddc,aak,bbk,cck + parameter(ui=(0.0d0,1.0d0)) + dimension tauv(jmx,kmx,nmx),alphav(jmx,kmx,nmx) c common/nray/nrayr,nrayth common/wrk/ywrk,ypwrk - common/psinv11/psinv11 + common/derip11/dpsi11dr,dpsi11dz,psinv11,bvx11,bvy11,bvz11 common/istep/istep + common/rwmax/rwmax common/ss/st + common/dnpar/dnpara + common/atjki/tauv,alphav c - rtimn=1.d+30 - rtimx=-1.d-30 + x4m=0.0d0 + x3ym=0.0d0 + x2y2m=0.0d0 + xy3m=0.0d0 + y4m=0.0d0 + x2zrm=0.0d0 + xyzrm=0.0d0 + y2zrm=0.0d0 + x2zwm=0.0d0 + xyzwm=0.0d0 + y2zwm=0.0d0 + z2wm=0.0d0 + z2rm=0.0d0 c - jd=1 - if(iproj.eq.0) jd=nrayr-1 - do j=1,nrayr,jd - kkk=nrayth - if(j.eq.1) kkk=1 - do k=1,kkk + do j=1,nrayr + kktx=nrayth + if(j.eq.1) kktx=1 + do k=1,kktx + zwj=(dble(j-1)*rwmax/dble(nrayr-1))**2+ + . 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) @@ -5703,39 +5723,172 @@ c xti= dx*csps1-dy*snps1 yti=(dx*snps1+dy*csps1)*csth1-dz*snth1 zti=(dx*snps1+dy*csps1)*snth1+dz*csth1 - rti=sqrt(xti**2+yti**2) + rti=sqrt(xti**2+yti**2) c dirxt= (dirx*csps1-diry*snps1)/dir diryt=((dirx*snps1+diry*csps1)*csth1-dirz*snth1)/dir dirzt=((dirx*snps1+diry*csps1)*snth1+dirz*csth1)/dir +c + bxti= bvx11*csps1-bvy11*snps1 + byti=(bvx11*snps1+bvy11*csps1)*csth1-bvz11*snth1 +c bzti=(bvx11*snps1+bvy11*csps1)*snth1+bvz11*csth1 +c + rr11=sqrt(ywrk(1,1,1)**2+ywrk(2,1,1)**2) + dpsx=dpsi11dr*ywrk(1,1,1)/rr11 + dpsy=dpsi11dr*ywrk(2,1,1)/rr11 + dpsz=dpsi11dz + xdpsi=dpsx*csps1-dpsy*snps1 + ydpsi=(dpsx*snps1+dpsy*csps1)*csth1-dpsz*snth1 + zdpsi=(dpsx*snps1+dpsy*csps1)*snth1+dpsz*csth1 + sq=sqrt(xdpsi**2+ydpsi**2+zdpsi**2) + if(sq.gt.0.0d0) then + xdpsi=dx*xdpsi/sq + ydpsi=dx*ydpsi/sq + zdpsi=dx*zdpsi/sq + end if c if(k.eq.1) then xti1=xti yti1=yti zti1=zti rti1=rti - end if + xdpsi1=xdpsi + ydpsi1=ydpsi + zdpsi1=zdpsi + end if +c + if(istep.eq.0) + . write(10,111) istep,j,k,xti,yti,zti,dirxt,diryt,dirzt,dir - if(.not.(iproj.eq.0.and.j.eq.1)) - . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 -c - if(rti.ge.rtimx.and.j.eq.nrayr) rtimx=rti - if(rti.le.rtimn.and.j.eq.nrayr) rtimn=rti + if((iproj.eq.1).or.(iproj.eq.0.and.j.eq.nrayr)) + . write(nfile,111) istep,j,k,xti,yti,zti,rti, + . xdpsi,ydpsi,zdpsi c + x4m=x4m+xti**4 + x3ym=x3ym+xti**3*yti + x2y2m=x2y2m+xti**2*yti**2 + xy3m=xy3m+xti*yti**3 + y4m=y4m+yti**4 + x2zrm=x2zrm+xti**2*zti + xyzrm=xyzrm+xti*yti*zti + y2zrm=y2zrm+yti**2*zti + x2zwm=x2zwm+xti**2*zwj + xyzwm=xyzwm+xti*yti*zwj + y2zwm=y2zwm+yti**2*zwj + z2wm=z2wm+zwj*zwj + z2rm=z2rm+zti*zti end do -c -c if(.not.(iproj.eq.0.and.j.eq.1)) -c . write(nfile,111) istep,j,k,xti,yti,zti,rti,psinv11 +c + if((iproj.eq.1.and.j.gt.1).or.(iproj.eq.0.and.j.eq.nrayr)) + . write(nfile,111) istep,j,k,xti1,yti1,zti1,rti1, + . xdpsi1,ydpsi1,zdpsi1 if(iproj.eq.1) write(nfile,*) ' ' end do -c +c write(nfile,*) ' ' c - write(12,99) istep,st,psinv11,rtimn,rtimx +c computation of the SI paraboloid +c + denw= x2y2m**3 + x4m*xy3m**2 + x3ym**2*y4m - + . x2y2m*(2*x3ym*xy3m + x4m*y4m) + aaw= -(x2y2m*xy3m*xyzwm) + x2y2m**2*y2zwm - + . x3ym*xy3m*y2zwm + x3ym*xyzwm*y4m + + . x2zwm*(xy3m**2 - x2y2m*y4m) + ccw= x2y2m**2*xyzwm + x4m*xy3m*y2zwm - + . x2y2m*(x2zwm*xy3m + x3ym*y2zwm) + x2zwm*x3ym*y4m - + . x4m*xyzwm*y4m + bbw= x2y2m**2*x2zwm - x2zwm*x3ym*xy3m + x4m*xy3m*xyzwm + + . x3ym**2*y2zwm - x2y2m*(x3ym*xyzwm + x4m*y2zwm) + aaw=aaw/denw + bbw=bbw/denw + ccw=ccw/denw + phiw = 0.5d0*atan(ccw/(aaw-bbw+1.0d-32)) + phiwdeg=phiw*180.d0/pi + aaaw = 0.5d0*(aaw+bbw + (aaw-bbw)/cos(2.0d0*phiw)) + bbbw = aaw+bbw - aaaw + errw= 2.0d0*aaw*bbw*x2y2m + ccw**2*x2y2m - 2.0d0*aaw*x2zwm + + - 2.0d0*aaw*ccw*x3ym + aaw**2*x4m + 2.0d0*bbw*ccw*xy3m - + - 2.0d0*ccw*xyzwm - 2.0d0*bbw*y2zwm + bbw**2*y4m + z2wm + errw=sqrt(abs(errw)/dble((nrayr-1)*nrayth)) + wcsi = 1.0d0/sqrt(aaaw) + weta = 1.0d0/sqrt(bbbw) +c +c computation of paraboloid Sr=z-(ax^2+cxy+by^2)=const +c + denr= x2y2m**3 + x4m*xy3m**2 + x3ym**2*y4m - + . x2y2m*(2*x3ym*xy3m + x4m*y4m) + aar= -(x2y2m*xy3m*xyzrm) + x2y2m**2*y2zrm - + . x3ym*xy3m*y2zrm + x3ym*xyzrm*y4m + + . x2zrm*(xy3m**2 - x2y2m*y4m) + ccr= x2y2m**2*xyzrm + x4m*xy3m*y2zrm - + . x2y2m*(x2zrm*xy3m + x3ym*y2zrm) + x2zrm*x3ym*y4m - + . x4m*xyzrm*y4m + bbr= x2y2m**2*x2zrm - x2zrm*x3ym*xy3m + x4m*xy3m*xyzrm + + . x3ym**2*y2zrm - x2y2m*(x3ym*xyzrm + x4m*y2zrm) + aar=aar/denr + bbr=bbr/denr + ccr=ccr/denr + phir = 0.5d0*atan(ccr/(aar-bbr+1.0d-32)) + phirdeg=phir*180.d0/pi + aaar = 0.5d0*(aar+bbr + (aar-bbr)/cos(2.0d0*phir)) + bbbr = aar+bbr - aaar + errr= 2.0d0*aar*bbr*x2y2m + ccr**2*x2y2m - 2.0d0*aar*x2zrm + + - 2.0d0*aar*ccr*x3ym + aar**2*x4m + 2.0d0*bbr*ccr*xy3m - + - 2.0d0*ccr*xyzrm - 2.0d0*bbr*y2zrm + bbr**2*y4m + z2rm + errr=sqrt(abs(errr)/dble((nrayr-1)*nrayth)) + rcicsi=-2.0d0*aaar + rcieta=-2.0d0*bbbr +c +c computation of Fourier Transform of exp[i k0 (Sr + i Si)] +c k spectrum +c + aar=-ak0*aar + bbr=-ak0*bbr + ccr=-ak0*ccr + aac=aar+ui*aaw + bbc=bbr+ui*bbw + ccc=ccr+ui*ccw + ddc=ccc*ccc-4.0d0*aac*bbc + aak=bbc/ddc + bbk=aac/ddc + cck=-ccc/ddc + akw=dimag(aak) + bkw=dimag(bbk) + ckw=dimag(cck) + dkpar2=4.0d0*(bxti**2*bkw+byti**2*akw-bxti*byti*ckw)/ + . (4.0d0*akw*bkw-ckw*ckw) + phik = 0.5d0*atan(ckw/(akw-bkw+1.0d-32)) + aakw = 0.5d0*(akw+bkw + (akw-bkw)/cos(2.0d0*phik)) + bbkw = akw+bkw - aakw + dk1 = 1.0d0/sqrt(aakw) + dk2 = 1.0d0/sqrt(bbkw) +c co = cos(phik) +c si = sin(phik) +c bxtir= co*bxti+si*byti +c bytir=-si*bxti+co*byti +c dkpar2=bxtir**2/aakw+bytir**2/bbkw + dkpar=sqrt(dkpar2) + dnpar=dkpar/ak0 + dnpara=0.5d0*dnpar +c +c spectral distribution of electric field E(k): +c exp[-(N//-N//0)^2/(DNPAR^2)] +c dnpar^2 = 2 (Delta N//)^2 , Maj !!! +c +c in vacuum dkpar^2 = 4/w0^2 +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 return - 99 format(i5,12(1x,e16.8e3)) + 99 format(i5,22(1x,e16.8e3)) 111 format(3i5,12(1x,e16.8e3)) end + c c c