example-003-generator.cpp

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


#define SIM_SPLIT     10
#define SIM_DT         1

#define MAX_PLOT_FIELD_NORM  100000
#define PLOT_FIELD_SCALE     5e-6
#define MESH_SIZE            40

#define CONNECTOR_POS    1.5
#define CONNECTOR_RADIUS  .2
#define WIRE_RADIUS       .1
#define RESISTANCE_RADIUS .2

#define NB_V_CONTOURS 20

#define VOLTAGE 10000

#define CONNECTION_STEP_1 20
#define CONNECTION_STEP_2 700

#define WIRE_RESISTIVITY .01
#define RESISTANCE       3

#define ANGLE_AMPL  20
#define DTHETA     .15

void dt(double& theta) {
  theta += DTHETA;
  if(theta > 180)
    theta -= 360;
}

double to_rad(double theta) {
  return theta*M_PI/180;
}

double azimuth(double theta, double angle_init) {
  return angle_init + ANGLE_AMPL*std::sin(to_rad(theta));
}


int main(int argc, char* argv[]) {
  if(argc != 3) {
    std::cout << "Usage : " << argv[0] << " <particle|field|both> <a|b|c...> | elec-plot.py" << std::endl
              << "        " << argv[0] << " <particle|field|both> <a|b|c...> | elec-plotV.py" << std::endl;
    return 0;
  }

  bool part  = true;
  bool field = true;
  std::string mode(argv[1]);
  std::string cond(argv[2]);
  double vmin = -(1.5)*VOLTAGE;
  double vmax =   .5*VOLTAGE;
  double alpha_0 = 0;
  double elevation = 0;

  
  if(mode == "particle") {
    part  = true;
    field = false;
  }
  else if(mode == "field") {
    part  = false;
    field = true;
  }

  elec::World world;
  double bound = CONNECTOR_POS - .08;
  elec::conductor::Rectangle ra1({-bound,-WIRE_RADIUS},{bound,WIRE_RADIUS},0);

  elec::conductor::Rectangle rb1({-bound,-WIRE_RADIUS},{-1,WIRE_RADIUS},WIRE_RESISTIVITY);
  elec::conductor::Rectangle rb2({     -.5,-WIRE_RADIUS},{bound,WIRE_RADIUS},WIRE_RESISTIVITY);
  elec::conductor::Rectangle rb3({     0,-RESISTANCE_RADIUS},{1,RESISTANCE_RADIUS},RESISTANCE);
  elec::conductor::Rectangle rb4({-1.1,-WIRE_RADIUS},{-.4,WIRE_RADIUS},WIRE_RESISTIVITY);


  elec::conductor::Rectangle rc1({-1.4,-0.1},{ 0.0, 0.1},WIRE_RESISTIVITY);
  elec::conductor::Rectangle rc2({ 1.0,-0.1},{ 1.4, 0.1},WIRE_RESISTIVITY);
  elec::conductor::Rectangle rc3({-0.1, 0.3},{ 1.2, 0.5},WIRE_RESISTIVITY);
  elec::conductor::Rectangle rc4({-0.1,-0.5},{ 1.2,-0.3},WIRE_RESISTIVITY);
  elec::conductor::Rectangle rc5({-0.1,-0.5},{ 0.1, 0.5},WIRE_RESISTIVITY);
  elec::conductor::Rectangle rc6({ 1.0,-0.5},{ 1.2, 0.5},WIRE_RESISTIVITY);

  elec::conductor::Rectangle rc7({-1.1,-0.2},{-0.4, 0.2},RESISTANCE);
  elec::conductor::Rectangle rc8({ 0.2, 0.2},{ 0.9, 0.6},RESISTANCE);
  elec::conductor::Rectangle rc9({ 0.2,-0.6},{ 0.9,-0.2},RESISTANCE);

  // This adds the two connectors of a generator in the world.
  world.set_generator({-CONNECTOR_POS,0},{CONNECTOR_POS,0},CONNECTOR_RADIUS,VOLTAGE);

  if(cond == "a") {}
  else if(cond == "b") {
    world += rb3;
    world += rb1;
    world += rb2;
    alpha_0 = 90;
    elevation = 20;
  }
  else if(cond == "c") {
    alpha_0 = 0;
    elevation = 50;
    world += rc7;
    world += rc8;
    world += rc9;

    world += rc1;
    world += rc2;
    world += rc3;
    world += rc4;
    world += rc5;
    world += rc6;
  }

  world.add_protons();
  world.add_electrons(world.nb_protons());
  world.end_of_particles();

  elec::Point min(-CONNECTOR_POS - CONNECTOR_RADIUS, -1);
  elec::Point max( CONNECTOR_POS + CONNECTOR_RADIUS,  1);
  std::vector<elec::Point> mesh;
  std::vector<elec::Point> Efield;
  std::vector<double>      Vfield;

  elec::mesh(min,max,MESH_SIZE,MESH_SIZE,std::back_inserter(mesh));

  auto plot = elec::plot::figure("Battery", min, max);

  Efield.clear(); Vfield.clear();
  if(field) {
    if(cond == "a")
      elec::E(world.particles().begin(), world.particles().end(),mesh.begin(), mesh.end(),std::back_inserter(Efield));
    elec::V(world.particles().begin(), world.particles().end(),mesh.begin(), mesh.end(),std::back_inserter(Vfield));
  }

  unsigned int nb_steps = 0;
  double theta = 0;
  if(cond == "a") 
    nb_steps = 400;
  else if(cond == "b") 
    nb_steps = 1100;
  else if(cond == "c") 
    nb_steps = 1700;
  for(unsigned int step = 0; step < nb_steps; ++step, dt(theta)) {
    if(step == CONNECTION_STEP_1 && cond == "a") {
      // New particles are added here, while the simulation is
      // running.
      world += ra1;
      world.add_protons();
      // world.nb_protons() returns the number of newly added protons.
      world.add_electrons(world.nb_protons());
      world.end_of_particles(); // Do not forget this.
    }
    if(step == CONNECTION_STEP_2 && cond == "b") {
      world += rb4;
      world.add_protons();
      world.add_electrons(world.nb_protons());
      world.end_of_particles();
    }

    plot.begin_frame("frame","jpg"); 
    plot.set_view(elevation,azimuth(theta,alpha_0));

    if(part) 
      plot.__plot_particles(world.particles().begin(), world.particles().end());
    if(field) {
      plot.__plot_potential(mesh.begin(), mesh.end(), Vfield.begin(), Vfield.end(), vmin, vmax, NB_V_CONTOURS);
      if(cond == "a")
        plot.__plot_field(mesh.begin(), mesh.end(), Efield.begin(), Efield.end(),PLOT_FIELD_SCALE,MAX_PLOT_FIELD_NORM);
    }
    plot.end_frame();
    
    world.sim(SIM_DT,SIM_SPLIT);

    if(field) {
      Efield.clear(); Vfield.clear();
      if(cond == "a")
        elec::E(world.particles().begin(), world.particles().end(),mesh.begin(), mesh.end(),std::back_inserter(Efield));
      elec::V(world.particles().begin(), world.particles().end(),mesh.begin(), mesh.end(),std::back_inserter(Vfield));
    }
  }

  
}