#include <iostream>
#include <cmath>
#include <limits>
#include "source/tessellation/VoronoiMesh.hpp"
#include "source/newtonian/two_dimensional/hdsim2d.hpp"
#include "source/newtonian/two_dimensional/interpolations/pcm2d.hpp"
#include "source/newtonian/two_dimensional/spatial_distributions/uniform2d.hpp"
#include "source/newtonian/common/ideal_gas.hpp"
#include "source/newtonian/two_dimensional/geometric_outer_boundaries/SquareBox.hpp"
#include "source/newtonian/two_dimensional/hydro_boundary_conditions/RigidWallHydro.hpp"
#include "source/newtonian/common/hllc.hpp"
#include "source/newtonian/two_dimensional/source_terms/zero_force.hpp"
#include "source/newtonian/two_dimensional/point_motions/eulerian.hpp"
#include "source/newtonian/two_dimensional/point_motions/lagrangian.hpp"
#include "source/newtonian/two_dimensional/point_motions/eulerian.hpp"
#include "source/newtonian/two_dimensional/point_motions/round_cells.hpp"
#include "source/newtonian/two_dimensional/diagnostics.hpp"
#include "source/newtonian/two_dimensional/source_terms/cylindrical_complementary.hpp"
#include "source/newtonian/two_dimensional/source_terms/CenterGravity.hpp"
#include "source/newtonian/two_dimensional/source_terms/SeveralSources.hpp"
#include "source/misc/simple_io.hpp"
#include "source/newtonian/test_2d/main_loop_2d.hpp"
#include "source/newtonian/test_2d/kill_switch.hpp"
#include "source/newtonian/two_dimensional/hdf5_diagnostics.hpp"
#include "source/misc/mesh_generator.hpp"
#include "source/tessellation/RoundGrid.hpp"
#include "source/misc/int2str.hpp"
#include "source/newtonian/two_dimensional/interpolations/linear_gauss_consistent.hpp"
#include "source/newtonian/test_2d/consecutive_snapshots.hpp"
#include <boost/foreach.hpp>
#include "source/newtonian/two_dimensional/hydro_boundary_conditions/FreeFlow.hpp"
#include "source/misc/utils.hpp"
#include "source/tessellation/shape_2d.hpp"
#include "source/newtonian/test_2d/piecewise.hpp"
#include "source/mpi/MeshPointsMPI.hpp"
#include "source/mpi/mpi_macro.hpp"
#include "source/mpi/ConstNumberPerProc.hpp"
#include "source/misc/hdf5_utils.hpp"
#include "source/newtonian/test_2d/contour.hpp"
#include "source/newtonian/two_dimensional/custom_evolutions/Ratchet.cpp"
#include "source/newtonian/two_dimensional/diagnostics.hpp"
#include <fenv.h>

using namespace std;
using namespace simulation2d;

namespace {

  class WriteConserved: public DiagnosticFunction
  {
  public:

    WriteConserved(string const& fname):
      tracer_(), cons_(), fname_(fname) {}

    void operator()(hdsim const& sim)
    {
      tracer_.push_back(total_tracer(sim,0));
      cons_.push_back(total_conserved(sim));
      time_.push_back(sim.GetTime());
    }

    ~WriteConserved(void)
    {
	ofstream f(fname_.c_str());
	for(size_t i=0;i<cons_.size();++i)
	  f << time_[i] << " "
	    << cons_[i].Mass << " "
	    << cons_[i].Momentum.x << " "
	    << cons_[i].Momentum.y << " "
	    << cons_[i].Energy << " "
	    << tracer_[i] << "\n";
	f.close();
    }

  private:
    mutable vector<double> tracer_;
    mutable vector<Conserved> cons_;
    mutable vector<double> time_;
    const string fname_;
  };

