-
Notifications
You must be signed in to change notification settings - Fork 0
/
simulateSpin2.m
82 lines (63 loc) · 2.42 KB
/
simulateSpin2.m
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
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
% Spin simulation 2
%
% mikael.mieskolainen@cern.ch, 2017
clear;
close all;
addpath ../../#matlabcodes/
addpath ./src/
%% Final state quantum numbers
s1 = 0; % Daughter 1 spin
s2 = 0; % Daughter 2 spin
T = 1; % Helicity coupling matrix (2*s1+1) x (2*s1+1) (user sets this up)
%% Plot pure spin states
% Number of events
N = 1e4;
% Spin
for J = 1:2
statistics = 1; % Bose
fig1 = figure;
fig2 = figure;
% Legend texts
legs = cell(2*J+1, 1);
xedge = linspace(-1,1,75); % cos(theta)
yedge = linspace(0,2*pi,75); % phi
xcenter = (xedge(1:end-1) + xedge(2:end))/2; % bin centers
ycenter = (yedge(1:end-1) + yedge(2:end))/2; % bin centers
% Simulate each pure states
for kk = 1:2*J+1
% Loop over diagonal [-J <= M <= J]
rho_i = zeros(2*J+1); rho_i(kk,kk) = 1;
[X,weights] = spindecay(rho_i, s1, s2, T, N, 0, 0);
M = J:-statistics:-J;
H = hist2w(X,weights,xedge,yedge);
% Plot different histograms
% 2D histogram
figure(fig1);
subplot(floor(sqrt(2*J+1) + 0.5), ceil(sqrt(2*J+1)), kk);
imagesc(xcenter, ycenter, H' / sum(H(:))); colormap(hot);% colorbar;
xlabel('$\cos(\theta)$','interpreter','latex');
ylabel('$\phi$ (rad)','interpreter','latex');
title(sprintf('$|%d,%d \\rangle$', J, M(kk)),'interpreter','latex');
axis square;
% 1D histogram
figure(fig2);
stephistedge(xedge, sum(H',1) ); hold on;
xlabel('$\cos(\theta)$','interpreter','latex');
title(sprintf('Pure states $|J,J_z \\rangle$'),'interpreter','latex');
legs{kk} = sprintf('$|%d,%d \\rangle$', J, M(kk));
end
figure(fig1);
filename = sprintf('./figs/J%d_pure_2D.pdf', J);
cmd = sprintf('print -dpdf %s', filename);
eval(cmd); pause(2); close(fig1);
system(sprintf('pdfcrop --margins ''10 10 10 10'' %s %s', filename, filename)); pause(2);
figure(fig2);
l = legend(legs);
set(l,'interpreter','latex','location','southeast');
axis tight; axis square;
set(gca,'yscale','log');
filename = sprintf('./figs/J%d_pure_1D.pdf', J);
cmd = sprintf('print -dpdf %s', filename);
eval(cmd); pause(2); close(fig2);
system(sprintf('pdfcrop --margins ''10 10 10 10'' %s %s', filename, filename)); pause(2);
end