Intensity plot gets cut up when attempting to propagate finite airy beam

[Romutulus]

Hi. I just started playing around with LightPipes recently. I was using an arbitrary finite airy beam as a test run and it's very clear that the output gets cut up to squares or dots if I attempt to use Fresnel or Forvard (even Forward seems to be the case from a single-time use). This gets worse as I use longer propagation distances. I am wondering if I am doing something wrong on the setup. Below is my code:

wavelength = 1310*nm
size = 30.0*mm
N = 512*4
w = 0.001
a = 0.001
X1 = np.linspace(-0.015,0.015, N)
Y1 = np.linspace(-0.015,0.015, N)
def xy_fin_airy(xx,yy):
    phi = special.airy(-xx/w)[0] * np.exp(a * (-xx+yy)) * special.airy(yy/w)[0]
    return phi
def airy_E(x_space, y_space):
    e_field = np.zeros((N,N), dtype=complex)
    for i in range(len(x_space)):
        for k in range(len(y_space)):
            e_field[i,k] = xy_fin_airy(x_space[i], y_space[k])
    return e_field
E = airy_E(X1,Y1)
F=Begin(size,wavelength,N)
F=SubIntensity(F, E)
E = airy_E(X1,Y1)
F=Begin(size,wavelength,N)
F=SubIntensity(F, E)

F3 = Fresnel(F,20*cm)
I3 = Intensity(0,F3)
plt.figure(figsize=(10,10))
plt.imshow(I3)

Output Intensity:

Initial Intensity:

I have attempted to write Fresnel diffraction manually before trying out LightPipes, but I had problems where I had weird edge intensities that do not exist. However, I did not get my intensity cut up like shown above. Additionally, papers I've read have results showing that Airy beams would shift diagonally during propagation in spatial domain, but it doesn't seem like it is the case here. If anyone has any insight on that part, that would be great, however it might also just be me not properly understanding the theory.

[ldoyle]

Hi Romutulus,

I finally had the chance to look at your problem, very interesting physics these Airy beams, never heard of before!
First general things that I tried (which did not solve the problem, but are always good to test):

• Increase the number of points N (e.g. double). If this changes the outcome, the field is not sufficiently finely sampled, often in the case of e.g. a lens/ wavefront curvature which becomes too steep.
• Test for edge effects by truncating the field to the center region, e.g. with F=CircAperture(F, size/4) before propagation. Forvard implies periodic boundary conditions (see Manual section of LightPipes) so there needs to be a (large) margin to the edges where the intensity should be 0.

Both did not seem to influence the problem, but I found the error in a misunderstanding of the SubIntensity command. SubIntensity should accept only the positive, real valued squared field amplitude. SubPhase has to be used to also insert the phase as angle in the complex plane. For your E field, the values are all real, i.e. the phase is flat, but half of the values are negative, so the phase is Pi and not 0. Check this by plotting plt.imshow(np.abs(E)) and plt.imshow(np.angle(E)).

By passing E to SubIntensity, the square root was taken to get from intensity to field amplitude, which lead to imaginary numbers for the negative values of E. The field F now has a phases of 0 and pi/2, which of course is not correct given your input field.

After fixing this, the propagation looks alright, even though I cannot judge if it equals the physically expected result.
I also sped up the calculation of the field by using numpy.meshgrid, see online tutorials for best practices using numpy etc.

import numpy as np
from scipy import special
import matplotlib.pyplot as plt
from LightPipes import *

wavelength = 1310*nm
size = 30.0*mm
N = 512*4
w = 0.001
a = 0.001
X1 = np.linspace(-0.015,0.015, N)
Y1 = np.linspace(-0.015,0.015, N)

def xy_fin_airy(xx,yy):
    phi = special.airy(-xx/w)[0] * np.exp(a * (-xx+yy)) * special.airy(yy/w)[0]
    return phi