  /*
    class Veil: public Acceleration
    {
    public:

    Veil(Acceleration& full):
    full_(full) {}

    Vector2D Calculate
    (Tessellation const &  tess,
    vector< Primitive > const &  cells,
    int  point,
    vector< Conserved > const &  fluxes,
    vector< Vector2D > const &  point_velocity,
    HydroBoundaryConditions const &  hbc,
    vector< vector< double > > const &  tracers,
    double  time,
    double  dt)
    {
    const Primitive& cell = cells[static_cast<size_t>(point)];
    if(abs(cell.Velocity.x)<std::numeric_limits<double>::epsilon()&&
    abs(cell.Velocity.y)<std::numeric_limits<double>::epsilon())
    return Vector2D(0,0);
    else
    return full_.Calculate(tess,cells,point,fluxes,point_velocity,hbc,tracers,
    time,dt);
    }

    private:
    Acceleration& full_;
    };*/

  double calc_mach_number(const Primitive& p)
  {
    return abs(p.Velocity)/p.SoundSpeed;
  }

  Vector2D find_pos_on_line(const Vector2D& p1,
			    double v1,
			    const Vector2D& p2,
			    double v2,
			    double vt)
  {
    return p1+((vt-v1)/(v2-v1))*(p2-p1);
  }

  /*
    vector<Vector2D> process_positions(const Index2Member<Vector2D>& gg)
    {
    vector<Vector2D> res(get_mpi_size());
    for(size_t i=0;i<res.size();++i)
    res[i] = gg(i*gg.getLength()/get_mpi_size());

    return res;
    }
  */

  class Constants
  {
  public:

    // Metric prefix
    const double kilo;
    const double centi;

    // Length / distance
    const double parsec;
    const double meter;

    // Time
    const double year;
    const double second;

    // Mass
    const double solar_mass;
    const double gram;

    // Energy
    const double erg;

    // Force
    const double newton;

    // Parameters
    const double adiabatic_index;
    const double gravitation_constant;
    const double rho_0;
    const double R_0;
    const double omega_in;
    const double inner_density_prefactor;
    const double R_b;
    const double omega_out;
    const double outer_density_prefactor;
    const double offset;
    const double supernova_energy;
    const double supernova_radius;
    const double supernova_volume;
    const double supernova_mass;
    const double supernova_density;
    const double supernova_pressure;
    const double black_hole_mass;
    const Vector2D lower_left;
    const Vector2D upper_right;

    Constants(void):
      kilo(1e3),
      centi(1e-2),
      parsec(1),
      meter(3.24e-17*parsec),
      year(1),
      second(year/3.16e7),
      solar_mass(1),
      gram(5.03e-34*solar_mass),
      erg(gram*pow(centi*meter/second,2)),
      newton(kilo*gram*meter/pow(second,2)),
      adiabatic_index(5./3.),
      gravitation_constant(6.673e-11*newton*pow(meter/(kilo*gram),2)),
      rho_0(2.2e-22*gram/pow(centi*meter,3)),
      R_0(0.04*parsec),
      omega_in(1),
      inner_density_prefactor(rho_0*pow(R_0,omega_in)),
      R_b(0.4*parsec),
      omega_out(3),
      outer_density_prefactor
      (inner_density_prefactor*pow(R_b,omega_out-omega_in)),
      offset(1*parsec),
      supernova_energy(1e51*erg),
      supernova_radius(0.2*parsec),
      supernova_volume((4.*M_PI/3)*pow(supernova_radius,3)),
      supernova_mass(5*solar_mass),
      supernova_density(supernova_mass/supernova_volume),
      supernova_pressure((adiabatic_index-1)*supernova_energy/supernova_volume),
      black_hole_mass(1e9*solar_mass),
      lower_left(parsec*Vector2D(0,-6)),
      upper_right(parsec*Vector2D(6,3)) {}
  };
}

namespace units
{
  // Time
  const double year = 1;
  //  const double second = year/3.16e7;

