Update the exit_wave_geometry to handle paralellogram detectors

CDTools/tools/initializers.py

@@ -88,14 +88,22 @@ def exit_wave_geometry(det_basis, det_shape, wavelength, distance, center=None,
 
     
     # Finally, generate the basis for the exit wave in real space
-    # I believe this calculation is incorrect for non-rectangular
-    # detectors, because the real space basis should be related to the
-    # dual of the original basis. Leaving this for now since
-    # non-rectangular detectors are not a pressing concern.
-    basis_dirs = det_basis / t.norm(det_basis, dim=0)
-    real_space_basis = basis_dirs * wavelength * distance / \
-        (full_shape.to(t.float32) * t.norm(det_basis,dim=0))
-    
+
+    # This method should work for a general parallelogram
+    # shaped detector
+    det_shape = det_basis * full_shape.to(t.float32)
+    pinv_basis = t.Tensor(np.linalg.pinv(det_shape).transpose()).to(t.float32)
+    real_space_basis = pinv_basis * wavelength * distance
+
+    # This is definitely correct, but less simple. Included here
+    # So future me can check that both versions are consistent.
+    #oop_dir = np.cross(det_basis[:,0],det_basis[:,1])
+    #oop_dir /= np.linalg.norm(oop_dir)
+    #full_basis = np.array([np.array(det_basis[:,0]),np.array(det_basis[:,1]),oop_dir]).transpose()
+    #inv_basis = t.Tensor(np.linalg.inv(full_basis)[:2,:].transpose()).to(t.float32)
+    #real_space_basis = inv_basis*wavelength * distance / \
+    #    full_shape.to(t.float32)
+
     # Finally, convert the shape back to a torch.Size
     full_shape = t.Size([dim * oversampling for dim in full_shape]) 
 

CDTools/tools/propagators.py

@@ -139,21 +139,21 @@ def generate_generalized_angular_spectrum_propagator(shape, basis, wavelength, p
     array of parallelograms. In addition, there is an assumed phase ramp
     applied to the wavefield before propagation, defined such that a feature
     with uniform phase will propagate along the direction of the 
-    defined propagation vector. This helps simplify 
+    defined propagation vector. This decision provides the best numerical
+    stabilit and allows for the simple setup of light fields copropagating
+    with the coordinate system.
 
 
     Parameters
     ----------
     shape : array
         The shape of the arrays to be propagated
-    spacing : array
+    basis : array
         The (2x3) set of basis vectors describing the array to be propagated
     wavelength : float
         The wavelength of light to simulate propagation of
     propagation_vector : array
         The displacement to propagate the wavefield along.
-    tilt : float
-        The tilt, in radians, of the plane that the wavefield is defined on
 
     Returns
     -------
@@ -161,10 +161,17 @@ def generate_generalized_angular_spectrum_propagator(shape, basis, wavelength, p
         A phase mask which accounts for the phase change that each plane wave will undergo.
     """
 
-    ki = 2 * np.pi * fftpack.fftfreq(shape[0],spacing[0])
-    kj = 2 * np.pi * fftpack.fftfreq(shape[1],spacing[1])
-    Kj, Ki = np.meshgrid(kj,ki)
+    
+    sides = basis * shape
+    fourier_basis = 2 * np.pi * np.linalg.pinv(sides).transpose()
+    print(fourier_basis)
 
+    ki = 2 * np.pi * fftpack.fftfreq(shape[0],basis[1,0])
+    kj = 2 * np.pi * fftpack.fftfreq(shape[1],basis[0,1])
+    print(ki[1]-ki[0])
+    print(kj[1]-kj[0])
+    Kj, Ki = np.meshgrid(kj,ki)
+    
     # Define this as complex so the square root properly gives
     # k>k0 components imaginary frequencies    
     k0 = np.complex128((2*np.pi/wavelength))

tests/tools/test_propagators.py

@@ -87,6 +87,51 @@ def test_near_field():
     # Again, 10^-3 is about all the accuracy we can expect
     assert np.max(np.abs(Emz-Emz_t)) < 1e-3 * np.max(np.abs(Emz))
 
+    
+def test_generalized_near_field():
+
+    # The strategy is to compare the propagation of a gaussian beam to
+    # the propagation in the paraxial approximation.
+    
+    x = (np.arange(800) - 400) * 1.5e-9
+    y = (np.arange(1200) - 600) * 1e-9
+    Ys,Xs = np.meshgrid(y,x)
+    Rs = np.sqrt(Xs**2+Ys**2)
+
+    wavelength = 3e-9 #nm
+    sigma = 20e-9 #nm
+    z = 1000e-9 #nm
+
+    k = 2 * np.pi / wavelength
+    w0 = np.sqrt(2)*sigma
+    zr = np.pi * w0**2 / wavelength
+    wz = w0 * np.sqrt(1 + (z / zr)**2)
+    Rz =  z * (1 + (zr / z)**2)
+    
+    E0 = np.exp(-Rs**2 / w0**2)
+
+    # The analytical expression for propagation of a gaussian beam in the
+    # paraxial approx
+    Ez = w0 / wz * np.exp(-Rs**2 / wz**2) * np.exp(-1j * k * ( z + Rs**2 / (2 * Rz)) + 1j * np.arctan(z / zr))
+
+    asp = propagators.generate_angular_spectrum_propagator(
+        E0.shape,(1.5e-9,1e-9),wavelength,z,dtype=t.float64)
+    
+    Ez_t = propagators.near_field(cmath.complex_to_torch(E0),asp)
+    Ez_t = cmath.torch_to_complex(Ez_t)
+    
+    # Check for at least 10^-3 relative accuracy in this scenario
+    assert np.max(np.abs(Ez-Ez_t)) < 1e-3 * np.max(np.abs(Ez))
+
+
+    Emz = np.conj(Ez)
+
+    Emz_t = propagators.inverse_near_field(cmath.complex_to_torch(E0),asp)
+    Emz_t = cmath.torch_to_complex(Emz_t)    
+
+    # Again, 10^-3 is about all the accuracy we can expect
+    assert np.max(np.abs(Emz-Emz_t)) < 1e-3 * np.max(np.abs(Emz))
+
 
 def test_inverse_near_field():