#include "epot_bicgstabsolver.hpp"
#include "geometry.hpp"
#include "func_solid.hpp"
#include "epot_efield.hpp"
#include "meshvectorfield.hpp"
#include "particledatabase.hpp"
#include "geomplotter.hpp"
#include "gtkplotter.hpp"
#include "ibsimu.hpp"
#include "error.hpp"


void simu( int argc, char **argv )
{
    std::string geom_fn = "geom.dat";
    std::ifstream is_geom( geom_fn.c_str() );
    if( !is_geom.good() )
        throw( Error( ERROR_LOCATION, (std::string)"couldn\'t open file \'" + geom_fn + "\'" ) );
    Geometry geom( is_geom );
    is_geom.close();
    geom.build_surface();

    EpotField epot( geom );
    MeshScalarField scharge( geom );
    MeshVectorField bfield;
    EpotEfield efield( epot );
    field_extrpl_e efldextrpl[6] = { FIELD_EXTRAPOLATE, FIELD_EXTRAPOLATE, 
				     FIELD_SYMMETRIC_POTENTIAL, FIELD_EXTRAPOLATE,
				     FIELD_EXTRAPOLATE, FIELD_EXTRAPOLATE };
    efield.set_extrapolation( efldextrpl );

    EpotBiCGSTABSolver solver( geom );

    ParticleDataBase3D pdb( geom );
    bool pmirror[6] = { false, false, true, false, false, false };
    pdb.set_mirror( pmirror );
    pdb.set_surface_collision( true );

    for( size_t i = 0; i < 5; i++ ) {
	solver.solve( epot, scharge );
	efield.recalculate();
	pdb.clear();
	pdb.add_cylindrical_beam_with_energy( 100000, 50.0, 1.0, 1.0, 
					      3.0e3, 0.0, 0.0, 
					      Vec3D(0,0,0),
					      Vec3D(1,0,0),
					      Vec3D(0,1,0),
					      0.012 );
	pdb.iterate_trajectories( scharge, efield, bfield );
    }

    GeomPlotter geomplotter( geom );
    geomplotter.set_size( 750, 750 );
    geomplotter.set_epot( &epot );
    geomplotter.set_particle_database( &pdb );
    geomplotter.set_view( VIEW_ZX, 0 );
    geomplotter.plot_png( "plot_zx.png" );

    if( true ) {
        MeshScalarField tdens( geom );
        pdb.build_trajectory_density_field( tdens );
        GTKPlotter plotter( &argc, &argv );
        plotter.set_geometry( &geom );
        plotter.set_epot( &epot );
        plotter.set_bfield( &bfield );
        plotter.set_efield( &efield );
        plotter.set_scharge( &scharge );
        plotter.set_trajdens( &tdens );
        plotter.set_particledatabase( &pdb );
        plotter.new_geometry_plot_window();
        plotter.run();
    }

}


int main( int argc, char **argv )
{
    try {
	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 );
}