  // parameters
  //  const double gravitation_constant = 6.673e-11*newton*pow(meter/(kilo*gram),2);
  //  const double black_hole_mass = 1e7*solar_mass;
  const double adiabatic_index = 5./3.;
  const double omega_in = 1;
  const double omega_out = 3;
  //  const double R_0 = 0.04*parsec;
  //  const double rho_0 = 2.2e-22*gram/pow(centi*meter,3);
  //  const double R_b = 0.4*parsec;
  //  const double inner_density_prefactor = rho_0*pow(R_0,omega_in);
  /*
    const double outer_density_prefactor = 
    inner_density_prefactor*pow(R_b,omega_out-omega_in);
  */
  //  const double supernova_radius = 0.2*parsec;
  //  const double supernova_energy = 1e51*erg;
  //  const double supernova_mass = 5*solar_mass;
  //  const double supernova_volume = (4.*M_PI/3.)*pow(supernova_radius,3);
  //const double supernova_density = supernova_mass / supernova_volume;
  //const double supernova_pressure = 
  //  (adiabatic_index-1)*(supernova_energy/supernova_volume);
  //const double ambient_pressure = supernova_pressure*1e-10;
  //  const double ambient_pressure = 1e-30;
  //  const double offset = 1*parsec;
  /*
    const Vector2D lower_left = parsec*Vector2D(1e-3,-6);
    const Vector2D upper_right = parsec*Vector2D(6,6);
  */
  //  const Vector2D lower_left = parsec*Vector2D(0,-6);
  //  const Vector2D upper_right = parsec*Vector2D(6,2);
}

using namespace units;

namespace {

  class HydrostaticPressure: public SpatialDistribution
  {
  public:

    HydrostaticPressure(double gm,
			double r0,
			double rho_0,
			double omega_1,
			double R_b,
			double omega_2):
      gm_(gm), r0_(r0), rho_0_(rho_0),
      omega_1_(omega_1), R_b_(R_b), omega_2_(omega_2) {}

    double operator()(const Vector2D& r) const
    {
      if(abs(r)<R_b_)
	return (gm_/r0_)*(rho_0_/(omega_1_+1))*pow(abs(r)/r0_,-omega_1_-1)-
	  (gm_/r0_)*(rho_0_/(omega_1_+1))*pow(R_b_/r0_,-omega_1_-1)+
	  (gm_/r0_)*(rho_0_/(omega_2_+1))*pow(R_b_/r0_,-omega_2_-1);
      else
	return (gm_/r0_)*(rho_0_/(omega_2_+1))*
	  pow(R_b_/r0_,-omega_1_-1)*
	  pow(abs(r)/R_b_,-omega_2_-1);
    }

  private:
    const double gm_;
    const double r0_;
    const double rho_0_;
    const double omega_1_;
    const double R_b_;
    const double omega_2_;
  };

  class BoundedRadialPowerLaw: public SpatialDistribution
  {
  public:

    BoundedRadialPowerLaw(const Vector2D& center,
			  double index,
			  double prefactor,
			  double lower_bound,
			  double upper_bound):
      center_(center),
      index_(index),
      prefactor_(prefactor),
      lower_bound_(lower_bound),
      upper_bound_(upper_bound) {}

    double operator()(const Vector2D& r) const
    {
      const double radius = max(min(abs(r-center_),
				    upper_bound_),
				lower_bound_);
      return prefactor_*pow(radius,index_);
    }

  private:
    const Vector2D center_;
    const double index_;
    const double prefactor_;
    const double lower_bound_;
    const double upper_bound_;
  };

  bool is_inside_rectangle(const Vector2D& point,
			   const Vector2D& lower_left,
			   const Vector2D& upper_right)
  {
    return ((point.x > lower_left.x) &&
	    (point.x < upper_right.x) &&
	    (point.y > lower_left.y) &&
	    (point.y < upper_right.y));
  }

  vector<Vector2D> rectangle_clip(const vector<Vector2D>& grid,
				  const Vector2D& lower_left,
				  const Vector2D& upper_right)
  {
    vector<Vector2D> res;
    for(size_t i=0, endp=grid.size();i<endp;++i){
      const Vector2D& point = grid[i];
      if(is_inside_rectangle(point,lower_left,upper_right))
	res.push_back(point);	 
    }
    return res;
  }

