.. note:: :class: sphx-glr-download-link-note Click :ref:`here ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_cookbook_SRWLIB_Example16.py: SRW Example #16 *************** Problem ------- The example was created by Timur Shaftan (BNL) for RadTrack project (https://github.com/radiasoft/radtrack). Adapted by Maksim Rakitin (BNL). The purpose of the example is to demonstrate good agreement of the SRW simulation of intensity distribution after diffraction on a circular aperture with an analytical approximation. The example requires SciPy library to perform comparison. Example Solution ---------------- .. rst-class:: sphx-glr-horizontal * .. image:: /cookbook/images/sphx_glr_SRWLIB_Example16_001.png :class: sphx-glr-multi-img * .. image:: /cookbook/images/sphx_glr_SRWLIB_Example16_002.png :class: sphx-glr-multi-img * .. image:: /cookbook/images/sphx_glr_SRWLIB_Example16_003.png :class: sphx-glr-multi-img * .. image:: /cookbook/images/sphx_glr_SRWLIB_Example16_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out Out: .. code-block:: none SRWLIB Python Example # 16: Calculation of intensity distribution due to diffraction on a circular aperture. Comparison with an analytical distribution... Maximum intensity before and after propagation: [1.5047206578732545e+24, 9.560987085080381e+22] Number of horizontal mesh points before and after propagation: [1120, 1152] Number of vertical mesh points before and after propagation: [1120, 1152] Wavefront horizontal start coordinates [mm] before and after propagation: [-0.01, -0.025075270515281814] Wavefront horizontal end coordinates [mm] before and after propagation: [0.01, 0.02507527051528182] Wavefront vertical start coordinates [mm] before and after propagation: [-0.01, -0.025075270515281814] Wavefront vertical end coordinates [mm] before and after propagation: [0.01, 0.02507527051528182] Vertical FWHM [mm] before and after propagation: [4.706168420219193, 3.381710448537803] Analytical FWHM [mm]: 3.4023930367156776 done | .. code-block:: default from __future__ import print_function # Python 2.7 compatibility import srwpy.uti_plot as uti_plot from srwpy.srwlib import * from srwpy.uti_math import fwhm print('SRWLIB Python Example # 16:') print('Calculation of intensity distribution due to diffraction on a circular aperture.') try: from scipy.special import jv scipy_imported = True print('Comparison with an analytical distribution...') except: scipy_imported = False print('Can not import scipy for comparison of the numerically calculated intensity distribution with an analytical approximation.') #************************************* Create examples directory if it does not exist example_folder = 'data_example_16' # example data sub-folder name if not os.path.isdir(example_folder): os.mkdir(example_folder) strIntOutFileName = 'ex16_res_int.dat' # file name for output SR intensity data before propagation strIntOutFileNameProp = 'ex16_res_int_prop_{}m.dat' # file name for output SR intensity data after propagation #************************************* Perform SRW calculations # Gaussian beam definition: GsnBm = SRWLGsnBm() GsnBm.x = 0 # Transverse Coordinates of Gaussian Beam Center at Waist [m] GsnBm.y = 0 GsnBm.z = 0.0 # Longitudinal Coordinate of Waist [m] GsnBm.xp = 0 # Average Angles of Gaussian Beam at Waist [rad] GsnBm.yp = 0 GsnBm.avgPhotEn = 0.5 # Photon Energy [eV] GsnBm.pulseEn = 1.0E7 # Energy per Pulse [J] - to be corrected GsnBm.repRate = 1 # Rep. Rate [Hz] - to be corrected GsnBm.polar = 6 # 0- Linear Horizontal / 1- Linear Vertical 2- Linear 45 degrees / 3- Linear 135 degrees / 4- Circular Right / 5- Circular / 6- Total GsnBm.sigX = 2.0E-3 # Horiz. RMS size at Waist [m] GsnBm.sigY = 2.0E-3 # Vert. RMS size at Waist [m] GsnBm.sigT = 1E-12 # Pulse duration [fs] (not used?) GsnBm.mx = 0 # Transverse Gauss-Hermite Mode Orders GsnBm.my = 0 # Wavefront definition: wfr = SRWLWfr() NEnergy = 1 # Number of points along Energy Nx = 300 # Number of points along X Ny = 300 # Number of points along Y wfr.allocate(NEnergy, Nx, Ny) # Numbers of points vs Photon Energy (1), Horizontal and Vertical Positions (dummy) wfr.mesh.zStart = 0.35 # Longitudinal Position [m] at which Electric Field has to be calculated, i.e. the position of the first optical element wfr.mesh.eStart = 0.5 # Initial Photon Energy [eV] wfr.mesh.eFin = 0.5 # Final Photon Energy [eV] firstHorAp = 2.0E-3 # First Aperture [m] firstVertAp = 2.0E-3 # [m] wfr.mesh.xStart = -0.01 # Initial Horizontal Position [m] wfr.mesh.xFin = 0.01 # Final Horizontal Position [m] wfr.mesh.yStart = -0.01 # Initial Vertical Position [m] wfr.mesh.yFin = 0.01 # Final Vertical Position [m] # Precision parameters for SR calculation: meth = 2 # SR calculation method: 0- "manual", 1- "auto-undulator", 2- "auto-wiggler" npTraj = 1 # number of points for trajectory calculation (not needed) relPrec = 0.01 # relative precision zStartInteg = 0.0 # longitudinal position to start integration (effective if < zEndInteg) zEndInteg = 0.0 # longitudinal position to finish integration (effective if > zStartInteg) useTermin = 1 # Use "terminating terms" (i.e. asymptotic expansions at zStartInteg and zEndInteg) or not (1 or 0 respectively) sampFactNxNyForProp = 1 # sampling factor for adjusting nx, ny (effective if > 0) arPrecPar = [meth, relPrec, zStartInteg, zEndInteg, npTraj, useTermin, sampFactNxNyForProp] # Calculating initial wavefront: srwl.CalcElecFieldGaussian(wfr, GsnBm, arPrecPar) meshIn = deepcopy(wfr.mesh) wfrIn = deepcopy(wfr) arIin = array('f', [0] * wfrIn.mesh.nx * wfrIn.mesh.ny) srwl.CalcIntFromElecField(arIin, wfrIn, 0, 0, 3, wfr.mesh.eStart, 0, 0) arIinY = array('f', [0] * wfrIn.mesh.ny) srwl.CalcIntFromElecField(arIinY, wfrIn, 0, 0, 2, wfrIn.mesh.eStart, 0, 0) # extracts intensity # Plotting initial wavefront: plotMeshInX = [1000 * wfrIn.mesh.xStart, 1000 * wfrIn.mesh.xFin, wfrIn.mesh.nx] plotMeshInY = [1000 * wfrIn.mesh.yStart, 1000 * wfrIn.mesh.yFin, wfrIn.mesh.ny] srwl_uti_save_intens_ascii(arIin, wfrIn.mesh, '{}/{}'.format(example_folder, strIntOutFileName), 0, ['Photon Energy', 'Horizontal Position', 'Vertical Position', ''], #_arUnits=['eV', 'm', 'm', 'ph/s/.1%bw/mm^2']) _arUnits=['eV', 'm', 'm', 'arb. units']) uti_plot.uti_plot2d1d(arIin, plotMeshInX, plotMeshInY, #labels=['Horizontal Position [mm]', 'Vertical Position [mm]', 'Intensity Before Propagation [a.u.]']) labels=['Horizontal Position', 'Vertical Position', 'Intensity Before Propagation'], units=['mm', 'mm', 'arb. units']) # Element definition: apertureSize = 0.00075 # Aperture radius, m defaultDriftLength = 1.0 # Drift length, m OpElement = [] ppOpElement = [] OpElement.append(SRWLOptA('c', 'a', apertureSize, apertureSize)) ppOpElement.append([0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0, 0, 0, 0]) arIinymax = [] for k, driftLength in enumerate([0.0, defaultDriftLength]): ppOpElement.append([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0, 0, 0]) if driftLength: OpElement.append(SRWLOptD(driftLength)) # Drift space ppOpElement.append([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0, 0, 0]) opBL = SRWLOptC(OpElement, ppOpElement) srwl.PropagElecField(wfr, opBL) # Propagate Electric Field polarization = 6 #0- Linear Horizontal / 1- Linear Vertical / 2- Linear 45 degrees / 3- Linear 135 degrees / 4- Circular Right / 5- Circular / 6- Total intensity = 0 #0- Single-e Int. / 1- Multi-e Int. / 2- Single-e Flux / 3- Multi-e Flux / 4- Single-e Rad. Phase / 5- Re Single-e E-field / 6- Im Single-e E-field dependArg = 3 #0- vs e / 1- vs x / 2- vs y / 3- vs x&y / 4- vs x&e / 5- vs y&e / 6- vs x&y&e # Calculating output wavefront: arIout = array('f', [0] * wfr.mesh.nx * wfr.mesh.ny) # "flat" array to take 2D intensity data arII = arIout arIE = array('f', [0] * wfr.mesh.nx * wfr.mesh.ny) srwl.CalcIntFromElecField(arII, wfr, polarization, intensity, dependArg, wfr.mesh.eStart, 0, 0) arI1y = array('f', [0] * wfr.mesh.ny) arRe = array('f', [0] * wfr.mesh.ny) arIm = array('f', [0] * wfr.mesh.ny) srwl.CalcIntFromElecField(arI1y, wfr, polarization, intensity, 2, wfr.mesh.eStart, 0, 0) # extracts intensity # Normalize intensities: arI1ymax = max(arI1y) arIinymax.append(max(arIinY)) for i in range(len(arI1y)): arI1y[i] /= arI1ymax for i in range(len(arIinY)): arIinY[i] /= arIinymax[-1] # Plotting output wavefront: plotNum = 1000 plotMeshx = [plotNum * wfr.mesh.xStart, plotNum * wfr.mesh.xFin, wfr.mesh.nx] plotMeshy = [plotNum * wfr.mesh.yStart, plotNum * wfr.mesh.yFin, wfr.mesh.ny] srwl_uti_save_intens_ascii(arII, wfr.mesh, '{}/{}'.format(example_folder, strIntOutFileNameProp.format(driftLength)), 0, ['Photon Energy', 'Horizontal Position', 'Vertical Position', ''], _arUnits=['eV', 'm', 'm', 'ph/s/.1%bw/mm^2']) uti_plot.uti_plot2d1d(arII, plotMeshx, plotMeshy, #labels=['Horizontal Position [mm]', 'Vertical Position [mm]', 'Intenisty at {}m After Aperture [a.u.]'.format(driftLength)]) labels=['Horizontal Position', 'Vertical Position', 'Intenisty at {} m After Aperture'.format(driftLength)], units=['mm', 'mm', 'arb. units']) srwl.CalcIntFromElecField(arRe, wfr, polarization, 5, 2, wfr.mesh.eStart, 0, 0) srwl.CalcIntFromElecField(arIm, wfr, polarization, 6, 2, wfr.mesh.eStart, 0, 0) def calc_fwhm(intensities, wavefront, shift=0.5, mesh=True, factor=1e3): if mesh: y = [] for i in range(wfrIn.mesh.ny): y.append((i - wavefront.mesh.ny / 2.0) / wavefront.mesh.ny * (wavefront.mesh.yFin - wavefront.mesh.yStart)) else: y = wavefront renormed_intensities = [] max_intensity = max(intensities) for i in intensities: renormed_intensities.append(float(i / max_intensity - shift)) return fwhm(y, renormed_intensities) * factor # in [mm] by default parameters = [ ['Maximum intensity before and after propagation', arIinymax[0], arI1ymax], ['Number of horizontal mesh points before and after propagation', wfrIn.mesh.nx, wfr.mesh.nx], ['Number of vertical mesh points before and after propagation', wfrIn.mesh.ny, wfr.mesh.ny], ['Wavefront horizontal start coordinates [mm] before and after propagation', wfrIn.mesh.xStart, wfr.mesh.xStart], ['Wavefront horizontal end coordinates [mm] before and after propagation', wfrIn.mesh.xFin, wfr.mesh.xFin], ['Wavefront vertical start coordinates [mm] before and after propagation', wfrIn.mesh.yStart, wfr.mesh.yStart], ['Wavefront vertical end coordinates [mm] before and after propagation', wfrIn.mesh.yFin, wfr.mesh.yFin], ['Vertical FWHM [mm] before and after propagation', calc_fwhm(arIinY, wfrIn), calc_fwhm(arI1y, wfr)], ] for i in range(len(parameters)): print('{}{}: [{}, {}]'.format(' ', *parameters[i])) #************************************* 2. Defining parameters for analytic calculation lam = 2.4796e-6 # 1.2398/0.5 eV numPointsIn = len(arIinY) numPointsOut = len(arI1y) meshSize = float(wfr.mesh.xFin) #************************************* 3. Computing intensity distribution as per Born & Wolf, Principles of Optics th = [] sIn = [] sOut = [] analyticalIntens = [] for i in range(numPointsOut): thx = 2.0 * (i - numPointsOut / 2.0 + 0.5) * meshSize / numPointsOut / driftLength th.append(thx) sOut.append(thx * driftLength * 1000) if scipy_imported: for i in range(numPointsOut): x = 3.1415 * apertureSize * sin(th[i]) / lam analyticalIntens.append((2 * jv(1, x) / x) ** 2) print('{}Analytical FWHM [mm]: {}'.format(' ', calc_fwhm(analyticalIntens, sOut, mesh=False, factor=1.0))) for i in range(numPointsIn): sIn.append(2000.0 * (i - numPointsIn / 2.0) * float(wfrIn.mesh.xFin) / numPointsIn) #************************************* 4. Plotting try: from matplotlib import pyplot as plt fig = plt.figure(figsize=(10, 7)) ax = fig.add_subplot(111) ax.plot(sIn, arIinY, '--g.', label='SRW (before aperture)') ax.plot(sOut, arI1y, '-r.', label='SRW ({}m after aperture)'.format(driftLength)) if analyticalIntens: ax.plot(sOut, analyticalIntens, '-b.', label='Analytical estimation') ax.legend() ax.set_xlabel('Vertical Position [mm]') ax.set_ylabel('Normalized Intensity') ax.set_title('Intensity Before and After Propagation\n(cut vs vertical position at x=0)') ax.grid() plt.savefig('{}/compare.png'.format(example_folder)) plt.show() except: pass print('done') .. rst-class:: sphx-glr-timing **Total running time of the script:** ( 0 minutes 9.909 seconds) .. _sphx_glr_download_cookbook_SRWLIB_Example16.py: .. only :: html .. container:: sphx-glr-footer :class: sphx-glr-footer-example .. container:: sphx-glr-download :download:`Download Python source code: SRWLIB_Example16.py ` .. container:: sphx-glr-download :download:`Download Jupyter notebook: SRWLIB_Example16.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_