#include <fstream>
#include <iomanip>
#include <limits>

#include <multimeshvectorfield.hpp>
#include <dxf_solid.hpp>
#include <gtkplotter.hpp>
#include <geomplotter.hpp>
#include <geometry.hpp>
#include <func_solid.hpp>
#include <epot_bicgstabsolver.hpp>
#include <epot_gssolver.hpp>
#include <epot_efield.hpp>
#include <error.hpp>
#include <ibsimu.hpp>
#include <mydxffile.hpp>
#include <trajectorydiagnostics.hpp>
#include <particledatabase.hpp>
#include <particlediagplotter.hpp>
#include <transformation.hpp>
#include <compmath.hpp>
#include <fielddiagplotter.hpp>
#include <random.hpp>



void simu( int argc, char **argv )
{
    std::ifstream is_geom( argv[1] );
    if( !is_geom.good() )
	throw( Error( ERROR_LOCATION, (std::string)"couldn\'t open file \'" + argv[1] + "\'" ) );
    Geometry geom( is_geom );
    is_geom.close();
 
    EpotField *epot = NULL;
    if( argc > 2 ) {
	std::ifstream is_epot( argv[2] );
	if( !is_epot.good() )
	    throw( Error( ERROR_LOCATION, (std::string)"couldn\'t open file \'" + argv[2] + "\'" ) );
	epot = new EpotField( is_epot, geom );
	is_epot.close();
    }

    ParticleDataBase3D *pdb = NULL;
    if( argc > 3 ) {
	std::ifstream is_pdb( argv[3] );
	if( !is_pdb.good() )
	    throw( Error( ERROR_LOCATION, (std::string)"couldn\'t open file \'" + argv[3] + "\'" ) );
	pdb = new ParticleDataBase3D( is_pdb, geom );
	is_pdb.close();
    }

    std::ofstream of("particles.txt");
    for( size_t i = 0; i < pdb->size(); i++ ) {
	Particle3D &p = pdb->particle(i);
	double x = p[1];
	double y = p[3];
	double z = p[5];
	double vx = p[2];
	double vy = p[4];
	double vz = p[6];

	double dx = p[1]-277.5e-3;
	double yp = vy/vx;
	double zp = vz/vx;
	if( x > geom.max(0)-geom.h() )
	    of << p.q()/CHARGE_E << " "
	       << p.IQ() << " "
	       << y-yp*dx << " "
	       << z-zp*dx << " "
	       << yp << " "
	       << zp << "\n";
    }



    
    // -----------------------------------
    /*
    GTKPlotter plotter( &argc, &argv );
    plotter.set_geometry( &geom );
    if( epot ) {
	plotter.set_epot( epot );
	EpotEfield *efield = new EpotEfield( *epot );
	plotter.set_efield( efield );
    }
    if( pdb ) {
	MeshScalarField *tdens = new MeshScalarField( geom );
	pdb->build_trajectory_density_field( *tdens );
	plotter.set_trajdens( tdens );
	plotter.set_particledatabase( pdb );
    }
    plotter.new_geometry_plot_window();
    plotter.run();
    */
}


int main( int argc, char **argv )
{
    if( argc <= 1 ) {
	std::cout << "Usage: gtkplot <geom> <epot> <pdb>\n";
	exit( 1 );
    }

    try {
	//ibsimu.set_message_output( outputFolder + "ibsimu.txt" );
	ibsimu.set_message_threshold( MSG_VERBOSE, 1 );
	ibsimu.set_thread_count( 4 );

	simu( argc, argv );
    } catch( Error e ) {
        e.print_error_message( ibsimu.message(0) );
        exit( 1 );
    }
    
    return( 0 );
}