  /*
    vector<double> calc_radius_list(void)
    {
    const double rmin = 1e-4;
    const double q = 1.01;
    const size_t n = 200;
    vector<double> res(n);
    for(size_t i=0;i<n;++i)
    res[i] = rmin*pow(q,i);
    write_number(0,"finished_calc_radius_list.txt");
    return res;
    }
  */

  /*
    vector<Vector2D> create_grid(Vector2D const& lower_left,
    Vector2D const& upper_right,
    double dx2x)
    {
    vector<Vector2D> res;
    for(double x = lower_left.x*(1+0.5*dx2x);
    x<upper_right.x; x*=(1+dx2x)){
    const double dx = x*dx2x;
    for(double y=lower_left.y+dx/2;
    y<upper_right.y; y += dx)
    res.push_back(Vector2D(x,y));
    }
    return res;
    }
  */

  /*
    vector<Vector2D> centered_hexagonal_grid(double r_min,
    double r_max)
    {
    const vector<double> r_list = arange(0,r_max,r_min);
    vector<Vector2D> res;
    for(size_t i=0;i<r_list.size();++i){
    const size_t angle_num = max<size_t>(6*i,1);
    vector<double> angle_list(angle_num,0);
    for(size_t j=0;j<angle_num;++j)
    angle_list.at(j) = 2*M_PI*(double)j/(double)(angle_num);
    for(size_t j=0;j<angle_num;++j)
    res.push_back(r_list.at(i)*Vector2D(cos(angle_list.at(j)),
    sin(angle_list.at(j))));
    }
    return res;
    }
  */

  vector<Vector2D> centered_logarithmic_spiral(double r_min,
					       double r_max,
					       double alpha,
					       const Vector2D& center)
  {
    const double theta_max = log(r_max/r_min)/alpha;
    const vector<double> theta_list = 
      arange(0,theta_max,2*M_PI*alpha);
  
    vector<double> r_list(theta_list.size(),0);
    for(size_t i=0;i<r_list.size();++i)
      r_list.at(i) = r_min*exp(alpha*theta_list.at(i));
  
    vector<Vector2D> res(r_list.size());
    for(size_t i=0;i<res.size();++i)
      res[i] = center+r_list[i]*Vector2D(cos(theta_list.at(i)),
					 sin(theta_list.at(i)));
    return res;
  }

  /*
    vector<Vector2D> complete_grid(double r_inner,
    double r_outer,
    double alpha)
    {
    const vector<Vector2D> inner = 
    centered_hexagonal_grid(r_inner*alpha*2*M_PI,
    r_inner);
    const vector<Vector2D> outer =
    centered_logarithmic_spiral(r_inner,
    r_outer,
    alpha);
    return join(inner, outer);
    }
  */

  template<class T> class VectorInitializer
  {
  public:

    VectorInitializer(const T& t):
      buf_(1,t) {}

    VectorInitializer& operator()(const T& t)
    {
      buf_.push_back(t);
      return *this;
    }

    vector<T> operator()(void)
    {
      return buf_;
    }

  private:
    vector<T> buf_;
  };
}

class SimData
{

public:

