Use numpy function in `interpolate.initialize_charge_planes`
Currently the restriction with the gauss cutoff for the charge planes is calculated via loops. Use numpy functions and features (np.argmax(a>=cutoff)
) to prevent those loops.
Check performance.
Consider future switch to CUDA accelerated initial conditions. Are these features supported?