Skip to content
Snippets Groups Projects
Commit 8bdf13e0 authored by Marc Pelizaeus's avatar Marc Pelizaeus
Browse files

Add ROOT macro for e/m

parent 27e2d7fc
No related branches found
No related tags found
No related merge requests found
simem.C 0 → 100644
void euler(double &x, double &y, double &vx, double &vy);
void expectedPosition(double t, double &x, double &y);
void predictor_corrector(double &x, double &y, double &vx, double &vy);
void runge_kutta(double &x, double &y, double &vx, double &vy);
double dt = 1E-12;
double em = 1.75882E11;
double B=1E-3;
double vinit = 1e6;
double expectedRadius = vinit/B/em;
void simem(TString algo, double time){
dt=time;
//TH1F* hr = new TH1F("hr","hr", 10, 0,10);
//TH1F* hv = new TH1F("hv","hv", 10, 0,10);
std::vector<double> times;
std::vector<double> radii;
std::vector<double> velocities;
TH2F* hxy = new TH2F("hxy", "hxy", 1000, -0.01,0.01, 1000, -0.02, 0.01);
double x(0), y(expectedRadius), vx(vinit), vy(0);
int i=0;
double vy_prev=-999.;
int count=0;
while(1){
//for(int i=0; i<1E4; i++){
if(algo=="ec") euler(x,y,vx,vy);
if(algo=="pc")predictor_corrector( x,y,vx,vy );
if(algo=="rk") runge_kutta( x,y,vx,vy );
if(vy_prev>0 && vy<0 ){
count++;
}
vy_prev = vy;
if(count==3) break;
hxy->Fill(x,y);
times.push_back(i*dt);
//double expected_x, expected_y;
//expectedPosition(i*dt, expected_x, expected_y);
double velocity = sqrt(vx*vx+vy*vy);
//hv->SetBinContent(binToFill, velocity/vinit);
velocities.push_back(velocity/vinit);
double radius = sqrt(x*x+y*y);
//hr->SetBinContent(binToFill, radius/expectedRadius);
radii.push_back(radius/expectedRadius);
i++;
}
cout << "Iterations " << i<< endl;
//hr->SetMarkerStyle(2);
//hv->SetMarkerStyle(2);
TGraph* gr = new TGraph(times.size(), times.data(), radii.data());
TCanvas* c1 = new TCanvas("c1","c1",900,450);
c1->Divide(2,2);
c1->cd(1);
gr->Draw("ALP");
c1->cd(2);
//hv->Draw("P");
c1->cd(3);
hxy->Draw();
}
void euler(double &x, double &y, double &vx, double &vy){
x=x+vx*dt;
y=y+vy*dt;
double vx_old = vx;
double vy_old = vy;
vx = vx_old + em * B * vy_old * dt;
vy = vy_old - em * B * vx_old * dt;
return;
}
void predictor_corrector(double &x, double &y, double &vx, double &vy){
double x_k1 = vx*dt;
double y_k1 = vy*dt;
double vx_k1 = em * B * vy * dt;
double vy_k1 = - em * B * vx * dt;
double x_pred = x + x_k1;
double y_pred = y + y_k1;
double vx_pred = vx + vx_k1;
double vy_pred = vy + vy_k1;
double x_k2 = vx_pred*dt;
double y_k2 = vy_pred*dt;
double vx_k2 = em * B * vy_pred * dt;
double vy_k2 = - em * B * vx_pred * dt;
x = x+0.5*(x_k1+x_k2);
y = y+0.5*(y_k1+y_k2);
vx = vx+0.5*(vx_k1+vx_k2);
vy = vy+0.5*(vy_k1+vy_k2);
return;
}
void runge_kutta(double &x, double &y, double &vx, double &vy){
double vx_k1 = vx;
double vy_k1 = vy;
double ax_k1 = em * B * vy;
double ay_k1 = - em * B * vx;
double vx_k2 = vx+0.5*ax_k1*dt ;
double vy_k2 = vy+0.5*ay_k1*dt;
double ax_k2 = em*B*vy_k2;
double ay_k2 = -em*B*vx_k2;
double vx_k3 = vx + 0.5*ax_k2*dt;
double vy_k3 = vy + 0.5*ay_k2*dt;
double ax_k3 = em*B*vy_k3;
double ay_k3 = -em*B*vx_k3;
double vx_k4 = vx + ax_k3*dt;
double vy_k4 = vy + ay_k3*dt;
double ax_k4 = em*B*vy_k4;
double ay_k4 = -em*B*vx_k4;
x = x+1./6. *(vx_k1+2.*vx_k2+2*vx_k3+vx_k4)*dt;
y = y+1./6. *(vy_k1+2.*vy_k2+2*vy_k3+vy_k4)*dt;
vx = vx+1./6. *(ax_k1+2.*ax_k2+2*ax_k3+ax_k4)*dt;
vy = vy+1./6. *(ay_k1+2.*ay_k2+2*ay_k3+ay_k4)*dt;
return;
}
void expectedPosition(double t, double &x, double &y){
double r = vinit/B/em;
double omega = vinit/r;
x = r*sin(omega*t);
y = r*cos(omega*t);
}
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment