main.cpp

Go to the documentation of this file.

 
#include <cmath>
#include <fstream>
#include <omp.h>
#include <string>
#include <sys/stat.h>
#include <vector>
 
#include "PenningTrap.hpp"
#include "utils.hpp"
 
#define PARTICLES 100
#define N 10000
#define CHARGE 1.
#define MASS 40. // unit: amu
 
Particle p1(CHARGE, MASS, vec_3d{20., 0., 20.}, vec_3d{0., 25., 0.});
Particle p2(CHARGE, MASS, vec_3d{25., 25., 0.}, vec_3d{0., 40., 5.});
 
void simulate_single_particle()
{
    DEBUG("Inside single particle sim");
    PenningTrap trap(std::vector<Particle>{p1});
 
    double time = 50.; // microseconds
 
    DEBUG("Write to dir");
    trap.write_simulation_to_dir("output/simulate_single_particle", time, N);
}
 
void simulate_two_particles()
{
    PenningTrap trap_no_interaction(std::vector<Particle>{p1, p2});
    PenningTrap trap_with_interaction(std::vector<Particle>{p1, p2});
 
    double time = 50.; // microseconds
 
    trap_no_interaction.write_simulation_to_dir(
        "output/simulate_2_particles/no_interaction", time, N, "rk4", false);
    trap_with_interaction.write_simulation_to_dir(
        "output/simulate_2_particles/with_interaction", time, N);
}
 
void simulate_single_particle_with_different_steps()
{
 
    double time = 50; // microseconds
 
    for (int i = 0; i < 4; i++) {
        int steps = 4000 * (i + 1);
        PenningTrap trap(std::vector<Particle>{p1});
        trap.write_simulation_to_dir("output/N_steps/RK4/" +
                                         std::to_string(steps) + "_steps",
                                     time, steps, "rk4", false);
    }
 
    for (int i = 0; i < 4; i++) {
        int steps = 4000 * (i + 1);
        PenningTrap trap(std::vector<Particle>{p1});
        trap.write_simulation_to_dir("output/N_steps/euler/" +
                                         std::to_string(steps) + "_steps",
                                     time, steps, "euler", false);
    }
}
 
void simulate_100_particles()
{
    PenningTrap trap((unsigned)100, T,
                     [](double t) {
                         return 25. * V / 1000. *
                                (1. + .4 * std::cos(1.5 * t));
                     },
                     500., 0);
 
    double time = 500.; // microseconds
 
    trap.write_simulation_to_dir("output/simulate_100_particles", time, N*5);
}
 
void simulate_100_particles_with_time_potential()
{
    double amplitudes[]{.1, .4, .7};
 
    double freq_start = .2;
    double freq_end = 2.5;
    double freq_increment = .02;
    size_t freq_iterations = (size_t) ((freq_end - freq_start) / freq_increment);
 
    std::string path = "output/time_dependent_potential/";
    mkpath(path);
 
    std::ofstream ofile;
 
    for (double f : amplitudes) {
        ofile.open(path + "f_" + std::to_string(f) + ".txt");
        #pragma omp parallel for ordered schedule(static, 1)
        for (size_t i=0; i < freq_iterations; i++) {
            double freq = freq_start + i*freq_increment;
            PenningTrap trap((unsigned)100, T,
                             [f, freq](double t) {
                                 return (25. * V / 1000.) *
                                        (1. + f * std::cos(freq * t));
                             },
                             500., 0.);
            double res = trap.fraction_of_particles_left(500., 40000, "rk4", true);
            #pragma omp ordered
            ofile << freq << "," << res << "\n";
        }
        ofile.close();
    }
}
 
int main()
{
 
    simulate_single_particle();
 
    simulate_two_particles();
 
    simulate_single_particle_with_different_steps();
 
    double start = omp_get_wtime();
 
    //simulate_100_particles();
 
    simulate_100_particles_with_time_potential();
 
    double end = omp_get_wtime();
 
    std::cout << "Time: " << end - start << " seconds" << std::endl;
 
    return 0;
}