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

#define COLOR
#define Y_IN
// #define Y_OUT

// Notation: kappa -> k

const double a = 10;
const double M = 1;
const double dt = 0.00001;

double phi(double r) { return 1 - (2*M/r) - r*r/(a*a); }

static inline double diff(double f(double), double t)
{
  return (f(t+dt) - f(t-dt))/(2*dt);
}

const double r1 = newton(phi,  2);
const double r2 = newton(phi, 10);

inline double dphi(double r, double k)
{
  return diff(phi, r)/(2*k);
}

inline double dy(double r, double k)
{
  return (1 - dphi(r,k)*dphi(r,k))/phi(r);
}

// for Yin
#ifdef Y_IN
const double k = 0.5*diff(phi, r1);

const double residue1 = sqrt(EPIX_EPSILON*dy(r1+EPIX_EPSILON, k));
const double residue2 = sqrt(EPIX_EPSILON*dy(r2-EPIX_EPSILON, k));

double dy(double r) 
{ 
  if ((r-r1)*(r2-r) == 0)
    return 0;
  else
    return sqrt(fabs(dy(r, k))) - 
      residue1/sqrt(fabs(r-r1)) - 
      residue2/sqrt(fabs(r-r2)); 
}

double Y(double r, double t)
{
  const int steps = 60;

  double x = r1;
  const double dx = (r-r1)/steps;

  // contribution from improper terms
  double sum = 2*(residue1*sqrt(fabs(r-r1)) - 
		  residue2*(sqrt(fabs(r2-r)) - sqrt(r2-r1)));

  // Simpson's rule
  for (int i=0; i <= steps; ++i, x += dx)
    sum +=  (1.0/6)*(dy(x) + 4*dy(x+0.5*dx) + dy(x+dx))*dx;

  return sum;
}
#endif // Y_IN

// for Yout
#ifdef Y_OUT
const double k = 0.5*diff(phi, r2);

const double residue1 = sqrt(EPIX_EPSILON*dy(r1+EPIX_EPSILON, k));
const double residue2 = sqrt(EPIX_EPSILON*dy(r2-EPIX_EPSILON, k));

double dy(double r) 
{ 
  if ((r-r1)*(r2-r) == 0)
    return 0;
  else
    return sqrt(fabs(dy(r, k))) - 
      residue1/sqrt(fabs(r-r1)) - 
      residue2/sqrt(fabs(r-r2)); 
}

double Y(double r, double t)
{
  const int steps = 60;

  double x = r2;
  const double dx = (r2-r)/steps;

  // contribution from improper terms
  double sum = 2*(residue1*sqrt(fabs(r2-r)) - 
		  residue2*(sqrt(fabs(r-r1)) - sqrt(r2-r1)));

  // Simpson's rule
  for (int i=0; i <= steps; ++i, x -= dx)
    sum +=  (1.0/6)*(dy(x) + 4*dy(x+0.5*dx) + dy(x+dx))*dx;

  return sum;
}

#endif

double X1p(double r, double t) { return  (1/k)*sqrt(phi(r))*cosh(k*t); }
double X1m(double r, double t) { return -(1/k)*sqrt(phi(r))*cosh(k*t); }
double T1p(double r, double t) { return  (1/k)*sqrt(phi(r))*sinh(k*t); }

double X2(double r, double t)  { return  (1/k)*sqrt(-phi(r))*sinh(k*t); }
double T2p(double r, double t) { return  (1/k)*sqrt(-phi(r))*cosh(k*t); }
double T2m(double r, double t) { return -(1/k)*sqrt(-phi(r))*cosh(k*t); }

triple Phi1(double s, double t) { return P(Y(s,t), X1p(s,t), T1p(s,t)); }
triple Phi2(double s, double t) { return P(Y(s,t), X2(s,t), T2p(s,t)); }
triple Phi3(double s, double t) { return P(Y(s,t), X1m(s,t), T1p(s,t)); }
triple Phi4(double s, double t) { return P(Y(s,t), X2(s,t), T2m(s,t)); }

const int coarse=32;
const int fine=64;
const double Y_MAX = 30;

#ifdef Y_IN
const double r2_my = 9;
#endif

#ifdef Y_OUT
const double r2_my = r2;
#endif

main() {
  bounding_box(P(-20,-20),P(20,20));
  unitlength("1pt");
  picture(P(360,360));

  //  viewpoint(6,2,5); // (Y,X,T)
  //  viewpoint(2.5,3,4); // frame on right
  //  viewpoint(3,2.5,4); // frame on left

  begin();

  clip_box(P(-10,-12,-10), P(Y_MAX,12,10));

  viewpoint(3,2.5,4);
  const double dist = 20;
  camera.range(3*dist);

  bold();

  degrees();
  // rin2
#ifdef COLOR
  rgb(0.6, 0, 0.6);
#endif
  clipmesh(Phi2, P(0,-20), P(r1,20), Net(coarse,coarse), Net(fine,fine));

  // rin4
#ifdef COLOR
  rgb(0.2, 0, 0.4);
#endif
  clipmesh(Phi4, P(0,-20), P(r1,20), Net(coarse,coarse), Net(fine,fine));

  plain();
  arrow(P(0,0,0), P(dist,0,0));
  label((1+dist)*E_1, "$Y$");

  arrow(P(0,0,0), P(0,dist,0));
  label((1+dist)*E_2, "$X$");

  arrow(P(0,0,0), P(0,0,dist));
  label((1+dist)*E_3, "$T$");

  bold();

  // draw bottom half
  clip_box(P(-10,-12,-10), P(Y_MAX,12,0));
  // rin3
#ifdef COLOR
  rgb(0.8, 0, 0.4);
#endif
  clipmesh(Phi3, P(r1,-20), P(r2_my,20), Net(coarse,coarse), Net(fine,fine));

  // rin1
#ifdef COLOR
  rgb(0.4, 0, 0.8);
#endif
  clipmesh(Phi1, P(r1,-20), P(r2_my,20), Net(coarse,coarse), Net(fine,fine));

  // erase top half
  clip_box(P(-10,-12,0), P(Y_MAX,12,10));
  white();
  pen(3);
  // rin3
  clipmesh(Phi3, P(r1,-20), P(r2_my,20), Net(coarse,coarse), Net(fine,fine));
  // rin1
  clipmesh(Phi1, P(r1,-20), P(r2_my,20), Net(coarse,coarse), Net(fine,fine));

  bold();
#ifndef COLOR
  black();
#endif

  // draw top half
  // rin3
#ifdef COLOR
  rgb(0.8, 0, 0.4);
#endif
  clipmesh(Phi3, P(r1,-20), P(r2_my,20), Net(coarse,coarse), Net(fine,fine));

  // rin1
#ifdef COLOR
  rgb(0.4, 0, 0.8);
#endif
  clipmesh(Phi1, P(r1,-20), P(r2_my,20), Net(coarse,coarse), Net(fine,fine));

  end();
}