fourier_transform function of GaussianSource should apply to discrete rather than analytic current amplitude

[oskooi]

As discussed in #1931, the fourier_transform member function of the GaussianSource class computes the Fourier transform of the analytic current amplitude G(t) rather than the discrete version as implemented in:

meep/src/meep.hpp

Lines 970 to 975 in 35b10e5

	// subclasses *can* override this method in order to specify the
	// current directly rather than as the derivative of dipole.
	// in that case you would probably ignore the dt argument.
	virtual std::complex<double> current(double time, double dt) const {
	return ((dipole(time + dt) - dipole(time)) / dt);
	}

meep/src/sources.cpp

Lines 97 to 116 in 35b10e5

	complex<double> gaussian_src_time::dipole(double time) const {
	double tt = time - peak_time;
	if (float(fabs(tt)) > cutoff) return 0.0;
	// correction factor so that current amplitude (= d(dipole)/dt) is
	// ~ 1 near the peak of the Gaussian.
	complex<double> amp = 1.0 / complex<double>(0, -2 * pi * freq);
	return exp(-tt * tt / (2 * width * width)) * polar(1.0, -2 * pi * freq * tt) * amp;
	}
	// (1/\sqrt{2*pi}) \int e^{i\omega t} G(t) dt
	// where G(t) is the *current* envelope, i.e. the time derivative
	// of the dipole envelope
	std::complex<double> gaussian_src_time::fourier_transform(const double f) {
	double omega = 2.0 * pi * f;
	double omega0 = 2.0 * pi * freq;
	double delta = (omega - omega0) * width;
	return width * polar(1.0, omega * peak_time) * exp(-0.5 * delta * delta);
	}

As it is currently set up, fourier_transform does not seem to have much practical use. #1931 discussed the possibility of using the fourier_transform feature as a replacement for the normalization run when computing the incident fields for scattering calculations. This can be useful for a number of applications including Tutorial/Diffraction Spectrum of a Binary Grating which involves an incident planewave source in a homogeneous medium.

(This is somewhat related to a similar feature request for the GaussianBeamSource in #1444.)

[stevengj]

Yes, seems reasonable.

In principle there is an efficient semi-analytical algorithm for this (you take the continuous Fourier transform of a Gaussian, and then add aliased copies of this, which will converge exponentially rapidly). But in any case you would want to implement the brute-force DTFT formula first before using the clever formula.