// Visualize the electric field of a dipole,
// and the motion of a test charge in that field.
// We write a file: test-dipole.png


// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program 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 program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "ephi.hpp"

int main (int argc, char *argv[])
{
  Statics statics;
  ElectroDynamics dyn(statics);
  int nobj(0);

#define TIME_TO_RUN .9e-3
#define PIXEL_INTERVAL 0.0000001
#define DIST 0.2

  dyn.setDTM (1, .01);   // speed things up a bit
  prec_t eps(.005);
  prec_t heavy(1e10*ELECTRON_MASS);
  prec_t fixed_q(50.*ELEMENTARY_CHARGE);

// Split each fixed charge so that half is above the plane
// and half is below the plane.
// This is a kludge, necessary because I do not know any
// better way to limit the max field, and control how detailed
// the colormap will be.
  dyn.inject_particle (heavy, fixed_q/2.,
        vect3d(-DIST, 0, eps), vect3d(1, 0, 0), 0);
  int left_pole = nobj++;
  dyn.inject_particle (heavy, fixed_q/2.,
        vect3d(-DIST, 0, -eps), vect3d(1, 0, 0), 0);
  nobj++;
  dyn.inject_particle (heavy, -fixed_q/2.,
        vect3d(DIST, 0, eps), vect3d(1, 0, 0), 0);
  nobj++;
  dyn.inject_particle (heavy, -fixed_q/2.,
        vect3d(DIST, 0, -eps), vect3d(1, 0, 0), 0);
  int right_pole = nobj++;

  fprintf (stderr, "Calculating...\n");
  prec_t pixworth = .0011;      // size of a pixel in analog SI units
  Screen screen(800, 600, pixworth*800, pixworth*600);
  Scene scene(dyn);
  scene.calc (screen, Scene::CALC_EFIELD);
  scene.set_coloring (0, true); // logarithmic, full colors
  scene.render_lic (screen);

  dyn.inject_particle (ELECTRON_MASS, -50.*ELEMENTARY_CHARGE,
        vect3d(0.02*DIST, -DIST, 0), vect3d(1, 0, 0), 0);
  int test_particle = nobj++;

  prec_t cp = 0;
  int count(0);
  vect3d const dx(pixworth, 0, 0);
  vect3d const dy(0, pixworth, 0);

  while ( true )
    {
      if ( dyn.get_elapsed_time() >= cp )
        {
          Screen::color tint = Screen::RED;
// Change color after each 100 timesteps:
          if ((count/100) & 1) tint = Screen::BLUE;
          vect3d where = dyn.get_pos(test_particle);
          screen.set_pixel (where, tint);
          screen.set_pixel (where + dx, tint);
          screen.set_pixel (where + dy, tint);
          cp += PIXEL_INTERVAL;
          count++;
        }

      if ( dyn.get_elapsed_time() >= TIME_TO_RUN )
          break;
      dyn.step ();
    }

  char buf[128];
  printf ("LEFT=%s\n", dyn.get_pos(left_pole).sprint(buf));
  printf ("RIGHT=%s\n", dyn.get_pos(right_pole).sprint(buf));
  printf ("TEST=%s\n", dyn.get_pos(test_particle).sprint(buf));

  screen.write ("test-dipole.png");
  return 0;
}