elecWorld.hpp

Go to the documentation of this file.

#pragma once

/*   This file is part of rl-lib
 *
 *   Copyright (C) 2010,  Herve FREZZA-BUET
 *
 *   Author : Herve Frezza-Buet
 *
 *   Contributor :
 *
 *   This library is free software; you can redistribute it and/or
 *   modify it under the terms of the GNU General Public
 *   License (GPL) as published by the Free Software Foundation; either
 *   version 3 of the License, or any later version.
 *   
 *   This library is distributed in the hope that it will be useful,
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 *   General Public License for more details.
 *   
 *   You should have received a copy of the GNU General Public
 *   License along with this library; if not, write to the Free Software
 *   Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 *
 *   Contact : Herve.Frezza-Buet@supelec.fr
 *
 */

#include <list>
#include <iterator>
#include <vector>
#include <map>
#include <algorithm>
#include <iterator>
#include <limits>
#include <elecConductor.hpp>
#include <elecParticle.hpp>
#include <cstddef>

namespace elec {
  class World {
  private:
    class Connector : public elec::conductor::Disk {
    public:
      Connector(const elec::Point& C,
           double R,
           double r) 
        : elec::conductor::Disk(C,R,r) {}
      virtual bool inside(const elec::Point& p) {
        double dist = elec::d(center,p);
        return .05 < dist && dist <= radius;
      }
    };

    std::vector<elec::conductor::Metal*> _conductors;
    std::list<elec::Particle> _particles; // Don't use a vector here, due to push_back reallocation.
    std::vector<elec::Particle*> _protons;
    std::vector<elec::Particle*> _electrons;
    elec::Point _min, _max;
    unsigned int nb_p, nb_e;
    unsigned int nb_part;
    unsigned int nb_cond;
    bool has_generator;
    Connector cathode, anode;
    double ddp;

    elec::conductor::Metal* owner(const elec::Point& p) {
      elec::conductor::Metal* res = nullptr;
      for(auto it = _conductors.begin();
          it != _conductors.end() && res == nullptr;
          ++it) 
        if((*it)->inside(p))
          res = *it;
      return res;
    }
    

  public:

    World() : _conductors(), _particles(), _protons(), _electrons(),
              nb_p(0), nb_e(0), nb_part(0), nb_cond(0), has_generator(false),
              cathode({0,0},0,1), anode({0,0},0,1),
              ddp(0) {
      _min = {std::numeric_limits<double>::max(),std::numeric_limits<double>::max()};
      _max = {std::numeric_limits<double>::lowest(),std::numeric_limits<double>::lowest()};
    }

    void set_generator(const elec::Point& cathode_pos,
                       const elec::Point& anode_pos,
                       double radius, double voltage) {
      has_generator  = true;
      cathode.radius = radius;
      anode.radius   = radius;
      cathode.center = cathode_pos;
      anode.center   = anode_pos;
      ddp            = voltage;
      
      *this += cathode;
      *this += anode;
    }

    const std::list<elec::Particle>& particles() {return _particles;}

    World& operator+=(elec::conductor::Metal& conductor) {
      conductor.idf = _conductors.size();
      _conductors.push_back(&conductor);
      auto bbox = conductor.bounding_box();
      if(bbox.first.x < _min.x) _min.x = bbox.first.x;
      if(bbox.first.y < _min.y) _min.y = bbox.first.y;
      if(bbox.second.x > _max.x) _max.x = bbox.second.x;
      if(bbox.second.y > _max.y) _max.y = bbox.second.y;
      return *this;
    }
    
    void add_electron(const elec::Point& pos) {
      _particles.push_back(elec::electron(pos));
      ++nb_e;
    }

    void add_proton(const elec::Point& pos) {
      _particles.push_back(elec::proton(pos));
      ++nb_p;
    }

    void add_protons() {
      double delta = 1/(double)elecPROTONS_PER_UNIT;
      elec::Point pos;
      for(pos.x = _min.x; pos.x <= _max.x; pos.x+=delta)
        for(pos.y = _min.y; pos.y <= _max.y; pos.y+=delta) 
          if(in_new_conductors(pos))
            add_proton(pos);
    }

    void add_electrons(unsigned int nb) {
      for(unsigned int i=0; i<nb; ++i) {
        elec::Point pos = elec::uniform(_min,_max);
        while(!in_new_conductors(pos))
          pos = elec::uniform(_min,_max);
        add_electron(pos);
      }
    }

    void end_of_particles() {
      auto it = _particles.begin();
      for(std::advance(it,nb_part); it != _particles.end(); ++it) {
        auto& p = *it;
        if(p.q > 0)
          _protons.push_back(&p);
        else
          _electrons.push_back(&p);
      }
      nb_cond  = _conductors.size();
      nb_part += nb_p + nb_e;
      
      nb_p = 0;
      nb_e = 0;
      // _min = {std::numeric_limits<double>::max(),std::numeric_limits<double>::max()};
      // _max = {std::numeric_limits<double>::lowest(),std::numeric_limits<double>::lowest()};
    }

    unsigned int nb_protons() {return nb_p;}


