-
Notifications
You must be signed in to change notification settings - Fork 0
/
guideSimLZ.m
184 lines (176 loc) · 6.01 KB
/
guideSimLZ.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
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
%% Guide Escape
% Here I wish to study the escape dynamics from the 2D trap generated by
% Decelerator Guiding potentials.
%
% Get fine guiding grid
%{
guide = efftrap3D('ppgg','ppgg',-180,0,180);
guide = mean(guide,3);
x = -.975:.025:.975;
x = x*1e-3;
sp = 1e-2;
xq = -1:sp:1;
xq = xq*1e-3;
[x, y] = meshgrid(x,x);
[xq, yq] = meshgrid(xq,xq);
guideq = interp2(x,y,guide,xq,yq,'spline');
% Get derivatives for acceleration
dVdxm = conv2(guideq,[-.5 0 .5],'same')/(sp*1e-3);
dvdx = griddedInterpolant(xq',yq',dVdxm,'linear','none');
%}
%% Run the odesolver
m = 2.82328e-26; %binding energy accounted for.
times = 0:1e-6:5e-3;
exited = @(t,y) isnan(y(1))-.5;
crossplane = @(t,y) sin(y(3)*pi/5e-3);
dealio = @(q) q{:};
escO = odeset('Events',@(t,y) dealio({[exited(t,y), crossplane(t,y)],[1 0],[0,0]}),...
'RelTol',1e-6,'AbsTol',1e-6);
[t, p, te, ye, ie] = ode45(@(t,p) [p(4);p(5);p(6);dvdx(p(2),p(1))/m;dvdx(p(1),p(2))/m;0],...
[0 5e-3],[0 .0008 .001 -7 4 800],escO);
length(t)
figure; plot(p(:,1),p(:,2))
hold on; plot(ye(:,1),ye(:,2),'ro')
length(te)
any(abs(ye(2:2:end,2))<5e-6)
any(abs(ye(1:2:end,1))<5e-6)
%% Run odesolver for collection of related initial positions, velocities
% See how S333/S40 varies with E-field and initial distribution in presence
% of spin-flip loss with different effective sizes.
tic
N = 20000;
temps = [0.25 0.35];
efields = 1:1:6;
allef = zeros(7,1);
alltimes = zeros(length(efields),length(temps),N);
allhops = false(length(efields),length(temps),N);
allexits = allhops;
all40 = alltimes;
all333= alltimes;
ratios = zeros(length(efields),length(temps));
for ef=1:length(efields)
for tt = 1:length(temps)
ri = sqrt(rand(N,1))*0.8e-3;
thi = rand(N,1)*2*pi;
k = 1.3806e-23;
T = temps(tt);
y0 = [ri.*cos(thi), ri.*sin(thi) , randn(N,1)*5e-3 , randn(N,2)*sqrt(k*T/m) , (randn(N,1)*sqrt(k*T/m) + 800)];
%y0 = [rand(N,2)*2e-3-1e-3 , randn(N,2)*sqrt(k*T/m)];
times = zeros(N,1);
hops = false(N,1);
exits = false(N,1);
fprintf('\n\n')
for i=1:N
fprintf('\b\b\b\b\b\b\b%7d',i)
e = efields(ef)/6;
[t, pw, te, ye, ie] = ode45(@(t,p) [p(4);p(5);p(6);e*dvdx(p(2),p(1))/m;e*dvdx(p(1),p(2))/m;0],...
[0 2.05e-3],y0(i,:),escO);
times(i) = t(end);
if any(ie==1)
times(i) = min(te(ie==1));
exits(i) = true;
end
te = te(ie==2);
ye = ye(ie==2,:);
yey = ye(2:2:end,2);
yex = ye(1:2:end,1);
isouty = abs(yey)<(7-efields(ef))*1e-6;
isoutx = abs(yex)<(7-efields(ef))*1e-6;
if any(isouty)
ter = te(2:2:end);
tm = min(ter(isouty));
times(i) = min(tm,times(i));
hops(i) = true;
end
if any(isoutx)
ter = te(1:2:end);
tm = min(ter(isoutx));
times(i) = min(tm,times(i));
hops(i) = true;
end
%pp(i,:) = p(end,:);
end
alltimes(ef,tt,:) = times;
allhops(ef,tt,:) = hops;
allexits(ef,tt,:) = exits;
all40(ef,tt,:) = times > .00025;
all333(ef,tt,:) = times > 2.04e-3;
ratios(ef,tt) = sum(all333(ef,tt,:))/sum(all40(ef,tt,:));
%figure;hist(times,50)
end
fprintf(' \n')
end
toc
figure;plot(efields,ratios(:,1)-ratios(:,2),'*')
%% Results Readme
% So far I have checked the escape dynamics for molecules loaded into the
% transverse trap created by DC guiding with +/-6.25 kV applied in a ++--
% manner to the rods, with different initial condition assumptions:
% 0.5 K, filled 2x2mm square: 1243 @ 40 stages, 1206 @ 333 stages.
% 0.3 K, filled 2x2mm square: 1317 @ 40 stages, 1273 @ 333 stages.
% The ratio of ratios gives a 0.4% effect, (experiment shows 10%).
%
% 0.5 K, filled 1mm diameter circle: 1244 @ 40 stages, 1229 @ 333 stages.
% 0.3 K, filled 1mm diameter circle: 2017 @ 40 stages, 1982 @ 333 stages.
% ratio of ratios: 0.5% effect.
%
% 1.0 K, filled 2x2mm square: 964 @ 40 stages, 940 @ 333 stages.
% 1.3 K, filled 2x2mm square: 839 @ 40 stages, 816 @ 333 stages.
% ratio of ratios: 5% effect.
%
% 1.0 K, offset 1mm circle: 770 @ 40 stages, 751 @ 333 stages.
% 1.3 K, offset 1mm circle: 633 @ 40 stages, 614 @ 333 stages.
% ratio of ratios: 0.5% effect.
%% Plot what's going on
figure(100)
hist(times,[2.5e-5:5e-5:2.025e-3])
hold on
times10 = repmat(times,1,10);
times10=times10(times10>1.5e-4 & times10<2.5e-4);
%hist(times10,[1.75e-4:5e-5:2.75e-4])
times100 = repmat(times,1,1000);
times100=times100(times100>2.5e-4 & times100<2e-3);
hist(times100,[2.75e-4:5e-5:1.975e-3])
%% Figure out how temperature is evolving.
KE = sum(pp(:,[3 4]).^2,2)*m*.5;
PE = interp2(xq,yq,guideq,pp(:,1),pp(:,2));
NN = ~isnan(pp(:,4));
T = mean(KE(NN))/k
V = mean(PE(NN))/k
E = T+V
% Ti=.05K, 1mm circle, 28, 40, 68 mK (T,V,E)
% Ti=0.3K, 1mm circle, 37, 56, 93 mK (T,V,E)
% Ti=0.5K, 1mm circle, 39, 58, 97 mK (T,V,E)
% Ti=1.0K, 1mm circle, 39, 60, 98 mK (T,V,E)
% Ti=.05K, 2mm square, 55, 59, 114 mK (T,V,E)
% Ti=0.3K, 2mm square, 59, 71, 130 mK (T,V,E)
% Ti=0.5K, 2mm square, 61, 71, 132 mK (T,V,E)
% Ti=1.0K, 2mm square, 58, 74, 132 mK (T,V,E)
% Ti=.05K, 1mm circle offset, 41, 50, 90 mK (T,V,E)
% Ti=0.3K, 1mm circle offset, 50, 62, 112 mK (T,V,E)
% Ti=0.5K, 1mm circle offset, 51, 64, 115 mK (T,V,E)
% Ti=1.0K, 1mm circle offset, 53, 64, 117 mK (T,V,E)
%% Put the above results in a figure
% First of all, for thermal equilibrium, we should have kT of thermal, 2kT
% of potential. Not surprisingly, we don't find this, since there's no
% mechanism for thermal equilibrium, we're just exchanging potential and
% kinetic by orbiting in the trap.
tp05 = [68 114 90];
tp3 = [93 130 112];
tp5 = [97 132 115];
tp10 =[98 132 117];
all = [tp05;tp3;tp5;tp10];
all = all(:,[1 3 2]);
x = [1 2 3];
figure; hold on;
colors = get(gca,'ColorOrder');
bar(all')
ll = legend('50 mK','300 mK','500 mK','1000 mK');
set(ll,'Location','NorthWest');
ylabel('Total Energy (mK)')
xlabel('Initialization Strategy')
set(gca,'XTickLabelMode','manual')
set(gca,'XTickLabel',{'','1mm circle','','offset circle','','2mm square',''})
set(gca,'FontSize',12);
title('Temperature after Guiding, Different Initializations')
grid on