Fix sign error in generalized near field and associated test

CDTools/tools/propagators.py

@@ -452,7 +452,7 @@ def generate_generalized_angular_spectrum_propagator(shape, basis, wavelength, o
     # input plane.
 
     # This may have a sign error - must be checked
-    phase_mask = np.exp(1j * np.tensordot(offset_vector,K_xyz,axes=1))
+    phase_mask = np.exp(-1j * np.tensordot(offset_vector,K_xyz,axes=1))
     
     # Next, we apply a shift to the k-space vectors which sets up
     # propagation such that a uniform phase object will propagate along the
@@ -474,7 +474,7 @@ def generate_generalized_angular_spectrum_propagator(shape, basis, wavelength, o
 
     # Only implement the shift if the flag is set to True
     if propagation_vector is not None:
-        
+            
         prop_perpendicular = np.dot(perpendicular_dir, propagation_vector)
         prop_parallel = propagation_vector - perpendicular_dir \
             * prop_perpendicular
@@ -487,12 +487,11 @@ def generate_generalized_angular_spectrum_propagator(shape, basis, wavelength, o
         else:  
             k_offset = prop_parallel * k0 / np.linalg.norm(propagation_vector)
 
-        K_xyz = K_xyz + k_offset[:,None,None] 
+        K_xyz = K_xyz - k_offset[:,None,None] 
 
         # There apparently is a sign correction that I need to apply
         #sign_correction = np.sign(np.dot(perpendicular_dir,propagation_vector))
         sign_correction = np.sign(np.dot(offset_vector,propagation_vector))
-    
 
         # we also need to remove the z-dependence on the phase
         # This time, though, the z-dependence actually has to do with

tests/tools/test_propagators.py

@@ -9,6 +9,7 @@ import torch as t
 import pytest
 import scipy.misc
 from scipy.fftpack import fftshift, ifftshift
+from scipy import stats
 from matplotlib import pyplot as plt
 
 
@@ -206,10 +207,13 @@ def test_generalized_near_field():
     # explicitly included
 
     basis= np.array([[0,-1.5e-9],[-1e-9,0],[0,0]])
-    x = (np.arange(901) - 450) * 1.5e-9
-    y = (np.arange(1200) - 600) * 1e-9
-    Xs_0,Ys_0 = np.meshgrid(x,y)
-    Zs_0 = np.zeros(Xs_0.shape)
+    i_vec,j_vec = np.arange(901) - 450 ,np.arange(1200)-600
+    Is, Js = np.meshgrid(i_vec,j_vec,indexing='ij')
+    Xs_0,Ys_0,Zs_0 = np.tensordot(basis,np.stack([Is,Js]),axes=1)
+    #x = (np.arange(901) - 450) * 1.5e-9
+    #y = (np.arange(1200) - 600) * 1e-9
+    #Xs_0,Ys_0 = np.meshgrid(x,y)
+    #Zs_0 = np.zeros(Xs_0.shape)
 
     Positions = np.stack([Xs_0,Ys_0,Zs_0])
     
@@ -311,6 +315,14 @@ def test_generalized_near_field():
                         z_dir, z_dir_large]
     purposes = ['standard']*3 + ['both-rot']*3 + ['shear-rot']*3 + ['backward']*2
     
+    #rot_mats = [Ry.transpose()]
+    #offset = np.cross(np.dot(Ry.transpose(),basis)[:,0],
+    #                  np.dot(Ry.transpose(),basis)[:,1])
+    #offset /= np.linalg.norm(offset) / 3e-6
+    #offset_vecs = [-offset]#[shear_offset]
+    #propagation_vecs = [z_dir]
+    #purposes=['meh']
+    
     for purpose,rot_mat,offset_vec, propagation_vec \
         in zip(purposes,rot_mats,offset_vecs,propagation_vecs):
         
@@ -319,7 +331,7 @@ def test_generalized_near_field():
         new_basis = np.dot(rot_mat, basis)
         Xs_prop, Ys_prop, Zs_prop = np.stack([Xs,Ys,Zs_0]) \
             + offset_vec[:,None,None]
-        
+
         print('Propagate Along',propagation_vec)
         
         if str(propagation_vec) == 'perp':
@@ -344,10 +356,13 @@ def test_generalized_near_field():
         Ez_t = cmath.torch_to_complex(Ez_t)
         # Check for at least 10^-3 relative accuracy in this scenario
         if not np.max(np.abs(Ez-Ez_t)) < 1e-3 * np.max(np.abs(Ez)):
+        #if True:
             plt.close('all')
+            plt.imshow(np.angle(E0))
+            plt.figure()
             plt.imshow(np.angle(Ez))
             plt.figure()
-            plt.imshow(np.angle(Ez_t))
+            plt.imshow(np.abs(Ez-Ez_t)/np.max(np.abs(Ez)))
             plt.show()
             
         assert np.max(np.abs(Ez-Ez_t)) < 1e-3 * np.max(np.abs(Ez))
@@ -361,7 +376,35 @@ def test_generalized_near_field():
             
         print('Test Successful')
 
+    # One final test, to see if any of a few arbitrary rotations will
+    # change the predicted propagation if everything else is kept
+    # constant
+    Rrands = [stats.ortho_group.rvs(3) for i in range(3)]
+    Xs,Ys,Zs_0 = np.tensordot(Rboth,Positions,axes=1)
+    new_basis = np.dot(rot_mat, basis)
+    offset_vec = shear_back_offset
+    propagation_vec = z_dir
+    
+    E0 = get_E(Xs,Ys,Zs_0, correct=True)
+    asp = propagators.generate_generalized_angular_spectrum_propagator(
+                    E0.shape, new_basis, wavelength,offset_vec,
+                    dtype=t.float64, propagation_vector=propagation_vec)
+    Ez_t = propagators.near_field(cmath.complex_to_torch(E0),asp)
+    Ez_t = cmath.torch_to_complex(Ez_t)
+        
 
+    for Rrand in Rrands:
+        Xs,Ys,Zs_0 = np.tensordot(Rrand, np.tensordot(Rboth,Positions,axes=1),axes=1)
+        rot_offset = np.dot(Rrand, offset_vec)
+        rot_basis = np.dot(Rrand, new_basis)
+        rot_prop = np.dot(Rrand, propagation_vec)
+        asp = propagators.generate_generalized_angular_spectrum_propagator(
+                    E0.shape, rot_basis, wavelength, rot_offset,
+                    dtype=t.float64, propagation_vector=rot_prop)
+        Ez_rot_t = propagators.near_field(cmath.complex_to_torch(E0),asp)
+        Ez_rot_t = cmath.torch_to_complex(Ez_rot_t)       
+        assert np.max(np.abs(Ez_t-Ez_rot_t)) < 1e-3 * np.max(np.abs(Ez_t))
+        
 
 def test_inverse_near_field():