From 15526b5a8a30b81107c8c9c23129057f98bbb757 Mon Sep 17 00:00:00 2001 From: Maksim Rakitin Date: Mon, 19 Nov 2018 00:08:41 -0500 Subject: [PATCH] DOC: added SRW examples #15 and #16 --- docs/source/_cookbook/SRWLIB_Example15.py | 263 ++++++++++++++++++ docs/source/_cookbook/SRWLIB_Example16.py | 235 ++++++++++++++++ .../_cookbook/data_example_15/.gitignore | 0 .../_cookbook/data_example_16/.gitignore | 0 4 files changed, 498 insertions(+) create mode 100644 docs/source/_cookbook/SRWLIB_Example15.py create mode 100644 docs/source/_cookbook/SRWLIB_Example16.py create mode 100644 docs/source/_cookbook/data_example_15/.gitignore create mode 100644 docs/source/_cookbook/data_example_16/.gitignore diff --git a/docs/source/_cookbook/SRWLIB_Example15.py b/docs/source/_cookbook/SRWLIB_Example15.py new file mode 100644 index 0000000..c67300a --- /dev/null +++ b/docs/source/_cookbook/SRWLIB_Example15.py @@ -0,0 +1,263 @@ +""" +SRW Example #15 +*************** + +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 propagation of a gaussian beam through a drift with an analytical +estimation. + +Example Solution +---------------- +""" + +from __future__ import print_function +import uti_plot +from srwlib import * +from uti_math import matr_prod, fwhm + +print('SRWLIB Python Example # 15:') +print('Calculating propagation of a gaussian beam through a drift and comparison with the analytical calculation.') + +#************************************* Create examples directory if it does not exist +example_folder = 'data_example_15' +if not os.path.isdir(example_folder): + os.mkdir(example_folder) + +#************************************* Auxiliary functions +def AnalyticEst(PhotonEnergy, WaistPozition, Waist, Dist): + """Perform analytical estimation""" + Lam = 1.24e-6 * PhotonEnergy + zR = 3.1415 * Waist ** 2 / Lam + wRMSan = [] + for l in range(len(Dist)): + wRMSan.append(1 * Waist * sqrt(1 + (Lam * (Dist[l] - WaistPozition) / 4 / 3.1415 / Waist ** 2) ** 2)) + return wRMSan + + +def BLMatrixMult(LensMatrX, LensMatrY, DriftMatr, DriftMatr0): + """Computes envelope in free space""" + InitDriftLenseX = matr_prod(LensMatrX, DriftMatr0) + tRMSfunX = matr_prod(DriftMatr, InitDriftLenseX) + InitDriftLenseY = matr_prod(LensMatrY, DriftMatr0) + tRMSfunY = matr_prod(DriftMatr, InitDriftLenseY) + return (tRMSfunX, tRMSfunY) + + +def Container(DriftLength, f_x, f_y): + """Container definition""" + OpElement = [] + ppOpElement = [] + OpElement.append(SRWLOptL(_Fx=f_x, _Fy=f_y, _x=0, _y=0)) + OpElement.append(SRWLOptD(DriftLength)) + ppOpElement.append([1, 1, 1.0, 0, 0, 1.0, 1.0, 1.0, 1.0]) # note that I use sel-adjust for Num grids + ppOpElement.append([1, 1, 1.0, 0, 0, 1.0, 1.0, 1.0, 1.0]) + OpElementContainer = [] + OpElementContainer.append(OpElement[0]) + OpElementContainer.append(OpElement[1]) + OpElementContainerProp = [] + OpElementContainerProp.append(ppOpElement[0]) + DriftMatr = [[1, DriftLength], [0, 1]] + OpElementContainerProp.append(ppOpElement[1]) + LensMatrX = [[1, 0], [-1.0 / f_x, 1]] + LensMatrY = [[1, 0], [-1.0 / f_y, 1]] + opBL = SRWLOptC(OpElementContainer, OpElementContainerProp) + return (opBL, LensMatrX, LensMatrY, DriftMatr) + + +def qParameter(PhotonEnergy, Waist, RadiusCurvature): + """Computing complex q parameter""" + Lam = 1.24e-6 * PhotonEnergy + qp = (1.0 + 0j) / complex(1 / RadiusCurvature, -Lam / 3.1415 / Waist ** 2) + return qp, Lam + +#************************************* 1. Defining Beam structure +GsnBm = SRWLGsnBm() # Gaussian Beam structure (just parameters) +GsnBm.x = 0 # Transverse Coordinates of Gaussian Beam Center at Waist [m] +GsnBm.y = 0 +GsnBm.z = 0 # Longitudinal Coordinate of Waist [m] +GsnBm.xp = 0 # Average Angles of Gaussian Beam at Waist [rad] +GsnBm.yp = 0 +GsnBm.avgPhotEn = 0.5 # 5000 #Photon Energy [eV] +GsnBm.pulseEn = 0.001 # Energy per Pulse [J] - to be corrected +GsnBm.repRate = 1 # Rep. Rate [Hz] - to be corrected +GsnBm.polar = 1 # 1- linear horizontal +GsnBm.sigX = 1e-03 # /2.35 #Horiz. RMS size at Waist [m] +GsnBm.sigY = 2e-03 # /2.35 #Vert. RMS size at Waist [m] +GsnBm.sigT = 10e-12 # Pulse duration [fs] (not used?) +GsnBm.mx = 0 # Transverse Gauss-Hermite Mode Orders +GsnBm.my = 0 + +#************************************* 2. Defining wavefront structure +wfr = SRWLWfr() # Initial Electric Field Wavefront +wfr.allocate(1, 100, 100) # Numbers of points vs Photon Energy (1), Horizontal and Vertical Positions (dummy) +wfr.mesh.zStart = 3.0 # Longitudinal Position [m] at which Electric Field has to be calculated, i.e. the position of the first optical element +wfr.mesh.eStart = GsnBm.avgPhotEn # Initial Photon Energy [eV] +wfr.mesh.eFin = GsnBm.avgPhotEn # Final Photon Energy [eV] +firstHorAp = 20.e-03 # First Aperture [m] +firstVertAp = 30.e-03 # [m] +wfr.mesh.xStart = -0.5 * firstHorAp # Initial Horizontal Position [m] +wfr.mesh.xFin = 0.5 * firstHorAp # Final Horizontal Position [m] +wfr.mesh.yStart = -0.5 * firstVertAp # Initial Vertical Position [m] +wfr.mesh.yFin = 0.5 * firstVertAp # Final Vertical Position [m] +DriftMatr0 = [[1, wfr.mesh.zStart], [0, 1]] + +#************************************* 3. Setting up propagation parameters +sampFactNxNyForProp = 5 # sampling factor for adjusting nx, ny (effective if > 0) +arPrecPar = [sampFactNxNyForProp] + +#************************************* 4. Defining optics properties +f_x = 3e+0 # focusing strength, m in X +f_y = 4e+0 # focusing strength, m in Y +StepSize = 0.3 # StepSize in meters along optical axis +InitialDist = 0 # Initial drift before start sampling RMS x/y after the lens +TotalLength = 6.0 # Total length after the lens +NumSteps = int((TotalLength - InitialDist) / StepSize) # Number of steps to sample RMS x/y after the lens + +#************************************* 5. Starting the cycle through the drift after the lens +xRMS = [] +yRMS = [] +s = [] +WRx = [] +WRy = [] +intensitiesToPlot = { # dict to store intensities and distances to plot at the end + 'intensity': [], + 'distance': [], + 'mesh': [], + 'mesh_x': [], + 'mesh_y': [] +} +# print('z xRMS yRMS mesh.nx mesh.ny xStart yStart') +# print('z xRMS yRMS nx ny xStart yStart') # OC +data_to_print = [] +#header = '{:3s} {:8s} {:8s} {:2s} {:2s} {:13s} {:13s}'.format('z', 'xRMS', 'yRMS', 'nx', 'ny', 'xStart', 'yStart') #MR27092016 +header = '{:3s} {:8s} {:8s} {:2s} {:2s} {:13s} {:13s}'.format('z', 'xFWHM', 'yFWHM', 'nx', 'ny', 'xStart', 'yStart') #OC18112017 +data_to_print.append(header) +print(header) +for j in range(NumSteps): + #********************************* Calculating Initial Wavefront and extracting Intensity: + srwl.CalcElecFieldGaussian(wfr, GsnBm, arPrecPar) + arI0 = array('f', [0] * wfr.mesh.nx * wfr.mesh.ny) # "flat" array to take 2D intensity data + srwl.CalcIntFromElecField(arI0, wfr, 6, 0, 3, wfr.mesh.eStart, 0, 0) # extracts intensity + wfrP = deepcopy(wfr) + + #********************************* Selecting radiation properties + Polar = 6 # 0- Linear Horizontal / 1- Linear Vertical 2- Linear 45 degrees / 3- Linear 135 degrees / 4- Circular Right / 5- Circular / 6- Total + Intens = 0 # 0=Single-e I/1=Multi-e I/2=Single-e F/3=Multi-e F/4=Single-e RadPhase/5=Re single-e Efield/6=Im single-e Efield + 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 + # plotNum = 1000 + + if InitialDist == 0.0: # plot initial intensity at 0 + intensitiesToPlot['intensity'].append(deepcopy(arI0)) + intensitiesToPlot['distance'].append(InitialDist) + intensitiesToPlot['mesh_x'].append(deepcopy([wfrP.mesh.xStart, wfrP.mesh.xFin, wfrP.mesh.nx])) + intensitiesToPlot['mesh'].append(deepcopy(wfrP.mesh)) + intensitiesToPlot['mesh_y'].append(deepcopy([wfrP.mesh.yStart, wfrP.mesh.yFin, wfrP.mesh.ny])) + + InitialDist = InitialDist + StepSize + (opBL, LensMatrX, LensMatrY, DriftMatr) = Container(InitialDist, f_x, f_y) + srwl.PropagElecField(wfrP, opBL) # Propagate E-field + + # plotMeshx = [plotNum * wfrP.mesh.xStart, plotNum * wfrP.mesh.xFin, wfrP.mesh.nx] + # plotMeshy = [plotNum * wfrP.mesh.yStart, plotNum * wfrP.mesh.yFin, wfrP.mesh.ny] + plotMeshx = [wfrP.mesh.xStart, wfrP.mesh.xFin, wfrP.mesh.nx] + plotMeshy = [wfrP.mesh.yStart, wfrP.mesh.yFin, wfrP.mesh.ny] + + #********************************* Extracting output wavefront + arII = array('f', [0] * wfrP.mesh.nx * wfrP.mesh.ny) # "flat" array to take 2D intensity data + arIE = array('f', [0] * wfrP.mesh.nx * wfrP.mesh.ny) + srwl.CalcIntFromElecField(arII, wfrP, Polar, Intens, DependArg, wfrP.mesh.eStart, 0, 0) + arIx = array('f', [0] * wfrP.mesh.nx) + srwl.CalcIntFromElecField(arIx, wfrP, 6, Intens, 1, wfrP.mesh.eStart, 0, 0) + arIy = array('f', [0] * wfrP.mesh.ny) + srwl.CalcIntFromElecField(arIy, wfrP, 6, Intens, 2, wfrP.mesh.eStart, 0, 0) + + if abs(InitialDist - 3.0) < 1e-10 or abs(InitialDist - 3.9) < 1e-10: # plot at these distances + intensitiesToPlot['intensity'].append(deepcopy(arII)) + # intensitiesToPlot['distance'].append(InitialDist) + intensitiesToPlot['distance'].append(round(InitialDist, 6)) # OC + intensitiesToPlot['mesh_x'].append(deepcopy(plotMeshx)) + intensitiesToPlot['mesh'].append(deepcopy(wfrP.mesh)) + intensitiesToPlot['mesh_y'].append(deepcopy(plotMeshy)) + + x = [] + y = [] + arIxmax = max(arIx) + arIxh = [] + arIymax = max(arIx) + arIyh = [] + for i in range(wfrP.mesh.nx): + x.append((i - wfrP.mesh.nx / 2.0) / wfrP.mesh.nx * (wfrP.mesh.xFin - wfrP.mesh.xStart)) + arIxh.append(float(arIx[i] / arIxmax - 0.5)) + for i in range(wfrP.mesh.ny): + y.append((i - wfrP.mesh.ny / 2.0) / wfrP.mesh.ny * (wfrP.mesh.yFin - wfrP.mesh.yStart)) + arIyh.append(float(arIy[i] / arIymax - 0.5)) + xRMS.append(fwhm(x, arIxh)) + yRMS.append(fwhm(y, arIyh)) + s.append(InitialDist) + (tRMSfunX, tRMSfunY) = BLMatrixMult(LensMatrX, LensMatrY, DriftMatr, DriftMatr0) + WRx.append(tRMSfunX) + WRy.append(tRMSfunY) + # print(InitialDist, xRMS[j], yRMS[j], wfrP.mesh.nx, wfrP.mesh.ny, wfrP.mesh.xStart, wfrP.mesh.yStart) + # print(round(InitialDist, 6), round(xRMS[j], 6), round(yRMS[j], 6), wfrP.mesh.nx, wfrP.mesh.ny, + # round(wfrP.mesh.xStart, 10), round(wfrP.mesh.yStart, 10)) # OC + data_row = '{:3.1f} {:8.6f} {:8.6f} {:2d} {:2d} {:12.10f} {:12.10f}'.format(InitialDist, xRMS[j], yRMS[j], wfrP.mesh.nx, + wfrP.mesh.ny, wfrP.mesh.xStart, wfrP.mesh.yStart) #MR27092016 + data_to_print.append(data_row) + print(data_row) + +#************************************* 6. Analytic calculations +xRMSan = AnalyticEst(GsnBm.avgPhotEn, wfr.mesh.zStart + s[0], GsnBm.sigX, s) +(qxP, Lam) = qParameter(GsnBm.avgPhotEn, GsnBm.sigX, wfr.mesh.zStart + s[0]) +qx0 = complex(0, 3.1415 / Lam * GsnBm.sigX ** 2) +qy0 = complex(0, 3.1415 / Lam * GsnBm.sigY ** 2) +Wthx = [] +Wthy = [] +for m in range(len(s)): + Wx = (WRx[m][0][0] * qx0 + WRx[m][0][1]) / (WRx[m][1][0] * qx0 + WRx[m][1][1]) # MatrixMultiply(WR,qxP) + RMSbx = sqrt(1.0 / (-1.0 / Wx).imag / 3.1415 * Lam) * 2.35 + Wthx.append(RMSbx) + Wy = (WRy[m][0][0] * qy0 + WRy[m][0][1]) / (WRy[m][1][0] * qy0 + WRy[m][1][1]) # MatrixMultiply(WR,qxP) + RMSby = sqrt(1.0 / (-1.0 / Wy).imag / 3.1415 * Lam) * 2.35 + Wthy.append(RMSby) + +#************************************* 7. Plotting +for i in range(len(intensitiesToPlot['intensity'])): + srwl_uti_save_intens_ascii(intensitiesToPlot['intensity'][i], intensitiesToPlot['mesh'][i], + '{}/intensity_{:.1f}m.dat'.format(example_folder, intensitiesToPlot['distance'][i]), + 0, ['Photon Energy', 'Horizontal Position', 'Vertical Position', ''], + _arUnits=['eV', 'm', 'm', 'ph/s/.1%bw/mm^2']) + uti_plot.uti_plot2d1d(intensitiesToPlot['intensity'][i], + intensitiesToPlot['mesh_x'][i], + intensitiesToPlot['mesh_y'][i], + x=0, y=0, + labels=['Horizontal Position', 'Vertical Position', + 'Intenisty at {} m'.format(intensitiesToPlot['distance'][i])], + units=['m', 'm', 'a.u.']) + +with open('{}/compare.dat'.format(example_folder), 'w') as f: + f.write('\n'.join(data_to_print)) + +try: + from matplotlib import pyplot as plt + fig = plt.figure() + ax = fig.add_subplot(111) + ax.plot(s, xRMS, '-r.', label="X envelope from SRW") + ax.plot(s, Wthx, '--ro', label="X envelope from analytical estimate") + ax.plot(s, yRMS, '-b.', label="Y envelope from SRW") + ax.plot(s, Wthy, '--bo', label="Y envelope from analytical estimate") + ax.legend() + ax.set_xlabel('Distance along beam line [m]') + #ax.set_ylabel('Horizontal and Vertical RMS beam sizes [m]') + ax.set_ylabel('Horizontal and Vertical FWHM beam sizes [m]') #OC18112017 + ax.set_title('Gaussian beam envelopes through a drift after a lens') + ax.grid() + plt.savefig('{}/compare.png'.format(example_folder)) + plt.show() +except: + pass + +print('Exit') diff --git a/docs/source/_cookbook/SRWLIB_Example16.py b/docs/source/_cookbook/SRWLIB_Example16.py new file mode 100644 index 0000000..c525861 --- /dev/null +++ b/docs/source/_cookbook/SRWLIB_Example16.py @@ -0,0 +1,235 @@ +""" +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 +---------------- +""" + +from __future__ import print_function # Python 2.7 compatibility +import uti_plot +from srwlib import * +from 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') diff --git a/docs/source/_cookbook/data_example_15/.gitignore b/docs/source/_cookbook/data_example_15/.gitignore new file mode 100644 index 0000000..e69de29 diff --git a/docs/source/_cookbook/data_example_16/.gitignore b/docs/source/_cookbook/data_example_16/.gitignore new file mode 100644 index 0000000..e69de29