Phase artifacts after creating laser from experimental data using ```FromArrayProfile```

[TimoEichner]

Under some circumstances, there seems to be a problem with the way the field -- and in particular the phase/wavefront -- is interpolated when creating a laser using experimental data and FromArrayProfile. This seems to be due to the way that the field is interpolated from the arrays axes/grid to the new grid of the Laserinstance, which is using a single interpolator for the whole complex envelope. This seems to cause problems in the interpolation of the phase, but on the amplitude and only be a problem if the grid of the array and the grid of the new Laserare different.

Below are the wavefront and amplitude for three different cases to show the problem:

Case1

Laser created using exactly the same grid as in the experimental data

file = np.load('data_file.npz')

t = file['time']
x = file['x']
y = file['y']

dimensions = "xyt" 

# Grid is exactly the same as in the file
lo = (x.min(),y.min(),t.min()) 
hi = (x.max(),y.max(),t.max())  
num_points = (x.size, y.size, t.size)

profile = FromArrayProfile(wavelength= file['center_wavelength'],
                           pol=(0,1), 
                           array=file['E_xyt'],
                           dim=dimensions,
                           axes={'x':x,
                                 'y':y,
                                 't':t},
                           axes_order=['x', 'y', 't'])

laser = Laser(dimensions, lo, hi, num_points, profile)

Case 2

A new custom grid is made and the field is interpolated to that field by FromArrayProfile

file = np.load('data_file.npz')

t = file['time']
x = file['x']
y = file['y']

dimensions = "xyt" 

# Custom grid is provided
lo = (-300e-6,-300e-6,-500e-15)
hi = (300e-6,300e-6,500e-15)
num_points = (200, 200, 300)


profile = FromArrayProfile(wavelength= file['center_wavelength'],
                           pol=(0,1), 
                           array=file['E_xyt'],
                           dim=dimensions,
                           axes={'x':x,
                                 'y':y,
                                 't':t},
                           axes_order=['x', 'y', 't'])

laser = Laser(dimensions, lo, hi, num_points, profile)

Case 3

Before creating the Laser, the grid is interpolated to the new custom grid using a different interpolator that does not interpolate the whole complex field envelope abs(E_field)*np.exp(1j*np.unwrap(np.angle(E_field), axis=-1)), but instead uses two interpolators: One for the real part and one for the imaginary part (shown further below).

file = np.load('data_file.npz')

t = file['time']
x = file['x']
y = file['y']

dimensions = "xyt" 

# Custom grid is provided
lo = (-300e-6,-300e-6,-500e-15) 
hi = (300e-6,300e-6,500e-15) 
num_points = (200, 200, 300) 

# Field is interpolated with function shown below
E_new, x_new, y_new, t_new = interpolate_grid(file['E_xyt'], 
                                                hi_old=(x.max(),y.max(),t.max()),
                                                 lo_old=(x.min(),y.min(),t.min()),
                                                 hi_new=hi,
                                                 lo_new=lo,
                                                 shape_new=num_points
                                                )

profile = FromArrayProfile(wavelength= file['center_wavelength'],#cwl,#
                           pol=(0,1), 
                           array=E_new,
                           dim=dimensions,
                           axes={'x':x_new,
                                 'y':x_new,
                                 't':t_new},
                           axes_order=['x', 'y', 't'])

laser = Laser(dimensions, lo, hi, num_points, profile)

with the interpolation function being the following:

from scipy.interpolate import RegularGridInterpolator

def interpolate_grid(E, lo_old, hi_old, lo_new, hi_new, shape_new=()):
    x_old = np.linspace(lo_old[0], hi_old[0], E.shape[0])
    y_old = np.linspace(lo_old[1], hi_old[1], E.shape[1])
    t_old = np.linspace(lo_old[2], hi_old[2], E.shape[2])

    x_new = np.linspace(lo_new[0], hi_new[0], shape_new[0])
    y_new = np.linspace(lo_new[1], hi_new[1], shape_new[1])
    t_new = np.linspace(lo_new[2], hi_new[2], shape_new[2])

    # build separate interpolators
    interp_real = RegularGridInterpolator((x_old, y_old, t_old),
                                            E.real,
                                            bounds_error=False,
                                            fill_value=0.0)
    
    interp_imag = RegularGridInterpolator((x_old, y_old, t_old),
                                            E.imag,
                                            bounds_error=False,
                                            fill_value=0.0)
    
    # create the interpolation grid
    X, Y, T = np.meshgrid(x_new, 
                            y_new, 
                            t_new, 
                            indexing='ij')
    
    points = np.stack((X.ravel(), Y.ravel(), T.ravel()), axis=-1)
    
    # interpolate real and imaginary parts separately
    E_real_interp = interp_real(points).reshape(X.shape)
    E_imag_interp = interp_imag(points).reshape(X.shape)
    
    # recombine
    E_new = E_real_interp + 1j * E_imag_interp
    return E_new, x_new, y_new, t_new

I would suggest to add this method of interpolation also to FromArrayProfile