- {
- “cells”: [
- {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“# Common PLEQUE tasksn”, “n”, “In this notebook, we cover the most common tasks involving a magnetic equilibrium. These are:n”, “n”, “- mapping along the magnetic surfacesn”, “- field line tracingn”, “- equilibrium visualisation using contour plotsn”, “- individual flux surface examinationn”, “- locating the separatrix position in a given profilen”, “- detector line of sight visualisation”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“%pylab inlinen”, “n”, “from pleque.io.readers import read_geqdskn”, “from pleque.utils.plotting import *n”, “from pleque.tests.utils import get_test_equilibria_filenames”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Load a testing equilibriumn”, “Several test equilibria come shipped with PLEQUE. Their location is:”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“gfiles = get_test_equilibria_filenames()n”, “gfiles”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“We store one of the text equilibria in the variable eq, an instance of the Equilibrium class.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“test_case_number = 5n”, “n”, “#Load equilibrium stored in the EQDSK formatn”, “eq = read_geqdsk(gfiles[test_case_number])n”, “n”, “#Plot basic overview of the equilibriumn”, “plt.figure()n”, “eq._plot_overview()n”, “n”, “#Plot X-pointsn”, “plot_extremes(eq, markeredgewidth=2)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Mapping along the magnetic surfacesn”, “n”, “First, a general plotting function plot_2d is defined.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“def plot_2d(R, Z, data, *args, title=None):n”, “ n”, “ #Define X and Y axis limits based on the vessel sizen”, “ rlim = [np.min(eq.first_wall.R), np.max(eq.first_wall.R)]n”, “ zlim = [np.min(eq.first_wall.Z), np.max(eq.first_wall.Z)]n”, “ size = rlim[1] - rlim[0]n”, “ rlim[0] -= size / 12n”, “ rlim[1] += size / 12n”, “ size = zlim[1] - zlim[0]n”, “ zlim[0] -= size / 12n”, “ zlim[1] += size / 12n”, “ n”, “ #Set up the figure: set axis limits, draw LCFS and first wall, write labelsn”, “ ax = plt.gca()n”, “ ax.set_xlim(rlim)n”, “ ax.set_ylim(zlim)n”, “ ax.plot(eq.lcfs.R, eq.lcfs.Z, color=’k’, ls=’–’, lw=2)n”, “ ax.plot(eq.first_wall.R, eq.first_wall.Z, ‘k-’, lw=2)n”, “ ax.set_xlabel(‘R [m]’)n”, “ ax.set_ylabel(‘Z [m]’)n”, “ ax.set_aspect(‘equal’)n”, “ if title is not None:n”, “ ax.set_title(title)n”, “ n”, “ #Finally, plot the desired quantityn”, “ cl = ax.contour(R, Z, data, *args)n”, “ n”, “ return cl”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Now we set up an $[R,Z]$ grid where these quantities are evaluated and plot the quantities.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“#Create an [R,Z] grid 200 by 300 pointsn”, “grid = eq.grid((200,300), dim=’size’)n”, “n”, “#Plot the poloidal flux and toroidal fluxn”, “plt.figure(figsize=(16,4))n”, “plt.subplot(131)n”, “plot_2d(grid.R, grid.Z, grid.psi, 20, title=r’$\psi$’)n”, “plt.subplot(132)n”, “plot_2d(grid.R, grid.Z, eq.tor_flux(grid), 100, title=’toroidal flux’)n”, “n”, “#Plot the poloidal magnetic field, toroidal magnetic field and the total magnetic fieldn”, “plt.figure(figsize=(16,4))n”, “plt.subplot(131)n”, “cl = plot_2d(grid.R, grid.Z, eq.B_pol(grid), 20, title=r’$B_\mathrm{p}$ [T]’)n”, “plt.colorbar(cl)n”, “plt.subplot(132)n”, “cl = plot_2d(grid.R, grid.Z, eq.B_tor(grid), 20, title=r’$B_\mathrm{t}$ [T]’)n”, “plt.colorbar(cl)n”, “plt.subplot(133)n”, “cl = plot_2d(grid.R, grid.Z, eq.B_abs(grid), 20, title=r’$|B|$ [T]’)n”, “plt.colorbar(cl)n”, “n”, “#Plot the total pressure, toroidal current density and poloidal current densityn”, “plt.figure(figsize=(16,4))n”, “plt.subplot(131)n”, “cl = plot_2d(grid.R, grid.Z, eq.pressure(grid)/1e3, np.linspace(0, 30, 21), title=r’$p$ [kPa]’)n”, “plt.colorbar(cl)n”, “plt.subplot(132)n”, “plot_2d(grid.R, grid.Z, eq.j_tor(grid), np.linspace(-5e6, 5e6, 30), title=r’$j_\phi$’)n”, “plt.subplot(133)n”, “plot_2d(grid.R, grid.Z, eq.j_pol(grid), np.linspace(0, 3e5, 21), title=r’$j_\theta$’)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Exploring flux surface propertiesn”, “n”, “With the eq._flux_surface(psi_n) function, one may study individual flux surfaces. In this section, we plot the $\psi_N=0.8$ flux surface and calculate its safety factor $q$, length in the poloidal direction, total 3D area, volume and toroidal current density.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“#Define the flux surface by its normalised poloidal fluxn”, “surf = eq._flux_surface(psi_n=0.8)[0]n”, “n”, “#Plot the flux surfacen”, “plt.figure()n”, “ax = gca()n”, “ax.plot(eq.lcfs.R, eq.lcfs.Z, color=’k’, ls=’–’, lw=2)n”, “ax.plot(eq.first_wall.R, eq.first_wall.Z, ‘k-’, lw=2)n”, “surf.plot(ls=’–‘)n”, “ax.set_xlabel(‘R [m]’)n”, “ax.set_ylabel(‘Z [m]’)n”, “ax.set_aspect(‘equal’)n”, “n”, “#Calculate several flux surface quantitiesn”, “print(‘Safety factor: %.2f’ % surf.eval_q[0])n”, “print(‘Length: %.2f m’ % surf.length)n”, “print(‘Area: %.4f m^2’ % surf.area)n”, “print(‘Volume: %.3f m^3’ % surf.volume)n”, “print(‘Toroidal current density: %.3f MA/m^2’ % (surf.tor_current/1e6))”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Profile mappingn”, “In experiment one often encounters the need to compare profiles which were measured at various locations in the tokamak. In this section, we show how such a profile may be mapped onto an arbitrary location and to the outer midplane.n”, “n”, “Our example profile is measured at the plasma top (plotted in red by the following code). It will be mapped to a chord located on the HFS (in violet) and to the outer midplane (not shown in the picture). Note that the outer midplane is defined by the O-point height, not with regard to the chamber ($Z=0$ etc.).”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Define the chord along which the profile was measured (in red)n”, “N = 200 #number of datapoints in the profilen”, “chord = eq.coordinates(R=0.6*np.ones(N), Z=np.linspace(0.3, 0., N))n”, “n”, “# Define the HFS chord where we wish to map the profile (in violet)n”, “chord_hfs = eq.coordinates(R=np.linspace(0.35, 0.6, 20), Z=-0.1*np.ones(20))n”, “n”, “# Plot both the chordsn”, “plt.figure()n”, “eq._plot_overview()n”, “chord.plot(lw=3, ls=’–’, color=’C3’, label=’measurement location’)n”, “chord_hfs.plot(lw=3, ls=’–’, color=’C4’, label=’HFS chord’)n”, “plt.legend(loc=3)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“The profile shape is defined somewhat arbitrarily using the error function erf.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“from scipy.special import erfn”, “n”, “# Define the profile valuesn”, “prof_func = lambda x, k1, xsep: k1/4 * (1 + erf((x-xsep)*20))*np.log((x+1)*1.2) - 4*np.exp(-(50*(x-1)**2))n”, “profile = prof_func(1 - chord.psi_n, 10, 0.15)n”, “n”, “# Plot the profile along the chord where it was measuredn”, “plt.figure()n”, “plt.plot(chord.Z, profile, color=’C3’)n”, “plt.xlabel(‘Z [m]’)n”, “plt.ylabel(‘profile value [a.u.]’)n”, “plt.tight_layout()”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“To begin the mapping, the profile is converted into a flux function by eq.fluxfuncs.add_flux_func(). The flux function is a spline, and therefore it can be evaluated at any $\psi_N$ coordinate covered by the original chord. This will allow its mapping to any other coordinate along the flux surfaces.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“eq.fluxfuncs.add_flux_func(‘test_profile’, profile, chord, spline_smooth=0)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“To evaluate the flux function along a chord, simply pass the chord (an instance of the Coordinates class) to the flux function. In the next figure the profile is mapped to the HFS cord.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Map the profile to the HFS cordn”, “plt.figure()n”, “plt.plot(chord_hfs.R, eq.fluxfuncs.test_profile(chord_hfs), ‘–’, color=’C4’)n”, “plt.xlabel(‘R [m]’)n”, “plt.ylabel(‘profile value [a.u.]’)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“For the outer midplane, no special chord need be specified. Every instance of the Coordinates class can automatically map its coordinates to the outer midplane. (Note that this doesn’t require a flux function to be specified. The conversion is performed in the coordinates only.)”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“#Map the profile to the outer midplanen”, “plt.figure()n”, “plt.plot(chord.r_mid, profile, color=’C1’)n”, “plt.xlabel(r’$R$ [m]’)n”, “plt.ylabel(‘profile value [a.u.]’)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Finally, the profile may be drawn along the entire poloidal cross section.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Assuming poloidal symmetry, plot the profile in the poloidal cross sectionn”, “plt.figure()n”, “ax = gca()n”, “ax.plot(eq.lcfs.R, eq.lcfs.Z, color=’k’, ls=’–’, lw=2)n”, “ax.plot(eq.first_wall.R, eq.first_wall.Z, ‘k-’, lw=2)n”, “grid = eq.grid()n”, “ax.pcolormesh(grid.R, grid.Z, eq.fluxfuncs.test_profile(grid))n”, “ax.set_xlabel(‘R [m]’)n”, “ax.set_ylabel(‘Z [m]’)n”, “ax.set_aspect(‘equal’)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Field line tracingn”, “n”, “Another task most commonly performed in edge plasma physics is field line tracing and length calculation. (In the core plasma, the field line length is defined as the parallel distance of one poloidal turn. In the SOL, it’s the so-called connection length.)n”, “n”, “To showcase field line tracing, first we define a set of five starting points, all located at the outer midplane ($Z=0$) with $R$ going from 0.55 m (core) to 0.76 m (SOL).”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Define the starting pointsn”, “N = 5n”, “Rs = np.linspace(0.57, 0.76, N, endpoint=True)n”, “Zs = np.zeros_like(Rs)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Next, the field lines beginning at these points are traced. The default tracing direction is direction=1, that is, following the direction of the toroidal magnetic field.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“traces = eq.trace_field_line(R=Rs, Z=Zs)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“To visualise the field lines, we plot them in top view, poloidal cross-section view and 3D view.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Define limiter as viewed from the topn”, “Ns = 100n”, “inner_lim = eq.coordinates(np.min(eq.first_wall.R)*np.ones(Ns), np.zeros(Ns), np.linspace(0, 2*np.pi, Ns))n”, “outer_lim = eq.coordinates(np.max(eq.first_wall.R)*np.ones(Ns), np.zeros(Ns), np.linspace(0, 2*np.pi, Ns))n”, “n”, “# Create a figuren”, “fig = plt.figure(figsize=(10,5))n”, “n”, “# Plot the top view of the field linesn”, “ax = plt.subplot(121)n”, “plt.plot(inner_lim.X, inner_lim.Y, ‘k-’, lw=4)n”, “plt.plot(outer_lim.X, outer_lim.Y, ‘k-’, lw=4)n”, “for fl in traces:n”, “ ax.plot(fl.X, fl.Y)n”, “ax.set_xlabel(‘$X$ [m]’)n”, “ax.set_ylabel(‘$Y$ [m]’)n”, “ax.set_aspect(‘equal’)n”, “n”, “# Plot the poloidal cross-section view of the field linesn”, “ax = plt.subplot(122)n”, “plt.plot(eq.first_wall.R, eq.first_wall.Z, ‘k-‘)n”, “plt.plot(eq.separatrix.R, eq.separatrix.Z, ‘C1–‘)n”, “for fl in traces:n”, “ plt.plot(fl.R, fl.Z)n”, “ax.set_xlabel(‘$R$ [m]’)n”, “ax.set_ylabel(‘$Z$ [m]’)n”, “ax.set_aspect(‘equal’)n”, “n”, “# Plot the 3D view of the field linesn”, “from mpl_toolkits.mplot3d import Axes3Dn”, “fig = plt.figure(figsize=(6,6))n”, “ax = fig.gca(projection=’3d’)n”, “for fl in traces:n”, “ ax.scatter(fl.X, fl.Y, fl.Z, s=0.3, marker=’.’)n”, “ax.set_xlabel(‘$X$ [m]’)n”, “ax.set_ylabel(‘$Y$ [m]’)n”, “ax.set_zlabel(‘$Z$ [m]’)n”, “#ax.set_aspect(‘equal’)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“The field line length may be calculated using its attribute length.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“traces[0].length”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“### Connection length profile in the SOLn”, “n”, “A subtask encountered e.g. in SOL transport analysis is to calculated the connection length $L$ of a number of SOL flux tubes. To this end, we define a couple more SOL field lines. Note that now the direction argument changes whether we trace to the HFS or LFS limiter/divertor. Also pay attention to the in_first_wall=True argument, which tells the field lines to terminate upon hitting the first wall. (Otherwise they would be terminated at the edge of a rectangle surrounding the vacuum vessel.)”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“Rsep = 0.7189 # You might want to change this when switching between different test equilibria.n”, “Rs_SOL = Rsep + 0.001*np.array([0, 0.2, 0.5, 0.7, 1, 1.5, 2.5, 4, 6, 9, 15, 20])n”, “Zs_SOL = np.zeros_like(Rs_SOL)n”, “n”, “SOL_traces = eq.trace_field_line(R=Rs_SOL, Z=Zs_SOL, direction=-1, in_first_wall=True)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Finally we calculate the connection length and plot its profile.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Calculate the connection lengthn”, “L_conn = np.array([SOL_traces[k].length for k in range(len(SOL_traces))])n”, “n”, “# Create a figuren”, “fig = plt.figure(figsize=(10,5))n”, “n”, “# Plot the poloidal cross-section view of the field linesn”, “ax = plt.subplot(121)n”, “ax.plot(eq.first_wall.R, eq.first_wall.Z, ‘k-‘)n”, “ax.plot(eq.separatrix.R, eq.separatrix.Z, ‘C1–‘)n”, “for fl in SOL_traces:n”, “ ax.plot(fl.R, fl.Z)n”, “ax.set_xlabel(‘R [m]’)n”, “ax.set_ylabel(‘Z [m]’)n”, “ax.set_aspect(‘equal’)n”, “n”, “# Plot the connection length profilen”, “ax = plt.subplot(122)n”, “ax.plot(Rs_SOL, L_conn, ‘ro’)n”, “ax.set_xlabel(‘R [m]’)n”, “ax.set_ylabel(‘L [m]’)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Equilibrium visualisation using contour plotsn”, “n”, “In this section PLEQUE is used to produce contour plots of the following quantities:n”, “n”, “- poloidal magnetic field flux $\psi$n”, “- toroidal magnetic field fluxn”, “- poloidal magnetic field $B_p$n”, “- toroidal magnetic field $B_t$n”, “- total magnetic field $|B|$n”, “- total pressure $p$n”, “- toroidal current density $j_\phi$n”, “- poloidal current density $j_\theta$n”, “n”, “First, a general plotting function plot_2d is defined.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“def plot_2d(R, Z, data, *args, title=None):n”, “ n”, “ # Define X and Y axis limits based on the vessel sizen”, “ rlim = [np.min(eq.first_wall.R), np.max(eq.first_wall.R)]n”, “ zlim = [np.min(eq.first_wall.Z), np.max(eq.first_wall.Z)]n”, “ size = rlim[1] - rlim[0]n”, “ rlim[0] -= size / 12n”, “ rlim[1] += size / 12n”, “ size = zlim[1] - zlim[0]n”, “ zlim[0] -= size / 12n”, “ zlim[1] += size / 12n”, “ n”, “ # Set up the figure: set axis limits, draw LCFS and first wall, write labelsn”, “ ax = plt.gca()n”, “ ax.set_xlim(rlim)n”, “ ax.set_ylim(zlim)n”, “ ax.plot(eq.lcfs.R, eq.lcfs.Z, color=’k’, ls=’–’, lw=2)n”, “ ax.plot(eq.first_wall.R, eq.first_wall.Z, ‘k-’, lw=2)n”, “ ax.set_xlabel(‘R [m]’)n”, “ ax.set_ylabel(‘Z [m]’)n”, “ ax.set_aspect(‘equal’)n”, “ if title is not None:n”, “ ax.set_title(title)n”, “ n”, “ # Finally, plot the desired quantityn”, “ cl = ax.contour(R, Z, data, *args)n”, “ n”, “ return cl”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Next we set up an $[R,Z]$ grid where these quantities are evaluated and plot the quantities.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Create an [R,Z] grid 200 by 300 pointsn”, “grid = eq.grid((200,300), dim=’size’)n”, “n”, “# Plot the poloidal flux and toroidal fluxn”, “plt.figure(figsize=(16,4))n”, “plt.subplot(131)n”, “plot_2d(grid.R, grid.Z, grid.psi, 20, title=r’$\psi$’)n”, “plt.subplot(132)n”, “plot_2d(grid.R, grid.Z, eq.tor_flux(grid), 100, title=’toroidal flux’)n”, “n”, “# Plot the poloidal magnetic field, toroidal magnetic field and the total magnetic fieldn”, “plt.figure(figsize=(16,4))n”, “plt.subplot(131)n”, “cl = plot_2d(grid.R, grid.Z, eq.B_pol(grid), 20, title=r’$B_\mathrm{p}$ [T]’)n”, “plt.colorbar(cl)n”, “plt.subplot(132)n”, “cl = plot_2d(grid.R, grid.Z, eq.B_tor(grid), 20, title=r’$B_\mathrm{t}$ [T]’)n”, “plt.colorbar(cl)n”, “plt.subplot(133)n”, “cl = plot_2d(grid.R, grid.Z, eq.B_abs(grid), 20, title=r’$|B|$ [T]’)n”, “plt.colorbar(cl)n”, “n”, “# Plot the total pressure, toroidal current density and poloidal current densityn”, “plt.figure(figsize=(16,4))n”, “plt.subplot(131)n”, “cl = plot_2d(grid.R, grid.Z, eq.pressure(grid)/1e3, np.linspace(0, 30, 21), title=r’$p$ [kPa]’)n”, “plt.colorbar(cl)n”, “plt.subplot(132)n”, “plot_2d(grid.R, grid.Z, eq.j_tor(grid), np.linspace(-5e6, 5e6, 30), title=r’$j_\phi$’)n”, “plt.subplot(133)n”, “plot_2d(grid.R, grid.Z, eq.j_pol(grid), np.linspace(0, 3e5, 21), title=r’$j_\theta$’)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Individual flux surface examinationn”, “n”, “With the eq._flux_surface(psi_n) function, one may study individual flux surfaces. In this section, we plot the $\psi_N=0.8$ flux surface and calculate its safety factor $q$, length in the poloidal direction, total 3D area, volume and toroidal current density.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Define the flux surface by its normalised poloidal fluxn”, “surf = eq._flux_surface(psi_n=0.8)[0]n”, “n”, “# Plot the flux surfacen”, “plt.figure()n”, “ax = gca()n”, “ax.plot(eq.lcfs.R, eq.lcfs.Z, color=’k’, ls=’–’, lw=2)n”, “ax.plot(eq.first_wall.R, eq.first_wall.Z, ‘k-’, lw=2)n”, “surf.plot(ls=’–‘)n”, “ax.set_xlabel(‘R [m]’)n”, “ax.set_ylabel(‘Z [m]’)n”, “ax.set_aspect(‘equal’)n”, “n”, “# Calculate several flux surface quantitiesn”, “print(‘Safety factor: %.2f’ % surf.eval_q[0])n”, “print(‘Length: %.2f m’ % surf.length)n”, “print(‘Area: %.4f m^2’ % surf.area)n”, “print(‘Volume: %.3f m^3’ % surf.volume)n”, “print(‘Toroidal current density: %.3f MA/m^2’ % (surf.tor_current/1e6))” “tr1 = eq.trace_field_line(r=eq.coordinates(psi_onq).r_mid[0], theta=0.09)[0]n”, “tr2 = eq.trace_field_line(tr1.R[-1], tr1.Z[-1], tr1.phi[-1])[0]n”, “tr3 = eq.trace_field_line(tr2.R[-1], tr2.Z[-1], tr2.phi[-1])[0]”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“## Locating the separatrix position in a given profilen”, “n”, “In experiment, one is often interested where the separatrix is along the chord of their measurement. In the following example the separatrix coordinates are calculated at the geometric outer midplane, that is, $Z=0$.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Define the measurement chord using two pointsn”, “chord = eq.coordinates(R=[0.6,0.8], Z=[0,0])n”, “n”, “# Calculate the intersection of the chord with the separatrix in 2Dn”, “intersection_point = chord.intersection(eq.lcfs, dim=2)n”, “n”, “# Plot the plasma with the intersection pointn”, “ax = plt.gca()n”, “eq.lcfs.plot()n”, “eq.first_wall.plot(c=’k’)n”, “chord.plot(color=’g’, marker=’x’)n”, “intersection_point.plot(marker=’o’, color=’r’)n”, “ax.set_aspect(‘equal’)n”, “ax.set_xlabel(‘$R$ [m]’)n”, “ax.set_ylabel(‘$Z$ [m]’)n”, “n”, “intersection_point.R”
]
}, {
“cell_type”: “markdown”, “metadata”: {
“collapsed”: true
}, “source”: [
“## Detector line of sight visualisationn”, “n”, “In this section, we demonstrate the flexibility of the Coordinates class by visualising a detector line of sight. Suppose we have a pixel detector at the position $[X, Y, Z] = [1.2 \, \mathrm{m}, 0 \, \mathrm{m}, -0.1 \, \mathrm{m}]$.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Define detector position [X, Y, Z]n”, “position = np.array((1.2, 0, -0.1))”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“The detector views the plasma mostly tangentially to the toroidal direction, but also sloping a little upward.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Define the line of sight direction (again along [X, Y, Z])n”, “direction = np.array((-1, 0.6, 0.2))n”, “n”, “# Norm the direction to unit lengthn”, “direction /= np.linalg.norm(direction)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“Now since the plasma geometry is curvilinear, the detector line of sight is not trivial. Luckily PLEQUE’s Coordinates class can easily express its stored coordinates both in the cartesian $[X,Y,Z]$ and the cylindrical $[R,Z,\phi]$ coordinate systems. In the following line, 20 points along the detector line of sight are calculated in 3D.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Calculate detector line of sight (LOS)n”, “LOS = eq.coordinates(position + direction[np.newaxis,:] * np.linspace(0, 2.0, 20)[:, np.newaxis],n”, “ coord_type=(‘X’, ‘Y’, ‘Z’)n”, “ )” “#Define flux surfaces where theta will be evaluatedn”, “psi_n = np.linspace(0, 1, 200)[1:-1]n”, “surfs = [eq._flux_surface(pn)[0] for pn in psi_n]n”, “n”, “#Define the flux surfaces which will show on the plotn”, “psi_n2 = np.linspace(0, 1, 7)[1:]n”, “surfs2 = [eq._flux_surface(pn)[0] for pn in psi_n2]n”, “n”, “#Define poloidal angles where theta isolines will be plottedn”, “thetas = np.linspace(0, 2*np.pi, 13, endpoint=False)n”, “n”, “#Prepare figuren”, “fig, axes = plt.subplots(1, 2, figsize=(10,6))n”, “ax1, ax2 = axesn”, “n”, “#Plot LCFS and several flux surfaces in both the plotsn”, “eq.lcfs.plot(ax = ax1, color = ‘k’, ls = ‘-’, lw=3)n”, “eq.lcfs.plot(ax = ax2, color = ‘k’, ls = ‘-’, lw=3)n”, “for s in surfs2:n”, “ s.plot(ax = ax1, color=’k’, lw = 1)n”, “ s.plot(ax = ax2, color=’k’, lw = 1)n”, “n”, “#Plot theta and theta_star isolinesn”, “for th in thetas:n”, “ # this is so ugly it has to implemented better as soon as possible (!)n”, “# print(th)n”, “ c = eq.coordinates(r = np.linspace(0, 0.4, 150), theta = np.ones(150)*th)n”, “ amin = np.argmin(np.abs(c.psi_n - 1))n”, “ r_lcfs = c.r[amin]n”, “ n”, “ psi_n = np.array([np.mean(s.psi_n) for s in surfs]) n”, “ c = eq.coordinates(r = np.linspace(0, r_lcfs, len(psi_n)), theta=np.ones(len(psi_n))*th)n”, “ c.plot(ax = ax1, color=’k’, lw=1)n”, “ n”, “ idxs = [np.argmin(np.abs(s.straight_fieldline_theta - th)) for s in surfs]n”, “ rs = [s.r[i] for s,i in zip(surfs,idxs)]n”, “ rs = np.hstack((0, rs))n”, “ thetas = [s.theta[i] for s,i in zip(surfs,idxs)]n”, “ thetas = np.hstack((0, thetas))n”, “ c = eq.coordinates(r = rs, theta = thetas)n”, “ c.plot(ax = ax2, color = ‘k’, lw=1)n”, “ n”, “# Make both the subplots prettyn”, “ax1.set_title(r’$\theta$’)n”, “ax1.set_aspect(‘equal’)n”, “ax1.set_xlabel(‘$R$ [m]’)n”, “ax1.set_ylabel(‘$Z$ [m]’)n”, “ax2.set_title(r’$\theta^*$’)n”, “ax2.set_aspect(‘equal’)n”, “ax2.set_xlabel(‘$R$ [m]’)n”, “ax2.set_ylabel(‘$Z$ [m]’)”
]
}, {
“cell_type”: “markdown”, “metadata”: {}, “source”: [
“To visualise the line of sight in top view $[X,Y]$ and poloidal cross-section view $[R,Z]$, we first define the limiter outline as viewed from the top. Then we proceed with the plotting.”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: [
“# Define the limiter outline as viewed from the topn”, “Ns = 100n”, “inner_lim = eq.coordinates(np.min(eq.first_wall.R)*np.ones(Ns), np.zeros(Ns), np.linspace(0, 2*np.pi, Ns))n”, “outer_lim = eq.coordinates(np.max(eq.first_wall.R)*np.ones(Ns), np.zeros(Ns), np.linspace(0, 2*np.pi, Ns))n”, “n”, “# Create a figuren”, “fig, axs = plt.subplots(1,2)n”, “n”, “# Plot the top view of the line of sightn”, “ax = axs[0]n”, “ax.plot(inner_lim.X, inner_lim.Y, ‘k-‘)n”, “ax.plot(outer_lim.X, outer_lim.Y, ‘k-‘)n”, “ax.plot(LOS.X, LOS.Y, ‘x–’, label=’Line of sight’)n”, “ax.plot(position[0], position[1], ‘d’, color=’C0’)n”, “ax.legend()n”, “ax.set_aspect(‘equal’)n”, “ax.set_xlabel(‘$X$ [m]’)n”, “ax.set_ylabel(‘$Y$ [m]’)n”, “n”, “# Plot the poloidal cross-section view of the line of sightn”, “ax = axs[1]n”, “ax.plot(eq.first_wall.R, eq.first_wall.Z, ‘k-‘)n”, “ax.plot(eq.lcfs.R, eq.lcfs.Z, ‘k–‘)n”, “ax.plot(LOS.R, LOS.Z, ‘x–‘)n”, “ax.plot(LOS.R[0], position[2], ‘d’, color=’C0’)n”, “ax.set_aspect(‘equal’)n”, “ax.set_xlabel(‘$R$ [m]’)n”, “ax.set_ylabel(‘$Z$ [m]’)”
]
}, {
“cell_type”: “code”, “execution_count”: null, “metadata”: {}, “outputs”: [], “source”: []
}
], “metadata”: {
- “kernelspec”: {
“display_name”: “Python 3”, “language”: “python”, “name”: “python3”
}, “language_info”: {
- “codemirror_mode”: {
“name”: “ipython”, “version”: 3
}, “file_extension”: “.py”, “mimetype”: “text/x-python”, “name”: “python”, “nbconvert_exporter”: “python”, “pygments_lexer”: “ipython3”, “version”: “3.7.0”
}
}, “nbformat”: 4, “nbformat_minor”: 1
}