karthik/Calculations
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

Cosmetic changes only, changed a few variable names, dynamical quantities for all atoms across the full simulation time are now returned instead of just the final values.

Browse Source
This commit is contained in:
Karthik 2021-07-11 14:43:36 +02:00
parent bc10c207f6
commit a41ae74cb7
1 changed files with 26 additions and 39 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

63
MOT Capture Process Simulation/@MOTSimulator/solver.m
View File

@ -1,54 +1,41 @@
function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPosition, InitialVelocity) function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelocity)
if this.Gravity if this.Gravity
g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration]; g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration];
else else
g = 0; g = 0;
end end
% Probability of Background Collisions ParticleDynamicalQuantities = zeros(int64(this.SimulationTime/this.TimeStep),6);
collision = rand(1); for i=1:int64(this.SimulationTime/this.TimeStep)
CollisionProbability = 1 - exp(-this.SimulationTime/this.CollisionTime);
if collision >= CollisionProbability || this.AtomicBeamCollision == false %flag for skipping the background collision ParticleDynamicalQuantities(i,1:3) = InitialPosition;
ParticleDynamicalQuantities(i,4:6) = InitialVelocity;
ParticleTrajectory = zeros(int64(this.SimulationTime/this.TimeStep),9); rt = InitialPosition;
vt = InitialVelocity;
for i=1:int64(this.SimulationTime/this.TimeStep) ga1 = this.calculateTotalAcceleration(rt,vt) + g;
gv1 = vt .* this.TimeStep;
rt = rt + 0.5 * gv1;
vt = vt + 0.5 * ga1 .* this.TimeStep;
ParticleTrajectory(i,1:3) = InitialPosition; ga2 = this.calculateTotalAcceleration(rt,vt) + g;
ParticleTrajectory(i,4:6) = InitialVelocity; gv2 = vt .* this.TimeStep;
rt = rt + 0.5 * gv2;
vt = vt + 0.5 *ga2 .* this.TimeStep;
rt = InitialPosition; ga3 = this.calculateTotalAcceleration(rt,vt) + g;
vt = InitialVelocity; gv3 = vt .* this.TimeStep;
rt = rt + 0.5 * gv3;
vt = vt + ga3 .* this.TimeStep;
ga1 = this.calculateTotalAcceleration(rt,vt) + g; ga4 = this.calculateTotalAcceleration(rt,vt) + g;
gv1 = vt .* this.TimeStep; gv4 = vt .* this.TimeStep;
rt = rt + 0.5 * gv1;
vt = vt + 0.5 * ga1 .* this.TimeStep;
ga2 = this.calculateTotalAcceleration(rt,vt) + g; InitialPosition = InitialPosition + (gv1+2*(gv2+gv3)+gv4)./6;
gv2 = vt .* this.TimeStep; InitialVelocity = InitialVelocity + this.TimeStep*(ga1+2*(ga2+ga3)+ga4)./6;
rt = rt + 0.5 * gv2; %Acceleration = (ga1+2*(ga2+ga3)+ ga4)./6;
vt = vt + 0.5 *ga2 .* this.TimeStep;
ga3 = this.calculateTotalAcceleration(rt,vt) + g;
gv3 = vt .* this.TimeStep;
rt = rt + 0.5 * gv3;
vt = vt + ga3 .* this.TimeStep;
ga4 = this.calculateTotalAcceleration(rt,vt) + g;
gv4 = vt .* this.TimeStep;
InitialPosition = InitialPosition + (gv1+2*(gv2+gv3)+gv4)./6;
InitialVelocity = InitialVelocity + this.TimeStep*(ga1+2*(ga2+ga3)+ga4)./6;
ParticleTrajectory(i,7:9) = (ga1+2*(ga2+ga3)+ ga4)./6;
end
FinalDynamicalQuantities = ParticleTrajectory(end,:);
else
ParticleTrajectory = zeros(int64(this.SimulationTime/this.TimeStep),9);
FinalDynamicalQuantities = -500*ones(1,9); % -500 is a flag for giving up this trajectory
end end
end end