1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
|
#include <iostream>
#include <cmath>
#include <fstream>
#include <vector>
using namespace std;
struct particle{
double x;
double y;
};
void readin_file(string& filename, vector <particle>& p, vector <particle>& v){
ifstream inFile;
inFile.open(filename);
if(inFile.fail()){ // Error test
cerr<< "Error opening file"<<endl;
exit(1);
}
double x,y,vx,vy;
while(inFile >> x >> y >> vx >> vy){
p.push_back({x,y});
v.push_back({vx,vy});
}
inFile.close();
}
void update_eta_force(vector <particle>& vel, double& eta_force, double& T){
cout<< vel.size()<<endl;
eta_force=-2*vel.size()*T;
cout << eta_force<<endl;
cout<<vel.size()<<endl;
for (int i=0;i<vel.size();i++){
eta_force+=(vel[i].x * vel[i].x) + (vel[i].y * vel[i].y);
}
}
int main(){
double dt=0.0005;
double R=pow(2,1./6);
double L=14;
double Ekin;
double Epot;
double T=1.5; //thermostat temperature
double temp; //instantaneous temperature
double eta=1;
double p_eta=1;
double p_eta_half;
double eta_force;
double Q_bar=100; //coupling constant
string initial_dat="initial.txt";
vector <particle> p,v;
readin_file(initial_dat,p,v);
vector <particle> F(p.size());
double Q=Q_bar*2*p.size();
update_forces(F,p,R,L,Epot);
update_eta_force(v,eta_force,T);
}
|