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

double f(double t, double a) { return 3*Sin(a*t)*exp(-0.5*t); }

// Completely fabricated data
// "Predicted"
inline double f1(double t) { return f(t,2); }
// "Observed"
inline double f2(double t) { return f1(t)+0.05*Cos(3*t); }
//"error"
inline double err(double t) { return 0.125*(1+fabs(f1(t))); }

// Number of error bars
const int samples=10;

void error_bar(P center, double error)
{
  P top    = center + P(0,error);
  P bottom = center - P(0,error);

  line(bottom,top);
  v_axis_tick(top);
  v_axis_tick(bottom);
}

// Half-size of bounding box
const int N=2;

main() {
  bounding_box(P(0,-N),P(2*N,N));
  unitlength("0.005in");
  picture(400*N,400*N);

  begin();

  double x_avg = 0.5*(x_min + x_max);
  // displacement for miscellaneous labels/objects
  const double delta=0.0625*x_size;

  // legend label locations
  P predicted=P(x_avg + delta, y_min + 2*delta);
  P observed =P(x_avg + delta, y_min +   delta);

  // graph paper
  grid(4*N,4*N);
  pen(0.15);
  grid(20*N,20*N);

  // "Predicted" curve
  rgb(0.2, 0.7, 0.8);
  bold();
  plot(f1, x_min, x_max, 120);

  // legend
  line(P(x_avg,y_min+1.9*delta), P(x_avg+delta,y_min+1.9*delta));

  black();
  masklabel(predicted, P(2,0), "Predicted", r);
  masklabel(observed,  P(2,0), "Observed" , r);

  // axes, names, and labels
  dart(P(x_min,0), P(x_max+delta,0));
  dart(P(0,y_min), P(0,y_max+delta));

  masklabel(P(x_max+delta,0), P(2,0), "Time (sec)", r);

  h_axis_masklabels(P(1,0), P(x_max,0), 3, P(0,4), t);
  v_axis_labels(y_size, P(-4, 0), l);

  plain();
  red();
  // legend
  error_bar(observed, 0.5*delta);

  // "error bars"
  double x=x_min;
  for (int i=1; i <= samples; ++i, x+=x_size/samples)
    error_bar(P(x, f2(x)), err(x));

  // legend
  circ(observed);

  // "data points"
  x=x_min;
  for (int i=1; i <= samples; ++i, x+=x_size/samples)
    circ(P(x,f2(x)));

  // figure title
  black();
  masklabel(P(x_avg, y_max), P(0,2), "A damped harmonic oscillator", t);

  degrees();
  label_angle(90);
  masklabel(P(0,0), P(-12,0), "Displacement (cm)",l);

  end();
}