  SimData(const Constants& c):
    pg_(Vector2D(0,0), Vector2D(0,1)),
    outer_(c.lower_left, c.upper_right),
    //    init_points_(cartesian_mesh(300,450,c.lower_left,c.upper_right)),
    //cgg_(11+2*150,7+2*300,lower_left, upper_right),
    init_points_(rectangle_clip
		 (centered_logarithmic_spiral(0.01,
					      abs(c.upper_right-c.lower_left),
					      0.002,
					      Vector2D(-0.2,-1)),
		  c.lower_left, c.upper_right)),
    /*
      init_points_(create_grid(lower_left,
      upper_right,
      0.1)),
    */
    tess_(),
    //proc_tess_(process_positions(Echo<Vector2D>(init_points_)),outer_),
    eos_(adiabatic_index),
    rs_(),
    hbc_(rs_),
    //interpm_(eos_,outer_,hbc_,true,false),
    interpm_(),
    raw_point_motion_(),
    point_motion_(raw_point_motion_,hbc_,
		  0.15, 0.02, false,0,&outer_),
    alt_point_motion_(),
    gravity_acc_(c.gravitation_constant*c.black_hole_mass,
		 0.1*c.parsec,
		 c.parsec*Vector2D(-0.05,0)),
  //    veil_(gravity_acc_),
    gravity_force_(gravity_acc_),
    geom_force_(pg_.getAxis()),
    force_(VectorInitializer<SourceTerm*>(&gravity_force_)
	   (&geom_force_)()),
    ce_(Ratchet::in),
    sim_(init_points_, //distribute_grid(proc_tess_,Echo<Vector2D>(init_points_)),
	 tess_,
	 //	 proc_tess_,
	 interpm_,
	 Piecewise
	 (Circle(Vector2D(0,-c.offset),c.supernova_radius),
	  Uniform2D(c.supernova_density),
	  Piecewise
	  (Circle(Vector2D(0,0),c.R_b),
	   BoundedRadialPowerLaw(Vector2D(0,0),
				 -omega_in,
				 c.inner_density_prefactor,
				 1e-6, 10),
	   BoundedRadialPowerLaw(Vector2D(0,0),
				 -omega_out,
				 c.outer_density_prefactor,
				 1e-6, 10))),
	 Piecewise(Circle(Vector2D(0,-c.offset),c.supernova_radius),
		   Uniform2D(c.supernova_pressure),
		   HydrostaticPressure(c.gravitation_constant*c.black_hole_mass,
				       c.R_0,
				       c.rho_0,
				       c.omega_in,
				       c.R_b,
				       c.omega_out)),
	 Uniform2D(0),
	 Uniform2D(0),
	 eos_,
	 rs_,
	 alt_point_motion_,
	 //	 point_motion_,
	 force_,
	 outer_,
	 hbc_) 
  {
    sim_.SetColdFlows(0.01,0.01);
    sim_.SetDensityFloor(1e-6,1e-23);
    sim_.changePhysicalGeometry(&pg_);
    sim_.custom_evolution_manager.addCustomEvolution(&ce_);
    for(int i=0;i<sim_.GetCellNo();++i){
      if(abs(sim_.GetMeshPoint(i))<0.02*c.parsec)
	sim_.custom_evolution_indices[static_cast<size_t>(i)] = 1;
    }
  }

  hdsim& getSim(void)
  {
    return sim_;
  }
  
private:
  const CylindricalSymmetry pg_;
  const SquareBox outer_;
  //  const CartesianGridGenerator cgg_;
  const vector<Vector2D> init_points_;
  VoronoiMesh tess_;
  VoronoiMesh proc_tess_;
  const IdealGas eos_;
  const Hllc rs_;
  const RigidWallHydro hbc_;
  //const FreeFlow hbc_;
  //  LinearGaussConsistent interpm_;
  PCM2D interpm_;
  Lagrangian raw_point_motion_;
  RoundCells point_motion_;
  Eulerian alt_point_motion_;
  CenterGravity gravity_acc_;
  //  Veil veil_;
  ConservativeForce gravity_force_;
  CylindricalComplementary geom_force_;
  SeveralSources force_;
  Ratchet ce_;
  hdsim sim_;
};

namespace {
  /*
    class HighestShockedPoint: public TerminationCondition
    {
    public:

    HighestShockedPoint(double position):
    position_(position) {}

    bool operator()(hdsim const& sim)
    {
    const double mn_thres = 0.01;
    double ys = -10;
    for(int i=0;i<sim.GetCellNo();++i){
    if(sim.GetCell(i).Velocity.y/sim.GetCell(i).SoundSpeed>mn_thres){
    ys = max(ys, sim.GetMeshPoint(i).y);
    }
    }
    cout << ys << endl;
    return ys < position_;
    }

    private:
    const double position_;
    };
  */
}