def airy_E_old(x_space, y_space):
    #NOTE in numpy, arrays are index[y, x], i.e. row major, so flipped here
    e_field = np.zeros((N,N), dtype=complex)
    for i in range(len(x_space)):
        for k in range(len(y_space)):
            e_field[i,k] = xy_fin_airy(x_space[i], y_space[k])
    return e_field

def airy_E(x_space, y_space):
    #speed up by vectorizing operations, numpy and scipy accept
    # multi-dimensional inputs
    xx, yy = np.meshgrid(x_space, y_space) #see Doc for example
    e_field = xy_fin_airy(xx, yy)
    #NOTE dtype of e_field is float, not complex, contains positive and negative
    # real values. If scipy.special.airy would return a complex, this would
    # automatically correctly be returned here
    return e_field

#NOTE E_field is now flipped from original code
E = airy_E(X1,Y1)
I_sub = np.abs(E)**2 #don't forget the square, can make all the difference sometimes!
Phase_sub = np.angle(E) #phase in radians

F=Begin(size,wavelength,N)
F=SubIntensity(F, I_sub)
F=SubPhase(F, Phase_sub)

F3 = Fresnel(F,20*cm)
I3 = Intensity(0,F3)
plt.figure(figsize=(10,10))
plt.imshow(I3)

PS: The image is now flipped compared to the old code, but I think this is correct since:

• airy function Ai gives oscillations for negative X/Y and falls off to 0 for positive values
• in the function definition xy_fin_airy, X sign is flipped and Y is unchanged
• the plot of imshow now shows oscillations in the top right quadrant, i.e. for positive X (since flipped in xy_fin_airy) and negative Y as expected

[ldoyle]

For the future, we should add some protection to SubIntensity to only accept positive, real valued arrays. In addition, we should add a command SubField to replace the complex field directly. For now, this can be done by directly accessing the fields F.field property.

[Romutulus]

Hi Idoyle. Thanks for helping me out. It seems to be working properly for me now (bar the shifting I mentioned earlier, but I suspect they were using the other coordinate system which scales based on propagation distance). I really do appreciate you taking the time to answer me. And yes, Airy beams do have pretty interesting physics. I find it weird that self-healing beams like these aren't really utilized to this day.

[FredvanGoor]

Hi Romutulus and Idoyle,
For me the Airy beam was new too. Very interesting and I want to add it to the LightPipes website. I studied the first article about it from G.A. Siviloglou et al. and wrote a script which propagates a laser beam through a spatial light modulator (phase-only) with a cubic phase profile followed by a Fourier transform lens. Indeed after the focus of the lens an Airy beam appeared.

Below is my script:

import numpy as np
from LightPipes import *
import matplotlib.pyplot as plt

'''
Generation of an Airy beam.
===========================
LightPipes for Python, Fred van Goor, 6-1-2022

A non-diffracting Airy beam can be generated by substituting a cubic phase distribution
in a laser beam. After a 2D Fourier transform with a positive lens the 
Airy beam will exist.
With the parameters below, the Airy beam exists from z = 20 cm to z= 150 cm.
Ref: 
G.A. Siviloglou, J. Broky, A. Dogariu and D.N. Christodoulides, PRL 99,213901(2007)
'''

wavelength = 780*nm #wavelength 
size = 25*mm #grid size
N = 4000 #grid dimension
alpha=130*np.pi #
w0=10*mm #beam waist laser
z=100*cm #propagation distance
f=55*cm #focal length of Fourier lens

#Generate input beam
F=Begin(size,wavelength,N)
F=GaussHermite(F,w0)
X,Y = F.mgrid_cartesian

#Cubic Phase Plate (CPP):
CPP=np.exp(1j*alpha*(X*X*X + Y*Y*Y)/wavelength)
F=SubPhase(F,CPP)
phase=Phase(F)

