From d88d4d683da9f6d1e750beb214fad6e94931f09d Mon Sep 17 00:00:00 2001 From: Karthik Chandrashekara Date: Tue, 29 Jun 2021 15:44:28 +0200 Subject: [PATCH] This is a class for simulating the capture process for Dy atoms in the 2-D and 3-D MOTs, including the dynamics from the push beam and slower beams. --- .../+Helper/ImageSelection.class | Bin 0 -> 1092 bytes .../+Helper/ImageSelection.java | 38 + .../+Helper/PhysicsConstants.m | 45 + .../+Helper/convertstruct2cell.m | 6 + .../+Helper/getFigureByTag.m | 191 ++++ MOT Capture Process Simulation/+Helper/ode5.m | 92 ++ .../+Helper/onenoteccdata.m | 55 ++ .../+Helper/parforNotifications.m | 148 ++++ .../+Helper/screencapture.m | 820 ++++++++++++++++++ .../plotPositionAndVelocitySampling.m | 56 ++ .../+Plotting/plotResultForTwoParameterScan.m | 65 ++ .../+Plotting/visualizeMagneticField.m | 70 ++ .../@MOTSimulator/MOTSimulator.m | 616 +++++++++++++ .../@MOTSimulator/accelerationDueToPushBeam.m | 26 + ...elerationDueToSpontaneousEmissionProcess.m | 15 + .../angularDistributionFunction.m | 37 + .../@MOTSimulator/calculateCaptureVelocity.m | 24 + .../calculateFreeMolecularRegimeFlux.m | 7 + .../@MOTSimulator/calculateLoadingRate.m | 67 ++ .../calculateLocalSaturationIntensity.m | 18 + .../calculateTotalAcceleration.m | 95 ++ .../@MOTSimulator/initialPositionSampling.m | 7 + .../@MOTSimulator/initialVelocitySampling.m | 47 + .../@MOTSimulator/magneticFieldForMOT.m | 14 + .../@MOTSimulator/reinitializeSimulator.m | 32 + .../@MOTSimulator/runSimulation.m | 66 ++ .../@MOTSimulator/setInitialConditions.m | 98 +++ .../@MOTSimulator/solver.m | 54 ++ .../velocityDistributionFunction.m | 5 + .../test_MOTSimulator.m | 100 +++ 30 files changed, 2914 insertions(+) create mode 100644 MOT Capture Process Simulation/+Helper/ImageSelection.class create mode 100644 MOT Capture Process Simulation/+Helper/ImageSelection.java create mode 100644 MOT Capture Process Simulation/+Helper/PhysicsConstants.m create mode 100644 MOT Capture Process Simulation/+Helper/convertstruct2cell.m create mode 100644 MOT Capture Process Simulation/+Helper/getFigureByTag.m create mode 100644 MOT Capture Process Simulation/+Helper/ode5.m create mode 100644 MOT Capture Process Simulation/+Helper/onenoteccdata.m create mode 100644 MOT Capture Process Simulation/+Helper/parforNotifications.m create mode 100644 MOT Capture Process Simulation/+Helper/screencapture.m create mode 100644 MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m create mode 100644 MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m create mode 100644 MOT Capture Process Simulation/+Plotting/visualizeMagneticField.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/angularDistributionFunction.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/initialPositionSampling.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/magneticFieldForMOT.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/runSimulation.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/solver.m create mode 100644 MOT Capture Process Simulation/@MOTSimulator/velocityDistributionFunction.m create mode 100644 MOT Capture Process Simulation/test_MOTSimulator.m diff --git a/MOT Capture Process Simulation/+Helper/ImageSelection.class b/MOT Capture Process Simulation/+Helper/ImageSelection.class new file mode 100644 index 0000000000000000000000000000000000000000..cfac41d4d552a6c83b184c77190049bc9ccc5c9a GIT binary patch literal 1092 zcma)5%Wl&^6g^`nabnywl(w{?v`zXT4`_e|3#eG|sDxC(A_b|ia*|2p%5{R{6uyOT zU;*V(2?XrjthEC>baV7;V=+==yzul^5H4_JLiqUj?<69oT_yd;PZbYY%wY3cKzIB%OV` zBLx=Y<}g#cH)yk2wjQZE8&jK(=LB~J3Z;LymY)eE?sr=xo!oXj`FOD3kp7O{a8;%w zgPmg`N{7I$5xUc4mZOQT?R9ET8hf%CP>}iXbyM~Nr|WUq*}r(B{a9ElmCxkEjMI;O zsSkR+t{=#j!pGa5D(|^Kdb8;s8>E+%1!lcF@SAeWQEOiaU93x&(kXaDJ&c7M7A#~j zX~DvTg$nWl*T=uvQ?GxbDOzo~yrQWJERV;P}Zs_nC9HJrT!tW|l&l&#*p}-H`1d-5?S03>np((?7CYjISJmVB^73MXbX5|Q? zQvC$&J#X}#F$3JKmam}YhwGwfEl+pH;EzIq5PO`xw0EqI z^2|}kJn@$>%SwVZrQ{;!7!~6JPoXL#jIpUOx5PNlO`^`?@oaNA`|WU6)W3=}=O{+S fyv}LrmrZ=1 && ischar(tag_name),'You must specify ``tag_name'' as a string.'); + + f_h = findobj('type','figure','tag',tag_name); + + if isempty(f_h) + f_h = figure('Tag',tag_name,varargin{:}); + + defaultNewFigProperties = {'Color','w','NumberTitle','off','Name',sprintf('Fig. %d',f_h.Number)}; + + varargin = [defaultNewFigProperties,varargin]; + else + f_h = f_h(1); + end + + if ~isempty(varargin) + set(f_h,varargin{:}); + end + + addCopyButton(f_h); + + if nargout > 0 + figure_handle = f_h; + else + set(groot,'CurrentFigure',f_h); + end + +end + +function addCopyButton(f_h) + + if(strcmp(f_h.ToolBar,'none')) + return + end + + tb = findall(f_h,'Type','uitoolbar'); + + pt = findall(tb, 'tag', 'Custom.CopyPlot' ); + if isempty(pt) + pt = uipushtool(tb); + else + pt = pt(1); + end + + cdata = zeros(16,16,3); + + % Evernote Logo +% cdata(:,:,1) =[255 NaN NaN NaN NaN 99 11 27 175 NaN NaN NaN NaN NaN NaN 255 +% NaN NaN NaN 251 93 14 0 0 0 66 70 106 210 NaN NaN NaN +% NaN NaN NaN 42 0 43 0 0 0 0 0 0 20 185 NaN NaN +% NaN 243 56 0 42 82 0 0 0 0 0 0 0 45 NaN NaN +% NaN 156 44 64 113 65 0 0 0 0 0 0 0 32 NaN NaN +% 136 9 26 28 11 0 0 0 0 0 0 0 0 10 188 NaN +% 132 0 0 0 0 0 0 0 0 0 136 175 16 0 133 NaN +% NaN 28 0 0 0 0 0 0 0 0 152 238 50 0 124 NaN +% NaN 58 0 0 0 0 0 0 0 0 0 9 0 0 71 NaN +% NaN 175 0 0 0 0 0 61 15 0 0 0 0 0 100 NaN +% NaN NaN 143 12 0 0 0 210 195 87 17 0 0 0 126 NaN +% NaN NaN NaN 183 118 50 150 NaN NaN 110 219 78 0 0 160 NaN +% NaN NaN NaN NaN NaN NaN NaN 191 0 35 NaN 150 0 23 NaN NaN +% NaN NaN NaN NaN NaN NaN NaN 124 0 172 NaN 81 0 93 NaN NaN +% 255 NaN NaN NaN NaN NaN NaN 183 0 0 0 0 51 228 NaN 245 +% 253 254 NaN NaN NaN NaN NaN NaN 156 63 45 100 NaN NaN 255 255]/255.; +% +% +% cdata(:,:,2) = [255 255 255 255 255 216 166 171 225 229 218 229 247 255 255 255 +% 255 255 255 255 201 166 159 157 167 188 189 200 243 255 255 255 +% 237 238 255 181 159 183 164 170 163 158 160 157 169 233 248 250 +% 224 235 188 140 182 195 161 168 168 168 168 169 147 186 244 240 +% 255 226 175 185 207 189 161 168 168 168 168 168 159 179 249 249 +% 227 172 172 179 172 163 169 168 168 170 163 155 160 173 231 237 +% 215 161 163 165 166 168 168 168 168 162 215 228 172 163 209 219 +% 248 178 159 168 168 168 168 168 168 159 220 249 185 158 208 222 +% 249 192 151 169 168 168 169 160 163 172 163 159 166 167 194 204 +% 246 229 155 157 168 169 159 188 174 154 162 167 166 166 202 214 +% 212 231 218 168 157 153 165 255 242 190 171 159 167 166 207 220 +% 218 203 251 243 206 181 230 210 208 207 242 196 154 168 223 232 +% 255 224 232 250 237 214 244 194 152 178 255 223 145 175 250 252 +% 255 255 244 239 222 213 240 214 149 228 254 199 136 203 244 232 +% 255 255 255 246 231 246 246 232 165 159 167 147 184 253 254 242 +% 253 254 255 255 254 255 255 255 231 183 178 199 249 255 255 255]/255.; +% +% +% cdata(:,:,3) = [255 255 255 255 255 117 38 50 187 211 170 190 234 255 255 255 +% 255 254 255 255 120 51 27 20 39 97 98 122 220 255 255 255 +% 238 252 246 73 22 71 37 49 35 20 24 18 49 196 231 231 +% 232 242 86 0 78 108 29 45 45 45 45 46 0 82 214 201 +% 255 175 63 85 139 98 27 45 45 45 45 45 23 72 233 231 +% 167 51 57 72 55 32 47 45 45 50 34 14 27 57 201 218 +% 154 30 33 38 39 45 45 45 45 31 157 188 53 34 153 180 +% 234 67 24 45 45 45 45 44 45 24 169 241 83 20 146 182 +% 241 99 4 48 45 45 47 28 35 53 32 26 39 44 104 127 +% 238 192 14 20 45 47 27 97 56 10 29 44 41 40 127 158 +% 214 253 169 37 20 16 34 218 207 105 55 23 42 40 147 182 +% 218 214 241 201 138 71 177 225 181 130 224 107 12 45 175 197 +% 255 233 202 218 212 132 230 196 27 61 255 172 0 64 240 242 +% 255 255 219 197 176 160 237 143 0 195 245 110 0 123 230 230 +% 255 255 255 227 197 241 244 202 36 24 39 0 81 228 242 245 +% 253 254 255 255 254 255 255 255 191 78 71 121 221 255 255 255]/255.; + + %OneNote logo + + cdata(:,:,1) =[255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 + 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 + 255 255 255 255 245 213 213 213 213 213 213 213 184 184 215 255 + 255 255 255 255 241 213 213 213 213 213 213 213 184 184 208 255 + 255 233 204 204 194 176 176 185 213 213 213 213 184 184 208 255 + 255 154 101 101 101 101 101 103 213 213 213 206 162 162 193 255 + 255 152 101 183 116 152 115 101 213 213 213 206 162 162 193 255 + 255 152 101 207 189 178 122 101 213 213 213 206 162 162 193 255 + 255 152 101 199 152 224 122 101 213 213 213 195 128 128 170 255 + 255 152 101 166 101 183 115 101 213 213 213 195 128 128 170 255 + 255 154 101 101 101 101 101 103 213 213 213 195 128 128 170 255 + 255 233 204 204 194 176 176 185 213 213 213 183 95 95 148 255 + 255 255 255 255 241 213 213 213 213 213 213 183 94 94 148 255 + 255 255 255 255 245 213 213 213 213 213 213 183 94 94 163 255 + 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 + 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255]/255.; + + + cdata(:,:,2) =[255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 + 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 + 255 255 255 255 219 112 110 110 110 110 110 134 84 84 158 255 + 255 255 255 255 207 110 110 110 110 110 110 134 84 84 141 255 + 255 222 178 178 146 81 81 88 110 110 110 134 84 84 141 255 + 255 102 23 23 23 23 23 24 110 110 110 125 58 58 123 255 + 255 100 23 147 46 100 44 23 110 110 110 125 58 58 123 255 + 255 100 23 183 156 139 55 23 110 110 110 125 58 58 123 255 + 255 100 23 170 99 208 55 23 110 110 110 119 38 38 109 255 + 255 100 23 121 23 146 44 23 110 110 110 119 38 38 109 255 + 255 102 23 23 23 23 23 24 110 110 110 119 38 38 109 255 + 255 222 178 178 146 81 81 88 110 110 110 118 37 37 109 255 + 255 255 255 255 207 110 110 110 110 110 110 118 37 37 110 255 + 255 255 255 255 219 112 110 110 110 110 110 118 37 37 131 255 + 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 + 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255]/255.; + + + cdata(:,:,3) =[255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 + 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 + 255 255 255 255 255 255 255 255 255 255 255 246 229 229 240 255 + 255 255 255 255 255 255 255 255 255 255 255 246 229 229 238 255 + 255 242 224 224 224 224 224 232 255 255 255 246 229 229 238 255 + 255 194 163 163 163 163 163 164 255 255 255 244 223 223 234 255 + 255 194 163 212 172 194 171 163 255 255 255 244 223 223 234 255 + 255 194 163 226 216 209 176 163 255 255 255 244 223 223 234 255 + 255 194 163 221 193 236 176 163 255 255 255 240 209 209 224 255 + 255 194 163 202 163 212 171 163 255 255 255 240 209 209 224 255 + 255 194 163 163 163 163 163 164 255 255 255 240 209 209 224 255 + 255 242 224 224 224 224 224 232 255 255 255 223 161 161 192 255 + 255 255 255 255 255 255 255 255 255 255 255 223 160 160 192 255 + 255 255 255 255 255 255 255 255 255 255 255 223 160 160 201 255 + 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 + 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255 255]/255.; + + + pt.Tag = 'Custom.CopyPlot'; + pt.CData = cdata; + pt.Separator = true; + pt.ClickedCallback = @copyToClipboard; + +end + +function copyToClipboard(~,~) + fig_h = get(get(gcbo,'Parent'),'Parent'); + if strcmp(fig_h.WindowStyle,'docked') + if ismac || ispc + matlab.graphics.internal.copyFigureHelper(fig_h); + else + %warning('Copy function to the clipboard only works if the figure is undocked.'); + Helper.screencapture(fig_h,[],'clipboard'); + end + else + pos = fig_h.Position; + Helper.screencapture(fig_h,[],'clipboard','position',[1,1,pos(3)-2,pos(4)]); + end +end + + + diff --git a/MOT Capture Process Simulation/+Helper/ode5.m b/MOT Capture Process Simulation/+Helper/ode5.m new file mode 100644 index 0000000..3eb003f --- /dev/null +++ b/MOT Capture Process Simulation/+Helper/ode5.m @@ -0,0 +1,92 @@ +function Y = ode5(odefun,tspan,y0,varargin) +%ODE5 Solve differential equations with a non-adaptive method of order 5. +% Y = ODE5(ODEFUN,TSPAN,Y0) with TSPAN = [T1, T2, T3, ... TN] integrates +% the system of differential equations y' = f(t,y) by stepping from T0 to +% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector. +% The vector Y0 is the initial conditions at T0. Each row in the solution +% array Y corresponds to a time specified in TSPAN. +% +% Y = ODE5(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters +% P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...). +% +% This is a non-adaptive solver. The step sequence is determined by TSPAN +% but the derivative function ODEFUN is evaluated multiple times per step. +% The solver implements the Dormand-Prince method of order 5 in a general +% framework of explicit Runge-Kutta methods. +% +% Example +% tspan = 0:0.1:20; +% y = ode5(@vdp1,tspan,[2 0]); +% plot(tspan,y(:,1)); +% solves the system y' = vdp1(t,y) with a constant step size of 0.1, +% and plots the first component of the solution. + +if ~isnumeric(tspan) + error('TSPAN should be a vector of integration steps.'); +end + +if ~isnumeric(y0) + error('Y0 should be a vector of initial conditions.'); +end + +h = diff(tspan); +if any(sign(h(1))*h <= 0) + error('Entries of TSPAN are not in order.') +end + +try + f0 = feval(odefun,tspan(1),y0,varargin{:}); +catch + msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr]; + error(msg); +end + +y0 = y0(:); % Make a column vector. +if ~isequal(size(y0),size(f0)) + error('Inconsistent sizes of Y0 and f(t0,y0).'); +end + +neq = length(y0); +N = length(tspan); +Y = zeros(neq,N); + +% Method coefficients -- Butcher's tableau +% +% C | A +% --+--- +% | B + +C = [1/5; 3/10; 4/5; 8/9; 1]; + +A = [ 1/5, 0, 0, 0, 0 + 3/40, 9/40, 0, 0, 0 + 44/45 -56/15, 32/9, 0, 0 + 19372/6561, -25360/2187, 64448/6561, -212/729, 0 + 9017/3168, -355/33, 46732/5247, 49/176, -5103/18656]; + +B = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]; + +% More convenient storage +A = A.'; +B = B(:); + +nstages = length(B); +F = zeros(neq,nstages); + +Y(:,1) = y0; +for i = 2:N + ti = tspan(i-1); + hi = h(i-1); + yi = Y(:,i-1); + + % General explicit Runge-Kutta framework + F(:,1) = feval(odefun,ti,yi,varargin{:}); + for stage = 2:nstages + tstage = ti + C(stage-1)*hi; + ystage = yi + F(:,1:stage-1)*(hi*A(1:stage-1,stage-1)); + F(:,stage) = feval(odefun,tstage,ystage,varargin{:}); + end + Y(:,i) = yi + F*(hi*B); + +end +Y = Y.'; diff --git a/MOT Capture Process Simulation/+Helper/onenoteccdata.m b/MOT Capture Process Simulation/+Helper/onenoteccdata.m new file mode 100644 index 0000000..5f3c177 --- /dev/null +++ b/MOT Capture Process Simulation/+Helper/onenoteccdata.m @@ -0,0 +1,55 @@ +cmap = zeros(16,16,3); + +cmap(:,:,1) = [0.0000 0.0118 0.4510 0.0039 0.2078 0.1569 0.4078 0.4431 0.4510 0.1922 0.4235 0.4196 0.2235 0.4235 0.4039 0.4392 + 0.4471 0.1647 0.4157 0.0000 0.0235 0.4353 0.0314 0.4314 0.0196 0.2392 0.0667 0.0392 0.4431 0.3804 0.2941 0.4275 + 0.3686 0.3608 0.2000 0.2824 0.3059 0.0549 0.1804 0.1882 0.4392 0.4314 0.3255 0.0078 0.0902 0.1961 0.4353 0.1412 + 0.2314 0.3647 0.0353 0.3804 0.1647 0.2431 0.1686 0.2745 0.2980 0.4235 0.3922 0.4157 0.2784 0.3333 0.2510 0.0588 + 0.1020 0.0745 0.2549 0.0471 0.1216 0.4000 0.3961 0.2627 0.1098 0.1725 0.3098 0.4314 0.3529 0.3412 0.0784 0.0824 + 0.4471 0.1490 0.1804 0.3529 0.2196 0.3137 0.3255 0.0941 0.0078 0.3294 0.3765 0.2706 0.0510 0.0157 0.4275 0.1176 + 0.1294 0.1333 0.1725 0.3451 0.2118 0.3843 0.1255 0.1569 0.2118 0.1608 0.0353 0.2039 0.1608 0.4510 1.0000 0.8000 + 0.9882 0.6510 0.9961 0.4549 0.4549 0.6824 0.7882 0.5686 0.5373 0.5490 0.7765 0.7137 0.8510 0.7176 0.5020 0.4902 + 0.8941 0.9020 0.4745 0.8980 0.9098 0.4824 0.6471 0.6353 0.9922 0.9647 0.6353 0.4588 0.9647 0.9020 0.4980 0.8118 + 0.5059 0.4941 0.9686 0.4863 0.5451 0.9725 0.8980 0.5451 0.5333 0.6824 0.4588 0.8196 0.8314 0.8980 0.8941 0.9961 + 0.5255 0.8392 0.9804 0.5216 0.8588 0.8078 0.5176 0.7647 0.5608 0.9725 0.9059 0.4627 0.9882 0.8275 0.7725 0.8745 + 0.8235 0.8431 0.7373 1.0000 0.5137 0.4706 0.4784 0.7412 0.8863 0.9373 0.5529 0.5804 0.4510 0.9255 0.8235 0.8667 + 0.7569 0.8824 0.5294 0.5176 0.5373 0.9569 0.5294 0.4824 0.5098 0.5137 0.5569 0.8471 0.5098 0.9490 0.8706 0.9412 + 0.4902 0.6000 0.6980 0.7882 0.5490 0.7216 0.6431 0.4824 0.5569 0.4667 0.6627 0.9922 0.7804 0.8039 0.6275 0.7333 + 0.5725 0.5647 0.8549 0.7529 0.6235 0.8784 0.5922 0.7294 0.6118 0.7922 0.7843 0.6667 0.9294 0.6902 0.6784 0.9176 + 0.6706 0.7490 0.7961 0.5882 0.8627 0.4627 0.6196 0.7059 0.6078 0.9765 0.6549 0.6863 0.5373 0.7098 0.7176 0.7765]; + +cmap(:,:,2) = [0.0000 0.0078 0.2157 0.0000 0.0980 0.0745 0.1922 0.2157 0.2157 0.0902 0.2000 0.1961 0.1059 0.2039 0.1882 0.2078 + 0.2078 0.0784 0.2000 0.0000 0.0118 0.2118 0.0157 0.2039 0.0078 0.1137 0.0314 0.0196 0.2118 0.1804 0.1373 0.2078 + 0.1765 0.1725 0.0941 0.1333 0.1451 0.0275 0.0863 0.0902 0.2078 0.2078 0.1529 0.0039 0.0431 0.0941 0.2039 0.0667 + 0.1098 0.1725 0.0157 0.1804 0.0784 0.1137 0.0824 0.1333 0.1412 0.2000 0.1882 0.2000 0.1333 0.1569 0.1176 0.0275 + 0.0471 0.0353 0.1216 0.0196 0.0588 0.1922 0.1882 0.1255 0.0510 0.0824 0.1451 0.2039 0.1686 0.1647 0.0392 0.0392 + 0.2157 0.0706 0.0863 0.1686 0.1020 0.1490 0.1529 0.0431 0.0039 0.1569 0.1804 0.1255 0.0235 0.0078 0.2000 0.0549 + 0.0627 0.0627 0.0824 0.1647 0.1020 0.1843 0.0588 0.0745 0.1020 0.0784 0.0157 0.0980 0.0784 0.2157 1.0000 0.7137 + 0.9843 0.4980 0.9961 0.2235 0.2196 0.5412 0.6980 0.3843 0.3373 0.3569 0.6824 0.5922 0.7843 0.6000 0.2902 0.2706 + 0.8510 0.8588 0.2471 0.8549 0.8667 0.2627 0.4980 0.4784 0.9843 0.9490 0.4745 0.2235 0.9451 0.8627 0.2824 0.7333 + 0.2941 0.2784 0.9529 0.2667 0.3490 0.9569 0.8510 0.3490 0.3333 0.5451 0.2275 0.7412 0.7608 0.8549 0.8471 0.9922 + 0.3255 0.7686 0.9725 0.3176 0.8000 0.7255 0.3098 0.6627 0.3725 0.9647 0.8627 0.2314 0.9804 0.7529 0.6745 0.8235 + 0.7451 0.7765 0.6235 0.9961 0.3020 0.2431 0.2510 0.6314 0.8392 0.9098 0.3608 0.4000 0.2196 0.8902 0.7490 0.8078 + 0.6549 0.8353 0.3294 0.3137 0.3412 0.9373 0.3255 0.2588 0.2980 0.3059 0.3686 0.7843 0.3020 0.9255 0.8157 0.9176 + 0.2745 0.4275 0.5686 0.6980 0.3569 0.6039 0.4863 0.2627 0.3647 0.2392 0.5137 0.9922 0.6863 0.7216 0.4706 0.6196 + 0.3882 0.3765 0.7882 0.6471 0.4588 0.8275 0.4157 0.6118 0.4431 0.7059 0.6902 0.5255 0.8980 0.5569 0.5412 0.8824 + 0.5333 0.6392 0.7098 0.4078 0.8039 0.2314 0.4549 0.5804 0.4392 0.9647 0.5059 0.5529 0.3373 0.5882 0.5961 0.6784]; + +cmap(:,:,3) = [0.0000 0.0157 0.4980 0.0039 0.2314 0.1725 0.4627 0.5020 0.5020 0.2196 0.4745 0.4706 0.2510 0.4784 0.4510 0.4980 + 0.4941 0.1882 0.4667 0.0000 0.0275 0.4941 0.0353 0.4902 0.0196 0.2667 0.0745 0.0471 0.4902 0.4314 0.3294 0.4784 + 0.4196 0.4000 0.2235 0.3216 0.3412 0.0627 0.2039 0.2118 0.4863 0.4863 0.3608 0.0078 0.1020 0.2196 0.4824 0.1569 + 0.2588 0.4118 0.0392 0.4235 0.1843 0.2745 0.1882 0.3059 0.3373 0.4784 0.4392 0.4627 0.3137 0.3765 0.2824 0.0667 + 0.1137 0.0824 0.2863 0.0510 0.1373 0.4510 0.4471 0.2941 0.1216 0.1961 0.3490 0.4824 0.3961 0.3804 0.0902 0.0941 + 0.4980 0.1647 0.2000 0.4000 0.2431 0.3529 0.3647 0.1059 0.0118 0.3686 0.4196 0.3020 0.0549 0.0196 0.4824 0.1294 + 0.1451 0.1529 0.1922 0.3882 0.2392 0.4353 0.1412 0.1765 0.2353 0.1804 0.0353 0.2275 0.1843 0.5059 1.0000 0.8196 + 0.9882 0.6863 0.9961 0.5098 0.5098 0.7137 0.8118 0.6118 0.5843 0.5922 0.8000 0.7412 0.8627 0.7451 0.5529 0.5412 + 0.9059 0.9137 0.5255 0.9098 0.9176 0.5333 0.6824 0.6706 0.9922 0.9686 0.6706 0.5098 0.9647 0.9137 0.5490 0.8314 + 0.5569 0.5451 0.9725 0.5373 0.5922 0.9725 0.9059 0.5882 0.5804 0.7137 0.5137 0.8353 0.8510 0.9059 0.9020 0.9961 + 0.5725 0.8549 0.9843 0.5725 0.8745 0.8275 0.5647 0.7882 0.6039 0.9765 0.9137 0.5176 0.9882 0.8431 0.7961 0.8863 + 0.8392 0.8588 0.7647 1.0000 0.5608 0.5216 0.5294 0.7686 0.8980 0.9412 0.6000 0.6235 0.5059 0.9333 0.8431 0.8784 + 0.7804 0.8941 0.5765 0.5686 0.5843 0.9608 0.5765 0.5333 0.5569 0.5647 0.6039 0.8627 0.5608 0.9569 0.8863 0.9490 + 0.5412 0.6392 0.7294 0.8078 0.5961 0.7490 0.6784 0.5373 0.6000 0.5216 0.6941 0.9922 0.8039 0.8235 0.6667 0.7608 + 0.6157 0.6078 0.8667 0.7765 0.6588 0.8902 0.6314 0.7569 0.6510 0.8157 0.8039 0.7020 0.9373 0.7216 0.7098 0.9255 + 0.7059 0.7725 0.8196 0.6314 0.8784 0.5137 0.6549 0.7373 0.6471 0.9804 0.6902 0.7176 0.5804 0.7412 0.7451 0.8000]; + +%% +[cdata, cmap] = imread('onenote.png'); \ No newline at end of file diff --git a/MOT Capture Process Simulation/+Helper/parforNotifications.m b/MOT Capture Process Simulation/+Helper/parforNotifications.m new file mode 100644 index 0000000..4ad3af4 --- /dev/null +++ b/MOT Capture Process Simulation/+Helper/parforNotifications.m @@ -0,0 +1,148 @@ +% Copyright (c) 2019 Andrea Alberti +% +% All rights reserved. +classdef parforNotifications < handle + properties + N; % number of iterations + text = 'Please wait ...'; % text to show + width = 50; + showWarning = true; + end + properties (GetAccess = public, SetAccess = private) + n; + end + properties (Access = private) + inProgress = false; + percent; + DataQueue; + usePercent; + Nstr; + NstrL; + lastComment; + end + methods + function this = parforNotifications() + this.DataQueue = parallel.pool.DataQueue; + afterEach(this.DataQueue, @this.updateStatus); + end + % Start progress bar + function PB_start(this,N,varargin) + assert(isscalar(N) && isnumeric(N) && N == floor(N) && N>0, 'Error: ''N'' must be a scalar positive integer.'); + + this.N = N; + + p = inputParser; + addParameter(p,'message','Please wait: '); + addParameter(p,'usePercentage',true); + + parse(p,varargin{:}); + + this.text = p.Results.message; + assert(ischar(this.text), 'Error: ''Message'' must be a string.'); + + this.usePercent = p.Results.usePercentage; + assert(isscalar(this.usePercent) && islogical(this.usePercent), 'Error: ''usePercentage'' must be a logical scalar.'); + + this.percent = 0; + this.n = 0; + this.lastComment = ''; + if this.usePercent + fprintf('%s [%s]: %3d%%\n',this.text, char(32*ones(1,this.width)),0); + else + this.Nstr = sprintf('%d',this.N); + this.NstrL = numel(this.Nstr); + fprintf('%s [%s]: %s/%s\n',this.text, char(32*ones(1,this.width)),[char(32*ones(1,this.NstrL-1)),'0'],this.Nstr); + end + + this.inProgress = true; + end + % Iterate progress bar + function PB_iterate(this,str) + if nargin == 1 + send(this.DataQueue,''); + else + send(this.DataQueue,str); + end + end + function warning(this,warn_id,msg) + if this.showWarning + msg = struct('Action','Warning','Id',warn_id,'Message',msg); + send(this.DataQueue,msg); + end + end + function PB_reprint(this) + p = round(100*this.n/this.N); + + this.percent = p; + + cursor_pos=1+round((this.width-1)*p/100); + + if p < 100 + sep_char = '|'; + else + sep_char = '.'; + end + + if this.usePercent + fprintf('%s [%s%s%s]: %3d%%\n', this.text, char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),p); + else + nstr=sprintf('%d',this.n); + fprintf('%s [%s%s%s]: %s/%s\n', this.text, char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),[char(32*ones(1,this.NstrL-numel(nstr))),nstr],this.Nstr); + end + end + function updateStatus(this,data) + + if ischar(data) + + this.n = this.n + 1; + + p = round(100*this.n/this.N); + + if p >= this.percent+1 || this.n == this.N + this.percent = p; + + cursor_pos=1+round((this.width-1)*p/100); + + if p < 100 + sep_char = '|'; + else + sep_char = '.'; + end + + if ~isempty(data) + comment = [' (',data,')']; + else + comment = ''; + end + + if this.usePercent + fprintf('%s%s%s%s]: %3d%%%s\n',char(8*ones(1,58+numel(this.lastComment))), char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),p,comment); + else + nstr=sprintf('%d',this.n); + fprintf('%s%s%s%s]: %s/%s%s\n',char(8*ones(1,55+2*numel(this.Nstr)+numel(this.lastComment))), char(46*ones(1,cursor_pos-1)), sep_char, char(32*ones(1,this.width-cursor_pos)),[char(32*ones(1,this.NstrL-numel(nstr))),nstr],this.Nstr,comment) + end + + this.lastComment = comment; + + + if p == 100 + this.inProgress = false; + end + end + + else + switch data.Action + case 'Warning' + warning(data.Id,[data.Message,newline]); + if this.inProgress + this.PB_reprint(); + end + end + + end + + end + end +end + + diff --git a/MOT Capture Process Simulation/+Helper/screencapture.m b/MOT Capture Process Simulation/+Helper/screencapture.m new file mode 100644 index 0000000..206312e --- /dev/null +++ b/MOT Capture Process Simulation/+Helper/screencapture.m @@ -0,0 +1,820 @@ +function imageData = screencapture(varargin) +% screencapture - get a screen-capture of a figure frame, component handle, or screen area rectangle +% +% ScreenCapture gets a screen-capture of any Matlab GUI handle (including desktop, +% figure, axes, image or uicontrol), or a specified area rectangle located relative to +% the specified handle. Screen area capture is possible by specifying the root (desktop) +% handle (=0). The output can be either to an image file or to a Matlab matrix (useful +% for displaying via imshow() or for further processing) or to the system clipboard. +% This utility also enables adding a toolbar button for easy interactive screen-capture. +% +% Syntax: +% imageData = screencapture(handle, position, target, 'PropName',PropValue, ...) +% +% Input Parameters: +% handle - optional handle to be used for screen-capture origin. +% If empty/unsupplied then current figure (gcf) will be used. +% position - optional position array in pixels: [x,y,width,height]. +% If empty/unsupplied then the handle's position vector will be used. +% If both handle and position are empty/unsupplied then the position +% will be retrieved via interactive mouse-selection. +% If handle is an image, then position is in data (not pixel) units, so the +% captured region remains the same after figure/axes resize (like imcrop) +% target - optional filename for storing the screen-capture, or the +% 'clipboard'/'printer' strings. +% If empty/unsupplied then no output to file will be done. +% The file format will be determined from the extension (JPG/PNG/...). +% Supported formats are those supported by the imwrite function. +% 'PropName',PropValue - +% optional list of property pairs (e.g., 'target','myImage.png','pos',[10,20,30,40],'handle',gca) +% PropNames may be abbreviated and are case-insensitive. +% PropNames may also be given in whichever order. +% Supported PropNames are: +% - 'handle' (default: gcf handle) +% - 'position' (default: gcf position array) +% - 'target' (default: '') +% - 'toolbar' (figure handle; default: gcf) +% this adds a screen-capture button to the figure's toolbar +% If this parameter is specified, then no screen-capture +% will take place and the returned imageData will be []. +% +% Output parameters: +% imageData - image data in a format acceptable by the imshow function +% If neither target nor imageData were specified, the user will be +% asked to interactively specify the output file. +% +% Examples: +% imageData = screencapture; % interactively select screen-capture rectangle +% imageData = screencapture(hListbox); % capture image of a uicontrol +% imageData = screencapture(0, [20,30,40,50]); % capture a small desktop region +% imageData = screencapture(gcf,[20,30,40,50]); % capture a small figure region +% imageData = screencapture(gca,[10,20,30,40]); % capture a small axes region +% imshow(imageData); % display the captured image in a matlab figure +% imwrite(imageData,'myImage.png'); % save the captured image to file +% img = imread('cameraman.tif'); +% hImg = imshow(img); +% screencapture(hImg,[60,35,140,80]); % capture a region of an image +% screencapture(gcf,[],'myFigure.jpg'); % capture the entire figure into file +% screencapture(gcf,[],'clipboard'); % capture the entire figure into clipboard +% screencapture(gcf,[],'printer'); % print the entire figure +% screencapture('handle',gcf,'target','myFigure.jpg'); % same as previous, save to file +% screencapture('handle',gcf,'target','clipboard'); % same as previous, copy to clipboard +% screencapture('handle',gcf,'target','printer'); % same as previous, send to printer +% screencapture('toolbar',gcf); % adds a screen-capture button to gcf's toolbar +% screencapture('toolbar',[],'target','sc.bmp'); % same with default output filename +% +% Technical description: +% http://UndocumentedMatlab.com/blog/screencapture-utility/ +% +% Bugs and suggestions: +% Please send to Yair Altman (altmany at gmail dot com) +% +% See also: +% imshow, imwrite, print +% +% Release history: +% 1.17 2016-05-16: Fix annoying warning about JavaFrame property becoming obsolete someday (yes, we know...) +% 1.16 2016-04-19: Fix for deployed application suggested by Dwight Bartholomew +% 1.10 2014-11-25: Added the 'print' target +% 1.9 2014-11-25: Fix for saving GIF files +% 1.8 2014-11-16: Fixes for R2014b +% 1.7 2014-04-28: Fixed bug when capturing interactive selection +% 1.6 2014-04-22: Only enable image formats when saving to an unspecified file via uiputfile +% 1.5 2013-04-18: Fixed bug in capture of non-square image; fixes for Win64 +% 1.4 2013-01-27: Fixed capture of Desktop (root); enabled rbbox anywhere on desktop (not necesarily in a Matlab figure); enabled output to clipboard (based on Jiro Doke's imclipboard utility); edge-case fixes; added Java compatibility check +% 1.3 2012-07-23: Capture current object (uicontrol/axes/figure) if w=h=0 (e.g., by clicking a single point); extra input args sanity checks; fix for docked windows and image axes; include axes labels & ticks by default when capturing axes; use data-units position vector when capturing images; many edge-case fixes +% 1.2 2011-01-16: another performance boost (thanks to Jan Simon); some compatibility fixes for Matlab 6.5 (untested) +% 1.1 2009-06-03: Handle missing output format; performance boost (thanks to Urs); fix minor root-handle bug; added toolbar button option +% 1.0 2009-06-02: First version posted on MathWorks File Exchange + +% License to use and modify this code is granted freely to all interested, as long as the original author is +% referenced and attributed as such. The original author maintains the right to be solely associated with this work. + +% Programmed and Copyright by Yair M. Altman: altmany(at)gmail.com +% $Revision: 1.17 $ $Date: 2016/05/16 17:59:36 $ + + % Ensure that java awt is enabled... + if ~usejava('awt') + error('YMA:screencapture:NeedAwt','ScreenCapture requires Java to run.'); + end + + % Ensure that our Java version supports the Robot class (requires JVM 1.3+) + try + robot = java.awt.Robot; %#ok + catch + uiwait(msgbox({['Your Matlab installation is so old that its Java engine (' version('-java') ... + ') does not have a java.awt.Robot class. '], ' ', ... + 'Without this class, taking a screen-capture is impossible.', ' ', ... + 'So, either install JVM 1.3 or higher, or use a newer Matlab release.'}, ... + 'ScreenCapture', 'warn')); + if nargout, imageData = []; end + return; + end + + % Process optional arguments + paramsStruct = processArgs(varargin{:}); + + % If toolbar button requested, add it and exit + if ~isempty(paramsStruct.toolbar) + + % Add the toolbar button + addToolbarButton(paramsStruct); + + % Return the figure to its pre-undocked state (when relevant) + redockFigureIfRelevant(paramsStruct); + + % Exit immediately (do NOT take a screen-capture) + if nargout, imageData = []; end + return; + end + + % Convert position from handle-relative to desktop Java-based pixels + [paramsStruct, msgStr] = convertPos(paramsStruct); + + % Capture the requested screen rectangle using java.awt.Robot + imgData = getScreenCaptureImageData(paramsStruct.position); + + % Return the figure to its pre-undocked state (when relevant) + redockFigureIfRelevant(paramsStruct); + + % Save image data in file or clipboard, if specified + if ~isempty(paramsStruct.target) + if strcmpi(paramsStruct.target,'clipboard') + if ~isempty(imgData) + imclipboard(imgData); + else + msgbox('No image area selected - not copying image to clipboard','ScreenCapture','warn'); + end + elseif strncmpi(paramsStruct.target,'print',5) % 'print' or 'printer' + if ~isempty(imgData) + hNewFig = figure('visible','off'); + imshow(imgData); + print(hNewFig); + delete(hNewFig); + else + msgbox('No image area selected - not printing screenshot','ScreenCapture','warn'); + end + else % real filename + if ~isempty(imgData) + imwrite(imgData,paramsStruct.target); + else + msgbox(['No image area selected - not saving image file ' paramsStruct.target],'ScreenCapture','warn'); + end + end + end + + % Return image raster data to user, if requested + if nargout + imageData = imgData; + + % If neither output formats was specified (neither target nor output data) + elseif isempty(paramsStruct.target) & ~isempty(imgData) %#ok ML6 + % Ask the user to specify a file + %error('YMA:screencapture:noOutput','No output specified for ScreenCapture: specify the output filename and/or output data'); + %format = '*.*'; + formats = imformats; + for idx = 1 : numel(formats) + ext = sprintf('*.%s;',formats(idx).ext{:}); + format(idx,1:2) = {ext(1:end-1), formats(idx).description}; %#ok + end + [filename,pathname] = uiputfile(format,'Save screen capture as'); + if ~isequal(filename,0) & ~isequal(pathname,0) %#ok Matlab6 compatibility + try + filename = fullfile(pathname,filename); + imwrite(imgData,filename); + catch % possibly a GIF file that requires indexed colors + [imgData,map] = rgb2ind(imgData,256); + imwrite(imgData,map,filename); + end + else + % TODO - copy to clipboard + end + end + + % Display msgStr, if relevant + if ~isempty(msgStr) + uiwait(msgbox(msgStr,'ScreenCapture')); + drawnow; pause(0.05); % time for the msgbox to disappear + end + + return; % debug breakpoint + +%% Process optional arguments +function paramsStruct = processArgs(varargin) + + % Get the properties in either direct or P-V format + [regParams, pvPairs] = parseparams(varargin); + + % Now process the optional P-V params + try + % Initialize + paramName = []; + paramsStruct = []; + paramsStruct.handle = []; + paramsStruct.position = []; + paramsStruct.target = ''; + paramsStruct.toolbar = []; + paramsStruct.wasDocked = 0; % no false available in ML6 + paramsStruct.wasInteractive = 0; % no false available in ML6 + + % Parse the regular (non-named) params in recption order + if ~isempty(regParams) & (isempty(regParams{1}) | ishandle(regParams{1}(1))) %#ok ML6 + paramsStruct.handle = regParams{1}; + regParams(1) = []; + end + if ~isempty(regParams) & isnumeric(regParams{1}) & (length(regParams{1}) == 4) %#ok ML6 + paramsStruct.position = regParams{1}; + regParams(1) = []; + end + if ~isempty(regParams) & ischar(regParams{1}) %#ok ML6 + paramsStruct.target = regParams{1}; + end + + % Parse the optional param PV pairs + supportedArgs = {'handle','position','target','toolbar'}; + while ~isempty(pvPairs) + + % Disregard empty propNames (may be due to users mis-interpretting the syntax help) + while ~isempty(pvPairs) & isempty(pvPairs{1}) %#ok ML6 + pvPairs(1) = []; + end + if isempty(pvPairs) + break; + end + + % Ensure basic format is valid + paramName = ''; + if ~ischar(pvPairs{1}) + error('YMA:screencapture:invalidProperty','Invalid property passed to ScreenCapture'); + elseif length(pvPairs) == 1 + if isempty(paramsStruct.target) + paramsStruct.target = pvPairs{1}; + break; + else + error('YMA:screencapture:noPropertyValue',['No value specified for property ''' pvPairs{1} '''']); + end + end + + % Process parameter values + paramName = pvPairs{1}; + if strcmpi(paramName,'filename') % backward compatibility + paramName = 'target'; + end + paramValue = pvPairs{2}; + pvPairs(1:2) = []; + idx = find(strncmpi(paramName,supportedArgs,length(paramName))); + if ~isempty(idx) + %paramsStruct.(lower(supportedArgs{idx(1)})) = paramValue; % incompatible with ML6 + paramsStruct = setfield(paramsStruct, lower(supportedArgs{idx(1)}), paramValue); %#ok ML6 + + % If 'toolbar' param specified, then it cannot be left empty - use gcf + if strncmpi(paramName,'toolbar',length(paramName)) & isempty(paramsStruct.toolbar) %#ok ML6 + paramsStruct.toolbar = getCurrentFig; + end + + elseif isempty(paramsStruct.target) + paramsStruct.target = paramName; + pvPairs = {paramValue, pvPairs{:}}; %#ok (more readable this way, although a bit less efficient...) + + else + supportedArgsStr = sprintf('''%s'',',supportedArgs{:}); + error('YMA:screencapture:invalidProperty','%s \n%s', ... + 'Invalid property passed to ScreenCapture', ... + ['Supported property names are: ' supportedArgsStr(1:end-1)]); + end + end % loop pvPairs + + catch + if ~isempty(paramName), paramName = [' ''' paramName '''']; end + error('YMA:screencapture:invalidProperty','Error setting ScreenCapture property %s:\n%s',paramName,lasterr); %#ok + end +%end % processArgs + +%% Convert position from handle-relative to desktop Java-based pixels +function [paramsStruct, msgStr] = convertPos(paramsStruct) + msgStr = ''; + try + % Get the screen-size for later use + screenSize = get(0,'ScreenSize'); + + % Get the containing figure's handle + hParent = paramsStruct.handle; + if isempty(paramsStruct.handle) + paramsStruct.hFigure = getCurrentFig; + hParent = paramsStruct.hFigure; + else + paramsStruct.hFigure = ancestor(paramsStruct.handle,'figure'); + end + + % To get the acurate pixel position, the figure window must be undocked + try + if strcmpi(get(paramsStruct.hFigure,'WindowStyle'),'docked') + set(paramsStruct.hFigure,'WindowStyle','normal'); + drawnow; pause(0.25); + paramsStruct.wasDocked = 1; % no true available in ML6 + end + catch + % never mind - ignore... + end + + % The figure (if specified) must be in focus + if ~isempty(paramsStruct.hFigure) & ishandle(paramsStruct.hFigure) %#ok ML6 + isFigureValid = 1; % no true available in ML6 + figure(paramsStruct.hFigure); + else + isFigureValid = 0; % no false available in ML6 + end + + % Flush all graphic events to ensure correct rendering + drawnow; pause(0.01); + + % No handle specified + wasPositionGiven = 1; % no true available in ML6 + if isempty(paramsStruct.handle) + + % Set default handle, if not supplied + paramsStruct.handle = paramsStruct.hFigure; + + % If position was not specified, get it interactively using RBBOX + if isempty(paramsStruct.position) + [paramsStruct.position, jFrameUsed, msgStr] = getInteractivePosition(paramsStruct.hFigure); %#ok jFrameUsed is unused + paramsStruct.wasInteractive = 1; % no true available in ML6 + wasPositionGiven = 0; % no false available in ML6 + end + + elseif ~ishandle(paramsStruct.handle) + % Handle was supplied - ensure it is a valid handle + error('YMA:screencapture:invalidHandle','Invalid handle passed to ScreenCapture'); + + elseif isempty(paramsStruct.position) + % Handle was supplied but position was not, so use the handle's position + paramsStruct.position = getPixelPos(paramsStruct.handle); + paramsStruct.position(1:2) = 0; + wasPositionGiven = 0; % no false available in ML6 + + elseif ~isnumeric(paramsStruct.position) | (length(paramsStruct.position) ~= 4) %#ok ML6 + % Both handle & position were supplied - ensure a valid pixel position vector + error('YMA:screencapture:invalidPosition','Invalid position vector passed to ScreenCapture: \nMust be a [x,y,w,h] numeric pixel array'); + end + + % Capture current object (uicontrol/axes/figure) if w=h=0 (single-click in interactive mode) + if paramsStruct.position(3)<=0 | paramsStruct.position(4)<=0 %#ok ML6 + %TODO - find a way to single-click another Matlab figure (the following does not work) + %paramsStruct.position = getPixelPos(ancestor(hittest,'figure')); + paramsStruct.position = getPixelPos(paramsStruct.handle); + paramsStruct.position(1:2) = 0; + paramsStruct.wasInteractive = 0; % no false available in ML6 + wasPositionGiven = 0; % no false available in ML6 + end + + % First get the parent handle's desktop-based Matlab pixel position + parentPos = [0,0,0,0]; + dX = 0; + dY = 0; + dW = 0; + dH = 0; + if ~isFigure(hParent) + % Get the reguested component's pixel position + parentPos = getPixelPos(hParent, 1); % no true available in ML6 + + % Axes position inaccuracy estimation + deltaX = 3; + deltaY = -1; + + % Fix for images + if isImage(hParent) % | (isAxes(hParent) & strcmpi(get(hParent,'YDir'),'reverse')) %#ok ML6 + + % Compensate for resized image axes + hAxes = get(hParent,'Parent'); + if all(get(hAxes,'DataAspectRatio')==1) % sanity check: this is the normal behavior + % Note 18/4/2013: the following fails for non-square images + %actualImgSize = min(parentPos(3:4)); + %dX = (parentPos(3) - actualImgSize) / 2; + %dY = (parentPos(4) - actualImgSize) / 2; + %parentPos(3:4) = actualImgSize; + + % The following should work for all types of images + actualImgSize = size(get(hParent,'CData')); + dX = (parentPos(3) - min(parentPos(3),actualImgSize(2))) / 2; + dY = (parentPos(4) - min(parentPos(4),actualImgSize(1))) / 2; + parentPos(3:4) = actualImgSize([2,1]); + %parentPos(3) = max(parentPos(3),actualImgSize(2)); + %parentPos(4) = max(parentPos(4),actualImgSize(1)); + end + + % Fix user-specified img positions (but not auto-inferred ones) + if wasPositionGiven + + % In images, use data units rather than pixel units + % Reverse the YDir + ymax = max(get(hParent,'YData')); + paramsStruct.position(2) = ymax - paramsStruct.position(2) - paramsStruct.position(4); + + % Note: it would be best to use hgconvertunits, but: + % ^^^^ (1) it fails on Matlab 6, and (2) it doesn't accept Data units + %paramsStruct.position = hgconvertunits(hFig, paramsStruct.position, 'Data', 'pixel', hParent); % fails! + xLims = get(hParent,'XData'); + yLims = get(hParent,'YData'); + xPixelsPerData = parentPos(3) / (diff(xLims) + 1); + yPixelsPerData = parentPos(4) / (diff(yLims) + 1); + paramsStruct.position(1) = round((paramsStruct.position(1)-xLims(1)) * xPixelsPerData); + paramsStruct.position(2) = round((paramsStruct.position(2)-yLims(1)) * yPixelsPerData + 2*dY); + paramsStruct.position(3) = round( paramsStruct.position(3) * xPixelsPerData); + paramsStruct.position(4) = round( paramsStruct.position(4) * yPixelsPerData); + + % Axes position inaccuracy estimation + if strcmpi(computer('arch'),'win64') + deltaX = 7; + deltaY = -7; + else + deltaX = 3; + deltaY = -3; + end + + else % axes/image position was auto-infered (entire image) + % Axes position inaccuracy estimation + if strcmpi(computer('arch'),'win64') + deltaX = 6; + deltaY = -6; + else + deltaX = 2; + deltaY = -2; + end + dW = -2*dX; + dH = -2*dY; + end + end + + %hFig = ancestor(hParent,'figure'); + hParent = paramsStruct.hFigure; + + elseif paramsStruct.wasInteractive % interactive figure rectangle + + % Compensate for 1px rbbox inaccuracies + deltaX = 2; + deltaY = -2; + + else % non-interactive figure + + % Compensate 4px figure boundaries = difference betweeen OuterPosition and Position + deltaX = -1; + deltaY = 1; + end + %disp(paramsStruct.position) % for debugging + + % Now get the pixel position relative to the monitor + figurePos = getPixelPos(hParent); + desktopPos = figurePos + parentPos; + + % Now convert to Java-based pixels based on screen size + % Note: multiple monitors are automatically handled correctly, since all + % ^^^^ Java positions are relative to the main monitor's top-left corner + javaX = desktopPos(1) + paramsStruct.position(1) + deltaX + dX; + javaY = screenSize(4) - desktopPos(2) - paramsStruct.position(2) - paramsStruct.position(4) + deltaY + dY; + width = paramsStruct.position(3) + dW; + height = paramsStruct.position(4) + dH; + paramsStruct.position = round([javaX, javaY, width, height]); + %paramsStruct.position + + % Ensure the figure is at the front so it can be screen-captured + if isFigureValid + figure(hParent); + drawnow; + pause(0.02); + end + catch + % Maybe root/desktop handle (root does not have a 'Position' prop so getPixelPos croaks + if isequal(double(hParent),0) % =root/desktop handle; handles case of hParent=[] + javaX = paramsStruct.position(1) - 1; + javaY = screenSize(4) - paramsStruct.position(2) - paramsStruct.position(4) - 1; + paramsStruct.position = [javaX, javaY, paramsStruct.position(3:4)]; + end + end +%end % convertPos + +%% Interactively get the requested capture rectangle +function [positionRect, jFrameUsed, msgStr] = getInteractivePosition(hFig) + msgStr = ''; + try + % First try the invisible-figure approach, in order to + % enable rbbox outside any existing figure boundaries + f = figure('units','pixel','pos',[-100,-100,10,10],'HitTest','off'); + drawnow; pause(0.01); + oldWarn = warning('off','MATLAB:HandleGraphics:ObsoletedProperty:JavaFrame'); + jf = get(handle(f),'JavaFrame'); + warning(oldWarn); + try + jWindow = jf.fFigureClient.getWindow; + catch + try + jWindow = jf.fHG1Client.getWindow; + catch + jWindow = jf.getFigurePanelContainer.getParent.getTopLevelAncestor; + end + end + com.sun.awt.AWTUtilities.setWindowOpacity(jWindow,0.05); %=nearly transparent (not fully so that mouse clicks are captured) + jWindow.setMaximized(1); % no true available in ML6 + jFrameUsed = 1; % no true available in ML6 + msg = {'Mouse-click and drag a bounding rectangle for screen-capture ' ... + ... %'or single-click any Matlab figure to capture the entire figure.' ... + }; + catch + % Something failed, so revert to a simple rbbox on a visible figure + try delete(f); drawnow; catch, end %Cleanup... + jFrameUsed = 0; % no false available in ML6 + msg = {'Mouse-click within any Matlab figure and then', ... + 'drag a bounding rectangle for screen-capture,', ... + 'or single-click to capture the entire figure'}; + end + uiwait(msgbox(msg,'ScreenCapture')); + + k = waitforbuttonpress; %#ok k is unused + %hFig = getCurrentFig; + %p1 = get(hFig,'CurrentPoint'); + positionRect = rbbox; + %p2 = get(hFig,'CurrentPoint'); + + if jFrameUsed + jFrameOrigin = getPixelPos(f); + delete(f); drawnow; + try + figOrigin = getPixelPos(hFig); + catch % empty/invalid hFig handle + figOrigin = [0,0,0,0]; + end + else + if isempty(hFig) + jFrameOrigin = getPixelPos(gcf); + else + jFrameOrigin = [0,0,0,0]; + end + figOrigin = [0,0,0,0]; + end + positionRect(1:2) = positionRect(1:2) + jFrameOrigin(1:2) - figOrigin(1:2); + + if prod(positionRect(3:4)) > 0 + msgStr = sprintf('%dx%d area captured',positionRect(3),positionRect(4)); + end +%end % getInteractivePosition + +%% Get current figure (even if its handle is hidden) +function hFig = getCurrentFig + oldState = get(0,'showHiddenHandles'); + set(0,'showHiddenHandles','on'); + hFig = get(0,'CurrentFigure'); + set(0,'showHiddenHandles',oldState); +%end % getCurrentFig + +%% Get ancestor figure - used for old Matlab versions that don't have a built-in ancestor() +function hObj = ancestor(hObj,type) + if ~isempty(hObj) & ishandle(hObj) %#ok for Matlab 6 compatibility + try + hObj = get(hObj,'Ancestor'); + catch + % never mind... + end + try + %if ~isa(handle(hObj),type) % this is best but always returns 0 in Matlab 6! + %if ~isprop(hObj,'type') | ~strcmpi(get(hObj,'type'),type) % no isprop() in ML6! + try + objType = get(hObj,'type'); + catch + objType = ''; + end + if ~strcmpi(objType,type) + try + parent = get(handle(hObj),'parent'); + catch + parent = hObj.getParent; % some objs have no 'Parent' prop, just this method... + end + if ~isempty(parent) % empty parent means root ancestor, so exit + hObj = ancestor(parent,type); + end + end + catch + % never mind... + end + end +%end % ancestor + +%% Get position of an HG object in specified units +function pos = getPos(hObj,field,units) + % Matlab 6 did not have hgconvertunits so use the old way... + oldUnits = get(hObj,'units'); + if strcmpi(oldUnits,units) % don't modify units unless we must! + pos = get(hObj,field); + else + set(hObj,'units',units); + pos = get(hObj,field); + set(hObj,'units',oldUnits); + end +%end % getPos + +%% Get pixel position of an HG object - for Matlab 6 compatibility +function pos = getPixelPos(hObj,varargin) + persistent originalObj + try + stk = dbstack; + if ~strcmp(stk(2).name,'getPixelPos') + originalObj = hObj; + end + + if isFigure(hObj) %| isAxes(hObj) + %try + pos = getPos(hObj,'OuterPosition','pixels'); + else %catch + % getpixelposition is unvectorized unfortunately! + pos = getpixelposition(hObj,varargin{:}); + + % add the axes labels/ticks if relevant (plus a tiny margin to fix 2px label/title inconsistencies) + if isAxes(hObj) & ~isImage(originalObj) %#ok ML6 + tightInsets = getPos(hObj,'TightInset','pixel'); + pos = pos + tightInsets.*[-1,-1,1,1] + [-1,1,1+tightInsets(1:2)]; + end + end + catch + try + % Matlab 6 did not have getpixelposition nor hgconvertunits so use the old way... + pos = getPos(hObj,'Position','pixels'); + catch + % Maybe the handle does not have a 'Position' prop (e.g., text/line/plot) - use its parent + pos = getPixelPos(get(hObj,'parent'),varargin{:}); + end + end + + % Handle the case of missing/invalid/empty HG handle + if isempty(pos) + pos = [0,0,0,0]; + end +%end % getPixelPos + +%% Adds a ScreenCapture toolbar button +function addToolbarButton(paramsStruct) + % Ensure we have a valid toolbar handle + hFig = ancestor(paramsStruct.toolbar,'figure'); + if isempty(hFig) + error('YMA:screencapture:badToolbar','the ''Toolbar'' parameter must contain a valid GUI handle'); + end + set(hFig,'ToolBar','figure'); + hToolbar = findall(hFig,'type','uitoolbar'); + if isempty(hToolbar) + error('YMA:screencapture:noToolbar','the ''Toolbar'' parameter must contain a figure handle possessing a valid toolbar'); + end + hToolbar = hToolbar(1); % just in case there are several toolbars... - use only the first + + % Prepare the camera icon + icon = ['3333333333333333'; ... + '3333333333333333'; ... + '3333300000333333'; ... + '3333065556033333'; ... + '3000000000000033'; ... + '3022222222222033'; ... + '3022220002222033'; ... + '3022203110222033'; ... + '3022201110222033'; ... + '3022204440222033'; ... + '3022220002222033'; ... + '3022222222222033'; ... + '3000000000000033'; ... + '3333333333333333'; ... + '3333333333333333'; ... + '3333333333333333']; + cm = [ 0 0 0; ... % black + 0 0.60 1; ... % light blue + 0.53 0.53 0.53; ... % light gray + NaN NaN NaN; ... % transparent + 0 0.73 0; ... % light green + 0.27 0.27 0.27; ... % gray + 0.13 0.13 0.13]; % dark gray + cdata = ind2rgb(uint8(icon-'0'),cm); + + % If the button does not already exit + hButton = findall(hToolbar,'Tag','ScreenCaptureButton'); + tooltip = 'Screen capture'; + if ~isempty(paramsStruct.target) + tooltip = [tooltip ' to ' paramsStruct.target]; + end + if isempty(hButton) + % Add the button with the icon to the figure's toolbar + hButton = uipushtool(hToolbar, 'CData',cdata, 'Tag','ScreenCaptureButton', 'TooltipString',tooltip, 'ClickedCallback',['screencapture(''' paramsStruct.target ''')']); %#ok unused + else + % Otherwise, simply update the existing button + set(hButton, 'CData',cdata, 'Tag','ScreenCaptureButton', 'TooltipString',tooltip, 'ClickedCallback',['screencapture(''' paramsStruct.target ''')']); + end +%end % addToolbarButton + +%% Java-get the actual screen-capture image data +function imgData = getScreenCaptureImageData(positionRect) + if isempty(positionRect) | all(positionRect==0) | positionRect(3)<=0 | positionRect(4)<=0 %#ok ML6 + imgData = []; + else + % Use java.awt.Robot to take a screen-capture of the specified screen area + rect = java.awt.Rectangle(positionRect(1), positionRect(2), positionRect(3), positionRect(4)); + robot = java.awt.Robot; + jImage = robot.createScreenCapture(rect); + + % Convert the resulting Java image to a Matlab image + % Adapted for a much-improved performance from: + % http://www.mathworks.com/support/solutions/data/1-2WPAYR.html + h = jImage.getHeight; + w = jImage.getWidth; + %imgData = zeros([h,w,3],'uint8'); + %pixelsData = uint8(jImage.getData.getPixels(0,0,w,h,[])); + %for i = 1 : h + % base = (i-1)*w*3+1; + % imgData(i,1:w,:) = deal(reshape(pixelsData(base:(base+3*w-1)),3,w)'); + %end + + % Performance further improved based on feedback from Urs Schwartz: + %pixelsData = reshape(typecast(jImage.getData.getDataStorage,'uint32'),w,h).'; + %imgData(:,:,3) = bitshift(bitand(pixelsData,256^1-1),-8*0); + %imgData(:,:,2) = bitshift(bitand(pixelsData,256^2-1),-8*1); + %imgData(:,:,1) = bitshift(bitand(pixelsData,256^3-1),-8*2); + + % Performance even further improved based on feedback from Jan Simon: + pixelsData = reshape(typecast(jImage.getData.getDataStorage, 'uint8'), 4, w, h); + imgData = cat(3, ... + transpose(reshape(pixelsData(3, :, :), w, h)), ... + transpose(reshape(pixelsData(2, :, :), w, h)), ... + transpose(reshape(pixelsData(1, :, :), w, h))); + end +%end % getInteractivePosition + +%% Return the figure to its pre-undocked state (when relevant) +function redockFigureIfRelevant(paramsStruct) + if paramsStruct.wasDocked + try + set(paramsStruct.hFigure,'WindowStyle','docked'); + %drawnow; + catch + % never mind - ignore... + end + end +%end % redockFigureIfRelevant + +%% Copy screen-capture to the system clipboard +% Adapted from http://www.mathworks.com/matlabcentral/fileexchange/28708-imclipboard/content/imclipboard.m +function imclipboard(imgData) + % Import necessary Java classes + import java.awt.Toolkit.* + import java.awt.image.BufferedImage + import java.awt.datatransfer.DataFlavor + + % Add the necessary Java class (ImageSelection) to the Java classpath + if ~exist('ImageSelection', 'class') + % Obtain the directory of the executable (or of the M-file if not deployed) + %javaaddpath(fileparts(which(mfilename)), '-end'); + if isdeployed % Stand-alone mode. + [status, result] = system('path'); %#ok + MatLabFilePath = char(regexpi(result, 'Path=(.*?);', 'tokens', 'once')); + else % MATLAB mode. + MatLabFilePath = fileparts(mfilename('fullpath')); + end + javaaddpath(MatLabFilePath, '-end'); + end + + % Get System Clipboard object (java.awt.Toolkit) + cb = getDefaultToolkit.getSystemClipboard; % can't use () in ML6! + + % Get image size + ht = size(imgData, 1); + wd = size(imgData, 2); + + % Convert to Blue-Green-Red format + imgData = imgData(:, :, [3 2 1]); + + % Convert to 3xWxH format + imgData = permute(imgData, [3, 2, 1]); + + % Append Alpha data (not used) + imgData = cat(1, imgData, 255*ones(1, wd, ht, 'uint8')); + + % Create image buffer + imBuffer = BufferedImage(wd, ht, BufferedImage.TYPE_INT_RGB); + imBuffer.setRGB(0, 0, wd, ht, typecast(imgData(:), 'int32'), 0, wd); + + % Create ImageSelection object + % % custom java class + imSelection = ImageSelection(imBuffer); + + % Set clipboard content to the image + cb.setContents(imSelection, []); +%end %imclipboard + +%% Is the provided handle a figure? +function flag = isFigure(hObj) + flag = isa(handle(hObj),'figure') | isa(hObj,'matlab.ui.Figure'); +%end %isFigure + +%% Is the provided handle an axes? +function flag = isAxes(hObj) + flag = isa(handle(hObj),'axes') | isa(hObj,'matlab.graphics.axis.Axes'); +%end %isFigure + +%% Is the provided handle an image? +function flag = isImage(hObj) + flag = isa(handle(hObj),'image') | isa(hObj,'matlab.graphics.primitive.Image'); +%end %isFigure + +%%%%%%%%%%%%%%%%%%%%%%%%%% TODO %%%%%%%%%%%%%%%%%%%%%%%%% +% find a way in interactive-mode to single-click another Matlab figure for screen-capture \ No newline at end of file diff --git a/MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m b/MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m new file mode 100644 index 0000000..0d2692a --- /dev/null +++ b/MOT Capture Process Simulation/+Plotting/plotPositionAndVelocitySampling.m @@ -0,0 +1,56 @@ +function plotPositionAndVelocitySampling(Simulator, NumberOfBins) + + f_h = Helper.getFigureByTag('RejectionSampling'); + set(groot,'CurrentFigure',f_h); + a_h = get(f_h, 'CurrentAxes'); + if ~isempty(get(a_h, 'Children')) + clf(f_h); + end + f_h.Name = 'Sampling'; + f_h.Units = 'pixels'; + + set(0,'units','pixels'); + screensize = get(0,'ScreenSize'); + f_h.Position = [[screensize(3)/7 screensize(4)/7] 1.357e+03 770]; + + initialPositions = Simulator.InitialPositions; + initialVelocities = Simulator.InitialVelocities; + + subplot(3,2,1) + histogram(initialPositions(:, 1)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','x-Component') + xlabel('Positions (mm)','FontSize', 14) + ylabel('Counts','FontSize', 14) + legend('FontSize', 14) + + subplot(3,2,3) + histogram(initialPositions(:, 2)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','y-Component') + xlabel('Positions (mm)','FontSize', 14) + ylabel('Counts','FontSize', 14) + legend('FontSize', 14) + + subplot(3,2,5) + histogram(initialPositions(:, 3)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','z-Component') + xlabel('Positions (mm)','FontSize', 14) + ylabel('Counts','FontSize', 14) + legend('FontSize', 14) + + subplot(3,2,2) + histogram(initialVelocities(:, 1),NumberOfBins, 'LineStyle', 'none', 'DisplayName','x-Component') + xlabel('Velocities (m/s)','FontSize', 14) + ylabel('Counts','FontSize', 14) + legend('FontSize', 14) + + subplot(3,2,4) + histogram(initialVelocities(:, 2)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','y-Component') + xlabel('Velocities (mm/s)','FontSize', 14) + ylabel('Counts','FontSize', 14) + legend('FontSize', 14) + + subplot(3,2,6) + histogram(initialVelocities(:, 3)*1e3,NumberOfBins, 'LineStyle', 'none', 'DisplayName','z-Component') + xlabel('Velocities (mm/s)','FontSize', 14) + ylabel('Counts','FontSize', 14) + legend('FontSize', 14) + + sgtitle('Rejection sampling for initial distributions','FontSize', 18) +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m b/MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m new file mode 100644 index 0000000..38ca88e --- /dev/null +++ b/MOT Capture Process Simulation/+Plotting/plotResultForTwoParameterScan.m @@ -0,0 +1,65 @@ +function plotResultForTwoParameterScan(XParameter, YParameter, ZQuantity, varargin) + + p = inputParser; + p.KeepUnmatched = true; + + addRequired(p, 'FirstParameterArray', @isvector) + addRequired(p, 'SecondParameterArray', @isvector) + addRequired(p, 'QuantityOfInterestArray', @ismatrix) + + addParameter(p, 'RescalingFactorForFirstParameter', 1, @isscalar) + addParameter(p, 'XLabelString', 'X parameter', @ischar) + addParameter(p, 'RescalingFactorForSecondParameter', 1, @isscalar) + addParameter(p, 'YLabelString', 'Y parameter', @ischar) + addParameter(p, 'ZLabelString', 'Z parameter', @ischar) + addParameter(p, 'TitleString', 'Two-Parameter Scan', @ischar) + + p.parse(XParameter, YParameter, ZQuantity, varargin{:}) + + XParameter = p.Results.FirstParameterArray; + RescalingFactorForXParameter = p.Results.RescalingFactorForFirstParameter; + XLabelString = p.Results.XLabelString; + YParameter = p.Results.SecondParameterArray; + RescalingFactorForYParameter = p.Results.RescalingFactorForSecondParameter; + YLabelString = p.Results.YLabelString; + ZQuantity = p.Results.QuantityOfInterestArray; + ZLabelString = p.Results.ZLabelString; + TitleString = p.Results.TitleString; + + f_h = Helper.getFigureByTag('Two-Parameter Scan'); + set(groot,'CurrentFigure',f_h); + a_h = get(f_h, 'CurrentAxes'); + if ~isempty(get(a_h, 'Children')) + clf(f_h); + end + f_h.Name = 'Two-Parameter Scan'; + 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]; + + RescaledXParameter = XParameter .* RescalingFactorForXParameter; + RescaledYParameter = YParameter .* RescalingFactorForYParameter; + + imagesc(RescaledXParameter, RescaledYParameter, ZQuantity(:,:)'); + + set(gca,'YDir','normal'); + + caxis([min(min(min(ZQuantity))) max(max(max(ZQuantity)))]); + + hXLabel = xlabel(XLabelString); + hYLabel = ylabel(YLabelString); + + shading flat; + c = colorbar; + c.Label.String= ZLabelString; + c.Label.FontSize = 14; + + hTitle = sgtitle(TitleString); + + set([hXLabel, hYLabel] , ... + 'FontSize' , 14 ); + set( hTitle , ... + 'FontSize' , 18 ); +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/+Plotting/visualizeMagneticField.m b/MOT Capture Process Simulation/+Plotting/visualizeMagneticField.m new file mode 100644 index 0000000..78f895a --- /dev/null +++ b/MOT Capture Process Simulation/+Plotting/visualizeMagneticField.m @@ -0,0 +1,70 @@ +function visualizeMagneticField(obj, x_range, y_range, z_range) + + f_h = Helper.getFigureByTag('VisualizeMagneticFieldFor2DMOT'); + set(groot,'CurrentFigure',f_h); + a_h = get(f_h, 'CurrentAxes'); + if ~isempty(get(a_h, 'Children')) + clf(f_h); + end + f_h.Name = 'Visualization'; + f_h.Units = 'pixels'; + + set(0,'units','pixels'); + screensize = get(0,'ScreenSize'); + f_h.Position = [[screensize(3)/3.5 screensize(4)/3.5] 820 645]; + + xmin = x_range(1); + xmax = x_range(2); + ymin = y_range(1); + ymax = y_range(2); + zmin = z_range(1); + zmax = z_range(2); + dx = (xmax-xmin)/8; + dy = (ymax-ymin)/8; + dz = (zmax-zmin)/8; + if dx ~= 0 + xm = xmin:dx:xmax; + else + xm = zeros(1,5); + end + + if dy ~= 0 + ym = ymin:dy:ymax; + else + ym = zeros(1,5); + end + + if dz ~= 0 + zm = zmin:dz:zmax; + else + zm = zeros(1,5); + end + [meshx,meshy,meshz] = meshgrid(xm,ym,zm); % construct data points + + switch obj.SimulationMode + case '2D' + alpha = obj.MagneticGradient; + Bx = @(x,y,z) alpha .* z; + By = @(x,y,z) 0 .* y; + Bz = @(x,y,z) alpha .* x; + Bx_val = Bx(meshx, meshy, meshz); + By_val = By(meshx, meshy, meshz); + Bz_val = Bz(meshx, meshy, meshz); + case '3D' + % Development in progress + end + + quiver3(meshx, meshy, meshz, Bx_val, By_val, Bz_val, 'Color', ' #6600ff'); + axis equal + + hXLabel = xlabel('x'); + hYLabel = ylabel('y'); + hZLabel = zlabel('z'); + + hTitle = sgtitle('Magnetic Field for 2-D MOT'); + + set([hXLabel, hYLabel, hZLabel] , ... + 'FontSize' , 14 ); + set( hTitle , ... + 'FontSize' , 18 ); +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m b/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m new file mode 100644 index 0000000..35b85ec --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/MOTSimulator.m @@ -0,0 +1,616 @@ +classdef MOTSimulator < handle & matlab.mixin.Copyable + + properties (Access = public) + SimulationMode; % MOT type + TimeStep; + SimulationTime; + NumberOfAtoms; + + InitialPositions; + InitialVelocities; + + NozzleLength; + NozzleRadius; + Beta; + ApertureCut; + OvenDistance; + OvenTemperature; + MagneticGradient; + NozzleExitDivergence; + MOTExitDivergence; + MOTDistance; + + BluePower; + BlueDetuning; + BlueBeamRadius; + BlueBeamWaist; + BlueWaveVector; + BlueSaturationIntensity; + + OrangePower; + OrangeDetuning; + OrangeBeamRadius; + OrangeBeamWaist; + OrangeWaveVector; + OrangeSaturationIntensity; + + CoolingBeamPower; + CoolingBeamWaveVector; + CoolingBeamLinewidth; + CoolingBeamDetuning; + CoolingBeamRadius; + CoolingBeamWaist; + CoolingBeamSaturationIntensity; + + SidebandPower; + SidebandDetuning; + SidebandBeamRadius; + SidebandBeamWaist; + SidebandBeamSaturationIntensity; + + PushBeamPower; + PushBeamWaveVector; + PushBeamLinewidth; + PushBeamDetuning; + PushBeamRadius; + PushBeamWaist; + PushBeamDistance; + DistanceBetweenPushBeamAnd3DMOTCenter; + PushBeamSaturationIntensity; + + ZeemanSlowerBeamPower; + ZeemanSlowerBeamDetuning; + ZeemanSlowerBeamRadius; + ZeemanSlowerBeamWaist; + ZeemanSlowerBeamSaturationIntensity; + + TotalPower; + LandegFactor; + MagneticSubLevel; + + CaptureVelocity; + VelocityCutoff; + CaptureFraction; + ClausingFactor; + AngularDistribution; + NormalizationConstantForAngularDistribution; + + %Flags + SpontaneousEmission; + Sideband; + ZeemanSlowerBeam; + Gravity; + + DebugMode; + DoSave; + end + + properties (SetAccess = private, GetAccess = public) + InitialParameters + end + + properties (Dependent, SetAccess = private) + CoolingBeamSaturationParameter; + SidebandSaturationParameter; + PushBeamSaturationParameter; + ZeemanSlowerBeamSaturationParameter; + OvenTemperatureinKelvin; + AverageVelocity; + AtomicBeamDensity; + MeanFreePath; + CollisionTime; + end + + methods + function s = MOTSimulator(varargin) + + p = inputParser; + p.KeepUnmatched = true; + addParameter(p, 'SimulationMode', '2D',... + @(x) any(strcmpi(x,{'2D','3D'}))); + addParameter(p, 'TimeStep', 10e-06,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'SimulationTime', 3e-03,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'SpontaneousEmission', false,... + @islogical); + addParameter(p, 'Sideband', false,... + @islogical); + addParameter(p, 'ZeemanSlowerBeam', false,... + @islogical); + addParameter(p, 'Gravity', false,... + @islogical); + addParameter(p, 'DebugMode', false,... + @islogical); + addParameter(p, 'SaveData', false,... + @islogical); + + p.parse(varargin{:}); + + s.SimulationMode = p.Results.SimulationMode; + + s.TimeStep = p.Results.TimeStep; + s.SimulationTime = p.Results.SimulationTime; + + s.SpontaneousEmission = p.Results.SpontaneousEmission; + s.Sideband = p.Results.Sideband; + s.ZeemanSlowerBeam = p.Results.ZeemanSlowerBeam; + s.Gravity = p.Results.Gravity; + + s.DebugMode = p.Results.DebugMode; + s.DoSave = p.Results.SaveData; + + s.reinitializeSimulator(); + + poolobj = gcp('nocreate'); % Check if pool is open + if isempty(poolobj) + parpool; + end + + end + end % - lifecycle + + methods + function set.TimeStep(this, val) + assert(val > 1e-06, 'Not time efficient to compute for time steps smaller than 1 microsecond!'); + this.TimeStep = val; + end + function ret = get.TimeStep(this) + ret = this.TimeStep; + end + function set.SimulationTime(this, val) + assert(val <= 5e-03, 'Not time efficient to compute for time spans longer than 5 milliseconds!'); + this.SimulationTime = val; + end + function ret = get.SimulationTime(this) + ret = this.SimulationTime; + end + function set.NumberOfAtoms(this, val) + assert(val <= 10000, 'Not time efficient to compute for atom numbers larger than 10,000!'); + this.NumberOfAtoms = val; + end + function ret = get.NumberOfAtoms(this) + ret = this.NumberOfAtoms; + end + + function set.InitialPositions(this,val) + this.InitialPositions = val; + end + function ret = get.InitialPositions(this) + ret = this.InitialPositions; + end + function set.InitialVelocities(this,val) + this.InitialVelocities = val; + end + function ret = get.InitialVelocities(this) + ret = this.InitialVelocities; + end + + function set.NozzleLength(this,val) + this.NozzleLength = val; + end + function ret = get.NozzleLength(this) + ret = this.NozzleLength; + end + function set.NozzleRadius(this,val) + this.NozzleRadius = val; + end + function ret = get.NozzleRadius(this) + ret = this.NozzleRadius; + end + function set.Beta(this,val) + this.Beta = val; + end + function ret = get.Beta(this) + ret = this.Beta; + end + function set.ApertureCut(this,val) + this.ApertureCut = val; + end + function ret = get.ApertureCut(this) + ret = this.ApertureCut; + end + function set.OvenDistance(this,val) + this.OvenDistance = val; + end + function ret = get.OvenDistance(this) + ret = this.OvenDistance; + end + function set.OvenTemperature(this,val) + this.OvenTemperature = val; + end + function ret = get.OvenTemperature(this) + ret = this.OvenTemperature; + end + function set.MagneticGradient(this,val) + this.MagneticGradient = val; + end + function ret = get.MagneticGradient(this) + ret = this.MagneticGradient; + end + function set.NozzleExitDivergence(this,val) + this.NozzleExitDivergence = val; + end + function ret = get.NozzleExitDivergence(this) + ret = this.NozzleExitDivergence; + end + function set.MOTExitDivergence(this,val) + this.MOTExitDivergence = val; + end + function ret = get.MOTExitDivergence(this) + ret = this.MOTExitDivergence; + end + function set.MOTDistance(this,val) + this.MOTDistance = val; + end + function ret = get.MOTDistance(this) + ret = this.MOTDistance; + end + + function set.BluePower(this,val) + this.BluePower = val; + end + function ret = get.BluePower(this) + ret = this.BluePower; + end + function set.BlueDetuning(this, val) + this.BlueDetuning = val; + end + function ret = get.BlueDetuning(this) + ret = this.BlueDetuning; + end + function set.BlueBeamRadius(this, val) + this.BlueBeamRadius = val; + end + function ret = get.BlueBeamRadius(this) + ret = this.BlueBeamRadius; + end + function set.BlueBeamWaist(this, val) + this.BlueBeamWaist = val; + end + function ret = get.BlueBeamWaist(this) + ret = this.BlueBeamWaist; + end + function set.BlueWaveVector(this, val) + this.BlueWaveVector = val; + end + function ret = get.BlueWaveVector(this) + ret = this.BlueWaveVector; + end + function set.BlueSaturationIntensity(this, val) + this.BlueSaturationIntensity = val; + end + function ret = get.BlueSaturationIntensity(this) + ret = this.BlueSaturationIntensity; + end + + function set.OrangePower(this,val) + this.OrangePower = val; + end + function ret = get.OrangePower(this) + ret = this.OrangePower; + end + function set.OrangeDetuning(this, val) + this.OrangeDetuning = val; + end + function ret = get.OrangeDetuning(this) + ret = this.OrangeDetuning; + end + function set.OrangeBeamRadius(this, val) + this.OrangeBeamRadius = val; + end + function ret = get.OrangeBeamRadius(this) + ret = this.OrangeBeamRadius; + end + function set.OrangeBeamWaist(this, val) + this.OrangeBeamWaist = val; + end + function ret = get.OrangeBeamWaist(this) + ret = this.OrangeBeamWaist; + end + function set.OrangeWaveVector(this, val) + this.OrangeWaveVector = val; + end + function ret = get.OrangeWaveVector(this) + ret = this.OrangeWaveVector; + end + function set.OrangeSaturationIntensity(this, val) + this.OrangeSaturationIntensity = val; + end + function ret = get.OrangeSaturationIntensity(this) + ret = this.OrangeSaturationIntensity; + end + + function set.CoolingBeamPower(this,val) + this.CoolingBeamPower = val; + end + function ret = get.CoolingBeamPower(this) + ret = this.CoolingBeamPower; + end + function set.CoolingBeamDetuning(this, val) + this.CoolingBeamDetuning = val; + end + function ret = get.CoolingBeamDetuning(this) + ret = this.CoolingBeamDetuning; + end + function set.CoolingBeamRadius(this, val) + this.CoolingBeamRadius = val; + end + function ret = get.CoolingBeamRadius(this) + ret = this.CoolingBeamRadius; + end + function set.CoolingBeamWaist(this, val) + this.CoolingBeamWaist = val; + end + function ret = get.CoolingBeamWaist(this) + ret = this.CoolingBeamWaist; + end + function set.CoolingBeamWaveVector(this, val) + this.CoolingBeamWaveVector = val; + end + function ret = get.CoolingBeamWaveVector(this) + ret = this.CoolingBeamWaveVector; + end + function set.CoolingBeamLinewidth(this, val) + this.CoolingBeamLinewidth = val; + end + function ret = get.CoolingBeamLinewidth(this) + ret = this.CoolingBeamLinewidth; + end + function set.CoolingBeamSaturationIntensity(this, val) + this.CoolingBeamSaturationIntensity = val; + end + function ret = get.CoolingBeamSaturationIntensity(this) + ret = this.CoolingBeamSaturationIntensity; + end + + function set.SidebandPower(this,val) + this.SidebandPower = val; + end + function ret = get.SidebandPower(this) + ret = this.SidebandPower; + end + function set.SidebandDetuning(this, val) + this.SidebandDetuning = val; + end + function ret = get.SidebandDetuning(this) + ret = this.SidebandDetuning; + end + function set.SidebandBeamRadius(this, val) + this.SidebandBeamRadius = val; + end + function ret = get.SidebandBeamRadius(this) + ret = this.SidebandBeamRadius; + end + function set.SidebandBeamWaist(this, val) + this.SidebandBeamWaist = val; + end + function ret = get.SidebandBeamWaist(this) + ret = this.SidebandBeamWaist; + end + function set.SidebandBeamSaturationIntensity(this, val) + this.SidebandBeamSaturationIntensity = val; + end + function ret = get.SidebandBeamSaturationIntensity(this) + ret = this.SidebandBeamSaturationIntensity; + end + + function set.PushBeamPower(this,val) + this.PushBeamPower = val; + end + function ret = get.PushBeamPower(this) + ret = this.PushBeamPower; + end + function set.PushBeamDetuning(this, val) + this.PushBeamDetuning = val; + end + function ret = get.PushBeamDetuning(this) + ret = this.PushBeamDetuning; + end + function set.PushBeamRadius(this, val) + this.PushBeamRadius = val; + end + function ret = get.PushBeamRadius(this) + ret = this.PushBeamRadius; + end + function set.PushBeamWaist(this, val) + this.PushBeamWaist = val; + end + function ret = get.PushBeamWaist(this) + ret = this.PushBeamWaist; + end + function set.PushBeamWaveVector(this, val) + this.PushBeamWaveVector = val; + end + function ret = get.PushBeamWaveVector(this) + ret = this.PushBeamWaveVector; + end + function set.PushBeamLinewidth(this, val) + this.PushBeamLinewidth = val; + end + function ret = get.PushBeamLinewidth(this) + ret = this.PushBeamLinewidth; + end + function set.PushBeamDistance(this, val) + this.PushBeamDistance = val; + end + function ret = get.PushBeamDistance(this) + ret = this.PushBeamDistance; + end + function set.DistanceBetweenPushBeamAnd3DMOTCenter(this, val) + this.DistanceBetweenPushBeamAnd3DMOTCenter = val; + end + function ret = get.DistanceBetweenPushBeamAnd3DMOTCenter(this) + ret = this.DistanceBetweenPushBeamAnd3DMOTCenter; + end + function set.PushBeamSaturationIntensity(this, val) + this.PushBeamSaturationIntensity = val; + end + function ret = get.PushBeamSaturationIntensity(this) + ret = this.PushBeamSaturationIntensity; + end + + function set.ZeemanSlowerBeamPower(this,val) + this.ZeemanSlowerBeamPower = val; + end + function ret = get.ZeemanSlowerBeamPower(this) + ret = this.ZeemanSlowerBeamPower; + end + function set.ZeemanSlowerBeamDetuning(this, val) + this.ZeemanSlowerBeamDetuning = val; + end + function ret = get.ZeemanSlowerBeamDetuning(this) + ret = this.ZeemanSlowerBeamDetuning; + end + function set.ZeemanSlowerBeamRadius(this, val) + this.ZeemanSlowerBeamRadius = val; + end + function ret = get.ZeemanSlowerBeamRadius(this) + ret = this.ZeemanSlowerBeamRadius; + end + function set.ZeemanSlowerBeamWaist(this, val) + this.ZeemanSlowerBeamWaist = val; + end + function ret = get.ZeemanSlowerBeamWaist(this) + ret = this.ZeemanSlowerBeamWaist; + end + function set.ZeemanSlowerBeamSaturationIntensity(this, val) + this.ZeemanSlowerBeamSaturationIntensity = val; + end + function ret = get.ZeemanSlowerBeamSaturationIntensity(this) + ret = this.ZeemanSlowerBeamSaturationIntensity; + end + + function set.TotalPower(this,val) + this.TotalPower = val; + end + function ret = get.TotalPower(this) + ret = this.TotalPower; + end + function set.LandegFactor(this,val) + this.LandegFactor = val; + end + function ret = get.LandegFactor(this) + ret = this.LandegFactor; + end + function set.MagneticSubLevel(this,val) + this.MagneticSubLevel = val; + end + function ret = get.MagneticSubLevel(this) + ret = this.MagneticSubLevel; + end + + function set.CaptureVelocity(this,val) + this.CaptureVelocity = val; + end + function ret = get.CaptureVelocity(this) + ret = this.CaptureVelocity; + end + function set.VelocityCutoff(this,val) + this.VelocityCutoff = val; + end + function ret = get.VelocityCutoff(this) + ret = this.VelocityCutoff; + end + function set.CaptureFraction(this,val) + this.CaptureFraction = val; + end + function ret = get.CaptureFraction(this) + ret = this.CaptureFraction; + end + function set.ClausingFactor(this,val) + this.ClausingFactor = val; + end + function ret = get.ClausingFactor(this) + ret = this.ClausingFactor; + end + function set.AngularDistribution(this,val) + this.AngularDistribution = val; + end + function ret = get.AngularDistribution(this) + ret = this.AngularDistribution; + end + function set.NormalizationConstantForAngularDistribution(this,val) + this.NormalizationConstantForAngularDistribution = val; + end + function ret = get.NormalizationConstantForAngularDistribution(this) + ret = this.NormalizationConstantForAngularDistribution; + end + + function set.DebugMode(this, val) + this.DebugMode = val; + end + function ret = get.DebugMode(this) + ret = this.DebugMode; + end + function set.DoSave(this, val) + this.DebugMode = val; + end + function ret = get.DoSave(this) + ret = this.DoSave; + end + end % - setters and getters + + methods + function ret = get.CoolingBeamSaturationParameter(this) + ret = 4* this.CoolingBeamPower/(pi*this.CoolingBeamWaist^2)/this.CoolingBeamSaturationIntensity/10; % two beams are reflected + end + + function ret = get.SidebandSaturationParameter(this) + ret = 4*this.SidebandPower/(pi*this.SidebandBeamWaist^2)/this.SidebandBeamSaturationIntensity/10; + end + + function ret = get.PushBeamSaturationParameter(this) + ret = this.PushBeamPower/(pi*this.PushBeamWaist^2)/this.PushBeamSaturationIntensity/10; + end + + function ret = get.ZeemanSlowerBeamSaturationParameter(this) + ret = this.ZeemanSlowerBeamPower/(pi*this.ZeemanSlowerBeamWaist^2)/this.ZeemanSlowerBeamSaturationIntensity/10; + end + + function ret = get.OvenTemperatureinKelvin(this) + ret = this.OvenTemperature + Helper.PhysicsConstants.ZeroKelvin; + end + + function ret = get.AverageVelocity(this) + %See Background collision probability section in Barbiero + ret = sqrt(((8*pi)/9) * ((Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)/Helper.PhysicsConstants.Dy164Mass)); + end + + function ret = get.AtomicBeamDensity(this) + %See Background collision probability section in Barbiero + ret = this.calculateFreeMolecularRegimeFlux() / (this.AverageVelocity * pi * (this.Beta*this.NozzleLength/2)^2); + end + + function ret = get.MeanFreePath(this) + % Cross section = pi ( 2 * Van-der-waals radius of Dy)^2; + % Van-der-waals radius of Dy = 281e-12 + %See Expected atomic flux section and Background collision probability section in Barbiero + ret = 1/(sqrt(2) * ( pi * (2*281e-12)^2) * this.AtomicBeamDensity); + end + + function ret = get.CollisionTime(this) + ret = 3 * this.MeanFreePath/this.AverageVelocity; %See Background collision probability section in Barbiero + end + + end % - getters for dependent properties + + methods(Access = protected) + function cp = copyElement(this) + % Shallow copy object + cp = copyElement@matlab.mixin.Copyable(this); + + % Forces the setter to redefine the function handles to the new copied object + % cp.potentialType = this.potentialType; + + pl = properties(this); + for k = 1:length(pl) + sc = superclasses(this.(pl{k})); + if any(contains(sc,{'matlab.mixin.Copyable'})) + cp.(pl{k}) = this.(pl{k}).copy(); + end + end + end + end + +end diff --git a/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m b/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m new file mode 100644 index 0000000..75286c6 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToPushBeam.m @@ -0,0 +1,26 @@ +function ret = accelerationDueToPushBeam(this, PositionVector, VelocityVector) + % is the distance between the chamber center and the cross point of push beam and z-axis (along the gravity) + WaveVectorEndPoint = [0, 1, this.DistanceBetweenPushBeamAnd3DMOTCenter/this.PushBeamDistance]; + WaveVectorEndPoint = WaveVectorEndPoint./norm(WaveVectorEndPoint); + Origin=[0,0,0]; + + SaturationIntensity = this.calculateLocalSaturationIntensity(this.PushBeamSaturationIntensity, PositionVector, Origin, WaveVectorEndPoint, this.PushBeamRadius, this.PushBeamWaist); + + DopplerShift = (VelocityVector * WaveVectorEndPoint(1:3)') * this.PushBeamWaveVector; + + Detuning = this.PushBeamDetuning - DopplerShift; + + a_push = (Helper.PhysicsConstants.PlanckConstantReduced * this.PushBeamWaveVector * WaveVectorEndPoint(1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.PushBeamLinewidth * 0.5) .* ... + (SaturationIntensity/(1 + SaturationIntensity + 4 * (Detuning./this.PushBeamLinewidth).^2)); + + if this.SpontaneousEmission + a_scatter = this.accelerationDueToSpontaneousEmissionProcess(SaturationIntensity, SaturationIntensity, Detuning, this.PushBeamLinewidth, this.PushBeamWaveVector); + else + a_scatter = [0,0,0]; + end + + a_total = a_push + a_scatter; + + ret = a_total(1:3); +end + diff --git a/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m b/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m new file mode 100644 index 0000000..d1ffd2b --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/accelerationDueToSpontaneousEmissionProcess.m @@ -0,0 +1,15 @@ +function ret = accelerationDueToSpontaneousEmissionProcess(this, SaturationIntensity, TotalSaturationIntensity, Detuning, Linewidth, WaveVector) + Vector = [2*rand(1)-1,2*rand(1)-1,2*rand(1)-1]; + Vector = Vector./norm(Vector); + + ScatteringRate = 0.5 * SaturationIntensity * Linewidth / (1 + TotalSaturationIntensity + 4 * (Detuning/Linewidth)^2); + NumberOfScatteringEvents = floor(ScatteringRate * this.TimeStep); + + if NumberOfScatteringEvents > 0 + ret = Vector.*((Helper.PhysicsConstants.PlanckConstantReduced * WaveVector) / ... + (Helper.PhysicsConstants.Dy164Mass * this.TimeStep)).* sqrt(NumberOfScatteringEvents); + else + ret = zeros(1,3); + end +end + diff --git a/MOT Capture Process Simulation/@MOTSimulator/angularDistributionFunction.m b/MOT Capture Process Simulation/@MOTSimulator/angularDistributionFunction.m new file mode 100644 index 0000000..b9c91ec --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/angularDistributionFunction.m @@ -0,0 +1,37 @@ +function ret = angularDistributionFunction(this, theta) + %This function calculate the angle distribution of atoms coming out + %from a single channel. + + KnudsenNumber = this.MeanFreePath/this.NozzleLength; + + alpha=0.5 - (3*this.Beta^2)^-1 * ... + ((1 - (2*this.Beta^3) + ((2*this.Beta^2) - 1) * sqrt(1+this.Beta^2)) / ... + (sqrt(1+this.Beta^2) - (this.Beta^2 * asinh((this.Beta^2)^-1)))); + + eta0 = alpha; + eta1 = 1 - alpha; + + delta = (eta0./sqrt(2*KnudsenNumber*(eta1-eta0)))./sqrt(cos(theta)); + + F = 2/sqrt(pi)* (1-eta1)/eta0 * delta.* exp( -(delta*eta1/eta0).^2 ); + q = this.Beta^-1 * tan(theta); + R = acos(q) - (q .* sqrt(1-q.^2)); + + if abs(q) >= 1 + t = linspace(0,1,10000); + S = sum(sqrt(1-t.^2).* ( erf(delta.*(1 + (t.*(eta1-eta0)./(q.*eta0)) ))-erf(delta)))*(t(2)-t(1)); + if S == 0 || isnan(S) + ret = eta0*cos(theta); + else + ret = eta0*cos(theta)+ 2/sqrt(pi)*eta0*cos(theta) * (exp(delta.^2)/delta) * S; + end + else + t = linspace(0,q,10000); + S = sum(sqrt(1-t.^2).* ( erf(delta.*(1 + (t.*(eta1-eta0)./(q.*eta0)) ))-erf(delta)))*(t(2)-t(1)); + if isnan(S) + S=0; + end + ret = 2/sqrt(pi)*eta0*cos(theta)*(exp(delta.^2)/delta) * (R./2*(erf(delta*eta1/eta0)-erf(delta)+F)+S)+eta0*cos(theta); + end +end + diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m b/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m new file mode 100644 index 0000000..837388e --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateCaptureVelocity.m @@ -0,0 +1,24 @@ +function ret = calculateCaptureVelocity(this, PositionVector, VelocityVector) + + VelocityVector = VelocityVector./norm(VelocityVector); + UpperLimit = 500; + LowerLimit = 0; + + for Index = 1:500 + InitialVelocity = 0.5 * (UpperLimit + LowerLimit) * VelocityVector; + [~, FinalDynamicalQuantities] = this.solver(PositionVector, InitialVelocity); + FinalPositionVector = FinalDynamicalQuantities(1:3); + if rssq(FinalPositionVector) <= this.OvenDistance + LowerLimit = 0.5 * (UpperLimit + LowerLimit); + else + UpperLimit = 0.5 * (UpperLimit + LowerLimit); + end + + if UpperLimit - LowerLimit < 1 + ret = InitialVelocity; + break; + end + end + clear Index +end + diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m b/MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m new file mode 100644 index 0000000..7a32584 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateFreeMolecularRegimeFlux.m @@ -0,0 +1,7 @@ +function ret = calculateFreeMolecularRegimeFlux(this) + %This function calculate the total flux of atoms coming out from a tube + %See Expected atomic flux section in Barbiero + Dy164VapourPressure = 133.322*exp(11.4103-2.3785e+04./(-219.4821+this.OvenTemperatureinKelvin)).*100; % Vapor Pressure Dysprosium for the given oven temperature + Dy164DensityinOven = Dy164VapourPressure/(Helper.PhysicsConstants.BoltzmannConstant*this.OvenTemperatureinKelvin); + ret =1/4 * Dy164DensityinOven * this.AverageVelocity * pi * this.NozzleRadius.^2; +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m b/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m new file mode 100644 index 0000000..3c2d225 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateLoadingRate.m @@ -0,0 +1,67 @@ +function [LoadingRate, StandardError] = calculateLoadingRate(this, CaptureFraction, ClausingFactor, FinalDynamicalQuantities) + + NumberOfLoadedAtoms = zeros(1, this.NumberOfAtoms); + AutocorrelationFunction = zeros(1, this.NumberOfAtoms); + + for i = 1:NumberOfLoadedAtoms + FinalPosition = FinalDynamicalQuantities(i,1:3); + DivergenceAngle = atan(sqrt((FinalPosition(1)^2+FinalPosition(3)^2)/(FinalPosition(2)^2))); + if (DivergenceAngle <= this.MOTExitDivergence) && (FinalPosition(2) >= 0) + if i == 1 + NumberOfLoadedAtoms(1) = 1; + else + NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1) + 1; + end + + else + if i == 1 + NumberOfLoadedAtoms(1) = 0; + else + NumberOfLoadedAtoms(i) = NumberOfLoadedAtoms(i-1); + end + end + end + + for i = 1:NumberOfLoadedAtoms-1 + MeanLoadingRate = 0; + MeanLoadingRateShifted = 0; + for j = 1:NumberOfLoadedAtoms-i + MeanLoadingRate = MeanLoadingRate + NumberOfLoadedAtoms(j)/j; + MeanLoadingRateShifted = MeanLoadingRateShifted + (NumberOfLoadedAtoms(i+j))/(i+j); + AutocorrelationFunction(i) = AutocorrelationFunction(i) + ((NumberOfLoadedAtoms(j)/j).*(NumberOfLoadedAtoms(i+j)/(i+j))); + end + AutocorrelationFunction(i) = ((NumberOfLoadedAtoms-i)^-1 * (AutocorrelationFunction(i)) - ((NumberOfLoadedAtoms-i)^-1 * MeanLoadingRate * MeanLoadingRateShifted)); + end + + if AutocorrelationFunction(1)~=0 % In case no atom loading + + AutocorrelationFunction = AutocorrelationFunction./AutocorrelationFunction(1); + x = linspace(1, NumberOfLoadedAtoms, NumberOfLoadedAtoms); + [FitObject,~] = fit(x',AutocorrelationFunction',"exp(-x/n)",'Startpoint', 100); + n = FitObject.n; % n is the autocorrelation factor + MeanLoadingRate = 0; + NumberOfBins = floor(NumberOfLoadedAtoms/(2*n+1)); + LoadingRateError = zeros(1,NumberOfBins); + BinNumberLimit = min(NumberOfBins-1,5); + for i = 1:NumberOfBins-BinNumberLimit + LoadingRateError(i) = NumberOfLoadedAtoms(NumberOfLoadedAtoms-ceil((i-1)*(2*n+1))) / ... + (NumberOfLoadedAtoms-ceil((i-1)*(2*n+1))); + MeanLoadingRate = MeanLoadingRate + LoadingRateError(i); + end + + MeanLoadingRate = MeanLoadingRate /(NumberOfBins-BinNumberLimit); + + StandardError = 0; + for i = 1:NumberOfBins-BinNumberLimit + StandardError = StandardError + (MeanLoadingRate-LoadingRateError(i))^2; + end + StandardError = sqrt(StandardError/(NumberOfBins-BinNumberLimit)); + + LoadingRate = (MeanLoadingRate * this.FreeMolecularRegimeFlux() * CaptureFraction * ClausingFactor); + + else + LoadingRate = 0; + StandardError = 0; + end + +end diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m b/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m new file mode 100644 index 0000000..5472d29 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateLocalSaturationIntensity.m @@ -0,0 +1,18 @@ +function ret = calculateLocalSaturationIntensity(~, PeakIntensity, PositionVector, WaveVectorOrigin, WaveVectorEndPoint, BeamRadius, BeamWaist) + WaveVector = WaveVectorEndPoint - WaveVectorOrigin; % Line + PositionVectorFromWaveVectorOrigin = PositionVector - WaveVectorOrigin; % Point = PositionVector + + %Height of parallelogram (Distance between point and line) = Area of parallelogram / Base + %One side of parallelogram = PositionVectorFromWaveVectorOrigin + %Base = Wavevector + %Area = One side of parallelogram X Base + %DistanceBetweenAtomAndLaserBeamAxis = norm(cross(PositionVectorFromWaveVectorOrigin, WaveVector))./ norm(WaveVector); % Slow + DistanceBetweenAtomAndLaserBeamAxis = norm((WaveVector*WaveVector')*PositionVectorFromWaveVectorOrigin-(WaveVector*PositionVectorFromWaveVectorOrigin')*WaveVector)./ ... + (WaveVector(1)^2+WaveVector(2)^2+WaveVector(3)^2); % Faster + + if DistanceBetweenAtomAndLaserBeamAxis <= BeamRadius + ret = PeakIntensity * exp(-2*DistanceBetweenAtomAndLaserBeamAxis^2 / BeamWaist^2); + else + ret = 0; + end +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m b/MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m new file mode 100644 index 0000000..b7da57c --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/calculateTotalAcceleration.m @@ -0,0 +1,95 @@ +function ret = calculateTotalAcceleration(this, PositionVector, VelocityVector) + + WaveVectorEndPoint = zeros(2,3); + WaveVectorEndPoint(1,:) = [1,0,1]; + WaveVectorEndPoint(1,:) = WaveVectorEndPoint(1,1:3)/norm(WaveVectorEndPoint(1,:)); + WaveVectorEndPoint(2,:) = [-1,0,1]; + WaveVectorEndPoint(2,:) = WaveVectorEndPoint(2,1:3)/norm(WaveVectorEndPoint(2,:)); + + Sigma = [1,-1]; + Origin = [0,0,0]; + + % Calculate the Saturation Intensity at the specified point along its Gaussian Profile + CoolingBeamLocalSaturationIntensity = [this.calculateLocalSaturationIntensity(1, PositionVector, Origin, WaveVectorEndPoint(1,:), this.CoolingBeamRadius, this.CoolingBeamWaist), ... + this.calculateLocalSaturationIntensity(1, PositionVector, Origin, WaveVectorEndPoint(2,:), this.CoolingBeamRadius, this.CoolingBeamWaist)]; + + SidebandLocalSaturationIntensity = [this.calculateLocalSaturationIntensity(1, PositionVector, Origin, WaveVectorEndPoint(1,:), this.SidebandBeamRadius, this.SidebandBeamWaist), ... + this.calculateLocalSaturationIntensity(1, PositionVector, Origin, WaveVectorEndPoint(2,:), this.SidebandBeamRadius, this.SidebandBeamWaist)]; + + + TotalAcceleration = zeros(1,3); + + Delta_Cooling = [0,0,0,0]; + Delta_Sideband = [0,0,0,0]; + + for i = 1:2 + + LocalMagneticField = this.magneticFieldForMOT(PositionVector); + + B = sign(LocalMagneticField(1:3) * WaveVectorEndPoint(i,1:3)') * LocalMagneticField(4); + + ZeemanShift = this.LandegFactor * this.MagneticSubLevel * Helper.PhysicsConstants.BohrMagneton * Helper.PhysicsConstants.PlanckConstantReduced * B; + + DopplerShift = (VelocityVector * WaveVectorEndPoint(i,1:3)') * this.CoolingBeamWaveVector; + + Delta_Cooling(i*2-1) = this.CoolingBeamDetuning + DopplerShift + ZeemanShift * Sigma(i); + Delta_Cooling(i*2) = this.CoolingBeamDetuning - DopplerShift - ZeemanShift * Sigma(i); + + if this.Sideband + Delta_Sideband(i*2-1) = this.SidebandDetuning + DopplerShift + ZeemanShift * Sigma(i); + Delta_Sideband(i*2) = this.SidebandDetuning - DopplerShift - ZeemanShift * Sigma(i); + end + end + + SaturationParameter = [0,0,0,0,0,0,0,0]; + + for i = 1:2 + SaturationParameter(2*i-1) = (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4*Delta_Cooling(2*i-1)^2 / this.CoolingBeamLinewidth^2); + SaturationParameter(2*i) = (0.25 * this.CoolingBeamSaturationParameter * CoolingBeamLocalSaturationIntensity(1)) /(1 + 4*Delta_Cooling(2*i)^2 / this.CoolingBeamLinewidth^2); + if this.Sideband + SaturationParameter(2*i-1+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(1)) /(1 + 4*Delta_Sideband(2*i-1)^2/ this.CoolingBeamLinewidth^2); + SaturationParameter(2*i+4) = (0.25 * this.SidebandSaturationParameter * SidebandLocalSaturationIntensity(2)) /(1 + 4*Delta_Sideband(2*i)^2 / this.CoolingBeamLinewidth^2); + end + end + + TotalSaturationParameter = sum(SaturationParameter); + + for i = 1:2 + + a_1 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ... + (SaturationParameter(2*i-1)/(1 + TotalSaturationParameter)); + a_2 = (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ... + (SaturationParameter(2*i) / (1 + TotalSaturationParameter)); + + if this.SpontaneousEmission + a_scattering = this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveVector) + ... + this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i), TotalSaturationParameter, Delta_Cooling(2*i) , this.CoolingBeamLinewidth, this.CoolingBeamWaveVector); + else + a_scattering = [0,0,0]; + end + + if this.Sideband + a_1 = a_1 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ... + (SaturationParameter(2*i-1+4)/(1 + TotalSaturationParameter)); + a_2 = a_2 + (Helper.PhysicsConstants.PlanckConstantReduced * this.CoolingBeamWaveVector * WaveVectorEndPoint(i,1:3)/Helper.PhysicsConstants.Dy164Mass).*(this.CoolingBeamLinewidth * 0.5) .* ... + (SaturationParameter(2*i+4)/(1 + TotalSaturationParameter)); + + if this.SpontaneousEmission + a_scattering = a_scattering + ... + this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i-1+4), TotalSaturationParameter, Delta_Cooling(2*i-1), this.CoolingBeamLinewidth, this.CoolingBeamWaveVector) + ... + this.accelerationDueToSpontaneousEmissionProcess(SaturationParameter(2*i+4), TotalSaturationParameter, Delta_Cooling(2*i) , this.CoolingBeamLinewidth, this.CoolingBeamWaveVector); + else + a_scattering = [0,0,0]; + end + end + + TotalAcceleration = TotalAcceleration + (a_2 - a_1) + a_scattering; + end + + if this.PushBeamRadius ~= 0 + TotalAcceleration = TotalAcceleration + this.accelerationDueToPushBeam(PositionVector, VelocityVector); + end + + ret = TotalAcceleration(1:3); + +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/initialPositionSampling.m b/MOT Capture Process Simulation/@MOTSimulator/initialPositionSampling.m new file mode 100644 index 0000000..85af394 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/initialPositionSampling.m @@ -0,0 +1,7 @@ +function ret = initialPositionSampling(this) + n = this.NumberOfAtoms; + phi = 2 * pi * rand(n,1); + rho = this.Beta * 0.5 * this.NozzleLength * sqrt(rand(n,1)); + ret = [-this.OvenDistance * ones(n,1), rho.*cos(phi), rho.*sin(phi)]; +end + diff --git a/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m b/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m new file mode 100644 index 0000000..129f17f --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/initialVelocitySampling.m @@ -0,0 +1,47 @@ +function ret = initialVelocitySampling(this, VelocityCutoff, AngularDistribution, NormalizationConstant) + n = this.NumberOfAtoms; + + ret = zeros(n,3); + SampledVelocityMagnitude =zeros(n,1); + SampledPolarAngle =zeros(n,1); + SampledAzimuthalAngle =zeros(n,1); + + MostProbableVelocity = sqrt((3 * Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperature) / Helper.PhysicsConstants.Dy164Mass); % For v * f(v) distribution + + if MostProbableVelocity > VelocityCutoff + MaximumVelocityAllowed = VelocityCutoff; + else + MaximumVelocityAllowed = MostProbableVelocity; + end + ProbabilityOfMaximumVelocityAllowed = this.velocityDistributionFunction(MaximumVelocityAllowed); + ProbabilityOfMaximumDivergenceAngleAllowed = NormalizationConstant * max(AngularDistribution); + + parfor i = 1:n + % Rejection Sampling of speed + y = ProbabilityOfMaximumVelocityAllowed * rand(1); + x = VelocityCutoff * rand(1); + + while y > this.velocityDistributionFunction(x) %As long as this loop condition is satisfied, reject the corresponding x value + y = ProbabilityOfMaximumVelocityAllowed * rand(1); + x = VelocityCutoff * rand(1); + end + SampledVelocityMagnitude(i) = x; % When loop condition is not satisfied, accept x value and store as sample + + % Rejection Sampling of polar angle + w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1); + z = this.MOTExitDivergence * rand(1); + + while w > (NormalizationConstant * this.angularDistributionFunction(z) * sin(z)) %As long as this loop condition is satisfied, reject the corresponding x value + w = ProbabilityOfMaximumDivergenceAngleAllowed * rand(1); + z = this.MOTExitDivergence * rand(1); + end + SampledPolarAngle(i) = z; %When loop condition is not satisfied, accept x value and store as sample + + % Sampling of azimuthal angle + SampledAzimuthalAngle(i)= 2 * pi * rand(1); + + ret(i,:)=[SampledVelocityMagnitude(i)*cos(SampledPolarAngle(i)), SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*cos(SampledAzimuthalAngle(i)), ... + SampledVelocityMagnitude(i)*sin(SampledPolarAngle(i))*sin(SampledAzimuthalAngle(i))]; + end +end + diff --git a/MOT Capture Process Simulation/@MOTSimulator/magneticFieldForMOT.m b/MOT Capture Process Simulation/@MOTSimulator/magneticFieldForMOT.m new file mode 100644 index 0000000..802c27b --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/magneticFieldForMOT.m @@ -0,0 +1,14 @@ +function ret = magneticFieldForMOT(this, r) + switch this.SimulationMode + case '2D' + ret = zeros(1,4); + alpha = this.MagneticGradient; + ret(1) = r(3)*alpha; + ret(2) = 0; + ret(3) = r(1)*alpha; + ret(4) = sqrt(ret(1)^2+ret(2)^2+ret(3)^2); + case '3D' + % Development in progress + + end +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m b/MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m new file mode 100644 index 0000000..82bdf49 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/reinitializeSimulator.m @@ -0,0 +1,32 @@ +function reinitializeSimulator(this) + %% PHYSICAL CONSTANTS + pc = Helper.PhysicsConstants; + %% SIMULATION PARAMETERS + this.NozzleLength = 60e-3; + this.NozzleRadius = 2.50e-3; + this.Beta = 2 * (this.NozzleRadius/this.NozzleLength); + this.ApertureCut = max(2.5e-3,this.NozzleRadius); + this.OvenDistance = ((25+12.5)*1e-3 + (this.NozzleRadius + this.ApertureCut)) / tan(15/360 * 2 * pi); + % Distance between the nozzle and the 2-D MOT chamber center + % 25 is the beam radius/sqrt(2) + % 12.5 is the radius of the oven + % 15 eg is the angle between the 2-D MOT chamber center and the nozzle + this.OvenTemperature = 1000; % Temperature in Celsius + this.MOTDistance = 320e-3; % Distance between the 2-D MOT the 3-D MOT + this.MagneticGradient = 0.425; % T/m + this.BlueWaveVector = 2*pi/pc.BlueWavelength; + this.BlueSaturationIntensity = 2*pi^2*pc.PlanckConstantReduced*pc.SpeedOfLight*pc.BlueLinewidth/3/(pc.BlueWavelength)^3/10; + this.OrangeWaveVector = 2*pi/pc.OrangeWavelength; + this.OrangeSaturationIntensity = 2*pi^2*pc.PlanckConstantReduced*pc.SpeedOfLight*pc.OrangeLinewidth/3/(pc.OrangeWavelength)^3/10; + this.BlueBeamRadius = min(0.035/2,sqrt(2)/2*this.OvenDistance); % Diameter of CF40 flange = 0.035 + Theta_Nozzle = atan((this.NozzleRadius+this.BlueBeamRadius*sqrt(2))/this.OvenDistance); % The angle of capture region towards the oven nozzle + Theta_Aperture = 15/360*2*pi; % The limitation angle of the second aperture in the oven + this.NozzleExitDivergence = min(Theta_Nozzle,Theta_Aperture); + this.MOTExitDivergence = 0.016; % The limitation angle between 2D-MOT and 3D-MOT + this.TotalPower = 0.4; + this.OrangeBeamRadius = 1.2e-03; + this.PushBeamRadius = 0; + this.PushBeamDistance = 0; + this.DistanceBetweenPushBeamAnd3DMOTCenter = 1; + this.ZeemanSlowerBeamRadius = 1; +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m b/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m new file mode 100644 index 0000000..6c4e0f5 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/runSimulation.m @@ -0,0 +1,66 @@ +function [LoadingRate, StandardError] = runSimulation(this) + n = this.NumberOfAtoms; + %% Calculate Background Collision Time --> Calculate Capture velocity --> Introduce velocity cutoff --> Calculate capture fraction + + this.CaptureVelocity = 1.05 * this.calculateCaptureVelocity([-this.OvenDistance,0,0], [1,0,0]); % Make 5% larger than the numerically obtained CV + this.VelocityCutoff = this.CaptureVelocity(1); %Should be the magnitude of the 3-D velocity vector but since here the obtained capture + %velocity is only along the x-axis, we take the first term which is the x-component of the velocity. + + VelocityDistributionFunction = @(velocity) sqrt(2 / pi) * (Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)) ... + * velocity.^3 .* exp(-velocity.^2 .* (Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ... + * this.OvenTemperatureinKelvin))); + + this.CaptureFraction = integral(VelocityDistributionFunction, 0, this.VelocityCutoff) / integral(VelocityDistributionFunction, 0, Inf); + + %% Calculate the Clausing Factor + + % Compute the angular distribution of the atomic beam + ThetaArray = linspace(0.0001, pi/2, 1000); + AngularDistribution = zeros(1,length(ThetaArray)); + parfor i = 1:length(ThetaArray) + AngularDistribution(i) = this.angularDistributionFunction(ThetaArray(i)); + end + + % Numerically integrate the angular distribution over the full solid angle + NormalizationConstant = 0; + for j = 1:length(ThetaArray) + if ThetaArray(j) <= this.NozzleExitDivergence + NormalizationConstant = NormalizationConstant + (2 * pi * sin(ThetaArray(j)) * AngularDistribution(j) * (ThetaArray(2)-ThetaArray(1))); + end + end + + this.ClausingFactor = (1/pi) * NormalizationConstant; %The complete intergration will give pi * ClausingFactor. + %Therefore, the Clausing Factor needs to be extracted from the + %result of the above integration by dividing out pi + + this.AngularDistribution = AngularDistribution; + this.NormalizationConstantForAngularDistribution = NormalizationConstant; + + %% + % - sampling the position distribution + this.InitialPositions = this.initialPositionSampling(); + % - sampling the velocity distribution + this.InitialVelocities = this.initialVelocitySampling(this.VelocityCutoff, this.AngularDistribution, this.NormalizationConstantForAngularDistribution); + + %% Solve ODE + progressbar = Helper.parforNotifications(); + progressbar.PB_start(n,'Message',['Simulating capture process for ' num2str(n,'%.0f') ' atoms:']); + + % calculate the final position of the atoms + FinalDynamicalQuantities = zeros(n,9); + Positions = this.InitialPositions; + Velocities = this.InitialVelocities; + parfor Index = 1:n + ret = this.solver(Positions(Index,:), Velocities(Index,:)); + FinalDynamicalQuantities(Index,:) = ret(2); + progressbar.PB_iterate(); + end + clear Index + %% Calculate the Loading Rate + [LoadingRate, StandardError] = this.calculateLoadingRate(this.CaptureFraction, this.ClausingFactor, FinalDynamicalQuantities); + %% Save + %if this.DoSave + + %end + +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m b/MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m new file mode 100644 index 0000000..64ffe39 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/setInitialConditions.m @@ -0,0 +1,98 @@ +function setInitialConditions(this, varargin) + + p = inputParser; + p.KeepUnmatched = true; + addParameter(p, 'NumberOfAtoms', 5000,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'BluePower', 200e-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'BlueDetuning', -1.92857*Helper.PhysicsConstants.BlueLinewidth,... + @(x) assert(isnumeric(x) && isscalar(x))); + addParameter(p, 'BlueBeamWaist', 10e-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'SidebandPower', 200e-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'SidebandDetuning', 0,... + @(x) assert(isnumeric(x) && isscalar(x))); + addParameter(p, 'SidebandBeamWaist', 12e-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'PushBeamPower', 10e-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'PushBeamDetuning', 0,... + @(x) assert(isnumeric(x) && isscalar(x))); + addParameter(p, 'PushBeamWaist', 0.81e-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'OrangePower', 70e-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'OrangeDetuning', 12e-3,... + @(x) assert(isnumeric(x) && isscalar(x))); + addParameter(p, 'OrangeBeamWaist', 12e-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'ZeemanSlowerBeamPower', 200e-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + addParameter(p, 'ZeemanSlowerBeamDetuning', -7*Helper.PhysicsConstants.BlueLinewidth,... + @(x) assert(isnumeric(x) && isscalar(x))); + addParameter(p, 'ZeemanSlowerBeamWaist', 7e-3,... + @(x) assert(isnumeric(x) && isscalar(x) && (x > 0))); + + p.parse(varargin{:}); + + this.NumberOfAtoms = p.Results.NumberOfAtoms; + this.BluePower = p.Results.BluePower; + this.BlueDetuning = p.Results.BlueDetuning; + this.BlueBeamWaist = p.Results.BlueBeamWaist; + this.SidebandPower = p.Results.SidebandPower; + this.SidebandDetuning = p.Results.SidebandDetuning; + this.SidebandBeamWaist = p.Results.SidebandBeamWaist; + this.PushBeamPower = p.Results.PushBeamPower; + this.PushBeamDetuning = p.Results.PushBeamDetuning; + this.PushBeamWaist = p.Results.PushBeamWaist; + this.OrangePower = p.Results.OrangePower; + this.OrangeDetuning = p.Results.OrangeDetuning; + this.OrangeBeamWaist = p.Results.OrangeBeamWaist; + this.ZeemanSlowerBeamPower = p.Results.ZeemanSlowerBeamPower; + this.ZeemanSlowerBeamDetuning = p.Results.ZeemanSlowerBeamDetuning; + this.ZeemanSlowerBeamWaist = p.Results.ZeemanSlowerBeamWaist; + + %% Set general parameters according to simulation mode + switch this.SimulationMode + case "2D" + this.CoolingBeamPower = this.BluePower; + this.CoolingBeamWaist = this.BlueBeamWaist; + this.CoolingBeamLinewidth = Helper.PhysicsConstants.BlueLinewidth; + this.CoolingBeamWaveVector = this.BlueWaveVector; + this.CoolingBeamDetuning = this.BlueDetuning; + this.CoolingBeamRadius = this.BlueBeamRadius; + this.CoolingBeamSaturationIntensity = this.BlueSaturationIntensity; + this.SidebandBeamRadius = this.BlueBeamRadius; + this.SidebandBeamSaturationIntensity = this.BlueSaturationIntensity; + this.PushBeamLinewidth = Helper.PhysicsConstants.OrangeLinewidth; + this.PushBeamWaveVector = this.OrangeWaveVector; + this.PushBeamDetuning = this.OrangeDetuning; + this.PushBeamSaturationIntensity = this.OrangeSaturationIntensity; + this.LandegFactor = Helper.PhysicsConstants.BlueLandegFactor; + this.MagneticSubLevel = 1; + case "3D" + % Development In progress + end + + + %% - store in struct + this.InitialParameters = struct; + this.InitialParameters.NumberOfAtoms = this.NumberOfAtoms; + this.InitialParameters.BluePower = this.BluePower; + this.InitialParameters.BlueDetuning = this.BlueDetuning; + this.InitialParameters.BlueBeamWaist = this.BlueBeamWaist; + this.InitialParameters.SidebandPower = this.SidebandPower; + this.InitialParameters.SidebandDetuning = this.SidebandDetuning; + this.InitialParameters.SidebandBeamWaist = this.SidebandBeamWaist; + this.InitialParameters.PushBeamPower = this.PushBeamPower; + this.InitialParameters.PushBeamDetuning = this.PushBeamDetuning; + this.InitialParameters.PushBeamWaist = this.PushBeamWaist; + this.InitialParameters.OrangePower = this.OrangePower; + this.InitialParameters.OrangeDetuning = this.OrangeDetuning; + this.InitialParameters.OrangeBeamWaist = this.OrangeBeamWaist; + this.InitialParameters.ZeemanSlowerBeamPower = this.ZeemanSlowerBeamPower; + this.InitialParameters.ZeemanSlowerBeamDetuning = this.ZeemanSlowerBeamDetuning; + this.InitialParameters.ZeemanSlowerBeamBeamWaist = this.ZeemanSlowerBeamWaist; +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/@MOTSimulator/solver.m b/MOT Capture Process Simulation/@MOTSimulator/solver.m new file mode 100644 index 0000000..b0c9dfd --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/solver.m @@ -0,0 +1,54 @@ +function [ParticleTrajectory, FinalDynamicalQuantities] = solver(this, InitialPosition, InitialVelocity) + if this.Gravity + g = [0,0,- -Helper.PhysicsConstants.GravitationalAcceleration]; + else + g = 0; + end + + % Probability of Background Collisions + collision = rand(1); + CollisionProbability = 1 - exp(-this.SimulationTime/this.CollisionTime); + + if collision >= CollisionProbability || this.CollisionTime == -500 % -500 is a flag for skipping the background collision + + ParticleTrajectory = zeros(int64(this.SimulationTime/this.TimeStep),9); + + for i=1:int64(this.SimulationTime/this.TimeStep) + + ParticleTrajectory(i,1:3) = InitialPosition; + ParticleTrajectory(i,4:6) = InitialVelocity; + + rt = InitialPosition; + vt = InitialVelocity; + + ga1 = this.calculateTotalAcceleration(rt,vt) + g; + gv1 = vt .* this.TimeStep; + rt = rt + 0.5 * gv1; + vt = vt + 0.5 * ga1 .* 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 + end +end + diff --git a/MOT Capture Process Simulation/@MOTSimulator/velocityDistributionFunction.m b/MOT Capture Process Simulation/@MOTSimulator/velocityDistributionFunction.m new file mode 100644 index 0000000..7894b06 --- /dev/null +++ b/MOT Capture Process Simulation/@MOTSimulator/velocityDistributionFunction.m @@ -0,0 +1,5 @@ +function ret = velocityDistributionFunction(this, velocity) + ret = sqrt(2 / pi) * (Helper.PhysicsConstants.Dy164Mass/(Helper.PhysicsConstants.BoltzmannConstant * this.OvenTemperatureinKelvin)) ... + * velocity^3 * exp(-velocity^2.*(Helper.PhysicsConstants.Dy164Mass / (2 * Helper.PhysicsConstants.BoltzmannConstant ... + * this.OvenTemperatureinKelvin))); +end \ No newline at end of file diff --git a/MOT Capture Process Simulation/test_MOTSimulator.m b/MOT Capture Process Simulation/test_MOTSimulator.m new file mode 100644 index 0000000..90a28fd --- /dev/null +++ b/MOT Capture Process Simulation/test_MOTSimulator.m @@ -0,0 +1,100 @@ +%% - Create solver object with specified options +clc +%% + +OptionsStruct = struct; +OptionsStruct.SimulationMode = '2D'; +OptionsStruct.TimeStep = 50e-06; +OptionsStruct.SimulationTime = 04e-03; +OptionsStruct.SpontaneousEmission = true; +OptionsStruct.Sideband = true; +OptionsStruct.ZeemanSlowerBeam = false; +OptionsStruct.Gravity = true; +OptionsStruct.DebugMode = false; +OptionsStruct.SaveData = false; +options = Helper.convertstruct2cell(OptionsStruct); + +Simulator = MOTSimulator(options{:}); + +clear OptionsStruct +%% - Set Initial Conditions: Run with default values +Simulator.setInitialConditions(); + +%% - Set Initial Conditions: Set manually +OptionsStruct = struct; +OptionsStruct.NumberOfAtoms = 5000; +OptionsStruct.BluePower = 0.2; % in W +OptionsStruct.BlueDetuning = -2 * Helper.PhysicsConstants.BlueLinewidth; % in Hz +OptionsStruct.BlueBeamWaist = 0.010; % in m +OptionsStruct.SidebandPower = 0.2; +OptionsStruct.SidebandDetuning = -3 * Helper.PhysicsConstants.BlueLinewidth; % in Hz +OptionsStruct.SidebandBeamWaist = 0.010; % in m +OptionsStruct.PushBeamPower = 0.010; % in W +OptionsStruct.PushBeamDetuning = 0; % in Hz +OptionsStruct.PushBeamWaist = 0.005; % in m + +options = Helper.convertstruct2cell(OptionsStruct); +Simulator.setInitialConditions(options{:}); +clear OptionsStruct +%% - Run Simulation +[LoadingRate, ~] = Simulator.runSimulation(); + +%% - Plot initial distribution +NumberOfBins = 100; +Plotting.plotPositionAndVelocitySampling(Simulator, NumberOfBins); + +%% - Plot Magnetic Field +XAxisRange = [-5 5]; +YAxisRange = [-5 5]; +ZAxisRange = [-5 5]; +Plotting.visualizeMagneticField(Simulator, XAxisRange, YAxisRange, ZAxisRange) + +%% - Scan parameters + %Use a for loop to change the parameter during set initial conditions + %Run simulation + +% TWO-PARAMETER SCAN + +NumberOfPointsForFirstParam = 10; %iterations of the simulation +NumberOfPointsForSecondParam = 10; +LoadingRateArray = zeros(NumberOfPointsForFirstParam, NumberOfPointsForSecondParam); + +% Scan Sideband Detuning and Power Ratio +DetuningArray = linspace(-0.5,-10, NumberOfPointsForFirstParam); +PowerArray = linspace(0.1,0.9, NumberOfPointsForSecondParam); + +tStart = tic; +for i=1:NumberOfPointsForFirstParam + OptionsStruct.SidebandDetuning = DetuningArray(i) * Helper.PhysicsConstants.BlueLinewidth; + for j=1:NumberOfPointsForSecondParam + OptionsStruct.BluePower = PowerArray(j) * Simulator.TotalPower; + OptionsStruct.SidebandPower = Simulator.TotalPower - OptionsStruct.BluePower; + options = Helper.convertstruct2cell(OptionsStruct); + Simulator.setInitialConditions(options{:}); + tic + [LoadingRate, ~] = Simulator.runSimulation(); + LoadingRateArray(i,j) = LoadingRate; + end +end +tEnd = toc(tStart); +fprintf('Total Computational Time: %0.1f seconds. \n', tEnd); + +%% - Plot results + +FirstParameterArray = DetuningArray * Helper.PhysicsConstants.BlueLinewidth; +SecondParameterArray = (Simulator.TotalPower - (PowerArray * Simulator.TotalPower)); +QuantityOfInterestArray = LoadingRateArray; + +OptionsStruct = struct; +OptionsStruct.RescalingFactorForFirstParameter = (Helper.PhysicsConstants.BlueLinewidth)^-1; +OptionsStruct.XLabelString = 'Detuning (\Delta/\Gamma)'; +OptionsStruct.RescalingFactorForSecondParameter = 1000; +OptionsStruct.YLabelString = 'Sideband Beam Power (mW)'; +OptionsStruct.ZLabelString = 'Loading rate (atoms/s)'; +OptionsStruct.TitleString = sprintf('Magnetic Gradient = %.0f (G/cm)', Simulator.MagneticGradient * 100); + +options = Helper.convertstruct2cell(OptionsStruct); + +Plotting.plotResultForTwoParameterScan(FirstParameterArray, SecondParameterArray, QuantityOfInterestArray, options{:}) + +clear OptionsStruct \ No newline at end of file