class CustomDiagnostic
{
public:

  CustomDiagnostic(double dt,
		   double t_start):
    dt_(dt),
    t_next_(t_start),
    counter_(0) {}

  void operator()(const hdsim& sim)
  {
    if(sim.GetTime()>t_next_){
      write_snapshot_to_hdf5(sim, "snapshot_"+int2str(counter_)+".h5");

      t_next_ += dt_;
      ++counter_;
    }
  }

private:
  const double dt_;
  double t_next_;
  int counter_;
};

namespace {
  /*
    vector<double> calc_decade_vector(double val, size_t decades)
    {
    vector<double> res;
    for(size_t i=0;i<decades;++i)
    res.push_back(val/pow(2.0,static_cast<double>(i)+1));
    return res;
    }
  */
}

namespace {

  /*
    double highest_shocked_point(const hdsim& sim)
    {
    const double mn_thres = 0.01;
    double ys = -10;
    for(int i=0;i<sim.GetCellNo();++i){
    if(sim.GetCell(i).Velocity.y/sim.GetCell(i).SoundSpeed>mn_thres){
    ys = max(ys,sim.GetMeshPoint(i).y);
    }
    }
    return ys;
    }
  */
}

namespace {
  /*
    class MyDiagnostic: public DiagnosticFunction
    {
    public:

    MyDiagnostic(double y_start,
    size_t decades):
    snapshot_times_(calc_decade_vector(y_start, decades)),
    y_next_(y_start), counter_(0), active_(true) {}

    void operator()(const hdsim& sim)
    {
    if(!active_)
    return;

    if(highest_shocked_point(sim)>y_next_){
    write_snapshot_to_hdf5(sim,"snapshot_"+int2str(counter_)+".h5");
    if(static_cast<size_t>(counter_)>snapshot_times_.size()-1)
    active_ = false;
    else{
    y_next_ = snapshot_times_.at(static_cast<size_t>(counter_));
    ++counter_;
    }
    }
    }

    private:
    const vector<double> snapshot_times_;
    mutable double y_next_;
    mutable int counter_;
    mutable bool active_;
    };
  */
}

namespace {
  class MultipleDiagnostics: public DiagnosticFunction
  {
  public:

    MultipleDiagnostics(const vector<DiagnosticFunction*>& diag_list):
      diag_list_(diag_list) {}

    void operator()(const hdsim& sim)
    {
      BOOST_FOREACH(DiagnosticFunction* df, diag_list_)
	{
	  (*df)(sim);
	}
    }

  private:
    const vector<DiagnosticFunction*> diag_list_;
  };
}

namespace {

  class ShockFrontDetector: public LocalContourCriterion
  {
  public:
    ShockFrontDetector(const Vector2D& center,
		       double mn_thres):
      center_(center), mn_thres_(mn_thres) {}

    std::pair<bool,Vector2D> operator()(const Edge& edge,
					const hdsim& sim) const
    {
      if(sim.getHydroBoundaryCondition().IsBoundary(edge,
						    sim.GetTessellation()))
	return std::pair<bool,Vector2D>(false,Vector2D(0,0));
      if(edge.neighbors.first>sim.GetCellNo() ||
	 edge.neighbors.second>sim.GetCellNo())
	return std::pair<bool,Vector2D>(false,Vector2D(0,0));
      const std::pair<double,double> mach_numbers
	(calc_mach_number(sim.GetCell(edge.neighbors.first)),
	 calc_mach_number(sim.GetCell(edge.neighbors.second)));
      if((mach_numbers.first-mn_thres_)*(mach_numbers.second-mn_thres_)>0 ||
	 (mach_numbers.second-mach_numbers.first)*
	 (abs(sim.GetMeshPoint(edge.neighbors.second)-center_)-
	  abs(sim.GetMeshPoint(edge.neighbors.first)-center_))>0)
	return std::pair<bool,Vector2D>(false,Vector2D(0,0));
      else
	return std::pair<bool,Vector2D>
	  (true,find_pos_on_line(sim.GetMeshPoint(edge.neighbors.first),
				 mach_numbers.first,
				 sim.GetMeshPoint(edge.neighbors.second),
				 mach_numbers.second,
				 mn_thres_));
						
    }

