Circularly polarized Alfvén waves

#define DBY (-pert*(sin(kx[0]*x+phi[0])+sin(kx[1]*x+phi[1])+sin(kx[2]*x+phi[2])))
#define DBZ (pert*(cos(kx[0]*x+phi[0])+cos(kx[1]*x+phi[1])+cos(kx[2]*x+phi[2])))

import numpy as np
from pywarpx import picmi
from space_analysis.simulation.warpx import HybridSimulation

constants = picmi.constants

def init_cp_alfven(
    ks,  #: wave numbers
    phases,  #: phases
    B0,
    A,  #: relative amplitude
    vA,
    n0,
    **kwargs,
):
    """
    Generate field with a wave propagating along the x axis at a large angle `theta` with respect to the background magnetic field lying in the x-z plane.

    The initial wave is an Alfven mode in which the magnetic field fluctuation points along the y and z axis and has a relative amplitude $A = \delta B_y / B_0$
    """

    wx_expression = "+".join([f"sin({k} * z + {phi})" for k, phi in zip(ks, phases)])
    wy_expression = "+".join([f"cos({k} * z + {phi})" for k, phi in zip(ks, phases)])


    Bz_expression = B0
    Bx_expression = f"{-A * B0} * ({wx_expression})"
    By_expression = f"{A * B0} * ({wy_expression})"

    pz_expression = 0
    px_expression = f"{A * vA} * ({wx_expression})"
    py_expression = f"{-A * vA} * ({wy_expression})"
    momentum_expressions = [px_expression, py_expression, pz_expression]

    field = picmi.AnalyticInitialField(
        Bx_expression=Bx_expression,
        By_expression=By_expression,
        Bz_expression=Bz_expression,
    )

    dist = picmi.AnalyticDistribution(
        density_expression=n0,
        momentum_expressions=momentum_expressions,
    )

    return field, dist

class CPAlfvenModes(HybridSimulation):
    # Applied field parameters
    dim: int = 1
    B0: float = 100 * 1e-9
    """Initial magnetic field strength (T)"""
    n0: float = 100 * 1e6
    """Initial plasma density (m^-3)"""

    A: float = 0.5  # relative amplitude

    ks: list[float] = None  # wave numbers
    phases: list[float] = None  # phases
    n_modes: int = 3  # number of modes

    # Spatial domain
    Lz_norm: float = 128
    Lx_norm: float = 64

    def setup_init_cond(self):
        """setup initial conditions"""
        self.ks = [
            (i + self.n_modes) * 2 * np.pi / self.Lz for i in range(self.n_modes)
        ]
        self.phases = [i * 2 * np.pi / self.n_modes for i in range(self.n_modes)]

        B_ext, dist = init_cp_alfven(
            ks=self.ks,
            phases=self.phases,
            B0=self.B0,
            A=self.A,
            vA=self.vA,
            n0=self.n0,
        )
        self._B_ext = B_ext
        self._dist = dist

        return self

simulation = CPAlfvenModes(
    diag_format='openpmd',
)

Numerical parameters:
    dt = 1.4e-04 s
    total steps = 25600

Initializing simulation with input parameters:
    Te = 24.588 eV
    n = 1.0e+02 cm^-3
    B0 = 100.00 nT
    M/m = 100

Plasma parameters:
    d_i = 5.3e+03 m
    t_ci = 3.6e-02 s
    v_ti = 2.1e+05 m/s
    vA = 9.3e+05 m/s
    vA/c = 0.0031021956193735692

simulation._sim.step()