/*
   * Load internal state from an archive.
   */
   void StructureFactorGrid::loadParameters(Serializable::IArchive &ar)
   {
      nAtomType_ = simulation().nAtomType();

      // Load parameter file parameters
      loadInterval(ar);
      loadOutputFileName(ar);
      loadParameter<int>(ar, "nMode", nMode_);
      modes_.allocate(nMode_, nAtomType_);
      loadDMatrix<double>(ar, "modes", modes_, nMode_, nAtomType_);
      loadParameter<int>(ar, "hMax", hMax_);
      loadParameter<LatticeSystem>(ar, "lattice", lattice_);

      // Load and broadcast other distributed members
      MpiLoader<Serializable::IArchive> loader(*this, ar);
      loader.load(nWave_);
      waveIntVectors_.allocate(nWave_);
      loader.load(waveIntVectors_, nWave_);
      loader.load(nStar_);
      starIds_.allocate(nStar_);
      loader.load(starIds_, nStar_);
      starSizes_.allocate(nStar_);
      loader.load(starSizes_, nStar_);
      loader.load(nSample_);

      if (simulation().domain().isMaster()) {
         structureFactors_.allocate(nWave_, nMode_);
         ar >> structureFactors_;
      }

    void TestCardiacSimulationArchiveDynamic() throw(Exception)
    {
#ifdef CHASTE_CAN_CHECKPOINT_DLLS
        // run a monodomain simulation
        {
            CardiacSimulation simulation("heart/test/data/xml/save_monodomain_dynamic.xml");
        }
        std::string foldername = "SaveMonodomainDynamic";

        // compare the files, using the CompareFilesViaHdf5DataReader() method
        TS_ASSERT( CompareFilesViaHdf5DataReader("heart/test/data/cardiac_simulations", "save_monodomain_dynamic", false,
                   foldername, "SimulationResults", true));

        FileFinder file(foldername + "_checkpoints/0.2ms/" + foldername + "_0.2ms/archive.arch.0",
                        RelativeTo::ChasteTestOutput);
        TS_ASSERT(file.Exists());

        /* If you want to update the .h5 results for this test for any reason, you need to stop the following lines adding to them.
         * So uncomment the assert(), run the test, and then do:
         cp /tmp/chaste/testoutput/SaveMonodomainDynamic/SimulationResults.h5 heart/test/data/cardiac_simulations/save_monodomain_dynamic.h5
         */
        //assert(0);

        //resume the simulation
        {
            CardiacSimulation simulation("heart/test/data/xml/resume_monodomain_dynamic.xml");
        }

        // compare the files, using the CompareFilesViaHdf5DataReader() method
        TS_ASSERT( CompareFilesViaHdf5DataReader("heart/test/data/cardiac_simulations", "resume_monodomain_dynamic", false,
                   foldername, "SimulationResults", true));

#endif // CHASTE_CAN_CHECKPOINT_DLLS
    }

bool PhysicsFixedFactory::create(const char *name, OscObject *object1, OscObject *object2)
{
    OscFixed *cons=NULL;
    cons = new OscFixedODE(simulation()->odeWorld(),
                           simulation()->odeSpace(),
                           name, m_parent, object1, object2);

    if (cons)
        return simulation()->add_constraint(*cons);
}

void solve()
{
  fc = find_pos( 'F' );
  cc = find_pos( 'C' );

  simulation( fc, f_p, f_path );
  simulation( cc, c_p, c_path );

  fout << process() << std::endl;
}

 void TestBiWithBath1dSmall() throw(Exception)
 {
     { CardiacSimulation simulation("heart/test/data/xml/bidomain_with_bath1d_small.xml"); }
     { CardiacSimulation simulation2("heart/test/data/xml/bidomain_with_bath1d_resume.xml"); }
     {
         // The default resume file specifies a simulation duration of zero.
         // In reality the user should edit the file to specify something sensible...
         FileFinder resume_xml("SaveBiWithBath1D_checkpoints/0.1ms/ResumeParameters.xml", RelativeTo::ChasteTestOutput);
         TS_ASSERT_THROWS_THIS(CardiacSimulation simulation(resume_xml.GetAbsolutePath()),
                               "The simulation duration must be positive, not -0.1");
     }
 }

bool HapticsSphereFactory::create(const char *name, float x, float y, float z)
{
    OscSphereCHAI *obj = new OscSphereCHAI(simulation()->world(),
                                           name, m_parent);

    if (!(obj && simulation()->add_object(*obj)))
        return false;

    obj->m_position.set(x, y, z);

    return true;
}

bool PhysicsPistonFactory::create(const char *name, OscObject *object1,
                                  OscObject *object2, double x, double y,
                                  double z, double ax, double ay, double az)
{
    OscPiston *cons=NULL;
    cons = new OscPistonODE(simulation()->odeWorld(),
                            simulation()->odeSpace(),
                            name, m_parent, object1, object2,
                            x, y, z, ax, ay, az);

    if (cons)
        return simulation()->add_constraint(*cons);
}

bool PhysicsFreeFactory::create(const char *name,
                                OscObject *object1, OscObject *object2)
{
    OscFree *cons=NULL;
    cons = new OscFreeODE(simulation()->odeWorld(),
                          simulation()->odeSpace(),
                          name, m_parent, object1, object2);

    cons->m_response->traceOn();

    if (cons)
        return simulation()->add_constraint(*cons);
}

bool PhysicsBallJointFactory::create(const char *name, OscObject *object1,
                                     OscObject *object2, double x, double y,
                                     double z)
{
    OscBallJoint *cons=NULL;
    cons = new OscBallJointODE(simulation()->odeWorld(),
                               simulation()->odeSpace(),
                               name, m_parent, object1, object2,
                               x, y, z);

    if (cons)
        return simulation()->add_constraint(*cons);
}

bool PhysicsPrismFactory::create(const char *name, float x, float y, float z)
{
    OscPrismODE *obj = new OscPrismODE(simulation()->odeWorld(),
                                         simulation()->odeSpace(),
                                         name, m_parent);

    if (!(obj && simulation()->add_object(*obj)))
            return false;

    obj->m_position.set(x, y, z);

    return true;
}

bool PhysicsHingeFactory::create(const char *name, OscObject *object1, OscObject *object2,
                                 double x, double y, double z, double ax, double ay, double az)
{
    OscHinge *cons=NULL;
    cons = new OscHingeODE(simulation()->odeWorld(),
                           simulation()->odeSpace(),
                           name, m_parent, object1, object2,
                           x, y, z, ax, ay, az);

    cons->m_response->traceOn();

    if (cons)
        return simulation()->add_constraint(*cons);
}

bool InterfaceSphereFactory::create(const char *name, float x, float y, float z)
{
    OscSphere *obj = new OscSphereInterface(NULL, name, m_parent);

    if (!(obj && simulation()->add_object(*obj)))
            return false;

    obj->m_position.set(x, y, z);
    obj->traceOn();

    simulation()->send(0, "/world/sphere/create", "sfff", name, x, y, z);

    return true;
}

PollenOnsetDateFromFile::PollenOnsetDateFromFile(QString fileName)
    : PollenOnsetDate(), data(simulation()->inputFilePath(fileName))
{
    DataGrid data(simulation()->inputFilePath(fileName));
    for (int i = 0; i < data.rowNumber(); ++i) {
        QDate date = stringToValue<QDate>(data.row(i).at(1));
        int year = date.year();
        if (!yearsFirstLine.contains(year)) {
            yearsFirstLine[year] = i;
        }
    }
    QList<int> years = yearsFirstLine.keys();
    lastYear = years.last();
}

bool InterfaceFixedFactory::create(const char *name, OscObject *object1, OscObject *object2)
{
    if (!object1) return false;

    OscFixed *cons = new OscFixedInterface(name, m_parent, object1, object2);
    if (!(cons && simulation()->add_constraint(*cons)))
            return false;

    cons->traceOn();

    simulation()->send(0, "/world/fixed/create", "sss",
                       name, object1->c_name(), object2?object2->c_name():"world");

    return true;
}

bool PhysicsHinge2Factory::create(const char *name, OscObject *object1,
                                  OscObject *object2, double x,
                                  double y, double z, double a1x,
                                  double a1y, double a1z, double a2x,
                                  double a2y, double a2z)
{
    OscHinge2 *cons=NULL;
    cons = new OscHinge2ODE(simulation()->odeWorld(),
                            simulation()->odeSpace(),
                            name, m_parent, object1, object2,
                            x, y, z, a1x, a1y, a1z, a2x, a2y, a2z);

    if (cons)
        return simulation()->add_constraint(*cons);
}

bool HapticsPrismFactory::create(const char *name, float x, float y, float z)
{
    printf("HapticsPrismFactory (%s) is creating a prism object called '%s'\n",
           m_parent->c_name(), name);

    OscPrismCHAI *obj = new OscPrismCHAI(simulation()->world(),
                                         name, m_parent);

    if (!(obj && simulation()->add_object(*obj)))
        return false;

    obj->m_position.set(x, y, z);

    return true;
}

OscCursorCHAI::OscCursorCHAI(cWorld *world, const char *name, OscBase *parent)
    : OscSphere(NULL, name, parent)
{
    // create the cursor object
    m_pCursor = new cMeta3dofPointer(world);
    world->addChild(m_pCursor);

    // User data points to the OscObject, used for identification
    // during object contact.
    m_pCursor->setUserData(this, 1);

    // replace the potential proxy algorithm with our own
    cGenericPointForceAlgo *old_proxy, *new_proxy;
    old_proxy = m_pCursor->m_pointForceAlgos[1];
    new_proxy = new cODEPotentialProxy(
        dynamic_cast<cPotentialFieldForceAlgo*>(old_proxy));
    m_pCursor->m_pointForceAlgos[1] = new_proxy;
    delete old_proxy;

    if (m_pCursor->initialize()) {
        m_bInitialized = false;
        printf("[%s] Could not initialize.\n", simulation()->type_str());
    } else {
        m_bInitialized = true;
        m_pCursor->start();

        printf("[%s] Using %s device.\n",
               simulation()->type_str(), device_str());
    }

    // rotate the cursor to match visual rotation
    m_pCursor->rotate(cVector3d(0,0,1),-90.0*M_PI/180.0);

    // make it a cursor tuned for a dynamic environment
    ((cProxyPointForceAlgo*)m_pCursor->m_pointForceAlgos[0])
    ->enableDynamicProxy(true);

    // this is necessary for the above rotation to take effect
    m_pCursor->computeGlobalPositions();

    // set up mass as zero to begin (transparent proxy)
    m_mass.set(0);

    // no extra force to begin with
    m_nExtraForceSteps = 0;

    m_pSpecial = new CHAIObject(this, m_pCursor, world);
}

int main(int argc, char **argv){
	 //printf( "usage: %s %s ", argv[1],argv[2] );
	 double mu,koff;
	 //int i;
	 //for (i=0;i<argc;i++){
		 //printf("%s\n",argv[i]);
	 //}
	mu=atof(argv[1]);
	koff=atof(argv[2]);
	//sprintf("%f %f %s %d",kon,koff,argv[3],argc);
	/*Initialize MPI*/
	MPI_Init(&argc,&argv);
	
	/*prints rank of MPI instance for debugging purposes*/
	int my_rank;
	MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
	//printf("outputs_metropolis_%1.1f_%1.1f/run%d.txt",kon,koff,my_rank);
	printf("outputs_%1.1f_%1.1f/run%d.txt",mu,koff,my_rank);
	printf("thread %d\n",my_rank);
	
	/*get final value of coverage, kinda only for debugging really*/
	
	simulation(mu,koff);
	
	/*finalize MPI instances*/
	MPI_Finalize();
	
	//fclose(fp);	
	//printf("final:  %llu\n",states);
    //printf("test\n");
    
    return 0;
}

int main(int argc, char *argv[])
{
    ExecutableSupport::StandardStartup(&argc, &argv);

    int exit_code = ExecutableSupport::EXIT_OK;

    try
    {
        if (argc<2)
        {
            ExecutableSupport::PrintError("Usage: Chaste parameters_file", true);
            exit_code = ExecutableSupport::EXIT_BAD_ARGUMENTS;
        }
        else
        {
            std::string xml_file_name(argv[1]);
            // Creates & runs the simulation
            CardiacSimulation simulation(xml_file_name, true);
        }
    }
    catch (const Exception& e)
    {
        ExecutableSupport::PrintError(e.GetMessage());
        exit_code = ExecutableSupport::EXIT_ERROR;
    }

    ExecutableSupport::FinalizePetsc();
    return exit_code;
}

int main(int argc, char *argv[])
{
	bool running = false;
	
	// Two verbosity levels
	myPrintf[QUIET] = quietprint;
	myPrintf[VERBOSE] = printf;
	
	// Before continuing, read the arguments
	running = readArguments( argc, argv );

	// If the arguments are correct, start a simulation
	while ( running ) {
		myPrintf[ verbosemode ] ( "main : running is true\n" );

		// Init graphics
		SDL_Surface * screen = NULL;
		screen = init_sdl();

		// Simulates
		simulation ( screen );

		// End of graphics
		SDL_FreeSurface(screen);
		quit_sdl();

		// Quit running mode
		running = false;
	};

	return 0;
}