  private:
    const Vector2D center_;
    const double mn_thres_;
  };

  /*
    class CellEdgesGetter: public Index2Member<Edge>
    {
    public:
    CellEdgesGetter(const Tessellation& tess,
    const int n):
    tess_(tess), edge_indices_(tess_.GetCellEdges(n)) {}

    size_t getLength(void) const
    {
    return edge_indices_.size();
    }

    Edge operator()(size_t i) const
    {
    return tess_.GetEdge(edge_indices_[i]);
    }

    private:
    const Tessellation& tess_;
    const vector<int> edge_indices_;
    };
  */

  /*
    class EnforceRigidWall: public Manipulate
    {
    public:

    EnforceRigidWall(void) {}

    void operator()(hdsim& sim)
    {
    const PhysicalGeometry& pg = sim.getPhysicalGeometry();
    const Tessellation& tess = sim.GetTessellation();
    for(size_t i=0;i<sim.GetAllCells().size();++i){
    if(tess.GetMeshPoint(i).x-lower_left.x<tess.GetWidth(i)){
    sim.getAllConserved()[i].Momentum.x = 0;
    const double volume = 
    pg.calcVolume(serial_generate(CellEdgesGetter(tess,i)));
    sim.GetAllCells()[i] = Conserved2Primitive
    (sim.getAllConserved()[i]/volume,sim.getEos());
    }
    }
    }
    };
  */

  void my_main_loop(hdsim& sim, const Constants& c)
  {
    //  sim.SetCfl(0.1);
    const double tf = 600*c.year;
    SafeTimeTermination term_cond_raw(tf,1e6);
    ConsecutiveSnapshots diag1(new ConstantTimeInterval(10*units::year),
			       new Rubric("snapshot_",".h5"));
    WriteTime diag2("time.txt");
    SequentialContour diag3(new ConstantTimeInterval(10*units::year),
			    new Rubric("contour_",".h5"),
			    new ShockFrontDetector(Vector2D(0,-c.offset),0.1));
    WriteConserved diag4("conserved.txt");
    vector<DiagnosticFunction*> diag_list;
    diag_list.push_back(&diag1);
    diag_list.push_back(&diag2);
    diag_list.push_back(&diag3);
    diag_list.push_back(&diag4);
    MultipleDiagnostics diag(diag_list);

    main_loop(sim, 
	      term_cond_raw,
	      &hdsim::TimeAdvance, 
	      &diag);
  }
}

namespace {
  void report_error(UniversalError const& eo)
  {
    cout << eo.GetErrorMessage() << endl;
    cout.precision(14);
    for(size_t i=0;i<eo.GetFields().size();++i)
      cout << eo.GetFields()[i] << " = "
	   << eo.GetValues()[i] << endl;
  }
}

int main(void)
{
  feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
  //  MPI_Init(NULL, NULL);

  try{

    Constants c;

    SimData sim_data(c);
    hdsim& sim = sim_data.getSim();

    HDF5Shortcut("units.h5")
      ("year",vector<double>(1,c.year))
      ("solar mass",vector<double>(1,c.solar_mass))
      ("parsec",vector<double>(1,c.parsec));

    write_snapshot_to_hdf5(sim,
			   "initial.h5");

    my_main_loop(sim,c);

    write_snapshot_to_hdf5(sim, 
			   "final.h5");

  }
  catch(const UniversalError& eo){
    report_error(eo);
    throw;
  }

  return 0;
}