From a41ae74cb757c9e7b156088c03f1497625fdce78 Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Sun, 11 Jul 2021 14:43:36 +0200 Subject: [PATCH] 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. --- .../@MOTSimulator/solver.m | 65 ++++++++----------- 1 file changed, 26 insertions(+), 39 deletions(-) diff --git a/MOT Capture Process Simulation/@MOTSimulator/solver.m b/MOT Capture Process Simulation/@MOTSimulator/solver.m index 49e5451..222baa8 100644 --- a/MOT Capture Process Simulation/@MOTSimulator/solver.m +++ b/MOT Capture Process Simulation/@MOTSimulator/solver.m @@ -1,54 +1,41 @@ -function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPosition, InitialVelocity) +function ParticleDynamicalQuantities = solver(this, InitialPosition, InitialVelocity) + if this.Gravity g = [0,0,-Helper.PhysicsConstants.GravitationalAcceleration]; else g = 0; end + + ParticleDynamicalQuantities = zeros(int64(this.SimulationTime/this.TimeStep),6); + for i=1:int64(this.SimulationTime/this.TimeStep) - % Probability of Background Collisions - collision = rand(1); - CollisionProbability = 1 - exp(-this.SimulationTime/this.CollisionTime); + ParticleDynamicalQuantities(i,1:3) = InitialPosition; + ParticleDynamicalQuantities(i,4:6) = InitialVelocity; - if collision >= CollisionProbability || this.AtomicBeamCollision == false %flag for skipping the background collision + rt = InitialPosition; + vt = InitialVelocity; - ParticleTrajectory = zeros(int64(this.SimulationTime/this.TimeStep),9); + ga1 = this.calculateTotalAcceleration(rt,vt) + g; + gv1 = vt .* this.TimeStep; + rt = rt + 0.5 * gv1; + vt = vt + 0.5 * ga1 .* this.TimeStep; - for i=1:int64(this.SimulationTime/this.TimeStep) - - ParticleTrajectory(i,1:3) = InitialPosition; - ParticleTrajectory(i,4:6) = InitialVelocity; + ga2 = this.calculateTotalAcceleration(rt,vt) + g; + gv2 = vt .* this.TimeStep; + rt = rt + 0.5 * gv2; + vt = vt + 0.5 *ga2 .* this.TimeStep; - rt = InitialPosition; - vt = InitialVelocity; + ga3 = this.calculateTotalAcceleration(rt,vt) + g; + gv3 = vt .* this.TimeStep; + rt = rt + 0.5 * gv3; + vt = vt + ga3 .* this.TimeStep; - ga1 = this.calculateTotalAcceleration(rt,vt) + g; - gv1 = vt .* this.TimeStep; - rt = rt + 0.5 * gv1; - vt = vt + 0.5 * ga1 .* this.TimeStep; + ga4 = this.calculateTotalAcceleration(rt,vt) + g; + gv4 = vt .* this.TimeStep; - ga2 = this.calculateTotalAcceleration(rt,vt) + g; - gv2 = vt .* this.TimeStep; - rt = rt + 0.5 * gv2; - 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 + InitialPosition = InitialPosition + (gv1+2*(gv2+gv3)+gv4)./6; + InitialVelocity = InitialVelocity + this.TimeStep*(ga1+2*(ga2+ga3)+ga4)./6; + %Acceleration = (ga1+2*(ga2+ga3)+ ga4)./6; end end