Page MenuHomec4science
Files F120053971

exercice8.cpp
No OneTemporary
Actions

File Metadata

Created
Tue, Jul 1, 15:29

exercice8.cpp
View Options

#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <cmath>
#include <complex> // Pour les nombres complexes
#include "ConfigFile.tpp"
#include "distrib.h"
#include "wavemap.h"
#include "vecNext.h"
using namespace std;
class Symulation
{
// symulation and output setup
const string output;
ofstream out_density;
ofstream out_state;
bool verbose;
const double tfin; // total time
double t; // current time
const double dt;
const double omega;
const double delta;
wavemap * system;
const size_t sampling;
enum outmode_t {
all = 0,
end = 1,
state = 2
};
const outmode_t outmode;
SVector<double, 2> detections; //DVector<double> detections;
size_t detection_index;
public:
Symulation(ConfigFile& config)
: output(config.get<string>("output")),
out_density(output + "_density.out"),
out_state(output + "_state.out"),
verbose(config.get<bool>("verbose")),
tfin(config.get<double>("tfin")),
t(0),
dt(config.get<double>("dt")),
omega(config.get<double>("omega")),
delta(config.get<double>("delta")),
sampling(config.get<int>("sampling")),
outmode(static_cast<outmode_t>(config.get<int>("outmode"))),
detections(0),
detection_index(0)
{
detections = config.get<SVector<double, 2>>("detections");
npdebug(detections)
npdebug(detections.size())
SVector<double, 2> bounds = config.get<SVector<double, 2>>("bounds");
size_t N = config.get<size_t>("N");
npdebug("Setting up system")
system = new wavemap(bounds, N, dt,
[&](double x) {
return V(x);
});
double x0 = config.get<double>("x0");
double s_norm = config.get<double>("sigma");
double n = config.get<double>("n");
dstb::gaussian init(x0, s_norm, n, bounds);
system->set(init);
cout.precision(15);
npdebug("Normalizing gaussian")
// normalize init conditions
double prob = sqrt(system->probability(bounds[0], bounds[1]));
for (size_t k = 0; k <= N; ++k)
(*system)[k] /= prob;
// Setup outputs
npdebug("Output stream setup")
out_density.precision(15);
out_state.precision(15);
ofstream potential_out(output + "_pot.out");
potential_out.precision(15);
for (size_t n = 0; n <= N; ++n)
potential_out << system->supp(n) << " " << V(system->supp(n)) << endl;
potential_out.close();
npdebug("Setup done")
}
~Symulation()
{
delete system;
out_density.close();
out_state.close();
}
double V(double x) const
{
return .5*omega*omega*min((x-delta)*(x-delta),(x+delta)*(x+delta));
}
bool check_detection() const
{
return detection_index < detections.size() && t >= detections[detection_index];
}
void reset_detection()
{
// set negative values to zero
for (size_t k = 0; k < system->N; ++k)
{
if (system->supp(k) < 0)
(*system)[k] = 0;
}
// normalize
double A = system->probability(system->bounds[0], system->bounds[1]);
for (size_t k = 0; k < system->N; ++k)
(*system)[k] /= sqrt(A);
++detection_index;
}
void write_state(ostream& out, double t)
{
out << t << " "
<< system->probability(system->bounds[0], 0) << " " // probabilite que la particule soit en x < 0
<< system->probability(0, system->bounds[1]) << " " // probabilite que la particule soit en x > 0
<< system->energy() << " " // Energie
<< system->position(1) << " " // Position moyenne
<< system->position(2) << " " // Position^2 moyenne
<< system->incert_pos() << " " // Incertitude position
<< system->momentum(1) << " " // Quantite de mouvement moyenne
<< system->momentum(2) << " "
<< system->incert_mom() << endl; // Incertitude momentum
}
void write_density(ostream& out, double t)
{
for (size_t n = 0; n <= system->N; ++n)
out << pow(abs((*system)[n]), 2) << " ";
out << endl;
}
void write_data(double t, bool dens = true)
{
if (dens) {
write_density(out_density, t);
if (verbose)
write_density(cerr, t);
}
write_state(out_state, t);
if (verbose)
write_state(cerr, t);
}
void start()
{
double stepsd = tfin / dt;
//size_t steps = ( abs(static_cast<double>(steps - static_cast<size_t>(stepsd + 1))) < (1e-2 * dt) ) ? stepsd + 1 : stepsd;
size_t steps = std::round(stepsd);
npdebug("dt = ", dt)
npdebug("stepsd = ", stepsd)
npdebug("steps = ", steps)
npdebug("detections size = ", detections.size())
for(size_t step = 1; step <= steps; ++step)
{
if (check_detection())
reset_detection();
switch (outmode)
{
case outmode_t::all:
if ((step-1) % sampling == 0)
write_data(t);
break;
case outmode_t::state:
if ((step-1) % sampling == 0)
write_data(t, false);
break;
}
//if (step % sampling == 0)
// npdebug("step = ", step);
system->evolve();
t += dt;
}
if (outmode == outmode_t::end)
write_data(t);
}
};
int main(int argc,char **argv)
{
string inputPath("configuration.in"); // Fichier d'input par defaut
if(argc>1) // Fichier d'input specifie par l'utilisateur ("./Exercice8 config_perso.in")
inputPath = argv[1];
ConfigFile configFile(inputPath); // Les parametres sont lus et stockes dans une "map" de strings.
for(int i(2); i<argc; ++i) // Input complementaires ("./Exercice8 config_perso.in input_scan=[valeur]")
configFile.process(argv[i]);
Symulation sym(configFile);
sym.start();
return 0;
}

Event Timeline

c4science · Help