plt.figure(num=1, figsize=(6, 6))
#Plot phase distribution behind CPP:
plt.subplot(2,1,1);plt.imshow(phase,cmap='jet');
plt.title('cubic phase'); plt.axis ('off');
#Zoom in to center part of the phase plate:
plt.xlim(N/2-600, N/2+600)
plt.ylim(N/2+600, N/2-600)

#Fourier transform lens:
F=Lens(F,f)
F=Fresnel(F,f)

#Propagate a distance z:
F=Fresnel(F,z)
I=Intensity(F)

#Find the coordinates of the maximum intensity:
result = np.where(I == np.amax(I))
Coordinates = list(zip(result[0], result[1]))
for Coord in Coordinates:
    Xmax=X[Coord[0],Coord[1]]; Ymax=Y[Coord[0],Coord[1]]
    print(Xmax/um,Ymax/um) #Xmax should be equal to Ymax

#Plot intensity distribution at a distance z:
plt.subplot(2,1,2);plt.imshow(I,cmap='jet');
txt = 'Intensity at z = ' + str(z/cm) + ' cm' + '\n'\
    + 'Deflection = {:.0f} microns'
plt.title(txt.format(Xmax/um))

plt.axis ('off');
#For zoom, remove comments below
plt.xlim(result[0]-200, result[0]+200)
plt.ylim(result[1]+200, result[1]-200)
plt.show()

I also saw deflection of the beam.
I agree with the comments of Idoyle, increasing the grid dimension should help ...

[Romutulus]

Hi FredvanGoor,

Adding Airy beam to the LightPipes website would be very cool. As for finding the acceleration of the beam, yeah I realized that they do exist after doing a max check on the intensity array like you have. It was small in the microns level so adding more points definitely did help. I guess I was expecting a more visible acceleration and so did not bother to actually sift through the array for confirmation. With that said, I am actually glad that the acceleration isn't as noticeable as I expected it to be, as that would be quite problematic in an experimental setting.

[ldoyle]

Hi Fred,

I was also thinking about trying to replicate the results of that paper myself, but you beat me to it :)
However it seems your propagation does not look quite like in the paper, so I checked the numbers and found 1 error: in your example, the SubPhase command suffers from the same problem as the SubIntensity problem identified in this issue, i.e. supplying a complex array instead of a real valued only leads to weird results.
If I simplify the phase calculation and use the numbers from the paper, I get very close to their results. The shift over 20cm is on the order of 300um, after 25-30cm the beam diffracts and no longer shows any Airy shape.

@Romutulus Instead of the trick of a Gauss beam + SLM + lens, you should also be able to insert directly your Airy formula and compare the two. Both should lead to similar shifts/"acceleration" for similar beam parameters I presume.

Best, Lenny

import numpy as np
from LightPipes import *
import matplotlib.pyplot as plt

'''
Generation of an Airy beam.
===========================
LightPipes for Python, Fred van Goor, 6-1-2022

A non-diffracting Airy beam can be generated by substituting a cubic phase distribution
in a laser beam. After a 2D Fourier transform with a positive lens the 
Airy beam will exist.
With the parameters below, the Airy beam exists from z = 0 cm to z= 25 cm.
Ref: 
G.A. Siviloglou, J. Broky, A. Dogariu and D.N. Christodoulides, PRL 99,213901(2007)
'''

#Grid definitions
wavelength = 488*nm #wavelength 
size = 20*mm #grid size
N = 1024 #grid dimension

#Gaussian beam definition, relation between w0 and FWHM(Intensity) is sqrt(2ln2)
w0=1/np.sqrt(2*np.log(2)) * 6.7*mm #beam waist laser

#SLM phase modulator / cubic phase plate:
# amplitude of modulation -20..20pi over 2cm (+-1cm)
alpha = 20*np.pi

#geometrical setup:
# focal length of lens for Fourier transform to generate Airy beam
f=120*cm
# propagation distance after lens focus for Airy beam packet
z=20*cm #propagation distance

