Field Analysis¶
%config InlineBackend.figure_format = 'retina' # Nicer plots
from math import pi
import h5py
import matplotlib.pyplot as plt
import numpy as np
import scipy
from scipy.constants import c, e, epsilon_0
from genesis.version4 import Genesis4, Write
h = scipy.constants.value("Planck constant in eV/Hz")
Z0 = pi * 119.916983 # V^2 / m
Create field data¶
%%time
G = Genesis4("data/basic4/cu_hxr.in", verbose=True)
G.input.main.time.sample = 200
G.input.main.track.zstop = 40
G.input.main.namelists.append(Write(field="end"))
G.nproc = 0 # auto-select
output = G.run()
Configured to run in: /tmp/tmpyj429inl Setting use_mpi = True because nproc = 2 Running Genesis4 in /tmp/tmpyj429inl /home/runner/miniconda3/envs/lume-genesis-dev/bin/mpirun -n 2 /home/runner/miniconda3/envs/lume-genesis-dev/bin/genesis4 -l hxr.lat genesis4.in --------------------------------------------- GENESIS - Version 4.6.6 has started... Compile info: Compiled by conda at 2024-06-10 16:25:29 [UTC] from Git Commit ID: 07004e4cb437003b4b6cf024122db18caee13722 Starting Time: Thu Aug 29 16:13:29 2024 MPI-Comm Size: 2 nodes Opened input file genesis4.in Parsing lattice file hxr.lat ... Setting up time window of 15.0027 microns with 544 sample points... Generating input radiation field for HARM = 1 ...
Adding profile with label: beamcurrent Adding profile with label: beamgamma Generating input particle distribution... Running Core Simulation... Time-dependent run with 544 slices for a time window of 15.0027 microns Initial analysis of electron beam and radiation field... Calculation: 0% done
Calculation: 10% done
Calculation: 20% done
Calculation: 30% done
Calculation: 40% done
Calculation: 50% done
Calculation: 60% done
Calculation: 70% done
Calculation: 80% done
Calculation: 90% done
Calculation terminated due to requested stop. Writing output file...
Core Simulation done. End of Track Writing field distribution to file: end.fld.h5 ... Program is terminating... Ending Time: Thu Aug 29 16:18:14 2024 Total Wall Clock Time: 284.99 seconds -------------------------------------
Success - execution took 285.13s.
CPU times: user 242 ms, sys: 167 ms, total: 408 ms Wall time: 4min 45s
G.path
'/tmp/tmpyj429inl'
G.plot("field_energy", y2=["field_xsize", "field_ysize"], ylim2=[0, 100e-6])
G.verbose = True
output.load_fields()
['end']
3D Field data¶
The full field data is stored as a 3D array of complex numbers DFL
in units of sqrt(W)
.
The relation of this and the electric field E
in V/m is:
E = DFL * sqrt(2*Z0) / Δ
Where Z0 = π * 119.9169832 V^2/W
exactly and Δ
is the grid spacing.
DFL = output.field3d["end"].dfl # sqrt(W)
DFL.dtype, DFL.shape
(dtype('complex128'), (101, 101, 544))
param = output.field3d["end"].param
param
FieldFileParams( gridpoints=101, gridsize=2e-06, wavelength=1.3789244869952112e-10, slicecount=544, slicespacing=2.7578489739904225e-08, )
Gather some convenient variables and arrays:
Ns = param.slicecount
ds = param.slicespacing
λ0 = param.wavelength
f0 = c / λ0
dt = ds / c
Δ = param.gridsize
s = np.arange(0, Ns) * ds
t = -s / c
Field power¶
The power array sums over the x and y components of the absolute square of the field data.
power0 = G.output["field_power"][-1, :] # W
power1 = np.sum(np.sum(np.abs(DFL) ** 2, axis=0), axis=0) # W
These are the same:
np.allclose(power0, power1)
True
fig, ax = plt.subplots()
ax.plot(s * 1e6, power0 / 1e9)
ax.set_xlabel("s (µm)")
ax.set_ylabel("power (GW)")
Text(0, 0.5, 'power (GW)')
Field energy¶
The total field energy is the integral:
energy1 = np.trapz(power1, dx=dt)
energy1 # J
/tmp/ipykernel_2578/1993435925.py:1: DeprecationWarning: `trapz` is deprecated. Use `trapezoid` instead, or one of the numerical integration functions in `scipy.integrate`. energy1 = np.trapz(power1, dx=dt)
np.float64(2.7945317690698246e-07)
np.sum(power1) * dt # J
np.float64(2.794531804544018e-07)
On-axis field intensity and phase¶
field = G.output.field
intensity0 = field.intensity_nearfield[-1, :] # W/m^2
phase0 = field.phase_nearfield[-1, :] # radian
icenter = param.gridpoints // 2
field1 = DFL[icenter, icenter, :] # sqrt(W)
phase1 = np.angle(field1) # radian
intensity1 = np.abs(field1**2) / Δ**2 # W/m^2
These are the same:
np.allclose(intensity0, intensity1)
True
np.allclose(phase0, phase1)
True
fig, ax = plt.subplots()
ax.plot(s * 1e6, intensity0)
ax.set_xlabel("s (µm)")
ax.set_ylabel(r"intensity (W/m$^2$)")
Text(0, 0.5, 'intensity (W/m$^2$)')
The same field can be reconstructed from these arrays
field0 = np.sqrt(intensity0) * np.exp(1j * phase0) * Δ # sqrt(W)
These are the same:
np.allclose(field0, field1)
True
Spectral fluence¶
The spectrum calculation takes some care with the FFT and units.
def spectrum_from_field(field, dt=1):
"""
Calculates the spectrum (fourier transformed field)
from a complex field array with spacing dt.
Parameters
----------
field: nd.array of shape (n,)
Complex field
dt: float
Spacing of the field data in some units (e.g. 's')
Returns
-------
freqs: nd.array of shape (n,)
Frequencies in reciprocal space with inverse units (e.g. 'Hz = 1/s')
spectrum: nd.array of shape (n,)
The fourier transformed field
"""
if len(field.shape) != 1:
raise ValueError("Only 1D arrays are currently supported")
spectrum = np.fft.fftshift(np.fft.fft(field)) * dt
ns = len(field)
freqs = np.fft.fftshift(np.fft.fftfreq(ns, dt))
return freqs, spectrum
freqs, spectrum = spectrum_from_field(field0, dt)
# Frequency spacing
# df = np.diff(freqs)[0] # Hz
df = 1 / (dt * len(field0))
df
19982569111045.027
Check Plancherel theorem
np.sum(np.abs(field0) ** 2) * dt # J
np.float64(5.181066945480859e-10)
np.sum(np.abs(spectrum) ** 2) * df # J
np.float64(5.181066945480859e-10)
freqs, spectrum = spectrum_from_field(field0, dt)
fig, ax = plt.subplots()
photon_energies = h * (freqs + f0) # eV
spectral_fluence = np.abs(spectrum**2) / Δ**2 / h # J/m^2/eV
ax.plot(photon_energies, spectral_fluence)
ax.set_xlabel("photon energy (eV)")
ax.set_ylabel(r"spectral fluence (J/m$^2$/eV)")
Text(0, 0.5, 'spectral fluence (J/m$^2$/eV)')
fig, ax = plt.subplots()
J_per_photon = h * (freqs + f0) * e # note h is in eV/Hz
ax.plot(photon_energies, spectral_fluence / J_per_photon * 1e-6)
ax.set_xlabel("photon energy (eV)")
ax.set_ylabel(r"spectral intensity (photons/mm$^2$/eV)")
Text(0, 0.5, 'spectral intensity (photons/mm$^2$/eV)')
3D FFT¶
WORK IN PROGRESS
field3d = DFL * np.sqrt(2 * Z0) / Δ # Electric field
# field3d = np.pad(field3d, 100, constant_values = 0) #Pad
nx, ny, nz = field3d.shape
dx, dy, dz = Δ, Δ, ds
nx, ny, nz, dx, dy, dz
# reciprocal spacings
dkx, dky, dkz = 1 / (nx * dx), 1 / (ny * dy), 1 / (nz * dz)
spectrum3d = np.fft.fftshift(np.fft.fftn(field3d)) * dx * dy * dz
# Check total energy integral
np.sum(np.abs(field3d) ** 2) * dx * dy * dz * epsilon_0 / 2
np.float64(2.794531798361843e-07)
np.sum(np.abs(spectrum3d) ** 2) * dkx * dky * dkz * epsilon_0 / 2
np.float64(2.794531798361847e-07)
kx = np.fft.fftshift(np.fft.fftfreq(nx, dx))
ky = np.fft.fftshift(np.fft.fftfreq(ny, dy))
kz = np.fft.fftshift(np.fft.fftfreq(nz, dz))
Kx, Ky, Kz = np.meshgrid(kx, ky, kz, indexing="ij")
kz0 = 2 * pi / λ0
Kz = Kz + kz0
K = np.sqrt(Kx**2 + Ky**2 + Kz**2)
k = K.flatten()
dEk = np.abs(spectrum3d.flatten()) ** 2 * dkx * dky * dkz * epsilon_0 / 2
np.sum(dEk)
np.float64(2.7945317983618417e-07)
x = k * c * h / (2 * pi) # eV
w = dEk
bins = 1000
hist, bin_edges = np.histogram(x, bins=bins, weights=w)
hist_x = bin_edges[:-1] + np.diff(bin_edges) / 2
hist_width = np.diff(bin_edges)
hist_y = hist / hist_width
# hist_y, hist_f, hist_prefix = nice_array(hist/hist_width)
fig, ax = plt.subplots()
ax.bar(hist_x, hist_y * 1e6, hist_width, color="gray")
ax2 = ax.twinx()
ax2.plot(hist_x, np.cumsum(hist) * 1e6)
ax2.set_ylabel("cumulative energy (µJ)")
ax2.set_ylim(0, None)
ax.set_xlabel("photon energy (eV)")
ax.set_ylabel("spectral energy (µJ/eV)")
Text(0, 0.5, 'spectral energy (µJ/eV)')
Field data formats¶
Genesis4 writes a custom field format. LUME-Genesis provides a reader for this, as well as a conversion tool to write in the openPMD-wavefront standard.
Read Field h5¶
G.output.load_fields()
['end']
end_field = G.output.field3d["end"]
DFL = end_field.dfl
PARAM = end_field.param
DFL.shape, DFL.dtype, PARAM
((101, 101, 544), dtype('complex128'), FieldFileParams( gridpoints=101, gridsize=2e-06, wavelength=1.3789244869952112e-10, slicecount=544, slicespacing=2.7578489739904225e-08, ))
Write Wavefront in openPMD-wavefront¶
end_field.write_openpmd_wavefront("genesis4_wavefront.h5", verbose=True)
Writing wavefront (dfl data) to file genesis4_wavefront.h5
# Read back
with h5py.File("genesis4_wavefront.h5", "r") as h5:
print(list(h5["data"]["000000"]["meshes"]))
g = h5["data/000000/meshes/electricField"]
print(dict(g.attrs))
# Get
E2 = h5["data/000000/meshes/electricField/x"][:]
['electricField'] {'axisLabels': array(['x', 'y', 'z'], dtype=object), 'geometry': 'cartesian', 'gridGlobalOffset': array([-1.00000000e-04, -1.00000000e-04, -7.50134921e-06]), 'gridSpacing': array([2.00000000e-06, 2.00000000e-06, 2.75784897e-08]), 'gridUnitDimension': array([[1, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0], [1, 0, 0, 0, 0, 0, 0]]), 'gridUnitSI': array([1., 1., 1.]), 'photonEnergy': np.float64(8991.3696944618), 'photonEnergyUnitDimension': array([ 2, 1, -2, 0, 0, 0, 0]), 'photonEnergyUnitSI': np.float64(1.602176634e-19), 'timeOffset': np.float64(0.0), 'unitDimension': array([ 1, 1, -3, -1, 0, 0, 0])}
# Data is the same
np.allclose(E2, DFL)
True
Plot¶
Simple plot
# Sum over y and compute the absolute square
dat2 = np.sum(np.abs(DFL) ** 2, axis=1)
# Plot the center
plt.figure(figsize=(12, 8))
plt.imshow(dat2[::, :], aspect="auto")
# plt.axis('off');plt.savefig('../../assets/field.png', bbox_inches='tight')
<matplotlib.image.AxesImage at 0x7f39d4fd2c30>
def plot_slice(i=0):
dat = np.angle(DFL[:, :, i])
dx = PARAM.gridsize * 1e6
plt.figure(figsize=(8, 8))
plt.xlabel("x (µm)")
plt.xlabel("y (µm)")
plt.title(f"Phase for slice {i}")
plt.imshow(dat.T, origin="lower", extent=[-dx, dx, -dx, dx])
plot_slice(i=100)
Interactive¶
try:
from ipywidgets import interact
except ImportError:
pass
else:
interact(plot_slice, i=(0, len(DFL[0, 0, :]) - 1, 1))
Cleanup¶
import os
os.remove("genesis4_wavefront.h5")