near to far field transformation for cylindrical coordinates

[oskooi]

With the demonstration of launching planewaves in cylindrical coordinates (#1085), a useful additional capability would be to compute the focusing properties of a rotationally-symmetric, diffractive lens via a near-to-far field transformation.

This feature can be validated/demonstrated by comparison with analytic solutions using three possible tests/examples:

1. an axicon
2. a Fresnel zone plate (analytic solution)
3. a circular aperture (or Poisson spot)

(thanks to Prof. Rajesh Menon at the U. Utah for the suggestions.)

[stevengj]

This is certainly possible. The difficulty is simply getting the Green's function that relates the solution at a point in the (r,z) plane to an arbitrary (x,y,z) point for the near-to-far calculation.

The simplest way is probably to use our existing green3d function and integrate it numerically over φ. (For a periodic integral like this, you can just use equally spaced points and keep doubling the number of points until it converges to the desired tolerance, since equally spaced points converge exponentially fast for smooth periodic integrands.) This is relatively expensive but should be easy to implement.

A more sophisticated approach would be to try and do the integral over φ analytically. I vaguely recall that this is possible but leads to some nasty special function (a hypergeometric function?), but I'm not confident in my recollection here. Even if we do this, however, we'll want the brute-force numeric-integration approach anyway in order to check it, so I would suggest implementing the numerical-integration approach first.

[oskooi]

To clarify the implementation strategy: this requires modifying dft_near2far::farfield_lowlevel such that every point in the DFT near fields object (below) is integrated over φ using the fixed angular dependence exp(imφ) by converting from cylindrical to 3d Cartesian coordinates and calling green3d? Any special tricks to improve convergence of the integral?

meep/src/near2far.cpp

Lines 268 to 285 in d068502

	LOOP_OVER_IVECS(f->fc->gv, f->is, f->ie, idx) {
	IVEC_LOOP_LOC(f->fc->gv, x0);
	x0 = f->S.transform(x0, f->sn) + rshift;
	vec xs(x0);
	for (int i0 = -periodic_n[0]; i0 <= periodic_n[0]; ++i0) {
	if (periodic_d[0] != NO_DIRECTION)
	xs.set_direction(periodic_d[0], x0.in_direction(periodic_d[0]) + i0 * period[0]);
	double phase0 = i0 * periodic_k[0];
	for (int i1 = -periodic_n[1]; i1 <= periodic_n[1]; ++i1) {
	if (periodic_d[1] != NO_DIRECTION)
	xs.set_direction(periodic_d[1], x0.in_direction(periodic_d[1]) + i1 * period[1]);
	double phase = phase0 + i1 * periodic_k[1];
	std::complex<double> cphase = std::polar(1.0, phase);
	green(EH6, x, freq, eps, mu, xs, c0, f->dft[Nfreq * idx_dft + i]);
	for (int j = 0; j < 6; ++j)
	EH[i * 6 + j] += EH6[j] * cphase;
	}
	}