#Generate input beam
F=Begin(size,wavelength,N)
F=GaussHermite(F,w0)
Y, X = F.mgrid_cartesian #observe X Y flip

#Cubic Phase Plate (CPP):
X_slm = X/(1*cm) #scaled to size of SLM (+-1cm each direction)
Y_slm = Y/(1*cm) #to make sure alpha is applied over correct distance/scale

#this term leads to strange intensity profile in SLM plane:
#CPP=np.exp(1j*alpha*(X*X*X + Y*Y*Y)/wavelength)
#instead, this seems to work as expected:
CPP = alpha * (X_slm**3 + Y_slm**3)
CPP += np.pi #plus Pi to match colorbar in publication
F=SubPhase(F,CPP)

fig, [axl, axr] = plt.subplots(1,2, sharex=True, sharey=True)
fig.suptitle('Intensity and phase at SLM')
axl.imshow(Intensity(F), cmap='jet')
axr.imshow(Phase(F), cmap='gray')

#Fourier transform lens:
F=Lens(F,f)
F=Fresnel(F,f)
I=Intensity(F)

#Find the coordinates of the maximum intensity in focal point (start of Airy beam):
(ymax,), (xmax,) = np.where(I == np.amax(I)) #weird unpack due to return value of where()
print('X of max at z=0: [um]', X[ymax, xmax]/um)
print('Y of max at z=0: [um]', Y[ymax, xmax]/um)

