elecParticle.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 <elecPoint.hpp>
#include <elecParams.hpp>
#include <cmath>
#include <utility>
#include <limits>

namespace elec {
  class Particle {
  public:
    Point pos,speed;
    double q; // -1, 1

    
    Particle() : pos(), speed(), q(0) {}
    Particle(const Point& p, double qq) 
      : pos(p), speed(0,0), q(qq) {}

    Particle(const Particle& cp) 
      : pos(cp.pos), speed(cp.speed), q(cp.q) {}
    Particle& operator=(const Particle& cp) {
      pos = cp.pos; 
      speed = cp.speed; 
      q = cp.q;
      return *this;
    }

  };

  inline std::ostream& operator<<(std::ostream& os, 
                                  const Particle& p) {
    if(p.q > 0)
      os << "p+@";
    else
      os << "e-@";
    os << p.pos << '/' << p.speed;
    return os;
  }
  inline Particle electron(const Point& pos) {
    return {pos,-1};
  }

  inline Particle proton(const Point& pos) {
    return {pos,1};
  }

  inline Point E(const Particle& p, const Point& at) {
    double dist = d2(p.pos,at);
    if(dist < elecEPS_E_RADIUS*elecEPS_E_RADIUS)
      return Point(0,0);
    dist = std::max(dist,elecMIN_E_RADIUS*elecMIN_E_RADIUS);
    Point d = (at-p.pos);
    return (*d)*(p.q*elecFIELD_COEF/dist);
  }

  template<typename Iter>
  Point E(const Iter& begin,
          const Iter& end,
          const Point& at) {
    Point e = {0,0};
    for(auto it = begin; it != end; ++it)
      e += E(*it,at);
    return e;
  }

  template<typename Iter, typename PartOf>
  Point E_(const Iter& begin,
          const Iter& end,
           const Point& at,
           const PartOf& part_of) {
    Point e = {0,0};
    for(auto it = begin; it != end; ++it)
      e += E(part_of(*it),at);
    return e;
  }

  template<typename PartIter, typename PosIter,
           typename OutputIt>
  double E(const PartIter& particle_begin,
           const PartIter& particle_end,
           const PosIter&  position_begin,
           const PosIter&  position_end,
           OutputIt out) {
    double max_norm = 0;
    for(auto it = position_begin; it != position_end; ++it) {
      auto e =  E(particle_begin,particle_end,*it);
      double tmp = e*e;
      if(tmp > max_norm)
        max_norm = tmp;
      *(out++) = e;
    }
    return std::sqrt(max_norm);
  }

  template<typename PartIter, typename PosIter,
           typename OutputIt, typename PartOf, typename PosOf>
  void E_(const PartIter&    particle_begin,
          const PartIter&    particle_end,
          const PosIter&     position_begin,
          const PosIter&     position_end,
          OutputIt out,
          const PartOf& part_of,
          const PosOf& pos_of) {
    for(auto it = position_begin; it != position_end; ++it)
      *(out++) = E_(particle_begin,particle_end,pos_of(*it),part_of);
  }
  inline double V(const Particle& p, const Point& at) {
    double dist = d(p.pos,at);
    if(dist < elecEPS_E_RADIUS)
      return 0;
    dist = std::max(dist,elecMIN_E_RADIUS);
    return p.q*elecFIELD_COEF/dist;
  }

  template<typename Iter>
  double V(const Iter& begin,
          const Iter& end,
          const Point& at) {
    double v = 0;
    for(auto it = begin; it != end; ++it)
      v += V(*it,at);
    return v;
  }

  template<typename PartIter, typename PosIter,
           typename OutputIt>
  std::pair<double,double> V(const PartIter& particle_begin,
         const PartIter& particle_end,
         const PosIter&  position_begin,
         const PosIter&  position_end,
         OutputIt out) {
    double min = std::numeric_limits<double>::max();
    double max = std::numeric_limits<double>::lowest();
    for(auto it = position_begin; it != position_end; ++it) {
      double v = V(particle_begin,particle_end,*it);
      if(v < min) min = v;
      if(v > max) max = v;
      *(out++) = v;
    }

    return {min,max};
  }
  
}