diff --git a/src/gray-externals.f90 b/src/gray-externals.f90 index 5cc4be0..640fa8c 100644 --- a/src/gray-externals.f90 +++ b/src/gray-externals.f90 @@ -421,16 +421,16 @@ implicit none write(4,*)' #sst R z phi psi rhot ne Te Btot '// & 'Nperp Npl ki alpha tau Pt dIds nh iohkw index_rt ddr' - write(8,*) ' #istep j k xt yt zt rt psin' - write(9,*) ' #istep j k xt yt zt rt psin' + write(8,*) ' #istep j k xt yt zt rt' + write(9,*) ' #istep j k xt yt zt rt' write(17,*) ' #sst Dr_Nr1 Di_Nr1' write(33,*) ' #i jk sst x y R z psi tauv Npl alpha index_rt' - write(12,*) ' #i sst psi w1 w2' + write(12,*) ' #sst w1 w2' write(7,*)'#Icd Pa Jphip dPdVp rhotj rhotjava rhotp rhotpav '// & 'drhotjava drhotpav ratjamx ratjbmx stmx psipol chipol index_rt '// & 'Jphimx dPdVmx drhotj drhotp' write(48,*) '#rhop rhot Jphi Jcdb dPdV Icdins Pins' - write(66,*) "# psipol0 chipol0 powrfl" +! write(66,*) "# psipol0 chipol0 powrfl" end subroutine prfile diff --git a/src/graycore.f90 b/src/graycore.f90 index 2765e51..7d3c460 100644 --- a/src/graycore.f90 +++ b/src/graycore.f90 @@ -372,6 +372,14 @@ contains snps=zero end if +! Gaussian beam: exp[-ik0 zt] exp[-i k0/2 S(xt,yt,zt)] +! xt,yt,zt, cartesian coordinate system with zt along the beamline and xt in the z = 0 plane +! S(xt,yt,zt) = S_real +i S_imag = Qxx(zt) xt^2 + Qyy(zt) yt^2 + 2 Qxy(zt) xt yt +! (csiw, etaw) and (csiR, etaR) intensity and phase ellipse, rotated by angle phiw and phiR +! S(xt,yt,zt) = csiR^2 / Rccsi +etaR^2 /Rceta - i (csiw^2 Wcsi +etaw^2 Weta) +! Rccsi,eta curvature radius at the launching point +! Wcsi,eta =2/(k0 wcsi,eta^2) with wcsi,eta^2 beam size at the launching point + phiwrad = phiw*degree phirrad = phir*degree csphiw = cos(phiwrad) @@ -383,21 +391,21 @@ contains wweta = two/(ak0*weta**2) if(phir/=phiw) then - sk = rcicsi + rcieta - sw = wwcsi + wweta - dk = rcicsi - rcieta - dw = wwcsi - wweta - ts = -(dk*sin(2*phirrad) - ui*dw*sin(2*phiwrad)) - tc = (dk*cos(2*phirrad) - ui*dw*cos(2*phiwrad)) - phic = half*catand(ts/tc) - ddd = dk*cos(2*(phirrad+phic)) - ui*dw*cos(2*(phiwrad+phic)) - sss = sk - ui*sw - qi1 = half*(sss + ddd) - qi2 = half*(sss - ddd) - rci1 = dble(qi1) - rci2 = dble(qi2) - ww1 =-dimag(qi1) - ww2 =-dimag(qi2) + sk = rcicsi + rcieta + sw = wwcsi + wweta + dk = rcicsi - rcieta + dw = wwcsi - wweta + ts = -(dk*sin(2*phirrad) - ui*dw*sin(2*phiwrad)) + tc = (dk*cos(2*phirrad) - ui*dw*cos(2*phiwrad)) + phic = half*catand(ts/tc) + ddd = dk*cos(2*(phirrad+phic)) - ui*dw*cos(2*(phiwrad+phic)) + sss = sk - ui*sw + qi1 = half*(sss + ddd) + qi2 = half*(sss - ddd) + rci1 = dble(qi1) + rci2 = dble(qi2) + ww1 = -dimag(qi1) + ww2 = -dimag(qi2) else rci1 = rcicsi rci2 = rcieta @@ -415,13 +423,13 @@ contains qqxx = qi1*cos(phic)**2 + qi2*sin(phic)**2 qqyy = qi1*sin(phic)**2 + qi2*cos(phic)**2 - qqxy = -(qi1 - qi2)*sin(2*phic) - wwxx = -dimag(qqxx) - wwyy = -dimag(qqyy) - wwxy = -half*dimag(qqxy) - rcixx = dble(qqxx) - rciyy = dble(qqyy) - rcixy = half* dble(qqxy) + qqxy = -(qi1 - qi2)*sin(phic)*cos(phic) + wwxx = -dimag(qqxx) + wwyy = -dimag(qqyy) + wwxy = -dimag(qqxy) + rcixx = dble(qqxx) + rciyy = dble(qqyy) + rcixy = dble(qqxy) dqi1 = -qi1**2 dqi2 = -qi2**2 @@ -429,20 +437,20 @@ contains d2qi2 = 2*qi2**3 dqqxx = dqi1*cos(phic)**2 + dqi2*sin(phic)**2 dqqyy = dqi1*sin(phic)**2 + dqi2*cos(phic)**2 - dqqxy = -(dqi1 - dqi2)*sin(2*phic) + dqqxy = -(dqi1 - dqi2)*sin(phic)*cos(phic) d2qqxx = d2qi1*cos(phic)**2 + d2qi2*sin(phic)**2 d2qqyy = d2qi1*sin(phic)**2 + d2qi2*cos(phic)**2 - d2qqxy = -(d2qi1 - d2qi2)*sin(2*phic) + d2qqxy = -(d2qi1 - d2qi2)*sin(phic)*cos(phic) - dwwxx = -dimag(dqqxx) - dwwyy = -dimag(dqqyy) - dwwxy = -half*dimag(dqqxy) - d2wwxx = -dimag(d2qqxx) - d2wwyy = -dimag(d2qqyy) - d2wwxy = -half*dimag(d2qqxy) - drcixx = dble(dqqxx) - drciyy = dble(dqqyy) - drcixy = half* dble(dqqxy) + dwwxx = -dimag(dqqxx) + dwwyy = -dimag(dqqyy) + dwwxy = -dimag(dqqxy) + d2wwxx = -dimag(d2qqxx) + d2wwyy = -dimag(d2qqyy) + d2wwxy = -dimag(d2qqxy) + drcixx = dble(dqqxx) + drciyy = dble(dqqyy) + drcixy = dble(dqqxy) if(nrayr > 1) then dr = rwmax/dble(nrayr-1) @@ -507,7 +515,7 @@ contains dy0t = dcsiw*snphiw + detaw*csphiw x0t = uj(j)*dx0t y0t = uj(j)*dy0t - z0t = -half*(rcixx*x0t**2 + rciyy*y0t**2 + 2*rcixy*x0t*y0t) + z0t = -half*(rcixx*x0t**2 + rciyy*y0t**2) + rcixy*x0t*y0t dx0 = x0t*csps + snps*(y0t*csth + z0t*snth) dy0 = -x0t*snps + csps*(y0t*csth + z0t*snth)