    void sim(double dt, unsigned int nb_splits) {
      std::random_shuffle(_electrons.begin(), _electrons.end());
      unsigned int nb = _electrons.size();
      unsigned int first = 0;
      unsigned int last;
      for(unsigned int i = 1; i <= nb_splits; ++i, first = last) {
        last = (unsigned int)(nb*i/(double)(nb_splits)+.5);
        sim_partial(dt,_electrons.begin() + first, _electrons.begin() + last);
      } 
      
      if(has_generator)
        update_generator();
    }

  private:

    bool in_new_conductors(const elec::Point& p) {
      auto cond_ptr = owner(p);
      return cond_ptr != nullptr 
        && cond_ptr->idf >= nb_cond;
    }

    std::vector<elec::Particle*> source_electrons;
    void update_generator() {
      // Let us compute the potentials at both sides.
      double V_plus  = elec::V(_particles.begin(), _particles.end(), cathode.center);
      double V_minus = elec::V(_particles.begin(), _particles.end(), anode.center);
      double DV = V_plus - V_minus;
      double EXTRA_V = ddp - DV;

      Connector* source;
      Connector* dest;
      double* V_source;
      double* V_dest;
      if(EXTRA_V > 0) {
        // we need to extract electrons from the cathode to the anode.
        source = &cathode;
        dest   = &anode;
        V_source = &V_plus;
        V_dest   = &V_minus;
      }
      else {
        // we need to extract electrons from the anode to the cathode.
        source = &anode;
        dest   = &cathode;
        V_source = &V_minus;
        V_dest   = &V_plus;
      }

      auto dest_bb = dest->bounding_box();

      // Let us first regulate the potential by transferring electrons
      // from source to destination.

      // We first need to identify the electrons in the source.
      source_electrons.clear();
      for(auto e_ptr : _electrons) 
        if(source->inside(e_ptr->pos))
          source_electrons.push_back(e_ptr);
      std::random_shuffle(source_electrons.begin(),source_electrons.end());

      double delta = std::fabs(ddp - (V_plus - V_minus));
      double delta_prev = delta+1;
      for(auto it = source_electrons.begin(); it != source_electrons.end() && (delta < delta_prev); ++it) {
        // Let us choose a random position in dest.
        elec::Point pos = dest->center;
        while(!(dest->inside(pos)))
          pos = elec::uniform(dest_bb.first,dest_bb.second);

        elec::Particle& e = *(*it);
        (*V_source) -= elec::V(e,source->center); // Remove from source
        (*V_dest)   -= elec::V(e,dest->center);
        e.pos   = pos;
        e.speed = 0;
        (*V_source) += elec::V(e,source->center); // add to dest;
        (*V_dest)   += elec::V(e,dest->center); 
        delta_prev = delta;
        delta = std::fabs(ddp - (V_plus - V_minus));
      }

      // Let us add the remaining electrons at dest
      while(delta < delta_prev) {
        elec::Point pos = dest->center;
        while(!(dest->inside(pos)))
          pos = elec::uniform(dest_bb.first,dest_bb.second);
        add_electron(pos);
        auto e = elec::electron(pos);
        (*V_source) += elec::V(e,source->center); // add to dest;
        (*V_dest)   += elec::V(e,dest->center); 
        delta_prev = delta;
        delta = std::fabs(ddp - (V_plus - V_minus));
      }
      end_of_particles();
    }


    void sim_partial(double dt,
                     std::vector<elec::Particle*>::iterator begin,
                     std::vector<elec::Particle*>::iterator end) {
      std::vector<elec::Point> electrons_field; // at considered electrons
      std::vector<elec::Point> protons_field;   // at considered electrons
      elec::E_(_electrons.begin(), _electrons.end(),
               begin,end,
               std::back_inserter(electrons_field),
               [](const elec::Particle* p) -> const elec::Particle& {return *p;},
               [](const elec::Particle* p) -> const elec::Point& {return p->pos;});
      elec::E_(_protons.begin(), _protons.end(),
               begin,end,
               std::back_inserter(protons_field),
               [](const elec::Particle* p) -> const elec::Particle& {return *p;},
               [](const elec::Particle* p) -> const elec::Point& {return p->pos;});
      auto pf = protons_field.begin();
      auto ef = electrons_field.begin();
      for(auto it = begin; it != end; ++it, ++pf, ++ef) {
        elec::Particle& electron = *(*(it));
        elec::Point force = -(*pf + *ef)*elecFORCE_COEF;

        auto conductor = owner(electron.pos);
        if(conductor != nullptr) {
          // Let us compute the friction.
          force += conductor->friction(electron);

          // Let us apply the dynamics law.
          electron.speed += force*dt/elecELECTROM_MASS;
          if(electron.speed*electron.speed > elecMAX_SPEED*elecMAX_SPEED) 
            electron.speed = (*electron.speed)*elecMAX_SPEED;
          elec::Point dp  = electron.speed*dt;
          electron.pos   += dp*dt;
        }
        else {
          // this is the Frenet Coordinates
          elec::Point T = *(-(*pf));
          elec::Point S = {-T.y, T.x};
          
          // We have to keep only the tangential component.
          elec::Point speed =  electron.speed + force*dt/elecELECTROM_MASS;
          electron.speed = S*(speed*S);
          
          // We put the particle back to the conductor (ignoring the speed computation).
          electron.pos   += T*elecISOLATION_COEF+electron.speed*dt;
        }
      }
    }
  };
}