# Penning Trap Simulation C++ simulation of charged particle dynamics in Penning traps with time-dependent electric fields. Includes numerical integration methods, analytical validation, and Python visualization tools. ## Project Structure ``` project3/ ├── src/ # C++ source files │ ├── main.cpp # Main simulation program │ ├── test.cpp # Analytical validation tests │ ├── particle.cpp # Particle class implementation │ └── penning_trap.cpp # Penning trap physics ├── include/ # Header files │ ├── particle.hpp # Particle class definitions │ └── penning_trap.hpp # Penning trap class ├── main.py # Python analysis and plotting ├── particles.in # Initial particle configuration ├── Makefile # Build configuration ├── tests/ # Test output (generated) ├── simulations/ # Simulation output (generated) └── figures/ # Generated plots (generated) ``` ## Requirements - **C++17** compiler (g++) - **Armadillo** library: `brew install armadillo` (macOS) or `sudo apt install libarmadillo-dev` (Ubuntu) - **Python** with numpy, matplotlib, pandas: `pip install numpy matplotlib pandas` ## Compilation ```bash make test # Build test program make main # Build main simulation ``` ## Running Simulations ### 1. Test Program (`./test`) Validates numerical methods against analytical solutions: ```bash ./test ``` - `n_steps`: Number of time steps (e.g., 4000, 8000, 16000, 32000) - `T_end`: Simulation time in microseconds (e.g., 50) - `method`: Integration method ("Euler" or "RK4") **Example:** ```bash ./test 4000 50 RK4 ``` ### 2. Main Program (`./main`) Two simulation modes: #### Frequency Analysis Mode Studies particle stability by sweeping the driving frequency and measuring how many particles remain confined in the trap. ```bash ./main f ``` - `omega_min/max`: Frequency range (rad/μs) - `step_size`: Frequency step size - `n_particles`: Number of particles (typically 100) - `amplitude`: Voltage oscillation amplitude (0.0-1.0) - `interactions`: 0 (disabled) or 1 (enabled) Coulomb interactions - `n_steps`: Integration steps (4000-32000) - `T_end`: Simulation time (microseconds) **Example:** ```bash ./main f 0.2 2.5 0.02 100 0.1 0 4000 500 ``` #### Motion Analysis Mode Tracks detailed particle trajectories over time, recording positions and velocities from initial conditions specified in `particles.in`. ```bash ./main m ``` - `amplitude`: Voltage oscillation amplitude (0.0-1.0) - `omega_V`: Angular frequency (rad/μs) - `interactions`: 0 (disabled) or 1 (enabled) Coulomb interactions - `n_steps`: Integration steps (4000-32000) - `T_end`: Simulation time (microseconds) **Example:** ```bash ./main m 0.1 2.0 0 4000 500 ``` ### 3. Automated Analysis (`python main.py`) Runs complete analysis pipeline with all tests and visualizations: ```bash python main.py ``` ## Input File Format `particles.in` specifies initial particle positions and velocities: ``` # x y z vx vy vz 20 0 20 0 25 0 25 25 0 0 40 5 ``` ## Output ### Generated Files - `tests/`: Numerical vs analytical comparison data - `simulations/`: Motion trajectories and frequency analysis - `figures/`: Plots of trajectories, phase space, error analysis ### Key Results - Particle motion with/without Coulomb interactions - Frequency-dependent stability analysis - Integration method accuracy comparison - Error convergence studies ## Reproducing Results 1. **Install dependencies** and compile: ```bash brew install armadillo # macOS pip install numpy matplotlib pandas make ``` 2. **Run complete analysis**: ```bash python main.py ``` 3. **View results**: Check `figures/` directory for generated plots ## Manual Execution Examples ```bash # Test numerical accuracy ./test 4000 50 RK4 ./test 8000 50 Euler # Frequency analysis without interactions ./main f 0.2 2.5 0.02 100 0.1 0 4000 500 # Motion analysis with interactions ./main m 0.1 2.0 1 4000 500 ``` ## Output Files and Data Analysis ### Simulation Output Structure #### Frequency Analysis Output ``` simulations/frequency_analysis/n4000_T500_f0.10_omegaV0.2-2.5_interactions0.txt ``` Contains frequency-dependent particle stability data: ``` # Frequency_analysis for amplitude f = 0.10, interactions = disabled and T_end = 500 microseconds omega_V fraction_left 2.000000e-01 0.95 4.000000e-01 0.92 ... ``` #### Motion Analysis Output ``` simulations/motion_analysis/n4000_T500_interactions0.txt ``` Time evolution of all particle coordinates: ``` # Motion analysis for amplitude f = 0.10, omega_V = 2.0, interactions = enabled, n_steps = 4000, T_end = 500 t x0 y0 z0 vx0 vy0 vz0 0.000000e+00 2.000000e+01 0.000000e+00 2.000000e+01 0.000000e+00 2.500000e+01 0.000000e+00 1.250000e-01 1.999688e+01 3.124844e+00 1.999375e+01 -2.499531e+00 2.499844e+01 -2.499687e-01 ... ``` #### Analytical Validation Output ``` tests/RK4_n50000_T50_numerical.txt # Numerical integration results tests/RK4_n50000_T50_analytical.txt # Exact analytical solutions tests/RK4_n50000_T50_error.txt # Error analysis with relative errors ``` ### Python Analysis Tools (`main.py`) #### Automated Analysis Pipeline ```python # Run comprehensive simulation campaign python main.py # Available functions in main.py: from pathlib import Path from main import ( run_tests, run_frequency_analysis, run_motion_analysis, plot_relative_errors, plot_z_motion, plot_frequency_analysis, plot_xy_motion_analysis, plot_phase_space_analysis ) # Run specific tests run_tests(['Euler', 'RK4'], [4000, 8000], 50.0) # Run frequency analysis run_frequency_analysis([(0.2, 2.5)], 0.02, 100, [0.1, 0.4], False, 4000, 500) # Generate plots plot_relative_errors(Path("tests"), Path("figures/errors"), ['RK4', 'Euler'], [4000, 8000], 50.0) ``` #### Generated Visualizations - **Error convergence plots**: Relative error vs. step size for different integration methods - **Particle trajectories**: Motion visualization - **Phase space analysis**: x-vx and z-vz phase portraits with interaction comparison - **Frequency response**: Particle stability across driving frequencies with resonance detection - **Comparative analysis**: Method accuracy, interaction effects, and performance comparisons - **Side-by-side plots**: Direct comparison of scenarios with/without Coulomb interactions ## Key Classes and Methods ### Particle Class (Base Class) ```cpp class Particle { public: Particle() = default; // Default constructor Particle(double q, double m, arma::vec3 r, arma::vec3 v); // Parameterized constructor // Accessors double charge() const; double mass() const; const arma::vec3& position() const; const arma::vec3& velocity() const; void print() const; }; ``` ### CalciumIon Class (Specialized Ca+ Ion) ```cpp class CalciumIon : public Particle { public: // Physical constants static constexpr double CHARGE = 1.0; // Elementary charge units static constexpr double MASS = 40.078; // Atomic mass units // Constructor (position and velocity required) CalciumIon(const arma::vec3& r, const arma::vec3& v); }; ``` ### PenningTrap Class ```cpp class PenningTrap { public: // Physical constants (publicly accessible) static constexpr double k_e = 1.38935333e+05; // Coulomb constant static constexpr double T = 9.64852558e+01; // Tesla conversion factor static constexpr double V = 9.64852558e+07; // Volt conversion factor // Constructors PenningTrap(); // Default parameters PenningTrap(double B0, double V0, double d); // Custom parameters // Particle management void add_particle(const Particle& p); void fill_random_particles(size_t n); // Add n random Ca+ ions const Particle& get_particle(int index = 0) const; size_t number_of_particles() const; // Time-dependent voltage control void set_time_dependent_voltage(double f, double omega_V); // Interaction control void set_interactions(bool enable); bool interactions_enabled() const; // Physics calculations (time-dependent) arma::vec3 external_E_field(const arma::vec3& r, double t = 0.0) const; arma::vec3 external_B_field(const arma::vec3& r) const; arma::vec3 force_particle(int i, int j) const; arma::vec3 total_force(int i, double t = 0.0) const; // Time evolution (time-aware) void evolve_RK4(double dt, double t = 0.0); void evolve_forward_Euler(double dt, double t = 0.0); // Accessors double magnetic_field_strength() const; double electric_potential() const; double characteristic_dimension() const; }; ``` ### Usage Examples #### Creating Particles ```cpp // Ca+ ions with position and velocity arma::vec3 r = {20.0, 0.0, 20.0}; // Position (μm) arma::vec3 v = {0.0, 25.0, 0.0}; // Velocity (μm/μs) CalciumIon ca_ion(r, v); // Generic particles Particle proton(1.0, 1.0, r, v); // Charge=1e, mass=1u ``` #### Setting Up Trap ```cpp PenningTrap trap; // Default: B0=1T, V0=25mV, d=500μm trap.add_particle(ca_ion); trap.set_time_dependent_voltage(0.1, 2.0); // f=0.1, ω=2.0 rad/μs trap.set_interactions(true); // Enable Coulomb forces ``` #### Time Evolution ```cpp double dt = 0.01, t = 0.0; for (int i = 0; i < 1000; ++i) { trap.evolve_RK4(dt, t); // Time-dependent evolution t += dt; } ``` ## License This project is developed for educational purposes as part of FYS4150 - Computational Physics. ## Authors - Vladislav Foss - Sebastian Koranda - Simen Norrud