example-001-basics.cpp

#include <elec.hpp>
#include <vector>
#include <iterator>
#include <cmath>

#define MAX_PLOT_FIELD_NORM 100
#define PLOT_FIELD_SCALE    .002
#define MESH_SIZE           40

#define V_MIN -10
#define V_MAX  20
#define NB_V_CONTOURS 20

#define NB_STALLED_FRAMES   20
#define NB_FRAMES_1        200
#define NB_FRAMES_2        100

#define MOTION_X0             .3
#define MOTION_Y0             -.5
#define MOTION_AMPLITUDE      2

#define ROTATION_VIEW .1

int main(int argc, char* argv[]) {

  if(argc != 2) {
    std::cout << "Usage : " << argv[0] << " run | elec-plot.py" << std::endl
              << "        " << argv[0] << " run | elec-plotV.py" << std::endl;
    return 0;
  }

  // Let us create 3 particles.
  std::vector<elec::Particle> particles;
  particles.push_back(elec::proton  ({-MOTION_X0, MOTION_Y0                    }));
  particles.push_back(elec::proton  ({ MOTION_X0, MOTION_Y0                    }));
  particles.push_back(elec::electron({         0, MOTION_Y0 + MOTION_AMPLITUDE }));
  
  // This is the plotting area.
  elec::Point min = {-1,-1};
  elec::Point max = { 1, 1};

  std::vector<elec::Point> mesh;    // This is where the E and V are computed for plotting.
  std::vector<elec::Point> Efield;  // There is one 2D value per mesh point.
  std::vector<double>      Vfield;  // There is one scalar value per mesh point.

  // Let us compute a regular mesh.
  elec::mesh(min,max,MESH_SIZE,MESH_SIZE,std::back_inserter(mesh));

  // This starts the plotting (warning, all plot functionality produce
  // an output on std::cout.... so use std::cerr for your own message
  // display.
  auto plot = elec::plot::figure("Electric Field", min, max);
  double a = 20;

  // We have to compute the electric field generated by the particles at all mesh points.
  Efield.clear();
  elec::E(particles.begin(), particles.end(),mesh.begin(), mesh.end(),std::back_inserter(Efield));
  // We do the same for the potential field.
  Vfield.clear();
  elec::V(particles.begin(), particles.end(),mesh.begin(), mesh.end(),std::back_inserter(Vfield));
  
  // .. and plot successive similar frames.
  for(unsigned int frame_id = 0; frame_id < NB_STALLED_FRAMES; ++frame_id, a+=ROTATION_VIEW) {
    plot.begin_frame("frame","jpg");
    plot.set_view(15,a); // This is used by 3D plots.
    plot.__plot_potential(mesh.begin(), mesh.end(), Vfield.begin(), Vfield.end(), V_MIN, V_MAX, NB_V_CONTOURS);
    plot.__plot_field(mesh.begin(), mesh.end(), Efield.begin(), Efield.end(),PLOT_FIELD_SCALE,MAX_PLOT_FIELD_NORM);
    plot.__plot_particles(particles.begin(), particles.end());
    plot.end_frame();
  }

  for(unsigned int frame_id = 0; frame_id < NB_FRAMES_1; ++frame_id, a+=ROTATION_VIEW) {
    // Let us reposition the last particle (the electron).
    particles[2].pos.y = MOTION_Y0 + MOTION_AMPLITUDE*(NB_FRAMES_1-1-frame_id)/(double)(NB_FRAMES_1-1);
    
    // We have to re-compute the electric field and potential.
    Efield.clear();
    Vfield.clear();
    elec::E(particles.begin(), particles.end(),mesh.begin(), mesh.end(),std::back_inserter(Efield));
    elec::V(particles.begin(), particles.end(),mesh.begin(), mesh.end(),std::back_inserter(Vfield));
    
    // Now, we can plot the figure
    plot.begin_frame("frame","jpg");
    plot.set_view(15,a);
    plot.__plot_potential(mesh.begin(), mesh.end(), Vfield.begin(), Vfield.end(), V_MIN, V_MAX, NB_V_CONTOURS);
    plot.__plot_field(mesh.begin(), mesh.end(), Efield.begin(), Efield.end(),PLOT_FIELD_SCALE,MAX_PLOT_FIELD_NORM);
    plot.__plot_particles(particles.begin(), particles.end());
    plot.end_frame();
  }

  for(unsigned int frame_id = 0; frame_id < NB_FRAMES_2; ++frame_id, a+=ROTATION_VIEW) {
    // Let us reposition the last particle (the electron).
    elec::Point O = {-MOTION_X0, MOTION_Y0};
    double r      = MOTION_X0;
    double theta  = M_PI_2*(frame_id)/(double)(NB_FRAMES_2-1);
    particles[2].pos = O + elec::Point(std::cos(theta),std::sin(theta))*r;
    
    Efield.clear();
    Vfield.clear();
    elec::E(particles.begin(), particles.end(),mesh.begin(), mesh.end(),std::back_inserter(Efield));
    elec::V(particles.begin(), particles.end(),mesh.begin(), mesh.end(),std::back_inserter(Vfield));
    
    // Now, we can plot the figure
    plot.begin_frame("frame","jpg");
    plot.set_view(15,a);
    plot.__plot_potential(mesh.begin(), mesh.end(), Vfield.begin(), Vfield.end(), V_MIN, V_MAX, NB_V_CONTOURS);
    plot.__plot_field(mesh.begin(), mesh.end(), Efield.begin(), Efield.end(),PLOT_FIELD_SCALE,MAX_PLOT_FIELD_NORM);
    plot.__plot_particles(particles.begin(), particles.end());
    plot.end_frame();
  }
  
  // We plot successive similar frames.
  for(unsigned int frame_id = 0; frame_id < NB_STALLED_FRAMES; ++frame_id, a+=ROTATION_VIEW) {
    plot.begin_frame("frame","jpg");
    plot.set_view(15,a);
    plot.__plot_potential(mesh.begin(), mesh.end(), Vfield.begin(), Vfield.end(), V_MIN, V_MAX, NB_V_CONTOURS);
    plot.__plot_field(mesh.begin(), mesh.end(), Efield.begin(), Efield.end(),PLOT_FIELD_SCALE,MAX_PLOT_FIELD_NORM);
    plot.__plot_particles(particles.begin(), particles.end());
    plot.end_frame();
  }

  return 0;
  
}