/* -*-ePiX-*- */
#include "epix.h"
using namespace ePiX;

// number of field lines per charge
const int N(13);

// charge magnitudes and locations
const int charge1(1);
const int charge2(1);

const P Q1( 1,0);
const P Q2(-1,0);

// inverse-square electric field from a point charge at the origin
P unit_charge(const P& arg) 
{ 
  return recip(arg|arg)*arg;
}

// electric field simulated by superimposing individual fields, then
// re-scaling so that field -> 0 near the charge location (improves
// plot quality:)
P E(double x, double y)
{ 
  P temp(x,y);
  // superposition of charge fields
  P E_temp(charge1*unit_charge(temp-Q1) + charge2*unit_charge(temp-Q2));

  // re-scale
  return (1.0/(E_temp|E_temp))*E_temp;
}

P potential(double x, double y)
{
  // J rotates a vector by 1/4 turn; parallel to equipotentials
  return J(E(x,y));
}

const double MAX(3);

int main()
{
  picture(P(-MAX,-MAX), P(MAX,MAX), "4x4in");

  begin();
  set_crop();
  degrees(); // change angle mode

  // plot field lines
  blue();
  for (int i=0; i < N; ++i)
    {
      // initial points trace a small circle about Q1 or Q2,
      ode_plot(E, Q1+polar(0.05, i*360.0/N), 10, 120);
      ode_plot(E, Q2-polar(0.05, i*360.0/N), 10, 120);

      // location of arrowhead
      P pt(flow(E, Q2-polar(0.05, i*360.0/N), 3, 12));
    }

  green();
  for (int i=-10; i < 10; ++i)
    {
      ode_plot(potential, Q1+polar(0.25*pow(0.8, i), 0), 2*M_PI, 120);
      ode_plot(potential, Q2-polar(0.25*pow(0.8, i), 0), 2*M_PI, 120);
    }

  dot_size(6);
  circ(Q1);
  circ(Q2);

  magenta();
  label(Q1, "$+$");
  label(Q2, "$+$");

  blue();
  for (int i=0; i < N; ++i)
    {
      P pt(flow(E, Q2-polar(0.05, i*360.0/N), 3.5, 12));
      arrow(pt, pt+0.01*E(pt.x1(), pt.x2()));
      arrow(-pt, -pt+0.01*E(-pt.x1(), -pt.x2()));
    }

  end();
}