| 622 | */ |
| 623 | |
| 624 | void gaussianbeam::get_fields(std::complex<double> *EH, const vec &x) const { |
| 625 | double n = sqrt(eps * mu); |
| 626 | double k = 2 * pi * freq * n; |
| 627 | double ZR = sqrt(mu / eps); |
| 628 | double z0 = k * w0 * w0 / 2; |
| 629 | double kz0 = k * z0; |
| 630 | vec xrel = x - x0; |
| 631 | |
| 632 | vec zhat = kdir / abs(kdir); |
| 633 | double rho = |
| 634 | abs(vec(zhat.y() * xrel.z() - zhat.z() * xrel.y(), zhat.z() * xrel.x() - zhat.x() * xrel.z(), |
| 635 | zhat.x() * xrel.y() - zhat.y() * xrel.x())); |
| 636 | double zhatdotxrel = zhat & xrel; |
| 637 | |
| 638 | bool UseRescaledFG = false; |
| 639 | |
| 640 | complex<double> zc = complex<double>(zhatdotxrel, -z0); |
| 641 | complex<double> Rsq = rho * rho + zc * zc; |
| 642 | complex<double> R = sqrt(Rsq); |
| 643 | complex<double> kR = k * R, kR2 = kR * kR, kR3 = kR2 * kR; |
| 644 | complex<double> f, g, fmgbRsq; |
| 645 | if (abs(kR) > 1e-4) { |
| 646 | complex<double> coskR, sinkR; |
| 647 | if (fabs(imag(kR)) > 30.0) { |
| 648 | UseRescaledFG = true; |
| 649 | complex<double> ExpI = exp(complex<double>(0, 1.0) * real(kR)); |
| 650 | complex<double> ExpPlus = exp(imag(kR) - kz0); |
| 651 | complex<double> ExpMinus = exp(-(imag(kR) + kz0)); |
| 652 | coskR = 0.5 * (ExpI * ExpMinus + conj(ExpI) * ExpPlus); |
| 653 | sinkR = -0.5 * complex<double>(0, 1.0) * (ExpI * ExpMinus - conj(ExpI) * ExpPlus); |
| 654 | } |
| 655 | else { |
| 656 | coskR = cos(kR); |
| 657 | sinkR = sin(kR); |
| 658 | } |
| 659 | f = -3.0 * (coskR / kR2 - sinkR / kR3); |
| 660 | g = 1.5 * (sinkR / kR + coskR / kR2 - sinkR / kR3); |
| 661 | fmgbRsq = (f - g) / Rsq; |
| 662 | } |
| 663 | else { |
| 664 | complex<double> kR4 = kR2 * kR2; |
| 665 | // Taylor series expansion for small R |
| 666 | f = kR4 / 280.0 - kR2 / 10.0 + 1.0; |
| 667 | g = 3.0 * kR4 / 280.0 - kR2 / 5.0 + 1.0; |
| 668 | fmgbRsq = (kR4 / 5040.0 - kR2 / 140.0 + 0.1) * (k * k); |
| 669 | } |
| 670 | complex<double> i2fk = 0.5 * complex<double>(0, 1) * f * k; |
| 671 | |
| 672 | complex<double> E[3], H[3]; |
| 673 | for (int j = 0; j < 3; ++j) { |
| 674 | E[j] = complex<double>(0, 0); |
| 675 | H[j] = complex<double>(0, 0); |
| 676 | } |
| 677 | |
| 678 | double rnorm = |
| 679 | sqrt(real(E0[0]) * real(E0[0]) + real(E0[1]) * real(E0[1]) + real(E0[2]) * real(E0[2])); |
| 680 | if (rnorm > 1e-13) { |
| 681 | vec xhat = vec(real(E0[0]), real(E0[1]), real(E0[2])) / rnorm; |