fig, [axl, axr] = plt.subplots(1,2, sharex=True, sharey=True)
fig.suptitle('Intensity and phase at focal point of lens (FT of input)')
axl.imshow(Intensity(F), cmap='jet')
zoom_px = 100
axl.set_xlim(N//2-zoom_px,N//2+zoom_px)
axl.set_ylim(N//2+zoom_px,N//2-zoom_px) #observe flip
axr.imshow(Phase(F), cmap='gray')

#Propagate a distance z:
F=Fresnel(F,z)
I=Intensity(F)

#Find the coordinates of the maximum intensity:
(ymax,), (xmax,) = np.where(I == np.amax(I)) #weird unpack due to return value of where()
print('X of max at z: [um]', X[ymax, xmax]/um)
print('Y of max at z: [um]', Y[ymax, xmax]/um)

fig, [axl, axr] = plt.subplots(1,2, sharex=True, sharey=True)
fig.suptitle('Intensity and phase after {:.1f}cm propagation'.format(z/cm))
axl.imshow(Intensity(F), cmap='jet')
zoom_px = 100
axl.set_xlim(N//2-zoom_px,N//2+zoom_px)
axl.set_ylim(N//2+zoom_px,N//2-zoom_px) #observe flip
axr.imshow(Phase(F), cmap='gray')

[FredvanGoor]

Hi Lenny and Romutulus,
You are right, I made a mistake with the phase!
I want to add the script to the website. Maybe it is possible to demonstrate the 'self-healing' feature of this beam as well?
Suggestions are, of course, very welcome.
Fred.

[Romutulus]

Hi Idolye,

I agree with your statement. I will probably get to confirming that in a bit. With that said, while playing around with propagation distances, I have found visible deflection (which in my case would be in mm) with distances of meters. However, with the parameters that I have meter-scale Fresnel propagation results are a bit iffy, so I will need to check that again with the other options.

Hi Fred,

On the note of demonstrating self-healing, I have played around with the simple trial demonstrated by Chu et al (https://doi.org/10.1103/PhysRevA.85.013815). Here I've just added a small circular blockage at the main lobe and looked at how the beam dealt with it:

from LightPipes import *
import matplotlib.pyplot as plt
import numpy as np
from scipy import special

def fres_prop(Field, dist, X_grid, Y_grid):
    F = Fresnel(Field, dist)
    I = Intensity(0, F)
    plt.figure(figsize=(10,10))
    plt.imshow(I, cmap='jet')
    plt.axis ('off')
    #Tracking max intensity position
    # result = np.where(I == np.amax(I))
    # Coordinates = list(zip(result[0], result[1]))
    # for Coord in Coordinates:
    #     Xmax=X_grid[Coord[0],Coord[1]]; Ymax=Y_grid[Coord[0],Coord[1]]
    #     print(Xmax/um,Ymax/um)
    return F

def xy_fin_airy(xx,yy):
    phi = special.airy(xx/w)[0] * np.exp(a * (xx+yy)) * special.airy(+yy/w)[0]
    return phi

def airy_E(x_space, y_space):
     #see Doc for example
    e_field = xy_fin_airy(x_space, y_space)
    return e_field

wavelength = 1310*nm
size = 30.0*mm
N = 256 * 8
w = 0.001
a = 0.001
X1 = np.linspace(-0.015,0.015, N)
Y1 = -np.linspace(-0.015,0.015, N)
XX, YY = np.meshgrid(X1, Y1)

E = airy_E(XX,YY)
I_sub = np.abs(E)**2 
Phase_sub = np.angle(E) 

# res1 = np.where(I_sub == np.amax(I_sub))
# Coordinates = list(zip(res1[0], res1[1]))
# for Coord in Coordinates:
#     Xmax=XX[Coord[0],Coord[1]]; Ymax=YY[Coord[0],Coord[1]]
#     print(Xmax/um,Ymax/um)

F = Begin(size,wavelength,N)
F = SubIntensity(F, I_sub)
F = SubPhase(F, Phase_sub)
F_start = CircScreen(F, 0.5*mm, -0.5*mm, 0.5*mm )
I_start = Intensity(2,F_start)
plt.figure(figsize=(10,10))
plt.imshow(I_start, cmap='jet')
plt.axis ('off')

F1 = fres_prop(F_start, 30*cm, XX, YY)
F2 = fres_prop(F1, 30*cm, XX, YY)
F3 = fres_prop(F2, 30*cm, XX, YY)
F4 = fres_prop(F3, 30*cm, XX, YY)
F5 = fres_prop(F4, 30*cm, XX, YY)
F6 = fres_prop(F5, 30*cm, XX, YY)
F7 = fres_prop(F6, 30*cm, XX, YY)

At z = 0

At z = 90cm

At z = 210cm

Additionally, I have not tried it yet, but I think dealing with smaller magnitude but many noise would be another good example (like a noisy channel or turbulence).

On another note, I did not expect to get this much conversation from a coding question haha. I'm kind of excited. I have not really been able to talk about general photonics since dropping out of my PhD program from stress exacerbated by covid.

[FredvanGoor]

Hi, I used the paper of T. Latychevskaia, D. Schachtler, and H.W. Fink, Applied Optics Vol. 55, Issue 22, pp. 6095-6101 (2016) in the simulation below. As you can see the LightPipes results are in good agreement with their figure 2e. However their simulated and theoretical peak intensities, shown in their figure 2f, do not agree with their experiment. Below I have plotted their experimental values in the LightPipes results. I believe that these results have a better agreement with their experimental values.

import numpy as np
from LightPipes import *
import matplotlib.pyplot as plt
from matplotlib.offsetbox import OffsetImage, AnnotationBbox
import matplotlib.image as mpimg
 
'''
Generation of an Airy beam.
===========================
LightPipes for Python, Fred van Goor, 11-1-2022

A non-diffracting Airy beam can be generated by substituting a cubic phase distribution
in a laser beam. After a 2D Fourier transform with a positive lens the 
Airy beam will exist.
With the parameters below, the Airy beam exists from z = 0 cm to z= 30 cm.
Ref: 
G.A. Siviloglou, J. Broky, A. Dogariu and D.N. Christodoulides, PRL 99,213901(2007)
T. Latychevskaia, D. Schachtler, and H.W. Fink, Applied Optics Vol. 55, Issue 22, pp. 6095-6101 (2016)
'''

wavelength = 650*nm #wavelength 
N=624
size=N*32.0*um
dN1=N//2-40
dN2=N//2+40
beta=117/m #
w0=10*mm #beam waist laser
f=80*cm #focal length of Fourier lens
k=2*np.pi/wavelength

#Generate input beam
F0=Begin(size,wavelength,N)
F0=GaussHermite(F0,w0)

#Cubic Phase Plate (CPP):
X,Y = F0.mgrid_cartesian
c = 2 * np.pi * beta
c3 = (c**3)/3
CPP=c3*(X**3 + Y**3) + np.pi
F0=SubPhase(F0,CPP)
phase=Phase(F0)

#initiate figure 1:
fig1, axs1 = plt.subplots(nrows=3, ncols=1,figsize=(5.0,6.0))
figure1 = mpimg.imread('figure1.jpg')
imagebox = OffsetImage(figure1, zoom=0.15)
ab = AnnotationBbox(imagebox, (0.4, 0.6),frameon=False)
s1=r'Generation of an Airy beam.'+'\n\n'
s2=r'References:'+ '\n'\
   r'(1) G.A. Siviloglou, J. Broky, A. Dogariu and D.N. Christodoulides, Phys. Rev. Lett, 99,213901(2007)'+'\n'\
   r'(2) T. Latychevskaia, D. Schachtler, and H.W. Fink, Appl. Optics, 55, 6095-6101(2016)'+'\n\n'
s3 = r'LightPipes for Python,' + '\n' + 'AiryBeam.py' + '\n\n'\
    r'Parameters from ref 2:' + '\n'\
    r'$\lambda = {:4.1f}$'.format(wavelength/nm) + r' $nm$' + '\n'\
    r'$size = {:4.2f}$'.format(size/mm) + r' $mm$' + '\n'\
    r'$N = {:4d}$'.format(N) + '\n'\
    r'$w_0 = {:4.2f}$'.format(w0/mm) + r' $mm$'+ '\n'\
    r'$f = {:4.2f}$'.format(f/cm) + r' $cm$'+ '\n'\
    r'$\beta = {:4.2f}$'.format(beta/m) + r' $m^{-1}$'+ '\n'\
    r'${\copyright}$ Fred van Goor, January 2022'+ '\n'
axs1[0].text(0.,1.3,s1,fontsize=15, verticalalignment='top')
axs1[0].text(0.,1.0,s2+s3,fontsize=5, verticalalignment='top');axs1[0].axis('off')
axs1[1].add_artist(ab);axs1[1].axis('off')
axs1[1].text(0.0,1.1,'figure  1 from reference 2',fontsize=5, verticalalignment='top')
s=r'Phase distribution SLM'
axs1[2].imshow(phase,cmap='gray'); axs1[2].axis('off'); axs1[2].set_title(s)

#Fourier transform lens:
F0=Lens(F0,f)

#Propagate a distance f:
F0=Fresnel(F0,f)

def AiryBeam(z):
	'''
	Propagates the Airy beam a distance z from the focal point..
	returns a tuple:
	Position of the maximum intensity,
	The intensity distribution at a distance z,
	The peak intensity of the Airy beam at z.
	'''
	F=Fresnel(F0,z)
	I=Intensity(F)
	#Find the coordinates of the maximum intensity:
	result = np.where(I == np.amax(I))
	Coordinates = list(zip(result[0], result[1]))
	for Coord in Coordinates:
		Xmax=X[Coord[0],Coord[1]]; Ymax=Y[Coord[0],Coord[1]]
		print(Xmax/mm,Ymax/mm) #Xmax should be equal to Ymax   
	return Xmax,Ymax,I,I[result][0]

def xmax(z):
	'''
	Returns:
	The theoretical deflection at z,
	The peak intensity of the Airy beam at a distance z,
	according to Latychevskaia et al.
	'''
	c=z*f/(f+z)
	zdivc=(f+z)/f
	xmax=(1.0188*(beta*wavelength*f)-c**2/(4*k**2*(beta*wavelength*f)**3))*zdivc
	return xmax, (f/(f+z))**2

zstart=0*cm
zend=31*cm
z=np.arange(zstart,zend,1*cm)
n=z.shape[0]
deflection=np.zeros(n)
Imax=np.zeros(n)

#initiate figure 2:
fig2, axs2 = plt.subplots(nrows=4, ncols=2,figsize=(5.0,6.0))
fig2.suptitle('simulation with data Latychevskaia et al.')
fig2.subplots_adjust(hspace=0.8)
for i in range(0,n):
	Xmax,Ymax,I,_Imax=AiryBeam(z[i])

	if z[i] == 0*cm:
		axs2[0,0].imshow(I,cmap='jet')
		s=r'$z = {:2.1f}$'.format(z[i]/cm) + r'$cm$'
		axs2[0,0].imshow(I,cmap='jet'); axs2[0,0].axis('off'); axs2[0,0].set_title(s)
		axs2[0,0].set_xlim(dN1,dN2)
		axs2[0,0].set_ylim(dN2,dN1)

	if z[i]== 10*cm:
		s=r'$z = {:2.1f}$'.format(z[i]/cm) + r'$cm$'
		axs2[1,0].imshow(I,cmap='jet'); axs2[1,0].axis('off'); axs2[1,0].set_title(s)
		axs2[1,0].set_xlim(dN1,dN2)
		axs2[1,0].set_ylim(dN2,dN1)

	if z[i]== 20*cm:
		s=r'$z = {:2.1f}$'.format(z[i]/cm) + r'$cm$'
		axs2[2,0].imshow(I,cmap='jet'); axs2[2,0].axis('off'); axs2[2,0].set_title(s)
		axs2[2,0].set_xlim(dN1,dN2)
		axs2[2,0].set_ylim(dN2,dN1)

	if z[i]== 30*cm:
		s=r'$z = {:2.1f}$'.format(z[i]/cm) + r'$cm$'
		axs2[3,0].imshow(I,cmap='jet'); axs2[3,0].axis('off'); axs2[3,0].set_title(s)
		axs2[3,0].set_xlim(dN1,dN2)
		axs2[3,0].set_ylim(dN2,dN1)

	deflection[i]=Xmax
	Imax[i]=_Imax

gs = axs2[0,1].get_gridspec()
for ax in axs2[0:,-1]:
	ax.remove()
ax1=fig2.add_subplot(gs[0:2,1])
ax1.plot(z/cm,-deflection/mm,'ro',z/cm,-xmax(z)[0]/mm)
ax1.set_xlabel('z/cm'); ax1.set_ylabel('deflection/mm')
ax1.legend(('simulated', 'theory ref 2'),
           shadow=True, loc=(0.1, 0.8), handlelength=1.5, fontsize=5)
           
ax2=fig2.add_subplot(gs[2:,1])           
ax2.plot(z/cm,Imax/Imax[0],'ro',z/cm,xmax(z)[1])
ax2.set_xlabel('z/cm'); ax2.set_ylabel('Peak intensity')

#Experimental data from ref 2, figure 2f:
#In pixels, as good as possible....
X_Exp=np.array([122,151,184,213,245,276,308,336,366,397,429,458,491,520,551,582])
Y_Exp=np.array([106,139,177,216,198,207,247,266,289,327,336,354,350,360,360,366])
Imax_Exp=1.0-1.0/(428-106)*(Y_Exp-106)
z_Exp=30*cm/(582-122)*(X_Exp-122)
ax2.plot(z_Exp/cm,Imax_Exp,'+')
ax2.legend(('simulated', 'theory ref 2', 'experiment ref 2'),
           shadow=True, loc=(0.45, 0.8), handlelength=1.5, fontsize=5)

plt.show()

Figure 2 from Latychevskaia et al.