{"id":695,"date":"2025-07-17T07:28:57","date_gmt":"2025-07-17T07:28:57","guid":{"rendered":"https:\/\/synchrotron-light.net\/wp\/?page_id=695"},"modified":"2025-07-20T05:13:26","modified_gmt":"2025-07-20T05:13:26","slug":"chapter-8-the-free-electron-laser","status":"publish","type":"page","link":"https:\/\/synchrotron-light.net\/wp\/?page_id=695","title":{"rendered":"Chapter 8: The free-electron laser"},"content":{"rendered":"

Below is a set of python codes associated with Chapter 8 of Daniele Pelliccia and David M. Paganin, “Synchrotron Light: A Physics Journey from Laboratory to Cosmos” (Oxford University Press, 2025).<\/p>\n

In order to run any of these python codes, you will need to include the following header.<\/p>\n

import numpy as np \r\nimport matplotlib.pyplot as plt\r\nfrom matplotlib.patches import Circle\r\nfrom matplotlib.patches import Ellipse\r\nfrom scipy.integrate import odeint\r\nimport math\r\n\r\nsynchrotron_style = {\r\n    'font.family': 'FreeSerif',\r\n    'font.size': 30,\r\n    'axes.linewidth': 2.0,\r\n    'axes.edgecolor': 'black',\r\n    'grid.color': '#5b5b5b',\r\n    'grid.linewidth': 1.2,\r\n}<\/pre>\n
Spatial modes<\/h2>\n
See Fig. 8.1.<\/p>\n
x = np.linspace(0, 2*np.pi, 1000)\r\nt1 = np.sin(x\/2)\r\nt2 = np.sin(x)\r\nt3 = np.sin(3*x\/2)\r\nt4 = np.sin(4*x\/2)\r\nt = [t1,t2,t3,t4]\r\ncolors = ['blue', 'red', 'green', 'gray']\r\nfig, ax = plt.subplots(4, 1,figsize=(25,10), sharex=True)\r\n\r\nfor i in range(1,11,1):\r\n    for j in range(4):\r\n        ax[j].plot(x, np.cos(np.pi*i\/11)*t[j], lw = 2, color=colors[j])\r\nfor j in range(4):    \r\n    ax[j].plot(x, t[j], lw = 2, color=colors[j])\r\n    ax[j].set_axis_off()\r\n\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
Planck’s law<\/h2>\n
See Fig. 8.2.<\/p>\n
h = 6.62607015e-34 # Planck's constant\r\nhbar = h\/(2*np.pi)\r\nkB = 1.380649e-23 # Boltzmann's constant\r\nc = 3.0e8 \r\n\r\nwith plt.style.context(('seaborn-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, axs = plt.subplots(1, 2,figsize=(16.5,6))\r\n    for T in [300, 1000, 3000]:\r\n        om = np.linspace(0.01*kB*T\/hbar, 20*kB*T\/hbar, 10000)\r\n        en_density = (hbar*om**3\/(np.pi**2 * c**3))\/( np.exp(hbar*om\/(kB*T)) -1)\r\n        axs[0].loglog(hbar*om\/(kB*T), en_density, lw = 4, label=\"T = \"+str(T)+\"K\")\r\n        axs[0].tick_params(axis='y', labelsize=24)\r\n        axs[0].tick_params(axis='x', labelsize=24, pad=5)\r\n        axs[0].set_ylabel('$\\\\langle \\\\mathcal{E}(\\\\omega)\\\\rangle$ ( $J s \/ m^{3})$', fontsize=28, labelpad = 15)\r\n        axs[0].set_xlabel('$\\\\hbar \\\\omega \/ (k_{B}T)$', fontsize=26)\r\n    axs[1].plot(hbar*om[:5000]\/(kB*T), en_density[:5000]\/en_density[:5000].max(), 'k', lw=4)  \r\n    axs[1].tick_params(axis='y', labelsize=24)\r\n    axs[1].tick_params(axis='x', labelsize=24)\r\n    axs[1].set_ylabel('Norm. energy density', fontsize=32)\r\n    axs[1].set_xlabel('$\\\\hbar \\\\omega \/ (k_{B}T)$', fontsize=26)\r\n\r\naxs[0].legend(fontsize=26, frameon=True, framealpha=1, handlelength=1, loc='lower left', bbox_to_anchor=(0.35, 1e-25))    \r\n\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
Ponderomotive phase<\/h2>\n
See Fig. 8.7.<\/p>\n
Phi_pos = -0.1\r\nlambda_u = 0.01\r\nk_u = 2*np.pi\/lambda_u\r\nc = 3.0e8\r\nomega_u = 2*k_u*c\r\nprint(1e-11*omega_u)\r\nN = 10\r\nt = np.linspace(0,N*2*np.pi\/omega_u,1000000)\r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax = plt.subplots(figsize=(12,5.8))\r\n    ax.plot(1e9*t, -np.ones_like(t)*np.sin(Phi_pos), lw = 4, label=\"$\\\\sin \\\\Phi^{+}$\")\r\n    ax.plot(1e9*t, -np.sin(Phi_pos-omega_u*t), lw = 4, label=\"$\\\\sin \\\\Phi^{-}$\")\r\n    ax.tick_params(axis='y', labelsize=26)\r\n    ax.tick_params(axis='x', labelsize=26)\r\n    ax.set_ylabel('Amplitude', fontsize=26)\r\n    ax.set_xlabel('$t$ (ns)', fontsize=26)\r\n    \r\n    ax.legend(fontsize=28, frameon=True, framealpha=1, handlelength=1)<\/pre>\n
\"\"<\/h1>\n
Numerical solution of the pendulum equation<\/h2>\n
See Fig. 8.8.<\/p>\n
def system(initial_conditions, t, Omega_sq):\r\n    # Set the initial conditions \r\n    Phi = initial_conditions[0]\r\n    eta = initial_conditions[1]\r\n    \r\n    # Write the coupled equations\r\n    d_Phi = eta\r\n    d_eta = -Omega_sq *np.sin(Phi)\r\n    \r\n    d_next_step = [d_Phi, d_eta]\r\n    \r\n    return d_next_step\r\n\r\nt = np.linspace(0,14.1,5000)\r\neta0 = np.arange(1.18,4.18, 0.2)\r\nOmega_sq=1\r\n\r\nPhi_sep = np.linspace(-np.pi,np.pi,256)\r\neta_sep = np.sqrt(2*Omega_sq*(np.cos(Phi_sep)+1))\r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax = plt.subplots(figsize=(10.3,10))\r\n    for i in eta0:\r\n\r\n        # Solve the system for different initial conditions\r\n        initial_conditions = [0,i]\r\n        final_conditions = odeint(system,initial_conditions, t, args=(Omega_sq,) )\r\n\r\n        ax.scatter((final_conditions[:,0]+np.pi)%(2*np.pi) - np.pi, final_conditions[:,1], s=3)\r\n        \r\n        ax.plot(Phi_sep, eta_sep, 'k--', lw=2)\r\n        ax.plot(Phi_sep, -eta_sep, 'k--', lw=2)\r\n        \r\n        ax.set_xticks([-np.pi, -np.pi\/2, 0, np.pi\/2, np.pi])\r\n        ax.set_xticklabels([\"$-\\pi$\", \"$-\\pi\/2$\", \"0\", \"$\\pi\/2$\", \"$\\pi$\"], fontsize=16)\r\n        \r\n        ax.set_xlabel(\"$\\\\Phi$\",fontsize=26)\r\n        ax.set_ylabel(\"$\\\\eta$\",fontsize=26)\r\n        ax.tick_params(axis='y', labelsize=30)\r\n        ax.tick_params(axis='x', labelsize=30)\r\nplt.show()\r\n<\/pre>\n
\"\"<\/p>\n
Phase-space trajectories<\/h2>\n
See Fig. 8.9.<\/p>\n
Phi = np.linspace(-np.pi,np.pi,65536)\r\nOmega_sq=1\r\nline_colour = [\"k\", \"b\", \"r\", \"g\", \"m\"]\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):  \r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax = plt.subplots(figsize=(10.2,10))\r\n    for i, C in enumerate([0,0.5,1,1.5,2]):\r\n        eta = np.sqrt(2*Omega_sq*(np.cos(Phi)+C))\r\n\r\n        ax.plot(Phi, eta, color = line_colour[i], lw=3)\r\n        ax.plot(Phi, -eta, color = line_colour[i], lw=3, label=\"C = \"+str(C))\r\n\r\n    ax.set_xticks([-np.pi, -np.pi\/2, 0, np.pi\/2, np.pi])\r\n    ax.set_xticklabels([\"$-\\pi$\", \"$-\\pi\/2$\", \"0\", \"$\\pi\/2$\", \"$\\pi$\"], fontsize=16)\r\n\r\n    ax.set_xlabel(\"$\\\\Phi$\",fontsize=26)\r\n    ax.set_ylabel(\"$\\\\eta$\",fontsize=26)\r\n    ax.tick_params(axis='y', labelsize=30)\r\n    ax.tick_params(axis='x', labelsize=30\r\n                  )\r\n\r\n    plt.legend(fontsize=26, frameon=True, framealpha=1, handlelength=1)\r\n\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
Plot of the gain function<\/h2>\n
See Fig. 8.10.<\/p>\n
det = np.linspace(-20,20,2000) #detuning parameter\r\ngain = (4\/det**3)*(1 - np.cos(det)- 0.5*det*np.sin(det))\r\nwith plt.style.context(('seaborn-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax = plt.subplots(figsize=(14,7.7))\r\n    \r\n    ax.plot(det, gain, lw=5)\r\n    \r\n    ax.set_xlabel(\"$\\\\tilde{D}$\",fontsize=26)\r\n    ax.set_ylabel(\"$g(\\\\tilde{D})$\",fontsize=26)\r\n    ax.tick_params(axis='y', labelsize=30)\r\n    ax.tick_params(axis='x', labelsize=30)\r\nplt.show()<\/pre>\n
\"\"<\/p>\n
Gain and line width: Madey’s theorem<\/h2>\n
See Fig. 8.11.<\/p>\n
det = np.linspace(-20,20,2000) #detuning parameter\r\ngain = (4\/det**3)*(1 - np.cos(det)- 0.5*det*np.sin(det))\r\nwith plt.style.context(('seaborn-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax = plt.subplots(figsize=(12.7,10))\r\n    \r\n    ax.plot(det, gain, 'k',lw=3, label = '$g(\\\\tilde{D})$')\r\n    ax.plot(det, (np.sin(det\/2) \/ (det\/2))**2, 'r', lw=3, label = 'sinc$^{2}(\\\\tilde{D}\/2)$+')\r\n    \r\n    ax.set_xlabel(\"$\\\\tilde{D}$\",fontsize=34)\r\n    ax.tick_params(axis='y', labelsize=34)\r\n    ax.tick_params(axis='x', labelsize=34)\r\n    \r\n    ax.legend(fontsize=36, loc=1, frameon=True, framealpha=1, handlelength=1)\r\n\r\n    plt.show()<\/pre>\n
\"\"<\/p>\n
Relativistic length contraction of electric field lines<\/h2>\n
See Fig. 8.12.<\/p>\n
First we define a function to calculate electric field components along two Cartesian coordinates.<\/p>\n
def E(q,r_0, x,y, gamma):\r\n\r\n    ''' Calculate components of Coulomb field'''    \r\n    eps_0 = 8.8541878188e-12\r\n    bet = np.sqrt(1-gamma**(-2))\r\n    r = np.hypot(x-r_0[0], y-r_0[1])\r\n    den = gamma**2 * r**3* (1-bet**2 * (y\/r)**2)**1.5\r\n    \r\n    E_x = q * (x - r_0[0]) \/ ((4*np.pi*eps_0)*den)\r\n    E_y = q * (y - r_0[1]) \/ ((4*np.pi*eps_0)*den)\r\n    \r\n    return E_x, E_y<\/pre>\n
Figure 8.12(b) can be generated using the code below. To generate Fig. 8.12(a), set the value of  \\gamma <\/span> to 1.<\/p>\n
# Define coordinates on a grid\r\nnx, ny = 128, 128\r\nlimit = 10\r\nx = np.linspace(-limit, limit, nx)\r\ny = np.linspace(-limit, limit, ny)\r\nX, Y = np.meshgrid(x, y)\r\n\r\n# Generate field\r\nq = 1e-9\r\nposition = [0,0]\r\ngamma=10\r\nEx, Ey = E(q, position,x=X,y=Y, gamma=gamma)\r\n\r\nwith plt.style.context(('default')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig = plt.figure(figsize=(10,10))\r\n    ax = fig.add_subplot(111)\r\n\r\n    # First generate the contour plot    \r\n    cp = ax.contourf(X, Y, np.log10(np.hypot(Ex, Ey)), levels=20, cmap=plt.cm.viridis)\r\n\r\n    # Then add field lines\r\n    for theta in np.linspace(0, 2*np.pi,30, endpoint=False):\r\n        ax.arrow(0, 0, 10*np.cos(theta)\/gamma, 10*np.sin(theta), head_width=0.2, length_includes_head=True, fc='k')\r\n\r\n    # Add a dot at the position of the charge\r\n    ax.add_artist(Circle(position, 0.02*limit, facecolor='red', edgecolor='black', zorder=10))\r\n\r\n    # Set limits and aspect ratio\r\n    ax.set_xlim(-limit,limit)\r\n    ax.set_ylim(-limit,limit)\r\n    ax.set_aspect('equal')\r\n\r\n    # Add a colorbar and scale its size to match the chart size\r\n    cbar = fig.colorbar(cp, ax=ax, fraction=0.046, pad=0.04)\r\n\r\n    cbar_ticks = np.arange(np.round(np.log10(np.hypot(Ex, Ey)).min()), np.round(np.log10(np.hypot(Ex, Ey)).max()),1)\r\n\r\n    cbar.ax.set_yticks(cbar_ticks)\r\n    cbar.ax.tick_params(axis='y', labelsize=34)\r\n\r\n    ax.set_xlabel('Parallel', fontsize=40)\r\n    ax.set_ylabel('Perpendicular', fontsize=40)\r\n    cbar.ax.set_ylabel('Log Field Amplitude', fontsize=40)\r\n\r\n    ax.set_xticks([-10,-5,0,5,10])\r\n    ax.set_yticks([-10,-5,0,5,10])\r\n    ax.tick_params(axis='x', labelsize=36)\r\n    ax.tick_params(axis='y', labelsize=36)\r\n    plt.show()<\/pre>\n
\"\"<\/p>\n
High-gain simulations<\/h2>\n
The code to generate the high-gain simulation results in Figs. 8.14, 8.15, and 8.19 is longer than what it would fit in a single snippet. The entire Python script can be downloaded here<\/a>.<\/p>\n
 <\/p>\n
Beam emittance<\/h2>\n
The following code generates the distribution in Fig. 8.18(b).<\/p>\n
# Generate a phase-sapce distribution of representative points\r\nn_p = 10000\r\nmu_x, sigma_x = 0, 0.1 # mean and standard deviation of position\r\nmu_xp, sigma_xp = 0, 0.01 # mean and standard deviation of angle\r\n\r\nsx1 = np.random.normal(mu_x, sigma_x, n_p)\r\nsxp = np.random.normal(mu_xp, sigma_xp, n_p)\r\nsx2 = sx1 + 10*np.tan(sxp)\r\n\r\ndef scatter_hist(x, y, ax, ax_histx, ax_histy):\r\n    \r\n    ax_histx.tick_params(axis=\"x\", labelbottom=False)\r\n    ax_histy.tick_params(axis=\"y\", labelleft=False)\r\n\r\n    # the scatter plot:\r\n    ax.scatter(x, y, c='darkslateblue', edgecolors='k', alpha=0.15)\r\n    \r\n    # now determine nice limits by hand:\r\n    nbins = 100\r\n    xmax = np.max(np.abs(x))\r\n    ymax = np.max(np.abs(y))\r\n    binwidthx = 2*xmax\/nbins\r\n    binwidthy = 2*ymax\/nbins              \r\n\r\n    binsx = np.arange(-xmax, xmax + binwidthx, binwidthx)\r\n    binsy = np.arange(-ymax, ymax + binwidthy, binwidthy)\r\n    ax_histx.hist(x, bins=binsx, color='darkslateblue')\r\n    ax_histy.hist(y, bins=binsy, color='darkslateblue', orientation='horizontal')\r\n\r\n    ax_histx.set_ylabel('Number of particles')\r\n    ax_histx.yaxis.label.set_fontsize(36)\r\n    ax_histx.tick_params(axis='y', labelsize=36)\r\n\r\n    ax_histy.set_xlabel('Number of particles')\r\n    ax_histy.xaxis.label.set_fontsize(36)\r\n    ax_histy.tick_params(axis='x', labelsize=36)\r\n    \r\n    ax.set_xlabel('Position (mm)')\r\n    ax.set_ylabel('Angle (rad)')\r\n    ax.xaxis.label.set_fontsize(36)\r\n    ax.tick_params(axis='x', labelsize=36)\r\n    ax.set_yticks([-0.03,0,0.03], [\"-0.03\",\"0.00\",\"0.03\"])\r\n    # ax.set_xticks([-0.3,0,0.3], [\"-0.3\",\"0.0\",\"0.3\"])\r\n    ax.yaxis.label.set_fontsize(36)\r\n    ax.tick_params(axis='y', labelsize=36)\r\n    \r\n# definitions for the axes\r\nleft, width = 0.1, 0.55\r\nbottom, height = 0.1, 0.55\r\nspacing = 0.02\r\n\r\nrect_scatter = [left, bottom, width, height]\r\nrect_histx = [left, bottom + height + spacing, width, 0.2]\r\nrect_histy = [left + width + spacing, bottom, 0.2, height]\r\n\r\nwith plt.style.context(('seaborn-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig = plt.figure(figsize=(12, 12))\r\n\r\n    ax = fig.add_axes(rect_scatter)\r\n    ax_histx = fig.add_axes(rect_histx, sharex=ax)\r\n    ax_histy = fig.add_axes(rect_histy, sharey=ax)\r\n\r\n    # use the previously defined function\r\n    scatter_hist(sx1, sxp, ax, ax_histx, ax_histy)\r\n    \r\n    t = np.linspace(0, 2*np.pi, 100)\r\n    ax.plot(sigma_x*np.cos(t) , sigma_xp*np.sin(t), color='darkred', lw=2, ls = '--')\r\n    ax.plot(3*sigma_x*np.cos(t) , 3*sigma_xp*np.sin(t), color='darkred', lw=2, ls = '--')\r\n\r\n    plt.show()<\/pre>\n
\"\"<\/p>\n
The following code generates the distribution in Fig. 8.18(c).<\/p>\n
# Generate the same representative points after propagation\r\n\r\n# definitions for the axes\r\nleft, width = 0.1, 0.55\r\nbottom, height = 0.1, 0.55\r\nspacing = 0.02\r\n\r\nrect_scatter = [left, bottom, width, height]\r\nrect_histx = [left, bottom + height + spacing, width, 0.2]\r\nrect_histy = [left + width + spacing, bottom, 0.2, height]\r\n\r\nwith plt.style.context(('seaborn-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig = plt.figure(figsize=(12, 12))\r\n\r\n    ax = fig.add_axes(rect_scatter)\r\n    ax_histx = fig.add_axes(rect_histx, sharex=ax)\r\n    ax_histy = fig.add_axes(rect_histy, sharey=ax)\r\n\r\n    # use the previously defined function\r\n    scatter_hist(sx2, sxp, ax, ax_histx, ax_histy)\r\n    \r\n    t = np.linspace(0, 2*np.pi, 100)\r\n    ax.plot(sigma_x*np.cos(t) + 10*sigma_xp*np.sin(t) , sigma_xp*np.sin(t), color='darkred', lw=2, ls = '--')\r\n    ax.plot(3*sigma_x*np.cos(t) + 30*sigma_xp*np.sin(t) , 3*sigma_xp*np.sin(t), color='darkred', lw=2, ls = '--')\r\n\r\n    plt.show()<\/pre>\n
\"\"<\/p>\n
Bakers’ transformation<\/h2>\n
Figure 8.21 is generated using this Python script<\/a>.<\/p>\n
 <\/p>\n
Separatrix<\/h2>\n
See Fig. 8.22.<\/p>\n
t = np.linspace(0,10.1,5000)\r\neta0 = [0.2]\r\nOmega_sq=1\r\n\r\nPhi_sep = np.linspace(-np.pi,np.pi,256)\r\neta_sep = np.sqrt(2*Omega_sq*(np.cos(Phi_sep)+1))\r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax = plt.subplots(figsize=(8,8))\r\n    for i in eta0:\r\n\r\n        ax.plot(Phi_sep, eta_sep, 'k--', lw=4)\r\n        ax.plot(Phi_sep, -eta_sep, 'k--', lw=4)\r\n        \r\n        ax.set_xticks([-np.pi, -np.pi\/2, 0, np.pi\/2, np.pi])\r\n        ax.set_xticklabels([\"$-\\pi$\", \"$-\\pi\/2$\", \"0\", \"$\\pi\/2$\", \"$\\pi$\"], fontsize=22)\r\n        \r\n        ax.set_xlabel(\"$\\\\Phi$\",fontsize=34)\r\n        ax.set_ylabel(\"$\\\\eta$\",fontsize=34)\r\n        ax.tick_params(axis='y', labelsize=36)\r\n        ax.tick_params(axis='x', labelsize=36)<\/pre>\n
\"\"<\/p>\n
Solution of eqn (8.57)<\/h2>\n
See Fig. 8.24.<\/p>\n
def pendulum(initial_conditions, t, Omega_sq):\r\n    # Set the initial conditions \r\n    Phi = initial_conditions[0]\r\n    eta = initial_conditions[1]\r\n    \r\n    # Write the coupled equations\r\n    d_Phi = eta\r\n    d_eta = -Omega_sq *np.sin(Phi)\r\n    \r\n    d_next_step = [d_Phi, d_eta]\r\n    \r\n    return d_next_step\r\n\r\nmu_0 = 1.25663706212e-6\r\nc = 3e8\r\ngamma_r = 1000\r\nK = 0.1\r\nl_u = 0.1 \r\nk_u = 2*np.pi\/l_u\r\nee = 1.602176634e-19\r\nme = 9.1093837015e-31\r\nlambda_L = 1e-9\r\nomega_L = 2*np.pi*c\/lambda_L\r\nE_0 = 5e8\r\n\r\nos = (ee*k_u*K*E_0)\/(2*me*c**2*gamma_r**2)\r\n\r\n\r\nt = np.linspace(0,555.5,5000000)\r\neps = 1.e-6\r\neta0 = [(1-eps)*2*np.sqrt(os)]\r\nprint(eta0)\r\nOmega_sq=os\r\n\r\ninitial_conditions = [3.138,0]\r\nfinal_conditions = odeint(pendulum,initial_conditions, t, args=(Omega_sq,) )\r\n\r\nPhi_sep = np.linspace(-np.pi,np.pi,256)\r\neta_sep = np.sqrt(2*Omega_sq*(np.cos(Phi_sep)+1))\r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax = plt.subplots(1,2, figsize=(16,7.8))\r\n    # for i in eta0:\r\n    #initial_conditions = [0,i]\r\n\r\n    ax[0].scatter((final_conditions[:,0]+np.pi)%(2*np.pi) - np.pi, final_conditions[:,1], s=3)\r\n    ax[0].plot(Phi_sep, eta_sep, 'k--', lw=2)\r\n    ax[0].plot(Phi_sep, -eta_sep, 'k--', lw=2)\r\n    \r\n    ax[1].scatter((final_conditions[:,0]+np.pi)%(2*np.pi) - np.pi, final_conditions[:,1], s=1)\r\n    \r\n    ax[1].plot(Phi_sep, eta_sep, 'k--', lw=2)\r\n    ax[1].plot(Phi_sep, -eta_sep, 'k--', lw=2)\r\n    ax[1].set_xlim(3.13,3.142)\r\n    ax[1].set_ylim(-0.0005,0.0005)\r\n    \r\n    ax[0].set_xticks([-np.pi, -np.pi\/2, 0, np.pi\/2, np.pi])\r\n    ax[0].set_xticklabels([\"$-\\pi$\", \"$-\\pi\/2$\", \"0\", \"$\\pi\/2$\", \"$\\pi$\"], fontsize=16)\r\n\r\n    ax[0].set_xlabel(\"$\\\\Phi$\",fontsize=30)\r\n    ax[0].set_ylabel(\"$\\\\eta$\",fontsize=30)\r\n    ax[0].tick_params(axis='y', labelsize=30)\r\n    ax[0].tick_params(axis='x', labelsize=30)\r\n    \r\n    ax[1].set_yticks([-0.00025, 0, 0.00025])\r\n    ax[1].set_yticklabels([\"-2.5\", \"0\", \"2.5\"], fontsize=30)\r\n    ax[1].set_xlabel(\"$\\\\Phi$\",fontsize=30)\r\n    ax[1].set_ylabel(\"$\\\\eta \\\\, (\\\\times 10^{-4})$\",fontsize=30)\r\n    ax[1].tick_params(axis='y', labelsize=30)\r\n    ax[1].tick_params(axis='x', labelsize=30)\r\n    \r\n    plt.tight_layout()<\/pre>\n
\"\"<\/p>\n
Laminar filamentation<\/h2>\n
See Fig. 8.23. The constants required in this code block are defined in the previous block.<\/p>\n
eps=1e-5\r\n\r\nt = np.linspace(0,50000.1,5000)\r\nn_p = 1000\r\neta_0 = 0.01*eps*np.random.randn(n_p) + (1-eps)*2*np.sqrt(os)\r\nphi_0 = 0.01*eps*np.random.randn(n_p)\r\n\r\neta = np.zeros((t.shape[0], n_p))\r\nphi = np.zeros((t.shape[0], n_p))\r\nfor i,j in enumerate(eta_0):\r\n    initial_conditions = [phi_0[i],j]\r\n    final_conditions = odeint(pendulum,initial_conditions, t, args=(os,) )\r\n    \r\n    phi[:,i] = final_conditions[:,0]\r\n    eta[:,i] = final_conditions[:,1] \r\n\r\nPhi_sep = np.linspace(-np.pi,np.pi,256)\r\neta_sep = np.sqrt(2*Omega_sq*(np.cos(Phi_sep)+1))\r\n\r\ntime_intervals = np.linspace(0,500000,65)#[:-2]\r\ncolour = 'red' #'darkslateblue'\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, axs = plt.subplots(4,1, figsize=(4,10))\r\n    axs[0].scatter(1e6*phi[0,:], eta[0,:]-eta[0,:].mean(), s=20, c=colour, ec='k')\r\n    axs[1].scatter(phi[0,:], eta[0,:], s=20, c=colour, ec='k')\r\n    axs[2].scatter(phi[2500,:], eta[2500,:], s=20, c=colour, ec='k')\r\n    axs[3].scatter(phi[-1,:], eta[-1:], s=20, c=colour, ec='k')\r\n    for idx, ax in enumerate(axs.flat):\r\n\r\n        if idx !=0:\r\n            ax.plot(Phi_sep, eta_sep, 'k--', lw=2)\r\n            ax.plot(Phi_sep, -eta_sep, 'k--', lw=2)    \r\n        \r\n        if idx ==0:\r\n            ax.set_xlabel(\"$\\\\Phi \\\\, (\\\\times 10^{-7})$\",fontsize=20)\r\n            ax.set_ylabel(\"$\\\\eta - \\\\langle \\\\eta \\\\rangle \\\\, (\\\\times 10^{-7})$\",fontsize=20)\r\n            ax.tick_params(axis='y', labelsize=22)\r\n            ax.tick_params(axis='x', labelsize=22)\r\n\r\n        else:    \r\n            ax.set_xlabel(\"$\\\\Phi$\",fontsize=20)\r\n            ax.set_ylabel(\"$\\\\eta$\",fontsize=20)\r\n            ax.tick_params(axis='y', labelsize=22)\r\n            ax.tick_params(axis='x', labelsize=22)\r\n        \r\nplt.tight_layout()<\/pre>\n
\"\"<\/h1>\n
Stochastic layer<\/h2>\n
See Fig. 8.25.<\/p>\n
def pendulum_perturbation(initial_conditions, t, Omega_sq, epsilon, omega_0):\r\n    # Set the initial conditions \r\n    Phi = initial_conditions[0]\r\n    eta = initial_conditions[1]\r\n    \r\n    # Write the coupled equations\r\n    d_Phi = eta\r\n    d_eta = -Omega_sq *np.sin(Phi) + epsilon*Omega_sq *np.sin(omega_0*t)\r\n    \r\n    d_next_step = [d_Phi, d_eta]\r\n    \r\n    return d_next_step\r\n\r\nmu_0 = 1.25663706212e-6\r\nc = 3e8\r\ngamma_r = 1000\r\nK = 0.1\r\nl_u = 0.1 \r\nk_u = 2*np.pi\/l_u\r\nee = 1.602176634e-19\r\nme = 9.1093837015e-31\r\nlambda_L = 1e-9\r\nomega_L = 2*np.pi*c\/lambda_L\r\nE_0 = 5e8\r\n\r\nos = (ee*k_u*K*E_0)\/(2*me*c**2*gamma_r**2)\r\nomega_0=0.9*os\r\n\r\nt = np.linspace(0,50000,5000000)\r\neps = 2.e-3\r\neta0 = [(1-eps)*2*np.sqrt(os), 0.5*(1-eps)*2*np.sqrt(os), 1.5*(1-eps)*2*np.sqrt(os)]#[0.11081005]  \r\nprint(eta0)\r\nOmega_sq=os\r\nepsilon=0.01\r\n\r\nPhi_sep = np.linspace(-np.pi,np.pi,256)\r\neta_sep = np.sqrt(2*Omega_sq*(np.cos(Phi_sep)+1))\r\n\r\nwith plt.style.context(('seaborn-v0_8-whitegrid')):\r\n    plt.rcParams.update(synchrotron_style)\r\n    fig, ax = plt.subplots(1,2, figsize=(16,8))\r\n    for i,j in enumerate(eta0):\r\n        \r\n        initial_conditions = [0,j]\r\n        final_conditions = odeint(pendulum_perturbation,initial_conditions, t, args=(Omega_sq,epsilon, omega_0) )\r\n\r\n        ax[0].scatter((final_conditions[:,0]+np.pi)%(2*np.pi) - np.pi, final_conditions[:,1], s=2, alpha=0.25)\r\n        ax[0].plot(Phi_sep, eta_sep, 'k--', lw=2)\r\n        ax[0].plot(Phi_sep, -eta_sep, 'k--', lw=2)\r\n        \r\n        ax[1].scatter((final_conditions[:,0]+np.pi)%(2*np.pi) - np.pi, final_conditions[:,1], s=2, alpha=0.25)\r\n        \r\n        ax[1].plot(Phi_sep, eta_sep, 'k--', lw=2)\r\n        ax[1].plot(Phi_sep, -eta_sep, 'k--', lw=2)\r\n        ax[1].set_xlim(3,3.15)\r\n        ax[1].set_ylim(-0.05,0.05)\r\n        \r\n        ax[0].set_xticks([-np.pi, -np.pi\/2, 0, np.pi\/2, np.pi])\r\n        ax[0].set_xticklabels([\"$-\\pi$\", \"$-\\pi\/2$\", \"0\", \"$\\pi\/2$\", \"$\\pi$\"], fontsize=16)\r\n\r\n        ax[0].set_xlabel(\"$\\\\Phi$\",fontsize=26)\r\n        ax[0].set_ylabel(\"$\\\\eta$\",fontsize=26)\r\n        ax[0].tick_params(axis='y', labelsize=28)\r\n        ax[0].tick_params(axis='x', labelsize=28)\r\n        \r\n        ax[1].set_xlabel(\"$\\\\Phi$\",fontsize=26)\r\n        ax[1].set_ylabel(\"$\\\\eta \\\\, (\\\\times 10^{-4})$\",fontsize=26)\r\n        ax[1].tick_params(axis='y', labelsize=28)\r\n        ax[1].tick_params(axis='x', labelsize=28)\r\n        \r\n        plt.tight_layout()<\/pre>\n
\"\"<\/p>\n
Turbulent chaotic evolution<\/h2>\n
See Fig. 8.26. The function “pendulum_perturbation” using in this block, has been defined in the previous block of code.<\/p>\n
t = np.linspace(0,50000,1000000)\r\nepsilon=0.01\r\neta0 = [(1-eps)*2*np.sqrt(os), 0.9999*(1-eps)*2*np.sqrt(os), 0.5*(1-eps)*2*np.sqrt(os)]\r\ncolor = ['#1f77b4', '#d62728', '#ff7f0e']\r\n\r\nomega_0=0.9*os\r\nwith plt.style.context(('seaborn-whitegrid')):\r\n    fig, ax = plt.subplots(figsize=(16,10))\r\n    for i,j in enumerate(eta0):\r\n        initial_conditions = [0,j]\r\n        final_conditions = odeint(pendulum_perturbation,initial_conditions, t, args=(Omega_sq,epsilon, omega_0) )\r\n\r\n        ax.plot(t[:],final_conditions[:,0], lw=3, color=color[i])\r\n    ax.set_xlabel(\"$t$ (s)\",fontsize=26)\r\n    ax.set_ylabel(\"$\\\\Phi (t)$\",fontsize=26)\r\n    ax.tick_params(axis='y', labelsize=20)\r\n    ax.tick_params(axis='x', labelsize=20)\r\nplt.show()<\/pre>\n
\"\"<\/p>\n","protected":false},"excerpt":{"rendered":"
Below is a set of python codes associated with Chapter 8 of Daniele Pelliccia and David M. Paganin, “Synchrotron Light: A Physics Journey from Laboratory to Cosmos” (Oxford University Press, 2025). In order to run any of these python codes, you will need to include the following header. import numpy as np import matplotlib.pyplot as […]<\/p>\n","protected":false},"author":1,"featured_media":0,"parent":79,"menu_order":0,"comment_status":"closed","ping_status":"closed","template":"","meta":{"_et_pb_use_builder":"","_et_pb_old_content":"","_et_gb_content_width":"","footnotes":""},"class_list":["post-695","page","type-page","status-publish","hentry"],"_links":{"self":[{"href":"https:\/\/synchrotron-light.net\/wp\/index.php?rest_route=\/wp\/v2\/pages\/695","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/synchrotron-light.net\/wp\/index.php?rest_route=\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/synchrotron-light.net\/wp\/index.php?rest_route=\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/synchrotron-light.net\/wp\/index.php?rest_route=\/wp\/v2\/users\/1"}],"replies":[{"embeddable":true,"href":"https:\/\/synchrotron-light.net\/wp\/index.php?rest_route=%2Fwp%2Fv2%2Fcomments&post=695"}],"version-history":[{"count":14,"href":"https:\/\/synchrotron-light.net\/wp\/index.php?rest_route=\/wp\/v2\/pages\/695\/revisions"}],"predecessor-version":[{"id":788,"href":"https:\/\/synchrotron-light.net\/wp\/index.php?rest_route=\/wp\/v2\/pages\/695\/revisions\/788"}],"up":[{"embeddable":true,"href":"https:\/\/synchrotron-light.net\/wp\/index.php?rest_route=\/wp\/v2\/pages\/79"}],"wp:attachment":[{"href":"https:\/\/synchrotron-light.net\/wp\/index.php?rest_route=%2Fwp%2Fv2%2Fmedia&parent=695"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}