82 lines
2.4 KiB
Matlab
82 lines
2.4 KiB
Matlab
function plotPhaseSpaceWithAccelerationField(obj, MinimumVelocity, MaximumVelocity, NumberOfBins, IncidentAtomDirection, IncidentAtomPosition)
|
|
|
|
f_h = Helper.getFigureByTag('Phase Space Plot');
|
|
set(groot,'CurrentFigure',f_h);
|
|
a_h = get(f_h, 'CurrentAxes');
|
|
if ~isempty(get(a_h, 'Children'))
|
|
clf(f_h);
|
|
end
|
|
f_h.Name = 'Phase Space Plot';
|
|
f_h.Units = 'pixels';
|
|
|
|
set(0,'units','pixels');
|
|
screensize = get(0,'ScreenSize');
|
|
f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 750 600];
|
|
|
|
N = obj.NumberOfAtoms;
|
|
L = obj.OvenDistance * 2;
|
|
Theta = IncidentAtomDirection;
|
|
z = IncidentAtomPosition;
|
|
T = obj.SimulationTime;
|
|
tau = obj.TimeStep;
|
|
|
|
[X,Y] = meshgrid(-L/2:L/NumberOfBins:L/2,-MaximumVelocity:2*MaximumVelocity/NumberOfBins:MaximumVelocity);
|
|
|
|
a=zeros(NumberOfBins+1,NumberOfBins+1,3);
|
|
|
|
obj.setInitialConditions();
|
|
|
|
for i=1:length(X)
|
|
for j=1:length(Y)
|
|
a(i,j,:) = obj.calculateTotalAcceleration([X(1,i), 0, z], [Y(j,1)*cos(Theta),0,Y(j,1)*sin(Theta)]);
|
|
end
|
|
end
|
|
for i=1:length(X)
|
|
for j=1:length(Y)
|
|
if isnan(a(i,j,1)) || isnan(a(i,j,2)) || isnan(a(i,j,3))
|
|
a(i,j,1)=0;
|
|
a(i,j,2)=0;
|
|
a(i,j,3)=0;
|
|
end
|
|
end
|
|
end
|
|
|
|
pcolor(X',Y',a(:,:,1))
|
|
hold on
|
|
col=colorbar;
|
|
col.Label.String='Aceleration (m/s^2)';
|
|
shading flat
|
|
|
|
%-------------------------------------------------------------------------
|
|
|
|
Y = linspace(MinimumVelocity, MaximumVelocity,N);
|
|
DynamicalQuantities = zeros(length(Y),int64(T/tau),6);
|
|
for i=1:length(Y)
|
|
x =-L/2;
|
|
vx = Y(i)*cos(Theta);
|
|
vz = Y(i)*sin(Theta);
|
|
r = [x,0,z];
|
|
v = [vx,0,vz];
|
|
DynamicalQuantities(i,:,:) = obj.solver(r, v);
|
|
end
|
|
|
|
hold on
|
|
|
|
for i=1:length(Y)
|
|
plot(DynamicalQuantities(i,:,1),DynamicalQuantities(i,:,4),'w','linewidth',1.3)
|
|
end
|
|
|
|
hold off
|
|
|
|
hXLabel = xlabel('Position: Along the x-axis (m)');
|
|
hYLabel = ylabel('Velocity (m/s)');
|
|
hTitle = sgtitle(sprintf("Magnetic gradient = %.2f T/m", obj.MagneticGradient));
|
|
set([hXLabel, hYLabel] , ...
|
|
'FontSize' , 14 );
|
|
set( hTitle , ...
|
|
'FontSize' , 18 );
|
|
|
|
Helper.bringFiguresWithTagInForeground();
|
|
end
|
|
|