Latest Calculations repo code.
This commit is contained in:
parent
294fe62c36
commit
461ca442de
5
.gitignore
vendored
Normal file
5
.gitignore
vendored
Normal file
@ -0,0 +1,5 @@
|
||||
2DMOT Simulation Code/Results
|
||||
AtomECS Simulation Code
|
||||
TimeSeriesAnalyzer/Time Series Data
|
||||
ULE Cavity Characterisitics/Data
|
||||
ULE Cavity Characterisitics/Figures
|
||||
BIN
GPE Solver/+Helper/ImageSelection.class
Normal file
BIN
GPE Solver/+Helper/ImageSelection.class
Normal file
Binary file not shown.
38
GPE Solver/+Helper/ImageSelection.java
Normal file
38
GPE Solver/+Helper/ImageSelection.java
Normal file
@ -0,0 +1,38 @@
|
||||
/*
|
||||
* Based on code snippet from
|
||||
* http://java.sun.com/developer/technicalArticles/releases/data/
|
||||
*
|
||||
* Copyright © 2008, 2010 Oracle and/or its affiliates. All rights reserved. Use is subject to license terms.
|
||||
*/
|
||||
|
||||
import java.awt.image.BufferedImage;
|
||||
import java.awt.datatransfer.*;
|
||||
|
||||
public class ImageSelection implements Transferable {
|
||||
|
||||
private static final DataFlavor flavors[] =
|
||||
{DataFlavor.imageFlavor};
|
||||
|
||||
private BufferedImage image;
|
||||
|
||||
public ImageSelection(BufferedImage image) {
|
||||
this.image = image;
|
||||
}
|
||||
|
||||
// Transferable
|
||||
public Object getTransferData(DataFlavor flavor) throws UnsupportedFlavorException {
|
||||
if (flavor.equals(flavors[0]) == false) {
|
||||
throw new UnsupportedFlavorException(flavor);
|
||||
}
|
||||
return image;
|
||||
}
|
||||
|
||||
public DataFlavor[] getTransferDataFlavors() {
|
||||
return flavors;
|
||||
}
|
||||
|
||||
public boolean isDataFlavorSupported(DataFlavor
|
||||
flavor) {
|
||||
return flavor.equals(flavors[0]);
|
||||
}
|
||||
}
|
||||
46
GPE Solver/+Helper/PhysicsConstants.m
Normal file
46
GPE Solver/+Helper/PhysicsConstants.m
Normal file
@ -0,0 +1,46 @@
|
||||
classdef PhysicsConstants < handle
|
||||
properties (Constant)
|
||||
% CODATA
|
||||
PlanckConstant=6.62607015E-34;
|
||||
PlanckConstantReduced=6.62607015E-34/(2*pi);
|
||||
FineStructureConstant=7.2973525698E-3;
|
||||
ElectronMass=9.10938291E-31;
|
||||
GravitationalConstant=6.67384E-11;
|
||||
ProtonMass=1.672621777E-27;
|
||||
AtomicMassUnit=1.66053878283E-27;
|
||||
BohrRadius=0.52917721092E-10;
|
||||
BohrMagneton=927.400968E-26;
|
||||
BoltzmannConstant=1.380649E-23;
|
||||
StandardGravityAcceleration=9.80665;
|
||||
SpeedOfLight=299792458;
|
||||
StefanBoltzmannConstant=5.670373E-8;
|
||||
ElectronCharge=1.602176634E-19;
|
||||
VacuumPermeability=1.25663706212E-6;
|
||||
DielectricConstant=8.8541878128E-12;
|
||||
ElectronGyromagneticFactor=-2.00231930436153;
|
||||
AvogadroConstant=6.02214076E23;
|
||||
ZeroKelvin = 273.15;
|
||||
GravitationalAcceleration = 9.80553;
|
||||
|
||||
% Dy specific constants
|
||||
Dy164Mass = 163.929174751*1.66053878283E-27;
|
||||
Dy164IsotopicAbundance = 0.2826;
|
||||
BlueWavelength = 421.291e-9;
|
||||
BlueLandegFactor = 1.22;
|
||||
BlueLifetime = 4.94e-9;
|
||||
BlueLinewidth = 1/4.94e-9;
|
||||
RedWavelength = 626.086e-9;
|
||||
RedLandegFactor = 1.29;
|
||||
RedLifetime = 1.2e-6;
|
||||
RedLinewidth = 1/1.2e-6;
|
||||
PushBeamWaveLength = 626.086e-9;
|
||||
PushBeamLifetime = 1.2e-6;
|
||||
PushBeamLinewidth = 1/1.2e-6;
|
||||
end
|
||||
|
||||
methods
|
||||
function pc = PhysicsConstants()
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
15
GPE Solver/+Helper/bringFiguresWithTagInForeground.m
Normal file
15
GPE Solver/+Helper/bringFiguresWithTagInForeground.m
Normal file
@ -0,0 +1,15 @@
|
||||
function output = bringFiguresWithTagInForeground()
|
||||
|
||||
figure_handles = findobj('type','figure');
|
||||
|
||||
for idx = 1:length(figure_handles)
|
||||
if ~isempty(figure_handles(idx).Tag)
|
||||
figure(figure_handles(idx));
|
||||
end
|
||||
end
|
||||
|
||||
if nargout > 0
|
||||
output = figure_handles;
|
||||
end
|
||||
|
||||
end
|
||||
10
GPE Solver/+Helper/calculateDistanceFromPointToLine.m
Normal file
10
GPE Solver/+Helper/calculateDistanceFromPointToLine.m
Normal file
@ -0,0 +1,10 @@
|
||||
function ret = calculateDistanceFromPointToLine(p0 , p1, p2)
|
||||
p01 = p0 - p1;
|
||||
p12 = p2 - p1;
|
||||
CrossProduct = [p01(2)*p12(3) - p01(3)*p12(2), p01(3)*p12(1) - p01(1)*p12(3), p01(1)*p12(2) - p01(2)*p12(1)];
|
||||
ret = norm(CrossProduct) / norm(p12);
|
||||
|
||||
%Height of parallelogram (Distance between point and line) = Area of parallelogram / Base
|
||||
%Area = One side of parallelogram X Base
|
||||
%ret = norm(cross(one side, base))./ norm(base);
|
||||
end
|
||||
6
GPE Solver/+Helper/convertstruct2cell.m
Normal file
6
GPE Solver/+Helper/convertstruct2cell.m
Normal file
@ -0,0 +1,6 @@
|
||||
function CellOut = convertstruct2cell(StructIn)
|
||||
% CellOut = Convertstruct2cell(StructIn)
|
||||
% converts a struct into a cell-matrix where the first column contains
|
||||
% the fieldnames and the second the contents
|
||||
CellOut = [fieldnames(StructIn) struct2cell(StructIn)]';
|
||||
end
|
||||
18
GPE Solver/+Helper/findAllZeroCrossings.m
Normal file
18
GPE Solver/+Helper/findAllZeroCrossings.m
Normal file
@ -0,0 +1,18 @@
|
||||
function ret = findAllZeroCrossings(x,y)
|
||||
% Finds all Zero-crossing of the function y = f(x)
|
||||
zci = @(v) find(v(:).*circshift(v(:), [-1 0]) <= 0); % Returns Approximate Zero-Crossing Indices Of Argument Vector
|
||||
zxidx = zci(y);
|
||||
if ~isempty(zxidx)
|
||||
for k1 = 1:numel(zxidx)
|
||||
idxrng = max([1 zxidx(k1)-1]):min([zxidx(k1)+1 numel(y)]);
|
||||
xrng = x(idxrng);
|
||||
yrng = y(idxrng);
|
||||
[yrng2, ~, jyrng] = unique(yrng); %yrng is a new array containing the unique values of yrng. jyrng contains the indices in yrng that correspond to the original vector. yrng = yrng2(jyrng)
|
||||
xrng2 = accumarray(jyrng, xrng, [], @mean); %This function creates a new array "xrng2" by applying the function "@mean" to all elements in "xrng" that have identical indices in "jyrng". Any elements with identical X values will have identical indices in jyrng. Thus, this function creates a new array by averaging values with identical X values in the original array.
|
||||
ret(k1) = interp1( yrng2(:), xrng2(:), 0, 'linear', 'extrap' );
|
||||
end
|
||||
else
|
||||
warning('No zero crossings found!')
|
||||
ret = nan;
|
||||
end
|
||||
end
|
||||
191
GPE Solver/+Helper/getFigureByTag.m
Normal file
191
GPE Solver/+Helper/getFigureByTag.m
Normal file
@ -0,0 +1,191 @@
|
||||
function figure_handle = getFigureByTag(tag_name, varargin)
|
||||
% figure_handle = getFigureByTag(tag_name, varargin)
|
||||
%
|
||||
% Example code:
|
||||
% f_h = getFigureByTag('survivalMeasurement','Name','Survival')
|
||||
%
|
||||
% clf(f_h);
|
||||
% a_h = gca(f_h);
|
||||
% xlim(a_h,[10,100]);
|
||||
% % custom position
|
||||
% f_h.Position = [4052.3 719.67 560 420];
|
||||
|
||||
assert(nargin>=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',[7,7,pos(3)-2,pos(4)]);
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
|
||||
92
GPE Solver/+Helper/ode5.m
Normal file
92
GPE Solver/+Helper/ode5.m
Normal file
@ -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.';
|
||||
55
GPE Solver/+Helper/onenoteccdata.m
Normal file
55
GPE Solver/+Helper/onenoteccdata.m
Normal file
@ -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');
|
||||
148
GPE Solver/+Helper/parforNotifications.m
Normal file
148
GPE Solver/+Helper/parforNotifications.m
Normal file
@ -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
|
||||
|
||||
|
||||
820
GPE Solver/+Helper/screencapture.m
Normal file
820
GPE Solver/+Helper/screencapture.m
Normal file
@ -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 <a href="http://www.mathworks.com/matlabcentral/fileexchange/authors/27420">MathWorks File Exchange</a>
|
||||
|
||||
% 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<NASGU>
|
||||
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<AGROW>
|
||||
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<LERR>
|
||||
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<ASGLU> 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<ASGLU>
|
||||
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
|
||||
50
GPE Solver/+Plotter/Analysis.m
Normal file
50
GPE Solver/+Plotter/Analysis.m
Normal file
@ -0,0 +1,50 @@
|
||||
set(0,'defaulttextInterpreter','latex')
|
||||
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
||||
format long
|
||||
|
||||
runIdx = 6;
|
||||
|
||||
load(sprintf('./Data/Run_%03i/psi_gs.mat',runIdx),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
|
||||
|
||||
x = Transf.x*Params.l0*1e6;
|
||||
y = Transf.y*Params.l0*1e6;
|
||||
z = Transf.z*Params.l0*1e6;
|
||||
%percentcomplete = linspace(0,1,Params.cut_off/200);
|
||||
|
||||
dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);
|
||||
%Plotting
|
||||
subplot(2,3,1)
|
||||
n = abs(psi).^2;
|
||||
nxz = squeeze(trapz(n*dy,2));
|
||||
nyz = squeeze(trapz(n*dx,1));
|
||||
nxy = squeeze(trapz(n*dz,3));
|
||||
|
||||
plotxz = pcolor(x,z,nxz');
|
||||
set(plotxz, 'EdgeColor', 'none');
|
||||
xlabel('$x$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
|
||||
|
||||
subplot(2,3,2)
|
||||
plotyz = pcolor(y,z,nyz');
|
||||
set(plotyz, 'EdgeColor', 'none');
|
||||
xlabel('$y$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
|
||||
|
||||
subplot(2,3,3)
|
||||
plotxy = pcolor(x,y,nxy');
|
||||
set(plotxy, 'EdgeColor', 'none');
|
||||
xlabel('$x$ [$\mu$m]'); ylabel('$y$ [$\mu$m]');
|
||||
|
||||
subplot(2,3,4)
|
||||
plot(-log10(Observ.residual),'-b')
|
||||
ylabel('$-\mathrm{log}_{10}(r)$'); xlabel('steps');
|
||||
|
||||
subplot(2,3,5)
|
||||
plot(Observ.EVec,'-b')
|
||||
ylabel('$E$'); xlabel('steps');
|
||||
|
||||
subplot(2,3,6)
|
||||
plot(Observ.mucVec,'-b')
|
||||
ylabel('$\mu$'); xlabel('steps');
|
||||
% xlim([0,1]); ylim([0,8]);
|
||||
% xlim([0,1]); ylim([0,8]);
|
||||
|
||||
Ecomp = energy_components(psi,Params,Transf,VDk,V);
|
||||
77
GPE Solver/+Plotter/MakeMovie.m
Normal file
77
GPE Solver/+Plotter/MakeMovie.m
Normal file
@ -0,0 +1,77 @@
|
||||
set(0,'defaulttextInterpreter','latex')
|
||||
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
||||
|
||||
RunIdx = 1;
|
||||
|
||||
FileDir = dir(sprintf('./Data/Run_%03i/TimeEvolution/*.mat',RunIdx));
|
||||
NumFiles = numel(FileDir);
|
||||
QuenchSettings = load(sprintf('./Data/Run_%03i/QuenchSettings',RunIdx),'Quench','Params','Transf','VDk','V');
|
||||
Transf = QuenchSettings.Transf; Params = QuenchSettings.Params;
|
||||
x = Transf.x; y = Transf.y; z = Transf.z;
|
||||
dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);
|
||||
|
||||
mkdir(sprintf('./Data/Run_%03i/Figures',RunIdx))
|
||||
outputVideo = VideoWriter(fullfile('./Data/Movie.avi'));
|
||||
outputVideo.FrameRate = 10;
|
||||
open(outputVideo)
|
||||
|
||||
figure(1);
|
||||
x0 = 800;
|
||||
y0 = 200;
|
||||
width = 800;
|
||||
height = 600;
|
||||
set(gcf,'position',[x0,y0,width,height])
|
||||
|
||||
EVecTemp = [];
|
||||
|
||||
for ii = 2:(NumFiles-1)
|
||||
load(sprintf('./Data/Run_%03i/TimeEvolution/psi_%i.mat',RunIdx,ii),'psi','muchem','T','Observ','t_idx');
|
||||
|
||||
%Plotting
|
||||
subplot(2,3,1)
|
||||
n = abs(psi).^2;
|
||||
nxz = squeeze(trapz(n*dy,2));
|
||||
nyz = squeeze(trapz(n*dx,1));
|
||||
nxy = squeeze(trapz(n*dz,3));
|
||||
|
||||
plotxz = pcolor(x,z,nxz'); shading interp
|
||||
set(plotxz, 'EdgeColor', 'none');
|
||||
xlabel('$x$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
|
||||
|
||||
subplot(2,3,2)
|
||||
plotyz = pcolor(y,z,nyz'); shading interp
|
||||
set(plotyz, 'EdgeColor', 'none');
|
||||
xlabel('$y$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
|
||||
|
||||
subplot(2,3,3)
|
||||
plotxy = pcolor(x,y,nxy'); shading interp
|
||||
set(plotxy, 'EdgeColor', 'none');
|
||||
xlabel('$x$ [$\mu$m]'); ylabel('$y$ [$\mu$m]');
|
||||
|
||||
subplot(2,3,4)
|
||||
plot(Observ.tVecPlot*1000/Params.w0,Observ.NormVec,'-b')
|
||||
ylabel('Normalization'); xlabel('$t$ [$m$s]');
|
||||
|
||||
subplot(2,3,5)
|
||||
plot(Observ.tVecPlot*1000/Params.w0,1-2*Observ.PCVec/pi,'-b')
|
||||
ylabel('Coherence'); xlabel('$t$ [$m$s]');
|
||||
ylim([0,1])
|
||||
|
||||
subplot(2,3,6)
|
||||
plot(Observ.tVecPlot*1000/Params.w0,Observ.EVec,'-b')
|
||||
ylabel('E'); xlabel('$t$ [$m$s]');
|
||||
|
||||
tVal = Observ.tVecPlot(end)*1000/Params.w0;
|
||||
sgtitle(sprintf('$\\mu =%.3f \\hbar\\omega_0$, $T=%.1f$nK, $t=%.1f$ms',muchem,T,tVal))
|
||||
|
||||
drawnow
|
||||
saveas(gcf,sprintf('./Data/Run_%03i/Figures/Image_%i.jpg',RunIdx,ii))
|
||||
img = imread(sprintf('./Data/Run_%03i/Figures/Image_%i.jpg',RunIdx,ii));
|
||||
writeVideo(outputVideo,img)
|
||||
% hold off;
|
||||
clf
|
||||
end
|
||||
|
||||
close(outputVideo)
|
||||
close(figure(1))
|
||||
delete(sprintf('./Data/Run_%03i/Figures/*.jpg',RunIdx)) % deleting images after movie is made
|
||||
58
GPE Solver/+Plotter/OrderParameter_m.m
Normal file
58
GPE Solver/+Plotter/OrderParameter_m.m
Normal file
@ -0,0 +1,58 @@
|
||||
function [m_Order] = OrderParameter_m(psi,Transf,Params,VDk,V,T,muchem)
|
||||
|
||||
NumRealiz = 100;
|
||||
|
||||
Mx = numel(Transf.x);
|
||||
My = numel(Transf.y);
|
||||
Mz = numel(Transf.z);
|
||||
|
||||
r = normrnd(0,1,size(psi));
|
||||
theta = rand(size(psi));
|
||||
noise = r.*exp(2*pi*1i*theta);
|
||||
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
Gamma = 1-1i*Params.gamma_S;
|
||||
dt = Params.dt;
|
||||
|
||||
avgpsi = 0;
|
||||
avgpsi2 = 0;
|
||||
|
||||
for jj = 1:NumRealiz
|
||||
%generate initial state
|
||||
xi = sqrt(2*Params.gamma_S*Params.kbol*T*10^(-9)*dt/(Params.hbar*Params.w0*Transf.dx*Transf.dy*Transf.dz));
|
||||
swapx = randi(length(Transf.x),1,length(Transf.x));
|
||||
swapy = randi(length(Transf.y),1,length(Transf.y));
|
||||
swapz = randi(length(Transf.z),1,length(Transf.z));
|
||||
psi_j = psi + xi * noise(swapx,swapy,swapz);
|
||||
|
||||
% --- % propagate forward in time 1 time step:
|
||||
%kin
|
||||
psi_j = fftn(psi_j);
|
||||
psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
|
||||
psi_j = ifftn(psi_j);
|
||||
|
||||
%DDI
|
||||
frho = fftn(abs(psi_j).^2);
|
||||
Phi = real(ifftn(frho.*VDk));
|
||||
|
||||
%Real-space
|
||||
psi_j = psi_j.*exp(-1i*Gamma*dt*(V + Params.gs*abs(psi_j).^2 + Params.gammaQF*abs(psi_j).^3 + Params.gdd*Phi - muchem));
|
||||
|
||||
%kin
|
||||
psi_j = fftn(psi_j);
|
||||
psi_j = psi_j.*exp(-0.5*1i*Gamma*dt*KEop);
|
||||
psi_j = ifftn(psi_j);
|
||||
|
||||
%Projection
|
||||
kcut = sqrt(2*Params.e_cut);
|
||||
K = (Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2)<kcut.^2;
|
||||
psi_j = ifftn(K.*fftn(psi_j));
|
||||
|
||||
% --- %
|
||||
avgpsi = avgpsi + abs(sum(psi_j(:)))/NumRealiz;
|
||||
avgpsi2 = avgpsi2 + sum(abs(psi_j(:)).^2)/NumRealiz;
|
||||
|
||||
end
|
||||
|
||||
m_Order = 1/sqrt(Mx*My*Mz)*avgpsi/sqrt(avgpsi2);
|
||||
end
|
||||
47
GPE Solver/+Plotter/runningplot.m
Normal file
47
GPE Solver/+Plotter/runningplot.m
Normal file
@ -0,0 +1,47 @@
|
||||
function RPlot = runningplot(psi,Params,Transf,Observ)
|
||||
set(0,'defaulttextInterpreter','latex')
|
||||
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
||||
|
||||
format long
|
||||
x = Transf.x*Params.l0*1e6;
|
||||
y = Transf.y*Params.l0*1e6;
|
||||
z = Transf.z*Params.l0*1e6;
|
||||
%percentcomplete = linspace(0,1,Params.cut_off/200);
|
||||
|
||||
dx = x(2)-x(1); dy = y(2)-y(1); dz = z(2)-z(1);
|
||||
|
||||
%Plotting
|
||||
|
||||
subplot(2,3,1)
|
||||
n = abs(psi).^2;
|
||||
nxz = squeeze(trapz(n*dy,2));
|
||||
nyz = squeeze(trapz(n*dx,1));
|
||||
nxy = squeeze(trapz(n*dz,3));
|
||||
|
||||
plotxz = pcolor(x,z,nxz');
|
||||
set(plotxz, 'EdgeColor', 'none');
|
||||
xlabel('$x$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
|
||||
|
||||
subplot(2,3,2)
|
||||
plotyz = pcolor(y,z,nyz');
|
||||
set(plotyz, 'EdgeColor', 'none');
|
||||
xlabel('$y$ [$\mu$m]'); ylabel('$z$ [$\mu$m]');
|
||||
|
||||
subplot(2,3,3)
|
||||
plotxy = pcolor(x,y,nxy');
|
||||
set(plotxy, 'EdgeColor', 'none');
|
||||
xlabel('$x$ [$\mu$m]'); ylabel('$y$ [$\mu$m]');
|
||||
|
||||
subplot(2,3,4)
|
||||
plot(-log10(Observ.residual),'-b')
|
||||
ylabel('$-\mathrm{log}_{10}(r)$'); xlabel('steps');
|
||||
|
||||
subplot(2,3,5)
|
||||
plot(Observ.EVec,'-b')
|
||||
ylabel('$E$'); xlabel('steps');
|
||||
|
||||
subplot(2,3,6)
|
||||
plot(Observ.mucVec,'-b')
|
||||
ylabel('$\mu$'); xlabel('steps');
|
||||
% xlim([0,1]); ylim([0,8]);
|
||||
% xlim([0,1]); ylim([0,8]);
|
||||
94
GPE Solver/+Simulator/Initialize.m
Normal file
94
GPE Solver/+Simulator/Initialize.m
Normal file
@ -0,0 +1,94 @@
|
||||
function [psi,V,VDk,gammaQF] = Initialize(Params,Transf)
|
||||
|
||||
format long
|
||||
X = Transf.X; Y = Transf.Y; Z = Transf.Z;
|
||||
m = Params.m;
|
||||
|
||||
Zcutoff = Params.Lz/2;
|
||||
|
||||
% == Potential == %
|
||||
V = 0.5*(Params.gx.*X.^2+Params.gy.*Y.^2+Params.gz*Z.^2);
|
||||
|
||||
% == Calulating the DDIs == %
|
||||
% For a cylindrical cutoff, we first construct a kr grid based on the 3D parameters using Bessel quadrature
|
||||
loadDDI = 1;
|
||||
|
||||
if loadDDI == 1
|
||||
VDk = load(sprintf('./Data/VDk_M.mat'));
|
||||
VDk = VDk.VDk;
|
||||
else
|
||||
Params.Lr = 0.5*min(Params.Lx,Params.Ly);
|
||||
Params.Nr = max(Params.Nx,Params.Ny);
|
||||
[TransfRad] = setup_space_radial(Params); %morder really doesn't matter
|
||||
VDk = VDcutoff(TransfRad.kr,TransfRad.kz,TransfRad.Rmax,Zcutoff);
|
||||
|
||||
disp('Calculated radial grid and cutoff')
|
||||
|
||||
% VDk = interp2(DDI.kz,DDI.kr,DDI.VDk,Transf.kz,Transf.kr,'spline');
|
||||
fullkr = [-flip(TransfRad.kr)',TransfRad.kr'];
|
||||
[KR,KZ] = ndgrid(fullkr,TransfRad.kz);
|
||||
|
||||
[KX3D,KY3D,KZ3D] = ndgrid(ifftshift(Transf.kx),ifftshift(Transf.ky),ifftshift(Transf.kz));
|
||||
KR3D = sqrt(KX3D.^2 + KY3D.^2);
|
||||
fullVDK = [flip(VDk',2),VDk']';
|
||||
VDk = interpn(KR,KZ,fullVDK,KR3D,KZ3D,'spline',-1/3); %Interpolating the radial VDk onto a new grid
|
||||
VDk = fftshift(VDk);
|
||||
save(sprintf('./Data/VDk_M.mat'),'VDk');
|
||||
end
|
||||
disp('Finished DDI')
|
||||
|
||||
|
||||
%XXXX
|
||||
% alph=acos((Transf.KZ*cos(Params.theta)+Transf.KX*sin(Params.theta))./sqrt(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2));
|
||||
% alph(1)=pi/2;
|
||||
% rhoc = max([abs(Transf.x),abs(Transf.y)]);
|
||||
% KR = sqrt(Transf.KX.^2+Transf.KY.^2);
|
||||
% func = @(n,u,v) v.^2./(u.^2+v.^2).*(v.*besselj(n,u).*besselk(n+1,v) - u.*besselj(n+1,u).*besselk(n,v));
|
||||
% % func2 = @(n,u) (8*besselj(n,u) - 2*u.*besselj(n+1,u))./u.^2;
|
||||
%
|
||||
% VDk = -0.5*func(0,KR*rhoc,abs(Transf.KZ)*rhoc) + (Transf.KX.^2./KR.^2 - 0.5).*func(2,KR*rhoc,abs(Transf.KZ)*rhoc);
|
||||
% VDk=(1/3)*(3*(cos(alph).^2)-1) - VDk;
|
||||
%
|
||||
% VDk(KR==0) = -1/3 + 1/2*abs(Transf.KZ(KR==0))*rhoc.*besselk(1,abs(Transf.KZ(KR==0))*rhoc);
|
||||
% VDk(Transf.KZ==0) = 1/6 + (Transf.KX(Transf.KZ==0).^2-Transf.KY(Transf.KZ==0).^2)./(KR(Transf.KZ==0).^2).*(1/2 - besselj(1,KR(Transf.KZ==0)*rhoc)./(KR(Transf.KZ==0)*rhoc));
|
||||
% VDk(1,1,1) = 1/6;
|
||||
%XXXX
|
||||
|
||||
|
||||
% % == Setting up the initial wavefunction == %
|
||||
loadstate = 0;
|
||||
|
||||
if loadstate == 1
|
||||
loadnumber = 1250;
|
||||
% load(sprintf('./Data/Seed/psi_%i.mat',loadnumber),'psi');
|
||||
load(sprintf('./Data/Run_004/psi_gs.mat'),'psi');
|
||||
|
||||
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
||||
|
||||
else
|
||||
ellx = sqrt(Params.hbar/(Params.m*Params.wx))/Params.l0;
|
||||
elly = sqrt(Params.hbar/(Params.m*Params.wy))/Params.l0;
|
||||
ellz = sqrt(Params.hbar/(Params.m*Params.wz))/Params.l0;
|
||||
|
||||
Rx = 4*sqrt(2)*ellx; Ry = 4*sqrt(2)*elly; Rz = sqrt(2)*ellz;
|
||||
X0 = 0.0*Transf.Xmax; Y0 = 0.0*Transf.Ymax; Z0 = 0*Transf.Zmax;
|
||||
|
||||
% psi = exp(-(X-X0).^2/Rx^2-(Y-Y0).^2/Ry^2-(Z-Z0).^2/Rz^2);
|
||||
% cur_norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
% psi = psi/sqrt(cur_norm);
|
||||
|
||||
psiz = exp(-(Z-Z0).^2/Rz^2)/sqrt(ellz*sqrt(pi));
|
||||
psi2d = load(sprintf('./Data/Seed/psi_2d_SS.mat'),'psiseed_2d'); psi2d = psi2d.psiseed_2d;
|
||||
psi = psiz.*repmat(psi2d,[1 1 length(Transf.z)]);
|
||||
|
||||
%add some noise
|
||||
r = normrnd(0,1,size(X));
|
||||
theta = rand(size(X));
|
||||
noise = r.*exp(2*pi*1i*theta);
|
||||
psi = psi + 0.00*noise;
|
||||
|
||||
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
Norm
|
||||
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
||||
end
|
||||
18
GPE Solver/+Simulator/PhaseCoherence.m
Normal file
18
GPE Solver/+Simulator/PhaseCoherence.m
Normal file
@ -0,0 +1,18 @@
|
||||
function [PhaseC] = PhaseCoherence(psi,Transf,Params)
|
||||
|
||||
norm = sum(sum(sum(abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = psi/sqrt(norm);
|
||||
|
||||
NumGlobalShifts = 800;
|
||||
betas = []; phishift = [];
|
||||
for jj = 1:NumGlobalShifts
|
||||
phishift(jj) = -pi + 2*pi*(jj-1)/NumGlobalShifts;
|
||||
betas(jj) = sum(sum(sum(abs(angle(psi*exp(-1i*phishift(jj)))).*abs(psi).^2)));
|
||||
end
|
||||
[minbeta,minidx] = min(betas);
|
||||
|
||||
psi = psi*exp(-1i*phishift(minidx));
|
||||
phi = angle(psi);
|
||||
|
||||
avgphi = sum(sum(sum(phi.*abs(psi).^2,1),2),3)*Transf.dx*Transf.dy*Transf.dz;
|
||||
PhaseC = sum(sum(sum(abs(angle(psi)-avgphi).*abs(psi).^2)))*Transf.dx*Transf.dy*Transf.dz;
|
||||
54
GPE Solver/+Simulator/VDcutoff.m
Normal file
54
GPE Solver/+Simulator/VDcutoff.m
Normal file
@ -0,0 +1,54 @@
|
||||
function VDk = VDcutoff(kr, kz, Rmax, Zmax, Nr)
|
||||
% makes the dipolar direct interaction matrix, size numel(kr) * numel(kz)
|
||||
% Rmax and Zmax are the interaction cutoffs. I use 4*(where the density goes to 10^-4 of its peak)
|
||||
% VDk needs to be multiplied by Cdd
|
||||
% approach is that of Lu, PRA 82, 023622 (2010)
|
||||
% blame Danny Baillie, 9 Aug 2011
|
||||
|
||||
% accuracy inputs for numerical integration
|
||||
if(nargin==4)
|
||||
Nr = 5e4;
|
||||
end
|
||||
Nz = 64;
|
||||
farRmultiple = 2000;
|
||||
|
||||
% analytical transform without cutoff
|
||||
[KR, KZ]=ndgrid(kr,kz);
|
||||
Ksq = KR.^2 + KZ.^2;
|
||||
cossq = KZ.^2./Ksq;
|
||||
VDk = cossq-1/3;
|
||||
|
||||
% analytical cutoff for slice 0<z<Zmax, 0<r<Inf Ronen, PRL 98, 030406 (2007)
|
||||
sinsq = 1 - cossq;
|
||||
VDk = VDk + exp(-Zmax*KR).*( sinsq .* cos(Zmax * KZ) - sqrt(sinsq.*cossq).*sin(Zmax * KZ) );
|
||||
|
||||
% midpoint grids for the integration over 0<z<Zmax, Rmax<r<farRmultiple*Rmax (i.e. starts at Rmax)
|
||||
dr=(farRmultiple-1)*Rmax/Nr;
|
||||
r = ((1:Nr)'-0.5)*dr+Rmax;
|
||||
dz=Zmax/Nz;
|
||||
z = ((1:Nz)-0.5)*dz;
|
||||
[R, Z] = ndgrid(r,z);
|
||||
|
||||
Rsq = R.^2 + Z.^2;
|
||||
|
||||
% real space interaction to be transformed
|
||||
igrandbase = (1 - 3*Z.^2./Rsq)./Rsq.^(3/2);
|
||||
|
||||
% do the Fourier transform numerically
|
||||
% prestore to ensure each besselj is only calculated once
|
||||
% cell is faster than (:,:,krn) slicing
|
||||
Nkr = numel(kr);
|
||||
besselr = cell(Nkr,1);
|
||||
for krn = 1:Nkr
|
||||
besselr{krn} = repmat(r.*besselj(0,kr(krn)*r),1,Nz);
|
||||
end
|
||||
for kzn = 1:numel(kz) % what goes wrong when kzn = 33?
|
||||
igrandbasez = repmat(cos(kz(kzn)*z),Nr,1) .* igrandbase;
|
||||
for krn = 1:Nkr
|
||||
igrand = igrandbasez.*besselr{krn};
|
||||
VDk(krn,kzn) = VDk(krn,kzn) - sum(igrand(:))*dz*dr;
|
||||
end
|
||||
end
|
||||
|
||||
% why are so few z values used?
|
||||
% are the z and kz values without the bounds intended?
|
||||
28
GPE Solver/+Simulator/chemicalpotential.m
Normal file
28
GPE Solver/+Simulator/chemicalpotential.m
Normal file
@ -0,0 +1,28 @@
|
||||
function muchem = chemicalpotential(psi,Params,Transf,VDk,V)
|
||||
|
||||
%Parameters
|
||||
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
|
||||
% DDIs
|
||||
frho=fftn(abs(psi).^2);
|
||||
Phi=real(ifftn(frho.*VDk));
|
||||
|
||||
Eddi = (Params.gdd*Phi.*abs(psi).^2);
|
||||
|
||||
%Kinetic energy
|
||||
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||
|
||||
%Potential energy
|
||||
Epot = V.*abs(psi).^2;
|
||||
|
||||
%Contact interactions
|
||||
Eint = Params.gs*abs(psi).^4;
|
||||
|
||||
%Quantum fluctuations
|
||||
Eqf = Params.gammaQF*abs(psi).^5;
|
||||
|
||||
%Total energy
|
||||
muchem = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz; %
|
||||
muchem = muchem / Params.N;
|
||||
35
GPE Solver/+Simulator/energy_components.m
Normal file
35
GPE Solver/+Simulator/energy_components.m
Normal file
@ -0,0 +1,35 @@
|
||||
function E = energy_components(psi,Params,Transf,VDk,V)
|
||||
|
||||
%Parameters
|
||||
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||
|
||||
% DDIs
|
||||
frho = fftn(abs(psi).^2);
|
||||
Phi = real(ifftn(frho.*VDk));
|
||||
|
||||
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
|
||||
E.Eddi = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
%Kinetic energy
|
||||
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
|
||||
|
||||
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||
E.Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||
|
||||
% Potential energy
|
||||
Epot = V.*abs(psi).^2;
|
||||
E.Epot = trapz(Epot(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
%Contact interactions
|
||||
Eint = 0.5*Params.gs*abs(psi).^4;
|
||||
E.Eint = trapz(Eint(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
%Quantum fluctuations
|
||||
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
|
||||
E.Eqf = trapz(Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
% plot(Transf.x,abs(psi(:,end/2,end/2+1)).^2)
|
||||
31
GPE Solver/+Simulator/energytotal.m
Normal file
31
GPE Solver/+Simulator/energytotal.m
Normal file
@ -0,0 +1,31 @@
|
||||
function E = energytotal(psi,Params,Transf,VDk,V)
|
||||
|
||||
%Parameters
|
||||
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
normfac = Params.Lx*Params.Ly*Params.Lz/numel(psi);
|
||||
|
||||
% DDIs
|
||||
frho = fftn(abs(psi).^2);
|
||||
Phi = real(ifftn(frho.*VDk));
|
||||
|
||||
Eddi = 0.5*Params.gdd*Phi.*abs(psi).^2;
|
||||
|
||||
% EddiTot = trapz(Eddi(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
|
||||
%Kinetic energy
|
||||
% psik = ifftshift(fftn(fftshift(psi)))*normfac;
|
||||
|
||||
Ekin = KEop.*abs(fftn(psi)*normfac).^2;
|
||||
Ekin = trapz(Ekin(:))*Transf.dkx*Transf.dky*Transf.dkz/(2*pi)^3;
|
||||
|
||||
% Potential energy
|
||||
Epot = V.*abs(psi).^2;
|
||||
|
||||
%Contact interactions
|
||||
Eint = 0.5*Params.gs*abs(psi).^4;
|
||||
|
||||
%Quantum fluctuations
|
||||
Eqf = 0.4*Params.gammaQF*abs(psi).^5;
|
||||
|
||||
E = Ekin + trapz(Epot(:) + Eint(:) + Eddi(:) + Eqf(:))*Transf.dx*Transf.dy*Transf.dz;
|
||||
24
GPE Solver/+Simulator/norm_resid.m
Normal file
24
GPE Solver/+Simulator/norm_resid.m
Normal file
@ -0,0 +1,24 @@
|
||||
function res = norm_resid(psi,Params,Transf,VDk,V,muchem)
|
||||
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
|
||||
% DDIs
|
||||
frho=fftn(abs(psi).^2);
|
||||
Phi=real(ifftn(frho.*VDk));
|
||||
|
||||
Eddi = Params.gdd*Phi.*psi;
|
||||
|
||||
%Kinetic energy
|
||||
Ekin = ifftn(KEop.*fftn(psi));
|
||||
|
||||
%Potential energy
|
||||
Epot = V.*psi;
|
||||
|
||||
%Contact interactions
|
||||
Eint = Params.gs*abs(psi).^2.*psi;
|
||||
|
||||
%Quantum fluctuations
|
||||
Eqf = Params.gammaQF*abs(psi).^3.*psi;
|
||||
|
||||
%Total energy
|
||||
res = trapz(abs(Ekin(:) + Epot(:) + Eint(:) + Eddi(:) + Eqf(:) - muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz)/trapz(abs(muchem*psi(:))*Transf.dx*Transf.dy*Transf.dz);
|
||||
99
GPE Solver/+Simulator/parameters.m
Normal file
99
GPE Solver/+Simulator/parameters.m
Normal file
@ -0,0 +1,99 @@
|
||||
function [Params] = parameters()
|
||||
%%--%% Parameters %%--%%
|
||||
%========= Simulation =========%
|
||||
pert=0; %0 = no perturbation during real-time, 1=perturbation
|
||||
%method=1; %0 = normal dipolar potential, 1=spherical cut-off, 2=cylindrical cut-off
|
||||
|
||||
%Energy tolerance
|
||||
Params.Etol=5e-10;
|
||||
Params.rtol = 1e-5;
|
||||
Params.theta=pi/2; %pi/2 dipoles along x, theta=0 dipoles along z
|
||||
Params.cut_off=2e6; %sometimes the imaginary time gets a little stuck
|
||||
%even though the solution is good, this just stops it going on forever
|
||||
|
||||
%========= Constants =========%
|
||||
hbar = 1.0545718e-34; %Planck constant [J.s]
|
||||
kbol = 1.38064852e-23; %Boltzmann Constant [J/K]
|
||||
mu0 = 1.25663706212e-6; %Vacuum Permeability [N/A^2] --
|
||||
muB = 9.274009994e-24; %Bohr Magneton [J/T]
|
||||
a0 = 5.2917721067e-11; %Bohr radius [m]
|
||||
m0 = 1.660539066e-27; %Atomic mass [kg]
|
||||
w0 = 2*pi*100; %angular frequency unit [s^-1]
|
||||
mu0factor = 0.3049584233607396;% =(m0/me)*pi*alpha^2 -- me=mass of electron, alpha=fine struct. const.
|
||||
% mu0=mu0factor *hbar^2*a0/(m0*muB^2)
|
||||
%=============================%
|
||||
|
||||
%Number of points in each direction
|
||||
Params.Nx = 128;
|
||||
Params.Ny = 128;
|
||||
Params.Nz = 96;
|
||||
|
||||
%Dimensions (in units of l0)
|
||||
Params.Lx = 40;
|
||||
Params.Ly = 40;
|
||||
Params.Lz = 20;
|
||||
|
||||
%Masses
|
||||
Params.m = 162*m0;
|
||||
l0 = sqrt(hbar/(Params.m*w0)); %Defining a harmonic oscillator length
|
||||
|
||||
%Atom numbers
|
||||
% Params.ppum = 2500; %particles per micron
|
||||
% Params.N = Params.Lz*Params.ppum*l0*1e6;
|
||||
|
||||
Params.N = 10^6;
|
||||
|
||||
%Dipole lengths (units of muB)
|
||||
Params.mu = 9.93*muB;
|
||||
|
||||
%scattering lengths
|
||||
Params.as = 86*a0;
|
||||
|
||||
%trapping frequencies
|
||||
Params.wx = 2*pi*125;
|
||||
Params.wy = 2*pi*125;
|
||||
Params.wz = 2*pi*250;
|
||||
|
||||
%Time step
|
||||
Params.dt = 0.0005;
|
||||
Params.mindt = 1e-6; %Minimum size for a time step using adaptive dt
|
||||
|
||||
%Stochastic GPE
|
||||
Params.gamma_S = 7.5*10^(-3); %gamma for the stochastic GPE
|
||||
|
||||
Params.muchem = 12.64*Params.wz/w0;
|
||||
%================ Parameters defined by those above ================%
|
||||
|
||||
% == Calculating quantum fluctuations == %
|
||||
eps_dd = Params.add/Params.as;
|
||||
if eps_dd == 0
|
||||
Q5 = 1;
|
||||
elseif eps_dd == 1
|
||||
Q5 = 3*sqrt(3)/2;
|
||||
else
|
||||
yeps = (1-eps_dd)/(3*eps_dd);
|
||||
Q5 = (3*eps_dd)^(5/2)*( (8+26*yeps+33*yeps^2)*sqrt(1+yeps) + 15*yeps^3*log((1+sqrt(1+yeps))/sqrt(yeps)) )/48;
|
||||
Q5 = real(Q5);
|
||||
end
|
||||
|
||||
Params.gammaQF = 128/3*sqrt(pi*(Params.as/l0)^5)*Q5;
|
||||
|
||||
|
||||
|
||||
%Contact interaction strength (units of l0/m)
|
||||
Params.gs = 4*pi*Params.as/l0;
|
||||
|
||||
%Dipole lengths
|
||||
Params.add = mu0*Params.mu^2*Params.m/(12*pi*hbar^2);
|
||||
|
||||
%DDI strength
|
||||
Params.gdd = 12*pi*Params.add/l0; %sometimes the 12 is a 4? --> depends on how Vdk (DDI) is defined
|
||||
|
||||
%Trap gamma
|
||||
Params.gx=(Params.wx/w0)^2;
|
||||
Params.gy=(Params.wy/w0)^2;
|
||||
Params.gz=(Params.wz/w0)^2;
|
||||
|
||||
%Loading the rest into Params
|
||||
Params.hbar = hbar; Params.kbol = kbol; Params.mu0 = mu0; Params.muB = muB; Params.a0 = a0;
|
||||
Params.w0 = w0; Params.l0 = l0;
|
||||
23
GPE Solver/+Simulator/running_file.m
Normal file
23
GPE Solver/+Simulator/running_file.m
Normal file
@ -0,0 +1,23 @@
|
||||
%-% Running file %-%
|
||||
clearvars
|
||||
|
||||
njob = 6;
|
||||
|
||||
mkdir(sprintf('./Data'))
|
||||
mkdir(sprintf('./Data/Run_%03i',njob))
|
||||
|
||||
%Obtain simulation parameters
|
||||
[Params] = parameters();
|
||||
|
||||
%Set up spatial grids and transforms
|
||||
[Transf] = setup_space(Params);
|
||||
|
||||
[psi,V,VDk] = Initialize(Params,Transf);
|
||||
|
||||
% --- Initialize
|
||||
Observ.EVec = []; Observ.NormVec = []; Observ.PCVec = []; Observ.tVecPlot = []; Observ.mucVec = [];
|
||||
t_idx = 1; %Start at t = 0;
|
||||
Observ.res_idx = 1;
|
||||
|
||||
|
||||
[psi] = ssfm_imag(psi,Params,Transf,VDk,V,njob,t_idx,Observ);
|
||||
32
GPE Solver/+Simulator/setup_space.m
Normal file
32
GPE Solver/+Simulator/setup_space.m
Normal file
@ -0,0 +1,32 @@
|
||||
function [Transf] = setup_space(Params)
|
||||
Transf.Xmax = 0.5*Params.Lx;
|
||||
Transf.Ymax = 0.5*Params.Ly;
|
||||
Transf.Zmax = 0.5*Params.Lz;
|
||||
|
||||
Nz = Params.Nz; Nx = Params.Nx; Ny = Params.Ny;
|
||||
|
||||
% Fourier grids
|
||||
x = linspace(-0.5*Params.Lx,0.5*Params.Lx-Params.Lx/Params.Nx,Params.Nx);
|
||||
Kmax = pi*Params.Nx/Params.Lx;
|
||||
kx = linspace(-Kmax,Kmax,Nx+1);
|
||||
kx = kx(1:end-1); dkx = kx(2)-kx(1);
|
||||
kx = fftshift(kx);
|
||||
|
||||
y = linspace(-0.5*Params.Ly,0.5*Params.Ly-Params.Ly/Params.Ny,Params.Ny);
|
||||
Kmax = pi*Params.Ny/Params.Ly;
|
||||
ky = linspace(-Kmax,Kmax,Ny+1);
|
||||
ky = ky(1:end-1); dky = ky(2)-ky(1);
|
||||
ky = fftshift(ky);
|
||||
|
||||
z = linspace(-0.5*Params.Lz,0.5*Params.Lz-Params.Lz/Params.Nz,Params.Nz);
|
||||
Kmax = pi*Params.Nz/Params.Lz;
|
||||
kz = linspace(-Kmax,Kmax,Nz+1);
|
||||
kz = kz(1:end-1); dkz = kz(2)-kz(1);
|
||||
kz = fftshift(kz);
|
||||
|
||||
[Transf.X,Transf.Y,Transf.Z]=ndgrid(x,y,z);
|
||||
[Transf.KX,Transf.KY,Transf.KZ]=ndgrid(kx,ky,kz);
|
||||
Transf.x = x; Transf.y = y; Transf.z = z;
|
||||
Transf.kx = kx; Transf.ky = ky; Transf.kz = kz;
|
||||
Transf.dx = x(2)-x(1); Transf.dy = y(2)-y(1); Transf.dz = z(2)-z(1);
|
||||
Transf.dkx = dkx; Transf.dky = dky; Transf.dkz = dkz;
|
||||
311
GPE Solver/+Simulator/setup_space_radial.m
Normal file
311
GPE Solver/+Simulator/setup_space_radial.m
Normal file
@ -0,0 +1,311 @@
|
||||
function [Transf] = setup_space_radial(Params,morder)
|
||||
Zmax = 0.5*Params.Lz;
|
||||
Rmax = Params.Lr;
|
||||
Nz = Params.Nz;
|
||||
Nr = Params.Nr;
|
||||
|
||||
if(nargin==1)
|
||||
morder=0; %only do Bessel J0
|
||||
end
|
||||
|
||||
% Fourier grids
|
||||
z=linspace(-Zmax,Zmax,Nz+1);
|
||||
z=z(1:end-1);
|
||||
dz=z(2)-z(1);
|
||||
Kmax=Nz*2*pi/(4*Zmax);
|
||||
kz=linspace(-Kmax,Kmax,Nz+1);
|
||||
kz=kz(1:end-1);
|
||||
|
||||
% Hankel grids and transform
|
||||
H = hankelmatrix(morder,Rmax,Nr);
|
||||
r=H.r(:);
|
||||
kr=H.kr(:);
|
||||
T = diag(H.J/H.kmax)*H.T*diag(Rmax./H.J)*dz*(2*pi);
|
||||
Tinv = diag(H.J./Rmax)*H.T'*diag(H.kmax./H.J)/dz/(2*pi);
|
||||
wr=H.wr;
|
||||
wk=H.wk;
|
||||
% H.T'*diag(H.J/H.vmax)*H.T*diag(Rmax./H.J)
|
||||
|
||||
[Transf.R,Transf.Z]=ndgrid(r,z);
|
||||
[Transf.KR,Transf.KZ]=ndgrid(kr,kz);
|
||||
Transf.T=T;
|
||||
Transf.Tinv=Tinv;
|
||||
Transf.r=r;
|
||||
Transf.kr=kr;
|
||||
Transf.z=z;
|
||||
Transf.kz=kz;
|
||||
Transf.wr=wr;
|
||||
Transf.wk=wk;
|
||||
Transf.Rmax=Rmax;
|
||||
Transf.Zmax=Zmax;
|
||||
Transf.dz=z(2)-z(1);
|
||||
Transf.dkz=kz(2)-kz(1);
|
||||
%b1=Transf;
|
||||
|
||||
function s_HT = hankelmatrix(order, rmax, Nr, eps_roots)
|
||||
%HANKEL_MATRIX: Generates data to use for Hankel Transforms
|
||||
%
|
||||
% s_HT = hankel_matrix(order, rmax, Nr, eps_roots)
|
||||
%
|
||||
% s_HT = Structure containing data to use for the pQDHT
|
||||
% order = Transform order
|
||||
% rmax = Radial extent of transform
|
||||
% Nr = Number of sample points
|
||||
% eps_roots = Error in estimation of roots of Bessel function (optional)
|
||||
%
|
||||
% s_HT:
|
||||
% order, rmax, Nr = As above
|
||||
% J_roots = Roots of the pth order Bessel fn.
|
||||
% J_roots_N1 = (N+1)th root
|
||||
% r = Radial co-ordinate vector
|
||||
% v = frequency co-ordinate vector
|
||||
% kr = Radial wave number co-ordinate vector
|
||||
% vmax = Limiting frequency
|
||||
% = roots_N1 / (2*pi*rmax)
|
||||
% S = rmax * 2*pi*vmax (S product)
|
||||
% T = Transform matrix
|
||||
% J = Scaling vector
|
||||
% = J_(order+1){roots}
|
||||
%
|
||||
% The algorithm used is that from:
|
||||
% "Computation of quasi-discrete Hankel transforms of the integer
|
||||
% order for propagating optical wave fields"
|
||||
% Manuel Guizar-Sicairos and Julio C. Guitierrez-Vega
|
||||
% J. Opt. Soc. Am. A 21(1) 53-58 (2004)
|
||||
%
|
||||
% The algorithm also calls the function:
|
||||
% zn = bessel_zeros(1, p, Nr+1, 1e-6),
|
||||
% where p and N are defined above, to calculate the roots of the bessel
|
||||
% function. This algorithm is taken from:
|
||||
% "An Algorithm with ALGOL 60 Program for the Computation of the
|
||||
% zeros of the Ordinary Bessel Functions and those of their
|
||||
% Derivatives".
|
||||
% N. M. Temme
|
||||
% Journal of Computational Physics, 32, 270-279 (1979)
|
||||
%
|
||||
% Example: Propagation of radial field
|
||||
%
|
||||
% % Note the use of matrix and element products / divisions
|
||||
% H = hankel_matrix(0, 1e-3, 512);
|
||||
% DR0 = 50e-6;
|
||||
% Ur0 = exp(-(H.r/DR0).^2);
|
||||
% Ukr0 = H.T * (Ur0./H.J);
|
||||
% k0 = 2*pi/800e-9;
|
||||
% kz = realsqrt((k0^2 - H.kr.^2).*(k0>H.kr));
|
||||
% z = (-5e-3:1e-5:5e-3);
|
||||
% Ukrz = (Ukr0*ones(1,length(z))).*exp(i*kz*z);
|
||||
% Urz = (H.T * Ukrz) .* (H.J * ones(1,length(z)));
|
||||
%
|
||||
% See also bessel_zeros, besselj
|
||||
|
||||
if (~exist('eps_roots', 'var')||isemtpy(eps_roots))
|
||||
s_HT.eps_roots = 1e-6;
|
||||
else
|
||||
s_HT.eps_roots = eps_roots;
|
||||
end
|
||||
|
||||
s_HT.order = order;
|
||||
s_HT.rmax = rmax;
|
||||
s_HT.Nr = Nr;
|
||||
|
||||
% Calculate N+1 roots:
|
||||
J_roots = bessel_zeros(1, s_HT.order, s_HT.Nr+1, s_HT.eps_roots);
|
||||
s_HT.J_roots = J_roots(1:end-1);
|
||||
s_HT.J_roots_N1 = J_roots(end);
|
||||
|
||||
% Calculate co-ordinate vectors
|
||||
s_HT.r = s_HT.J_roots * s_HT.rmax / s_HT.J_roots_N1;
|
||||
s_HT.v = s_HT.J_roots / (2*pi * s_HT.rmax);
|
||||
s_HT.kr = 2*pi * s_HT.v;
|
||||
s_HT.kmax = s_HT.J_roots_N1 / (s_HT.rmax);
|
||||
s_HT.vmax = s_HT.J_roots_N1 / (2*pi * s_HT.rmax);
|
||||
s_HT.S = s_HT.J_roots_N1;
|
||||
|
||||
% Calculate hankel matrix and vectors
|
||||
% I use (p=order) and (p1=order+1)
|
||||
Jp = besselj(s_HT.order, (s_HT.J_roots) * (s_HT.J_roots.') / s_HT.S);
|
||||
Jp1 = abs(besselj(s_HT.order+1, s_HT.J_roots));
|
||||
s_HT.T = 2*Jp./(Jp1 * (Jp1.') * s_HT.S);
|
||||
s_HT.J = Jp1;
|
||||
s_HT.wr=2./((s_HT.kmax)^2*abs(Jp1).^2);
|
||||
s_HT.wk=2./((s_HT.rmax)^2*abs(Jp1).^2);
|
||||
|
||||
return
|
||||
|
||||
|
||||
|
||||
|
||||
function z = bessel_zeros(d, a, n, e)
|
||||
%BESSEL_ZEROS: Finds the first n zeros of a bessel function
|
||||
%
|
||||
% z = bessel_zeros(d, a, n, e)
|
||||
%
|
||||
% z = zeros of the bessel function
|
||||
% d = Bessel function type:
|
||||
% 1: Ja
|
||||
% 2: Ya
|
||||
% 3: Ja'
|
||||
% 4: Ya'
|
||||
% a = Bessel order (a>=0)
|
||||
% n = Number of zeros to find
|
||||
% e = Relative error in root
|
||||
%
|
||||
% This function uses the routine described in:
|
||||
% "An Algorithm with ALGOL 60 Program for the Computation of the
|
||||
% zeros of the Ordinary Bessel Functions and those of their
|
||||
% Derivatives".
|
||||
% N. M. Temme
|
||||
% Journal of Computational Physics, 32, 270-279 (1979)
|
||||
|
||||
z = zeros(n, 1);
|
||||
aa = a^2;
|
||||
mu = 4*aa;
|
||||
mu2 = mu^2;
|
||||
mu3 = mu^3;
|
||||
mu4 = mu^4;
|
||||
|
||||
if (d<3)
|
||||
p = 7*mu - 31;
|
||||
p0 = mu - 1;
|
||||
if ((1+p)==p)
|
||||
p1 = 0;
|
||||
q1 = 0;
|
||||
else
|
||||
p1 = 4*(253*mu2 - 3722*mu+17869)*p0/(15*p);
|
||||
q1 = 1.6*(83*mu2 - 982*mu + 3779)/p;
|
||||
end
|
||||
else
|
||||
p = 7*mu2 + 82*mu - 9;
|
||||
p0 = mu + 3;
|
||||
if ((p+1)==1)
|
||||
p1 = 0;
|
||||
q1 = 0;
|
||||
else
|
||||
p1 = (4048*mu4 + 131264*mu3 - 221984*mu2 - 417600*mu + 1012176)/(60*p);
|
||||
q1 = 1.6*(83*mu3 + 2075*mu2 - 3039*mu + 3537)/p;
|
||||
end
|
||||
end
|
||||
|
||||
if (d==1)|(d==4)
|
||||
t = .25;
|
||||
else
|
||||
t = .75;
|
||||
end
|
||||
tt = 4*t;
|
||||
|
||||
if (d<3)
|
||||
pp1 = 5/48;
|
||||
qq1 = -5/36;
|
||||
else
|
||||
pp1 = -7/48;
|
||||
qq1 = 35/288;
|
||||
end
|
||||
|
||||
y = .375*pi;
|
||||
if (a>=3)
|
||||
bb = a^(-2/3);
|
||||
else
|
||||
bb = 1;
|
||||
end
|
||||
a1 = 3*a - 8;
|
||||
% psi = (.5*a + .25)*pi;
|
||||
|
||||
for s=1:n
|
||||
if ((a==0)&(s==1)&(d==3))
|
||||
x = 0;
|
||||
j = 0;
|
||||
else
|
||||
if (s>=a1)
|
||||
b = (s + .5*a - t)*pi;
|
||||
c = .015625/(b^2);
|
||||
x = b - .125*(p0 - p1*c)/(b*(1 - q1*c));
|
||||
else
|
||||
if (s==1)
|
||||
switch (d)
|
||||
case (1)
|
||||
x = -2.33811;
|
||||
case (2)
|
||||
x = -1.17371;
|
||||
case (3)
|
||||
x = -1.01879;
|
||||
otherwise
|
||||
x = -2.29444;
|
||||
end
|
||||
else
|
||||
x = y*(4*s - tt);
|
||||
v = x^(-2);
|
||||
x = -x^(2/3) * (1 + v*(pp1 + qq1*v));
|
||||
end
|
||||
u = x*bb;
|
||||
v = fi(2/3 * (-u)^1.5);
|
||||
w = 1/cos(v);
|
||||
xx = 1 - w^2;
|
||||
c = sqrt(u/xx);
|
||||
if (d<3)
|
||||
x = w*(a + c*(-5/u - c*(6 - 10/xx))/(48*a*u));
|
||||
else
|
||||
x = w*(a + c*(7/u + c*(18 - 14/xx))/(48*a*u));
|
||||
end
|
||||
end
|
||||
j = 0;
|
||||
|
||||
while ((j==0)|((j<5)&(abs(w/x)>e)))
|
||||
xx = x^2;
|
||||
x4 = x^4;
|
||||
a2 = aa - xx;
|
||||
r0 = bessr(d, a, x);
|
||||
j = j+1;
|
||||
if (d<3)
|
||||
u = r0;
|
||||
w = 6*x*(2*a + 1);
|
||||
p = (1 - 4*a2)/w;
|
||||
q = (4*(xx-mu) - 2 - 12*a)/w;
|
||||
else
|
||||
u = -xx*r0/a2;
|
||||
v = 2*x*a2/(3*(aa+xx));
|
||||
w = 64*a2^3;
|
||||
q = 2*v*(1 + mu2 + 32*mu*xx + 48*x4)/w;
|
||||
p = v*(1 + (40*mu*xx + 48*x4 - mu2)/w);
|
||||
end
|
||||
w = u*(1 + p*r0)/(1 + q*r0);
|
||||
x = x+w;
|
||||
end
|
||||
z(s) = x;
|
||||
end
|
||||
end
|
||||
|
||||
|
||||
function FI = fi(y)
|
||||
c1 = 1.570796;
|
||||
if (~y)
|
||||
FI = 0;
|
||||
elseif (y>1e5)
|
||||
FI = c1;
|
||||
else
|
||||
if (y<1)
|
||||
p = (3*y)^(1/3);
|
||||
pp = p^2;
|
||||
p = p*(1 + pp*(pp*(27 - 2*pp) - 210)/1575);
|
||||
else
|
||||
p = 1/(y + c1);
|
||||
pp = p^2;
|
||||
p = c1 - p*(1 + pp*(2310 + pp*(3003 + pp*(4818 + pp*(8591 + pp*16328))))/3465);
|
||||
end
|
||||
pp = (y+p)^2;
|
||||
r = (p - atan(p+y))/pp;
|
||||
FI = p - (1+pp)*r*(1 + r/(p+y));
|
||||
end
|
||||
return
|
||||
|
||||
function Jr = bessr(d, a, x)
|
||||
switch (d)
|
||||
case (1)
|
||||
Jr = besselj(a, x)./besselj(a+1, x);
|
||||
case (2)
|
||||
Jr = bessely(a, x)./bessely(a+1, x);
|
||||
case (3)
|
||||
Jr = a./x - besselj(a+1, x)./besselj(a, x);
|
||||
otherwise
|
||||
Jr = a./x - bessely(a+1, x)./bessely(a, x);
|
||||
end
|
||||
return
|
||||
98
GPE Solver/+Simulator/ssfm_imag.m
Normal file
98
GPE Solver/+Simulator/ssfm_imag.m
Normal file
@ -0,0 +1,98 @@
|
||||
function [psi] = ssfm_imag(psi,Params,Transf,VDk,V,njob,t_idx,Observ)
|
||||
|
||||
set(0,'defaulttextInterpreter','latex')
|
||||
set(groot, 'defaultAxesTickLabelInterpreter','latex'); set(groot, 'defaultLegendInterpreter','latex');
|
||||
|
||||
dt=-1j*abs(Params.dt);
|
||||
|
||||
KEop= 0.5*(Transf.KX.^2+Transf.KY.^2+Transf.KZ.^2);
|
||||
Observ.residual = 1; Observ.res = 1;
|
||||
|
||||
figure(1)
|
||||
muchem = chemicalpotential(psi,Params,Transf,VDk,V);
|
||||
runningplot(psi,Params,Transf,Observ)
|
||||
drawnow
|
||||
|
||||
AdaptIdx = 0;
|
||||
|
||||
while t_idx < Params.cut_off
|
||||
%kin
|
||||
psi = fftn(psi);
|
||||
psi = psi.*exp(-0.5*1i*dt*KEop);
|
||||
psi = ifftn(psi);
|
||||
|
||||
%DDI
|
||||
frho = fftn(abs(psi).^2);
|
||||
Phi = real(ifftn(frho.*VDk));
|
||||
|
||||
%Real-space
|
||||
psi = psi.*exp(-1i*dt*(V + Params.gs*abs(psi).^2 + Params.gammaQF*abs(psi).^3 + Params.gdd*Phi - muchem));
|
||||
|
||||
%kin
|
||||
psi = fftn(psi);
|
||||
psi = psi.*exp(-0.5*1i*dt*KEop);
|
||||
psi = ifftn(psi);
|
||||
|
||||
%Renorm
|
||||
Norm = trapz(abs(psi(:)).^2)*Transf.dx*Transf.dy*Transf.dz;
|
||||
psi = sqrt(Params.N)*psi/sqrt(Norm);
|
||||
|
||||
muchem = chemicalpotential(psi,Params,Transf,VDk,V);
|
||||
|
||||
%Plotting loop
|
||||
if mod(t_idx,1000) == 0
|
||||
|
||||
%Change in Energy
|
||||
E = energytotal(psi,Params,Transf,VDk,V);
|
||||
E = E/Norm;
|
||||
Observ.EVec = [Observ.EVec E];
|
||||
|
||||
%Chemical potential
|
||||
Observ.mucVec = [Observ.mucVec muchem];
|
||||
|
||||
%Normalized residuals
|
||||
res = norm_resid(psi,Params,Transf,VDk,V,muchem);
|
||||
Observ.residual = [Observ.residual res];
|
||||
|
||||
Observ.res_idx = Observ.res_idx + 1;
|
||||
figure(1)
|
||||
runningplot(psi,Params,Transf,Observ)
|
||||
drawnow
|
||||
|
||||
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
|
||||
|
||||
%Adaptive time step -- Careful, this can quickly get out of control
|
||||
relres = abs(Observ.residual(Observ.res_idx)-Observ.residual(Observ.res_idx-1))/Observ.residual(Observ.res_idx);
|
||||
if relres <1e-5
|
||||
if AdaptIdx > 4 && abs(dt) > Params.mindt
|
||||
dt = dt / 2;
|
||||
fprintf('Time step changed to '); disp(dt);
|
||||
AdaptIdx = 0;
|
||||
elseif AdaptIdx > 4 && abs(dt) < Params.mindt
|
||||
break
|
||||
else
|
||||
AdaptIdx = AdaptIdx + 1;
|
||||
end
|
||||
else
|
||||
AdaptIdx = 0;
|
||||
end
|
||||
end
|
||||
if any(isnan(psi(:)))
|
||||
disp('Idiot.')
|
||||
break
|
||||
end
|
||||
t_idx=t_idx+1;
|
||||
end
|
||||
|
||||
|
||||
%Change in Energy
|
||||
E = energytotal(psi,Params,Transf,VDk,V)
|
||||
E = E/Norm;
|
||||
Observ.EVec = [Observ.EVec E];
|
||||
|
||||
% Phase coherence
|
||||
[PhaseC] = PhaseCoherence(psi,Transf,Params);
|
||||
Observ.PCVec = [Observ.PCVec PhaseC];
|
||||
|
||||
Observ.res_idx = Observ.res_idx + 1;
|
||||
save(sprintf('./Data/Run_%03i/psi_gs.mat',njob),'psi','muchem','Observ','t_idx','Transf','Params','VDk','V');
|
||||
233
GaussianBeamABCD/BeamPropagation.py
Normal file
233
GaussianBeamABCD/BeamPropagation.py
Normal file
@ -0,0 +1,233 @@
|
||||
"""
|
||||
@author: Adam Newton Wright, https://github.com/adamnewtonwright/GaussianBeamPropagation
|
||||
"""
|
||||
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import sympy as sym
|
||||
from sympy import oo
|
||||
|
||||
|
||||
# Input Ray parameter, i.e. height and angle
|
||||
def ray(y,theta):
|
||||
'''
|
||||
Parameters
|
||||
----------
|
||||
y : float or integer or sympy symbol in meters
|
||||
The vertical height of a ray.
|
||||
theta : float or integer in radians
|
||||
The angle of divergence of the ray.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mat : 2x1 matrix
|
||||
[
|
||||
[y],
|
||||
[teta]
|
||||
]
|
||||
|
||||
'''
|
||||
|
||||
mat = np.array([[y],[theta]])
|
||||
return mat
|
||||
|
||||
# Ray Transfer Matrix for ideal lens with focal length f
|
||||
def lens(f):
|
||||
'''
|
||||
Parameters
|
||||
----------
|
||||
f : float or integer or sympy symbol in meters
|
||||
Thin lens focal length in meters
|
||||
|
||||
Returns
|
||||
-------
|
||||
mat : 2x2 matrix
|
||||
[
|
||||
[ 1, 0],
|
||||
[-1/f, 1]
|
||||
]
|
||||
|
||||
'''
|
||||
|
||||
mat = np.array([[1,0], [-1/f, 1]])
|
||||
return mat
|
||||
|
||||
# Ray Transfer Matrix for propagation of distance d
|
||||
def prop(d):
|
||||
'''
|
||||
Parameters
|
||||
----------
|
||||
d : float or integer or sympy symbol
|
||||
Distance light is propagating along the z-axis.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mat: 2x2 matrix
|
||||
[
|
||||
[1, d],
|
||||
[0, 1]
|
||||
]
|
||||
|
||||
'''
|
||||
mat = np.array([[1,d], [0,1]])
|
||||
return mat
|
||||
|
||||
# multiplying the matrices together. mat1 is the last matrix the light interacts with
|
||||
def mult(mat1,*argv):
|
||||
'''
|
||||
Parameters
|
||||
----------
|
||||
mat1 : 2x2 ABCD matrix
|
||||
Last matrix light interacts with.
|
||||
*argv : 2x2 ABCD matrices
|
||||
From left to right, the matrices should be entered such that the leftmost matrix interacts
|
||||
with light temporally after the rightmost matrix.
|
||||
|
||||
Returns
|
||||
-------
|
||||
Mat : 2x2 matrix
|
||||
The ABCd matrix describing the whole optical system.
|
||||
|
||||
'''
|
||||
|
||||
Mat = mat1
|
||||
for arg in argv:
|
||||
Mat = np.dot(Mat, arg)
|
||||
return Mat
|
||||
|
||||
# Adding Gaussian beam parameters
|
||||
def Zr(wo, lam):
|
||||
'''
|
||||
Parameters
|
||||
----------
|
||||
wo : float, integer, or symbol
|
||||
Beam waist radius in meters.
|
||||
lam : float, integer, or symbol
|
||||
Wavelength of light in meters.
|
||||
|
||||
Returns
|
||||
-------
|
||||
zr : float, int, symbols
|
||||
Rayleigh range for given beam waist and wavelength.
|
||||
|
||||
'''
|
||||
|
||||
zr = np.pi * wo**2 / lam
|
||||
return zr
|
||||
|
||||
def W0(zr, lam):
|
||||
'''
|
||||
Parameters
|
||||
----------
|
||||
zr : float, integer, symbol
|
||||
Rayleigh range in meters
|
||||
lam : float, integer, symbol
|
||||
Wavelength of light in meters
|
||||
|
||||
Returns
|
||||
-------
|
||||
w0 : float, integer, symbol
|
||||
Beam waist radius in meters
|
||||
|
||||
'''
|
||||
|
||||
w0 = np.sqrt(lam * zr / np.pi)
|
||||
return w0
|
||||
|
||||
# Remember, there should be an i in front of zr
|
||||
# but this complicates the calculations, so we usually just let z = 0
|
||||
# and don't explicitly deal with the i, but still do the math accordingly
|
||||
#def q0_func(z,zr):
|
||||
# qz = z + zr
|
||||
# return qz
|
||||
|
||||
def q1_func(z, w0, lam, mat):
|
||||
'''
|
||||
Parameters
|
||||
----------
|
||||
z : float, int, symbol
|
||||
Position of the beam waist in meters.
|
||||
w0 : float, int, symbol
|
||||
Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
|
||||
lam : float, int, symbol
|
||||
Wavelength of light in meters.
|
||||
mat : float, int, symbol
|
||||
The ABCD 2x2 matrix describing the optical system.
|
||||
|
||||
Returns
|
||||
-------
|
||||
z: float, int, symbol
|
||||
Position of the beam waist after the optical system
|
||||
zr: float, int, symbol
|
||||
Rayleigh range of the beam after the optical system
|
||||
'''
|
||||
|
||||
A = mat[0][0]
|
||||
B = mat[0][1]
|
||||
C = mat[1][0]
|
||||
D = mat[1][1]
|
||||
zr = Zr(w0, lam)
|
||||
real = (A*C*(z**2 + zr**2) + z*(A*D + B*C) + B*D) / (C**2*(z**2 + zr**2) + 2*C*D*z + D**2)
|
||||
imag = (zr * (A*D - B*C)) / (C**2*(z**2 + zr**2) + 2*C*D*z + D**2)
|
||||
z = real
|
||||
zr = imag
|
||||
return z, zr
|
||||
|
||||
def q1_inv_func(z, w0, lam, mat):
|
||||
'''
|
||||
Parameters
|
||||
----------
|
||||
z : float, int, symbol
|
||||
Position of the beam waist in meters.
|
||||
w0 : float, int, symbol
|
||||
Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
|
||||
lam : float, int, symbol
|
||||
Wavelength of light in meters.
|
||||
mat : float, int, symbol
|
||||
The ABCD 2x2 matrix describing the optical system.
|
||||
|
||||
Returns
|
||||
-------
|
||||
R : float, int, symbol
|
||||
Radius of curvature of the wavefront in meters.
|
||||
w : float, int, symbol
|
||||
Radius of the beam in meters.
|
||||
|
||||
'''
|
||||
A = mat[0][0]
|
||||
B = mat[0][1]
|
||||
C = mat[1,0]
|
||||
D = mat[1][1]
|
||||
zr = Zr(w0, lam)
|
||||
real = (A*C*(z**2 + zr**2) + z*(A*D + B*C) + B*D) / (A**2*(z**2 + zr**2) + 2*A*B*z + B**2)
|
||||
imag = -zr * (A*D-B*C) / (A**2 *(z**2 + zr**2) + 2*A*B*z + B**2)
|
||||
R = 1/real
|
||||
w = (-lam / imag / np.pi)**.5
|
||||
return R, w
|
||||
|
||||
def plot(func, var, rang = np.arange(0,3,.01)):
|
||||
'''
|
||||
Parameters
|
||||
----------
|
||||
func : Sympy function of one variable
|
||||
Sympy function defining the beam width after the last optical element.
|
||||
var : sympy variable
|
||||
Variable in func that will be plotted.
|
||||
rang : numpy array
|
||||
Array of the values along the optical axis to be plotted
|
||||
|
||||
Returns
|
||||
-------
|
||||
plot : matplotlib graph
|
||||
Graph of the beam width of var
|
||||
|
||||
|
||||
'''
|
||||
func = sym.lambdify(var, func)
|
||||
plt.figure()
|
||||
plt.plot(rang, func(rang), color = 'b')
|
||||
plt.plot(rang, -func(rang), color = 'b')
|
||||
plt.grid()
|
||||
plt.xlabel('Optic Axis (m)')
|
||||
plt.ylabel('Beam size (m)')
|
||||
plt.show()
|
||||
70
GaussianBeamABCD/DMD_Setup.py
Normal file
70
GaussianBeamABCD/DMD_Setup.py
Normal file
@ -0,0 +1,70 @@
|
||||
import BeamPropagation as bp # This is the script that handles the propagation
|
||||
import sympy as sym # For Symbolic examples
|
||||
import numpy as np # Handling of lists and for plotting
|
||||
import matplotlib.pyplot as plt # Plotting
|
||||
|
||||
|
||||
"""A Gaussian beam can be defined by it's (radial) waist wo, it's Rayleigh range zR, and the location of its waist zO"""
|
||||
w0 = 5.2E-3 # 1mm beam waist
|
||||
lam = 532E-9 # wavelength of 355 nm (UV)
|
||||
zR = bp.Zr(w0, lam) # Rayleigh range in m
|
||||
z0 = 0 # location of waist in m
|
||||
|
||||
"""Define first 4-f optical system using matrices"""
|
||||
d1, d2, d3, f1, f2 = sym.symbols('d1 d2 d3 f1 f2')
|
||||
M = bp.mult(bp.prop(d3),bp.lens(f2),bp.prop(d2), bp.lens(f1), bp.prop(d1))
|
||||
|
||||
"""Use script to do all the ABCD and q-parameter math, and return the waist and radius of curvature functions"""
|
||||
R, w = bp.q1_inv_func(0, w0, lam, M)
|
||||
|
||||
"""Substitute and extract the required separation between lenses of first 4-f system"""
|
||||
distance_between_dmd_first_lens = 250E-3
|
||||
first_focal_length = 250.9E-3
|
||||
second_focal_length = 50E-3
|
||||
demag = 1/5
|
||||
target_w0 = demag*w0
|
||||
w = w.subs(f1, first_focal_length).subs(f2, second_focal_length).subs(d1, distance_between_dmd_first_lens).subs(d3,0)
|
||||
eq = sym.solve(w - target_w0, d2)[0]
|
||||
distance_between_lens_of_first_4f = list(eq.atoms())[0]
|
||||
print('Distance between lenses of first 4-f system = {} mm'.format(distance_between_lens_of_first_4f * 1E3))
|
||||
|
||||
# Sanity check
|
||||
# expansion_factor = w.subs(d2,distance_between_lens_of_first_4f).subs(d3,0)/ w0
|
||||
# print('beam is w = {:.2f} x w0'.format(expansion_factor))
|
||||
|
||||
# """Plot beam propagation up to 3 m after the first 4-f system"""
|
||||
# M = bp.mult(bp.prop(d3),bp.lens(second_focal_length),bp.prop(distance_between_lens_of_first_4f), bp.lens(first_focal_length), bp.prop(distance_between_dmd_first_lens))
|
||||
# R, w = bp.q1_inv_func(0, w0, lam, M)
|
||||
# bp.plot(w,d3, rang = np.arange(0,0.050,.0005))
|
||||
|
||||
|
||||
"""Define the full optical system of two 4-f setups using matrices"""
|
||||
d1, d2, d3, d4, f1, f2, f3 = sym.symbols('d1 d2 d3 d4 f1 f2 f3')
|
||||
M = bp.mult(bp.prop(d4),bp.lens(f3), bp.prop(d3),bp.lens(f2),bp.prop(d2), bp.lens(f1), bp.prop(d1))
|
||||
|
||||
"""Use script to do all the ABCD and q-parameter math, and return the waist and radius of curvature functions"""
|
||||
R, w = bp.q1_inv_func(0, w0, lam, M)
|
||||
|
||||
# """Find the focal length of lens required after the first 4-f system to have a collimated beam, given a certain separation between the first 4-f system and this lens"""
|
||||
distance_between_4fs = 550E-3
|
||||
R_coll = R.subs(d1,distance_between_dmd_first_lens).subs(d2,distance_between_lens_of_first_4f).subs(d3,distance_between_4fs).subs(d4,0).subs(f1,first_focal_length).subs(f2,second_focal_length)
|
||||
f3_coll = sym.solve(1/R_coll,f3)[0]
|
||||
third_focal_length = list(f3_coll.atoms())[0]
|
||||
print('For a fixed separation between first 4-f and third lens of {:.3f} mm, f3 = {:.3f} mm for a collimated beam'.format(distance_between_4fs* 1E3, third_focal_length * 1E3))
|
||||
|
||||
# # """Plot beam propagation up to 3 m after the first 4-f system"""
|
||||
# M = bp.mult(bp.prop(d4),bp.lens(third_focal_length),bp.prop(distance_between_4fs), bp.lens(second_focal_length), bp.prop(distance_between_lens_of_first_4f), bp.lens(first_focal_length), bp.prop(distance_between_dmd_first_lens))
|
||||
# R, w = bp.q1_inv_func(0, w0, lam, M)
|
||||
# bp.plot(w,d4, rang = np.arange(0,0.050,.0005))
|
||||
|
||||
third_focal_length = 501.8E-3
|
||||
R_coll = R.subs(d1,distance_between_dmd_first_lens).subs(d2,distance_between_lens_of_first_4f).subs(d4,0).subs(f1,first_focal_length).subs(f2,second_focal_length).subs(f3,third_focal_length)
|
||||
d3_coll = sym.solve(1/R_coll,d3)[1]
|
||||
distance_between_4fs = list(d3_coll.atoms())[0]
|
||||
print('For a fixed third focal length of {:.3f} mm, d3 = {:.3f} mm, for a collimated beam'.format(third_focal_length* 1E3, distance_between_4fs * 1E3))
|
||||
|
||||
# """Plot beam propagation up to 3 m after the first 4-f system"""
|
||||
# M = bp.mult(bp.prop(d4),bp.lens(third_focal_length),bp.prop(distance_between_4fs), bp.lens(second_focal_length), bp.prop(distance_between_lens_of_first_4f), bp.lens(first_focal_length), bp.prop(distance_between_dmd_first_lens))
|
||||
# R, w = bp.q1_inv_func(0, w0, lam, M)
|
||||
# bp.plot(w,d4, rang = np.arange(0,0.050,.0005))
|
||||
|
||||
589
GaussianBeamABCD/Example.ipynb
Normal file
589
GaussianBeamABCD/Example.ipynb
Normal file
File diff suppressed because one or more lines are too long
359
GaussianBeamABCD/README.md
Normal file
359
GaussianBeamABCD/README.md
Normal file
@ -0,0 +1,359 @@
|
||||
# Gaussian Beam Propagation
|
||||
|
||||
## Import files
|
||||
|
||||
|
||||
```python
|
||||
import BeamProp_Script as bs # This is the script that handles the propagation
|
||||
import sympy as sym # For Symbolic examples
|
||||
import numpy as np # Handling of lists and for plotting
|
||||
import matplotlib.pyplot as plt # Plotting
|
||||
```
|
||||
|
||||
### Let's show what BeamProp_Script has
|
||||
|
||||
|
||||
```python
|
||||
help(bs)
|
||||
```
|
||||
|
||||
Help on module BeamProp_Script:
|
||||
|
||||
NAME
|
||||
BeamProp_Script - Created on Wed Feb 19 15:51:54 2020
|
||||
|
||||
DESCRIPTION
|
||||
@author: wrighta
|
||||
|
||||
FUNCTIONS
|
||||
W0(zr, lam)
|
||||
Parameters
|
||||
----------
|
||||
zr : float, integer, symbol
|
||||
Rayleigh range in meters
|
||||
lam : float, integer, symbol
|
||||
Wavelength of light in meters
|
||||
|
||||
Returns
|
||||
-------
|
||||
w0 : float, integer, symbol
|
||||
Beam waist radius in meters
|
||||
|
||||
Zr(wo, lam)
|
||||
Parameters
|
||||
----------
|
||||
wo : float, integer, or symbol
|
||||
Beam waist radius in meters.
|
||||
lam : float, integer, or symbol
|
||||
Wavelength of light in meters.
|
||||
|
||||
Returns
|
||||
-------
|
||||
zr : float, int, symbols
|
||||
Rayleigh range for given beam waist and wavelength.
|
||||
|
||||
lens(f)
|
||||
Parameters
|
||||
----------
|
||||
f : float or integer or sympy symbol in meters
|
||||
Thin lens focal length in meters
|
||||
|
||||
Returns
|
||||
-------
|
||||
mat : 2x2 matrix
|
||||
[
|
||||
[ 1, 0],
|
||||
[-1/f, 1]
|
||||
]
|
||||
|
||||
mult(mat1, *argv)
|
||||
Parameters
|
||||
----------
|
||||
mat1 : 2x2 ABCD matrix
|
||||
Last matrix light interacts with.
|
||||
*argv : 2x2 ABCD matrices
|
||||
From left to right, the matrices should be entered such that the leftmost matrix interacts
|
||||
with light temporally after the rightmost matrix.
|
||||
|
||||
Returns
|
||||
-------
|
||||
Mat : 2x2 matrix
|
||||
The ABCd matrix describing the whole optical system.
|
||||
|
||||
plot(func, var, rang=array([0. , 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1 ,
|
||||
0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2 , 0.21,
|
||||
0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3 , 0.31, 0.32,
|
||||
0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4 , 0.41, 0.42, 0.43,
|
||||
0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5 , 0.51, 0.52, 0.53, 0.54,
|
||||
0.55, 0.56, 0.57, 0.58, 0.59, 0.6 , 0.61, 0.62, 0.63, 0.64, 0.65,
|
||||
0.66, 0.67, 0.68, 0.69, 0.7 , 0.71, 0.72, 0.73, 0.74, 0.75, 0.76,
|
||||
0.77, 0.78, 0.79, 0.8 , 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87,
|
||||
0.88, 0.89, 0.9 , 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98,
|
||||
0.99, 1. , 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09,
|
||||
1.1 , 1.11, 1.12, 1.13, 1.14, 1.15, 1.16, 1.17, 1.18, 1.19, 1.2 ,
|
||||
1.21, 1.22, 1.23, 1.24, 1.25, 1.26, 1.27, 1.28, 1.29, 1.3 , 1.31,
|
||||
1.32, 1.33, 1.34, 1.35, 1.36, 1.37, 1.38, 1.39, 1.4 , 1.41, 1.42,
|
||||
1.43, 1.44, 1.45, 1.46, 1.47, 1.48, 1.49, 1.5 , 1.51, 1.52, 1.53,
|
||||
1.54, 1.55, 1.56, 1.57, 1.58, 1.59, 1.6 , 1.61, 1.62, 1.63, 1.64,
|
||||
1.65, 1.66, 1.67, 1.68, 1.69, 1.7 , 1.71, 1.72, 1.73, 1.74, 1.75,
|
||||
1.76, 1.77, 1.78, 1.79, 1.8 , 1.81, 1.82, 1.83, 1.84, 1.85, 1.86,
|
||||
1.87, 1.88, 1.89, 1.9 , 1.91, 1.92, 1.93, 1.94, 1.95, 1.96, 1.97,
|
||||
1.98, 1.99, 2. , 2.01, 2.02, 2.03, 2.04, 2.05, 2.06, 2.07, 2.08,
|
||||
2.09, 2.1 , 2.11, 2.12, 2.13, 2.14, 2.15, 2.16, 2.17, 2.18, 2.19,
|
||||
2.2 , 2.21, 2.22, 2.23, 2.24, 2.25, 2.26, 2.27, 2.28, 2.29, 2.3 ,
|
||||
2.31, 2.32, 2.33, 2.34, 2.35, 2.36, 2.37, 2.38, 2.39, 2.4 , 2.41,
|
||||
2.42, 2.43, 2.44, 2.45, 2.46, 2.47, 2.48, 2.49, 2.5 , 2.51, 2.52,
|
||||
2.53, 2.54, 2.55, 2.56, 2.57, 2.58, 2.59, 2.6 , 2.61, 2.62, 2.63,
|
||||
2.64, 2.65, 2.66, 2.67, 2.68, 2.69, 2.7 , 2.71, 2.72, 2.73, 2.74,
|
||||
2.75, 2.76, 2.77, 2.78, 2.79, 2.8 , 2.81, 2.82, 2.83, 2.84, 2.85,
|
||||
2.86, 2.87, 2.88, 2.89, 2.9 , 2.91, 2.92, 2.93, 2.94, 2.95, 2.96,
|
||||
2.97, 2.98, 2.99]))
|
||||
Parameters
|
||||
----------
|
||||
func : Sympy function of one variable
|
||||
Sympy function defining the beam width after the last optical element.
|
||||
var : sympy variable
|
||||
Variable in func that will be plotted.
|
||||
rang : numpy array
|
||||
Array of the values along the optical axis to be plotted
|
||||
|
||||
Returns
|
||||
-------
|
||||
plot : matplotlib graph
|
||||
Graph of the beam width of var
|
||||
|
||||
prop(d)
|
||||
Parameters
|
||||
----------
|
||||
d : float or integer or sympy symbol
|
||||
Distance light is propagating along the z-axis.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mat: 2x2 matrix
|
||||
[
|
||||
[1, d],
|
||||
[0, 1]
|
||||
]
|
||||
|
||||
q1_func(z, w0, lam, mat)
|
||||
Parameters
|
||||
----------
|
||||
z : float, int, symbol
|
||||
Position of the beam waist in meters.
|
||||
w0 : float, int, symbol
|
||||
Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
|
||||
lam : float, int, symbol
|
||||
Wavelength of light in meters.
|
||||
mat : float, int, symbol
|
||||
The ABCD 2x2 matrix describing the optical system.
|
||||
|
||||
Returns
|
||||
-------
|
||||
z: float, int, symbol
|
||||
Position of the beam waist after the optical system
|
||||
zr: float, int, symbol
|
||||
Rayleigh range of the beam after the optical system
|
||||
|
||||
q1_inv_func(z, w0, lam, mat)
|
||||
Parameters
|
||||
----------
|
||||
z : float, int, symbol
|
||||
Position of the beam waist in meters.
|
||||
w0 : float, int, symbol
|
||||
Radial waist size in meters (of the embedded Gaussian, i.e. W0/M).
|
||||
lam : float, int, symbol
|
||||
Wavelength of light in meters.
|
||||
mat : float, int, symbol
|
||||
The ABCD 2x2 matrix describing the optical system.
|
||||
|
||||
Returns
|
||||
-------
|
||||
R : float, int, symbol
|
||||
Radius of curvature of the wavefront in meters.
|
||||
w : float, int, symbol
|
||||
Radius of the beam in meters.
|
||||
|
||||
ray(y, theta)
|
||||
Parameters
|
||||
----------
|
||||
y : float or integer or sympy symbol in meters
|
||||
The vertical height of a ray.
|
||||
theta : float or integer in radians
|
||||
The angle of divergence of the ray.
|
||||
|
||||
Returns
|
||||
-------
|
||||
mat : 2x1 matrix
|
||||
[
|
||||
[y],
|
||||
[teta]
|
||||
]
|
||||
|
||||
DATA
|
||||
oo = oo
|
||||
|
||||
FILE
|
||||
c:\users\wrighta\documents\beamprop\beamprop_script.py
|
||||
|
||||
|
||||
|
||||
|
||||
## Let's first see how we define a beam and how we can visualize it propagating.
|
||||
|
||||
### A Gaussian beam can be defined by it's (radial) waist, $w_0$, it's Rayleigh range, $z_R = \frac{\pi * w_0^2}{\lambda}$, and the location of its waist, $z_0$.
|
||||
|
||||
|
||||
```python
|
||||
w0 = 1E-3 # 1mm beam waist
|
||||
lam = 355E-9 # wavelength of 355 nm (UV)
|
||||
zR = bs.Zr(w0, lam) # Rayleigh range in m
|
||||
z0 = 0 # location of waist in m
|
||||
```
|
||||
|
||||
### We now want to define our "optical system" using matrices. For this first example, we will just use a free space propagation matrix, and let the beam propagate a distance $d$ which we will define using a symbol.
|
||||
|
||||
|
||||
```python
|
||||
d = sym.symbols('d')
|
||||
M = bs.prop(d)
|
||||
```
|
||||
|
||||
### We now use the bs script to do all the ABCD and q-parameter math, and return the waist and radius of curvature functions
|
||||
|
||||
|
||||
```python
|
||||
R, w = bs.q1_inv_func(0, w0, lam, M)
|
||||
```
|
||||
|
||||
|
||||
```python
|
||||
print('w = {}'.format(w))
|
||||
```
|
||||
|
||||
w = 0.001*(0.0127690021685256*d**2 + 1)**0.5
|
||||
|
||||
|
||||
### And as simple as that, we have a function for our waist. Let's plot it and see what it looks like
|
||||
|
||||
|
||||
```python
|
||||
bs.plot(w, d, rang = np.arange(0,10))
|
||||
```
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
### Let's show what happens when a beam travels through a lens. We use the "mult" function to multiply multiple ABCD matrices together.
|
||||
|
||||
|
||||
```python
|
||||
w0 = 1E-3 # 1mm beam waist
|
||||
lam = 355E-9 # wavelength of 355 nm (UV)
|
||||
zR = bs.Zr(w0, lam) # Rayleigh range in m
|
||||
z0 = 0 # location of waist in m
|
||||
|
||||
d = sym.symbols('d')
|
||||
M = bs.mult(bs.prop(d), bs.lens(.5), bs.prop(1))
|
||||
|
||||
R, w = bs.q1_inv_func(0, w0, lam, M)
|
||||
|
||||
bs.plot(w, d, rang = np.arange(0,1,.01))
|
||||
```
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
### Lets look at how to expand and collimate a beam with a two lens system
|
||||
|
||||
|
||||
```python
|
||||
w0 = 1E-3 # 1mm beam waist
|
||||
lam = 355E-9 # wavelength of 355 nm (UV)
|
||||
zR = bs.Zr(w0, lam) # Rayleigh range in m
|
||||
z0 = 0 # location of waist in m
|
||||
|
||||
d1, d2, d3, f1, f2 = sym.symbols('d1 d2 d3 f1 f2')
|
||||
|
||||
M = bs.mult(bs.prop(d3),bs.lens(f2),bs.prop(d2), bs.lens(f1), bs.prop(d1))
|
||||
|
||||
R, w = bs.q1_inv_func(0, w0, lam, M)
|
||||
```
|
||||
|
||||
### For example, lets say the beam travels 1 m before hitting the first lens, and we want the beam to be 5x w0 after coming out of the second lens. We substitute d1 for 1 meter, since the beam propagates 1 meter, and we substitute d3 for 0, since we only care about the beam size right at the second lens. This gives us a relation between f1 and d2 (the separation between the lenses).
|
||||
|
||||
|
||||
```python
|
||||
w = w.subs(d1,1).subs(d3,0)
|
||||
f1_eq = sym.solve(w - 5*w0, f1)[0]
|
||||
print('f = {}'.format(f1_eq))
|
||||
```
|
||||
|
||||
f = 1.0084642216545e+15*d2*(1.12051580183833e+27*d2 - 4.41556446152598e+29*sqrt(1 - 0.000504320418227052*d2**2) + 8.88733242867719e+28)/(1.13000009595246e+42*d2**2 + 2.26000019190491e+42*d2 - 2.12276362486616e+45)
|
||||
|
||||
|
||||
#### Suppose we wanted the distance between the lenses to be 1 meter, we could find what f1 we need.
|
||||
|
||||
|
||||
```python
|
||||
print('f1 = {:.2f} m, for a lens separation of 1 meter'.format(f1_eq.subs(d2, 1)))
|
||||
```
|
||||
|
||||
f1 = 0.17 m, for a lens separation of 1 meter
|
||||
|
||||
|
||||
### Now we need to collimate the beam. Lets still assume the beam propagates 1 m, and f1 = .17 m.
|
||||
|
||||
There are a couple different ways to think about collimation. One is that the beam size doesn't change over a long distance. The other is that the radius of curvature is infinite (i.e. a plane wave). Lets us the latter interpretation. Thus, we want to find the focal length f2 that makes R infinite, or that makes 1/R =0.
|
||||
|
||||
|
||||
```python
|
||||
R_coll = R.subs(d1,1).subs(d2,1).subs(f1,.17).subs(d3,0)
|
||||
f2_coll = sym.solve(1/R_coll,f2)[0]
|
||||
print('f2 = {:.2f}, for a collimated beam, 5x the original waist, after propagating 1m to the first lens of f1 = .17m, and propagating another 1m to the second lens'.format(f2_coll))
|
||||
```
|
||||
|
||||
f2 = 0.83, for a collimated beam, 5x the original waist, after propagating 1m to the first lens of f1 = .17m, and propagating another 1m to the second lens
|
||||
|
||||
|
||||
### Lets plot the beam profile after the second lens, and see if it is collimated.
|
||||
|
||||
|
||||
```python
|
||||
M = bs.mult(bs.prop(d3),bs.lens(.83),bs.prop(1), bs.lens(.17), bs.prop(1))
|
||||
|
||||
R, w = bs.q1_inv_func(0, w0, lam, M)
|
||||
|
||||
bs.plot(w,d3)
|
||||
```
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
### Looks very collimated. Lets check the beam size (to make sure its 5* w0) and check the collimation
|
||||
|
||||
|
||||
```python
|
||||
expansion_factor = w.subs(d3,0)/ w0
|
||||
print('beam is w = {:.2f} x w0'.format(expansion_factor))
|
||||
```
|
||||
|
||||
beam is w = 4.90 x w0
|
||||
|
||||
|
||||
|
||||
```python
|
||||
beam_size_change = (w.subs(d3,10) - w.subs(d3,0)) / w.subs(d3,0) * 100
|
||||
print('Over 10 m after second lens, beam changes by {:.0f}%'.format(beam_size_change))
|
||||
```
|
||||
|
||||
Over 10 m after second lens, beam changes by 1%
|
||||
|
||||
|
||||
|
||||
```python
|
||||
|
||||
```
|
||||
552
TimeSeriesAnalyzer/TimeSeriesAnalysis.ipynb
Normal file
552
TimeSeriesAnalyzer/TimeSeriesAnalysis.ipynb
Normal file
File diff suppressed because one or more lines are too long
333
TimeSeriesAnalyzer/TimeSeriesAnalyzer.py
Normal file
333
TimeSeriesAnalyzer/TimeSeriesAnalyzer.py
Normal file
@ -0,0 +1,333 @@
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.gridspec as gridspec
|
||||
from scipy.signal import find_peaks, resample, detrend
|
||||
import csv
|
||||
|
||||
"""
|
||||
NOTES
|
||||
|
||||
When you compute the Fourier transform of a signal, you obtain a set of complex-valued coefficients corresponding to different frequency components present in the signal. These coefficients represent the amplitude and phase of sinusoidal components at specific frequencies.
|
||||
|
||||
The frequency range covered by the Fourier transform output is divided into discrete frequency bins, each representing a specific frequency component. The width of these bins depends on the sampling rate and the length of the input signal. In a typical implementation, the frequency bins are evenly spaced.
|
||||
|
||||
In practical terms, the power spectrum bins correspond to the frequency components at which the power spectral density (PSD) or magnitude squared of the Fourier coefficients are evaluated. These bins are used to represent the distribution of signal power across different frequency components, providing insights into the frequency content of the signal.
|
||||
|
||||
The bin width, also called the resolution bandwidth, is simply sampling rate / Total number of samples = (1/dt)/N
|
||||
|
||||
N can be re-written in terms of t, dt. Putting this in to the expression for Sxx = 2 * |fft|^2/(Resolution Bandwidth * Noise Power Bandwidth) We get the expression used here. Note that the fft computed above is scaled by N which results eventually in the factor dt^2/T.
|
||||
|
||||
The noise power bandwidth is typically 1 Hz if no windowing/tapering function is used.
|
||||
|
||||
|
||||
Compute the broadband noise level in Vrms2/Hz by summing all the power spectrum bins, excluding any peaks and the DC component, and dividing the sum by the equivalent noise bandwidth of the window
|
||||
|
||||
The equivalent noise bandwidth (ENBW) of a window in the context of a power spectrum refers to a measure of the effective bandwidth of the window function applied to the signal before taking the Fourier transform.
|
||||
|
||||
When you compute the power spectrum of a signal using a windowing function (e.g., Hamming window, Hann window, etc.), the window modifies the original signal by tapering its edges. This tapering reduces the spectral leakage and improves frequency resolution but also introduces a smoothing effect, which can affect the estimation of the signal's power at different frequencies.
|
||||
|
||||
The equivalent noise bandwidth provides a way to quantify the effective bandwidth of the window function in terms of its impact on noise power. It represents the width of a rectangular filter that would have the same noise power as the windowed signal.
|
||||
|
||||
In practical terms, when calculating the power spectrum of a signal using a window, the ENBW is used to adjust the power spectrum to account for the smoothing effect of the window. Dividing the sum of the power spectrum bins by the ENBW yields an estimate of the noise power per unit frequency bandwidth.
|
||||
|
||||
ENBW is often used in the context of noise measurements or signal processing applications where accurate estimation of noise power is important. It helps ensure that the power spectrum accurately reflects the true power distribution of the signal, accounting for the effects of windowing.
|
||||
|
||||
|
||||
The noise floor refers to the minimum level of signal that can be reliably distinguished from the background noise. It represents the lowest amplitude of a signal that can be detected or measured with reasonable accuracy.
|
||||
|
||||
The noise floor is often defined as the RMS (Root Mean Square) value of the background noise in a given frequency band.
|
||||
|
||||
To compute the noise floor value from the power spectral density (PSD) values, you typically need to analyze the portion of the PSD that corresponds to the background noise.
|
||||
|
||||
"""
|
||||
|
||||
def compute_autocorrelation(data):
|
||||
|
||||
print("Calculating autocorrelation...")
|
||||
yunbiased = data-np.mean(data)
|
||||
ynorm = np.sum(yunbiased**2)
|
||||
acor = np.correlate(yunbiased, yunbiased, "same")/ynorm
|
||||
return acor
|
||||
|
||||
def pre_process_data(data, new_sampling_rate):
|
||||
# Resample the time series
|
||||
# Define the new sampling rate and calculate the new time values
|
||||
n_points = int((len(data[:, 0]) * (data[:, 0][1] - data[:, 0][0])) * new_sampling_rate)
|
||||
t = np.linspace(data[:, 0][0], data[:, 0][-1], n_points)
|
||||
x_array = resample(data[:, 1], n_points)
|
||||
x = detrend(x_array - x_array.mean())
|
||||
|
||||
return np.column_stack((t, x))
|
||||
|
||||
def smoothen_data(data, window_size):
|
||||
# Smooth data by doing a moving average
|
||||
return np.convolve(data, np.ones(window_size, dtype=int)/window_size, mode='valid')
|
||||
|
||||
def compute_psd(data, new_sampling_rate, window_size = 21):
|
||||
"""
|
||||
A power spectral density (PSD) takes the amplitude of the FFT, multiplies it by its complex conjugate and normalizes it to the frequency bin width.
|
||||
|
||||
"""
|
||||
|
||||
processed_data = pre_process_data(data, new_sampling_rate)
|
||||
t, x = processed_data[:, 0], processed_data[:, 1]
|
||||
|
||||
dt = t[1] - t[0] # Define the sampling interval
|
||||
N = x.shape[0] # Define the total number of data points
|
||||
T = N * dt # Define the total duration of the data
|
||||
|
||||
# Calculate fft
|
||||
print("Calculating power spectrum...")
|
||||
fft_ts = np.fft.fft(x) # Compute Fourier transform of x
|
||||
Sxx = 2 * (dt ** 2 / T) * (fft_ts * fft_ts.conj()) # Compute spectrum
|
||||
Sxx = Sxx[:int(len(x) / 2)] # Ignore negative frequencies, we have accounted for this by the scaling factor of 2 in the previous step
|
||||
|
||||
return processed_data, smoothen_data(Sxx.real, window_size)
|
||||
|
||||
def compute_RIN(time, voltages, Sxx_smooth):
|
||||
|
||||
dt = time[1] - time[0] # Define the sampling interval
|
||||
N = voltages.shape[0] # Define the total number of data points
|
||||
T = N * dt # Define the total duration of the data
|
||||
df = 1 / T.max()
|
||||
|
||||
# Compute the average power
|
||||
average_P = np.mean(np.squared(voltages))
|
||||
|
||||
# Calculate the RIN
|
||||
RIN_Sxx_smooth = 10 * np.log10(Sxx_smooth / (average_P * df))
|
||||
|
||||
return RIN_Sxx_smooth
|
||||
|
||||
def find_noise_peaks(psd, faxis, freq_range, threshold):
|
||||
"""
|
||||
Compute the peak power in the specified frequency range.
|
||||
|
||||
Parameters:
|
||||
psd_values: array-like
|
||||
Power spectral density values.
|
||||
faxis: array-like
|
||||
Frequencies corresponding to the PSD values.
|
||||
freq_range: tuple
|
||||
Tuple containing the start and end frequencies of the range of interest.
|
||||
threshold: scalar
|
||||
Threshold for peak heights
|
||||
|
||||
Returns:
|
||||
float: Peak power in the specified frequency range.
|
||||
"""
|
||||
start_freq, end_freq = freq_range
|
||||
idx_start = np.argmax(faxis >= start_freq)
|
||||
idx_end = np.argmax(faxis >= end_freq)
|
||||
sliced_psd = psd[idx_start:idx_end]
|
||||
sliced_faxis = faxis[idx_start:idx_end]
|
||||
|
||||
peak_indices, _ = find_peaks(sliced_psd, height=threshold)
|
||||
peak_powers = 10 * np.log10(sliced_psd[peak_indices])
|
||||
peak_frequencies = np.around(sliced_faxis[peak_indices], 2)
|
||||
|
||||
return peak_powers, peak_frequencies
|
||||
|
||||
def compute_noise_level(psd, resolution_bandwidth, exclude_peaks=False, faxis=None, freq_range=None, threshold=None):
|
||||
"""
|
||||
Compute the noise level from a power spectral density (PSD).
|
||||
|
||||
Parameters:
|
||||
psd: array-like
|
||||
One-sided power spectral density.
|
||||
resolution_bandwidth: float
|
||||
Bin width
|
||||
Returns:
|
||||
float: Noise level (Vrms^2).
|
||||
"""
|
||||
|
||||
noise_level = None
|
||||
# Exclude peaks from the sum
|
||||
if exclude_peaks and threshold is not None:
|
||||
|
||||
threshold = 10**(threshold/10)
|
||||
|
||||
if freq_range is None:
|
||||
peak_indices, _ = find_peaks(psd, height=threshold)
|
||||
noise_level = resolution_bandwidth * np.sum([psd[i] for i in range(len(psd)) if i not in peak_indices])
|
||||
else:
|
||||
start_freq, end_freq = freq_range
|
||||
idx_start = np.argmax(faxis >= start_freq)
|
||||
idx_end = np.argmax(faxis >= end_freq)
|
||||
sliced_psd = psd[idx_start:idx_end]
|
||||
peak_indices, _ = find_peaks(sliced_psd, height=threshold)
|
||||
noise_level = resolution_bandwidth * np.sum([sliced_psd[i] for i in range(len(sliced_psd)) if i not in peak_indices])
|
||||
else:
|
||||
|
||||
if freq_range is None:
|
||||
noise_level = resolution_bandwidth * np.sum([psd[i] for i in range(len(psd))])
|
||||
else:
|
||||
start_freq, end_freq = freq_range
|
||||
idx_start = np.argmax(faxis >= start_freq)
|
||||
idx_end = np.argmax(faxis >= end_freq)
|
||||
sliced_psd = psd[idx_start:idx_end]
|
||||
noise_level = resolution_bandwidth * np.sum([sliced_psd[i] for i in range(len(sliced_psd))])
|
||||
|
||||
return noise_level
|
||||
|
||||
def extract_data(filepath):
|
||||
|
||||
# Open the CSV file
|
||||
with open(filepath, newline='') as csvfile:
|
||||
# Skip the first line (header)
|
||||
next(csvfile)
|
||||
|
||||
# Read the CSV file using csv.reader
|
||||
reader = csv.reader(csvfile)
|
||||
|
||||
# Read the headers from the second line
|
||||
next_header = next(reader)
|
||||
string_number = next_header[-1]
|
||||
|
||||
try:
|
||||
time_step = int(string_number)
|
||||
except ValueError:
|
||||
try:
|
||||
time_step = float(string_number)
|
||||
except ValueError:
|
||||
print("The string does not represent a valid number.")
|
||||
|
||||
# Initialize lists to store the first and second values
|
||||
first_column = []
|
||||
second_column = []
|
||||
|
||||
# Iterate over each row in the CSV file
|
||||
for row in reader:
|
||||
# Extract the first and second values from the row and convert to float
|
||||
first_value = float(row[0])
|
||||
second_value = float(row[1])
|
||||
|
||||
# Append the values to their respective lists
|
||||
first_column.append(first_value)
|
||||
second_column.append(second_value)
|
||||
|
||||
# Convert the lists into numpy arrays
|
||||
time_array = np.arange(0, len(first_column)*time_step, time_step)
|
||||
voltage_array = np.array(second_column)
|
||||
|
||||
# Stack the arrays horizontally to form a single 2D array
|
||||
data_array = np.column_stack((time_array, voltage_array))
|
||||
|
||||
return data_array
|
||||
|
||||
def plot_analysis(data, data_bkg, Sxx, Sxx_bkg, data_str, bkg_str, peak_find_threshold, window_size = 21, plot_only_psd = True):
|
||||
|
||||
|
||||
time, voltages = data[:, 0], data[:, 1]
|
||||
time_bkg, voltages_bkg = data_bkg[:, 0], data_bkg[:, 1]
|
||||
|
||||
dt = time[1] - time[0] # Define the sampling interval
|
||||
N = voltages.shape[0] # Define the total number of data points
|
||||
T = N * dt # Define the total duration of the data
|
||||
df = 1 / T.max()
|
||||
|
||||
fNQ = 1 / dt / 2 # Determine Nyquist frequency
|
||||
faxis = smoothen_data(np.linspace(0,fNQ,N//2), window_size) # Construct frequency axis
|
||||
|
||||
""" Noise levels in units of Vrms^2/Hz"""
|
||||
# resolution_bandwidth = (1/dt)/N
|
||||
# broadband_noise_level = compute_noise_level(Sxx, resolution_bandwidth) # Integrates across PSD from DC to Nyquist frequency, gives result in in units of Vrms^2/Hz
|
||||
# noise_floor = np.mean(Sxx_bkg)
|
||||
|
||||
freq_range = (50, max(faxis))
|
||||
threshold = 10**(peak_find_threshold/10)
|
||||
peak_powers, peak_frequencies = find_noise_peaks(Sxx, faxis, freq_range, threshold)
|
||||
|
||||
if plot_only_psd:
|
||||
|
||||
plt.figure(figsize=(12, 8))
|
||||
|
||||
# Plot Power Spectrum in dB
|
||||
plt.semilogx(faxis, 10 * np.log10(Sxx_bkg), color='orange', linewidth=0.5, label = bkg_str)
|
||||
plt.semilogx(faxis, 10 * np.log10(Sxx), color='green', linewidth=2, label = data_str)
|
||||
|
||||
# plt.axhline(y=10 * np.log10(broadband_noise_level), color='red', linewidth=2, linestyle='--', label=f'Broadband cumulative noise level: {10 * np.log10(broadband_noise_level):.1f} dB')
|
||||
# plt.axhline(y=10 * np.log10(noise_floor), color='blue', linewidth=2, linestyle='--', label=f'Broadband noise floor: {10 * np.log10(noise_floor):.1f} dB')
|
||||
|
||||
plt.plot(peak_frequencies, peak_powers, 'o', markerfacecolor='none', markeredgecolor='r', markersize=10) # Plot power against frequency as hollow circles
|
||||
for freq, power in zip(peak_frequencies, peak_powers):
|
||||
plt.text(freq, power, str(freq)+' Hz', verticalalignment='bottom', horizontalalignment='right') # Add text next to each circle indicating the frequency
|
||||
|
||||
plt.grid(True, which="both", linestyle='-', linewidth=0.5, color='gray') # Thin lines for non-decade grid
|
||||
plt.grid(True, which="both", linestyle=':', linewidth=1, color='gray', axis='x') # Thick lines for decade grid
|
||||
|
||||
# Calculate the x-axis values for multiples of 10
|
||||
x_multiples_of_10 = [10**i for i in range(int(np.log10(min(faxis[faxis > 0]))), int(np.log10(max(faxis[faxis > 0]))) + 1)]
|
||||
# Add thick lines for multiples of 10
|
||||
for val in x_multiples_of_10:
|
||||
plt.axvline(x=val, color='black', linestyle='-', linewidth=2) # Thick lines for multiples of 10
|
||||
|
||||
f_sig_idx = np.argmax(Sxx)
|
||||
# SNR_f = 10 * np.log10(Sxx[f_sig_idx] / np.sum(np.delete(Sxx, f_sig_idx)))
|
||||
# SNR_f = 10 * np.log10(Sxx[f_sig_idx] / noise_floor)
|
||||
|
||||
plt.xlim([min(faxis), max(faxis)])
|
||||
# plt.ylim([-100, 10])
|
||||
plt.legend(loc = 3, fontsize=12)
|
||||
plt.xlabel('Frequency [Hz]', fontsize=14)
|
||||
plt.ylabel('Power Spectral Density [dB/Hz]', fontsize=14)
|
||||
# plt.title('SNR= %.2f dB' % (SNR_f), fontsize=14)
|
||||
|
||||
# Adjust layout
|
||||
plt.tight_layout()
|
||||
|
||||
# Show plot
|
||||
plt.show()
|
||||
|
||||
else:
|
||||
# Create subplots
|
||||
plt.figure(figsize=(12, 8))
|
||||
gs = gridspec.GridSpec(2, 3, width_ratios=[1, 1, 1], height_ratios=[1, 1])
|
||||
|
||||
# Plot 1: Time vs Voltage
|
||||
axs1 = plt.subplot(gs[0, 0:])
|
||||
axs1.plot(time_bkg, voltages_bkg, marker='o', color='orange', linewidth=0.5, ms=1, label = bkg_str)
|
||||
axs1.plot(time, voltages, marker='o', color='green', linewidth=0.5, ms=1, label = data_str)
|
||||
axs1.set_ylim([-0.5, 0.5])
|
||||
axs1.set_xlabel('Time (s)', fontsize=14)
|
||||
axs1.set_ylabel('Voltage (V)', fontsize=14)
|
||||
axs1.legend(loc = 1, fontsize=12)
|
||||
axs1.autoscale(tight=True)
|
||||
axs1.grid(True)
|
||||
|
||||
# Plot 2: Power Spectrum in dB
|
||||
axs2 = plt.subplot(gs[1, 0:])
|
||||
axs2.semilogx(faxis, 10 * np.log10(Sxx_bkg), color='orange', linewidth=0.5, label = bkg_str)
|
||||
axs2.semilogx(faxis, 10 * np.log10(Sxx), color='green', linewidth=2, label = data_str)
|
||||
|
||||
# axs2.axhline(y=10 * np.log10(broadband_noise_level), color='red', linewidth=2, linestyle='--', label=f'Broadband cumulative noise level: {10 * np.log10(broadband_noise_level):.1f} dB')
|
||||
# axs2.axhline(y=10 * np.log10(noise_floor), color='blue', linewidth=2, linestyle='--', label=f'Broadband noise floor: {10 * np.log10(noise_floor):.1f} dB')
|
||||
|
||||
axs2.plot(peak_frequencies, peak_powers, 'o', markerfacecolor='none', markeredgecolor='r', markersize=10) # Plot power against frequency as hollow circles
|
||||
for freq, power in zip(peak_frequencies, peak_powers):
|
||||
axs2.text(freq, power, str(freq)+' Hz', verticalalignment='bottom', horizontalalignment='right') # Add text next to each circle indicating the frequency
|
||||
|
||||
axs2.grid(True, which="both", linestyle='-', linewidth=0.5, color='gray') # Thin lines for non-decade grid
|
||||
axs2.grid(True, which="both", linestyle=':', linewidth=1, color='gray', axis='x') # Thick lines for decade grid
|
||||
# Calculate the x-axis values for multiples of 10
|
||||
x_multiples_of_10 = [10**i for i in range(int(np.log10(min(faxis[faxis > 0]))), int(np.log10(max(faxis[faxis > 0]))) + 1)]
|
||||
# Add thick lines for multiples of 10
|
||||
for val in x_multiples_of_10:
|
||||
axs2.axvline(x=val, color='black', linestyle='-', linewidth=2) # Thick lines for multiples of 10
|
||||
|
||||
f_sig_idx = np.argmax(Sxx)
|
||||
# SNR_f = 10 * np.log10(Sxx[f_sig_idx] / np.sum(np.delete(Sxx, f_sig_idx)))
|
||||
# SNR_f = 10 * np.log10(Sxx[f_sig_idx] / noise_floor)
|
||||
|
||||
axs2.set_xlim([min(faxis), max(faxis)])
|
||||
# axs2.set_ylim([-100, 10])
|
||||
axs2.legend(loc = 3, fontsize=12)
|
||||
axs2.set_xlabel('Frequency [Hz]', fontsize=14)
|
||||
axs2.set_ylabel('Power Spectral Density [dB/Hz]', fontsize=14)
|
||||
# axs2.set_title('SNR= %.2f dB' % (SNR_f), fontsize=14)
|
||||
|
||||
# Adjust layout
|
||||
plt.tight_layout()
|
||||
|
||||
# Show plot
|
||||
plt.show()
|
||||
@ -0,0 +1,50 @@
|
||||
%read CSV file
|
||||
filename = 'C:\Users\Karthik\Documents\Git Repos\ULE Cavity Characterisitics\ReflectivityCurve_ULECavity.csv';
|
||||
delimiter = ',';
|
||||
startRow = 1;
|
||||
formatSpec = '%f%f';
|
||||
fileID = fopen(filename,'r');
|
||||
dataset = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'EmptyValue', NaN, 'ReturnOnError', false, 'EndOfLine', '\r\n');
|
||||
fclose(fileID);
|
||||
|
||||
wavelengths = dataset{1};
|
||||
reflectance = dataset{2};
|
||||
|
||||
figure(1);
|
||||
clf;
|
||||
|
||||
xq = 549:1:1100;
|
||||
vq1 = interp1(wavelengths(30:end),reflectance(30:end),xq);
|
||||
plot(wavelengths,reflectance,'o',xq, vq1,':.');
|
||||
%scatter(wavelengths, reflectance, 'LineWidth', 2)
|
||||
|
||||
R_626 = vq1(xq==626);
|
||||
R_842 = vq1(xq==842);
|
||||
F_626 = pi * sqrt(R_626) / (1 - R_626);
|
||||
F_842 = pi * sqrt(R_842) / (1 - R_842);
|
||||
hold on
|
||||
plot(626, R_626, 'o', 'MarkerSize', 15, 'LineWidth', 3)
|
||||
line([626 626], [0 2],'Color','red','LineStyle','--')
|
||||
line([500 1100], [R_626 R_626],'Color','red','LineStyle','--')
|
||||
text(630, R_626 - 0.1, sprintf('626, %.3f', R_626))
|
||||
text(630, R_626 - 0.15, sprintf('F = %.3f', F_626))
|
||||
%annotation('arrow', [0.293 0.293], [0.11 0.3]);
|
||||
|
||||
plot(842, R_842, 'o', 'MarkerSize', 15, 'LineWidth', 3)
|
||||
line([842 842], [0 2],'Color','red','LineStyle','--')
|
||||
line([500 1100], [R_842 R_842],'Color','red','LineStyle','--')
|
||||
text(846, R_842 - 0.1, sprintf('842, %.3f', R_842))
|
||||
text(846, R_842 - 0.15, sprintf('F = %.3f', F_842))
|
||||
%annotation('arrow', [0.571 0.571], [0.11 0.2]);
|
||||
|
||||
hold off
|
||||
|
||||
|
||||
hXLabel = xlabel('Wavelength (nm)');
|
||||
hYLabel = ylabel('?? Mirror Reflectivity R(%) ??');
|
||||
hTitle = sgtitle('??Reflectivity Curve of the ULE Cavity Mirrors??');
|
||||
set([hXLabel, hYLabel] , ...
|
||||
'FontSize' , 14 );
|
||||
set( hTitle , ...
|
||||
'FontSize' , 18 );
|
||||
grid on
|
||||
@ -0,0 +1,94 @@
|
||||
%read CSV file
|
||||
filename = 'C:\Users\Karthik\Documents\Git Repos\ULE Cavity Characterisitics\TransmissionCurve_ULECavity.csv';
|
||||
delimiter = ',';
|
||||
startRow = 1;
|
||||
formatSpec = '%f%f';
|
||||
fileID = fopen(filename,'r');
|
||||
dataset = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'EmptyValue', NaN, 'ReturnOnError', false, 'EndOfLine', '\r\n');
|
||||
fclose(fileID);
|
||||
|
||||
wavelengths = dataset{1};
|
||||
transmission = dataset{2};
|
||||
|
||||
f_h = Helper.getFigureByTag('CavityCharacteristics');
|
||||
set(groot,'CurrentFigure',f_h);
|
||||
a_h = get(f_h, 'CurrentAxes');
|
||||
if ~isempty(get(a_h, 'Children'))
|
||||
clf(f_h);
|
||||
end
|
||||
f_h.Name = 'Wavelength dependence';
|
||||
f_h.Units = 'pixels';
|
||||
|
||||
set(0,'units','pixels');
|
||||
screensize = get(0,'ScreenSize');
|
||||
f_h.Position = [[screensize(3)/30 screensize(4)/4] 1812 429];
|
||||
clf;
|
||||
|
||||
W = 600:1:850;
|
||||
T = interp1(wavelengths(51:149), transmission(51:149), W);
|
||||
|
||||
subplot(1,3,1)
|
||||
plot(wavelengths,transmission,'o', 'MarkerSize', 5);
|
||||
hold on
|
||||
plot(W, T,':.');
|
||||
|
||||
FSR = 1.5e+9;
|
||||
|
||||
T_626 = T(W==626) * 1e-2;
|
||||
R_626 = 1 - T_626;
|
||||
T_842 = T(W==842) * 1e-2;
|
||||
R_842 = 1 - T_842;
|
||||
F_626 = pi * sqrt(R_626) / (1 - R_626);
|
||||
F_842 = pi * sqrt(R_842) / (1 - R_842);
|
||||
L_626 = FSR / F_626;
|
||||
L_842 = FSR / F_842;
|
||||
|
||||
R = 1 - (T * 1e-2);
|
||||
F = pi .* sqrt(R) ./ (1 - R);
|
||||
L = FSR ./ F;
|
||||
|
||||
% plot(626, T_626 * 1e2, 'o', 'MarkerSize', 15, 'LineWidth', 3)
|
||||
% line([626 626], [0 0.06],'Color',[0.9586 0.7372 0.2537],'LineStyle','--')
|
||||
% line([500 1000], [T_626 T_626] * 1e2,'Color',[0.9586 0.7372 0.2537],'LineStyle','--')
|
||||
text(630, T_626 * 1e2 - 0.002, sprintf('R @ 626 = %.5f', R_626), 'FontSize' , 10)
|
||||
% text(630, T_626 * 1e2 - 0.006, sprintf('F = %.3f', F_626), 'FontSize' , 10)
|
||||
% annotation('arrow', [0.293 0.293], [0.11 0.3]);
|
||||
|
||||
% plot(842, T_842 * 1e2, 'o', 'MarkerSize', 15, 'LineWidth', 3)
|
||||
% line([842 842], [0 0.06],'Color','red','LineStyle','--')
|
||||
% line([500 1000], [T_842 T_842] * 1e2,'Color','red','LineStyle','--')
|
||||
text(825, T_842 * 1e2 + 0.006, sprintf('R @ 842 = %.5f', R_842), 'FontSize' , 10)
|
||||
% text(846, T_842 * 1e2 + 0.002, sprintf('F = %.3f', F_842), 'FontSize' , 10)
|
||||
% annotation('arrow', [0.571 0.571], [0.11 0.2]);
|
||||
|
||||
hold off
|
||||
|
||||
hXLabel = xlabel('Wavelength (nm)');
|
||||
hYLabel = ylabel('Mirror Transmission T(%)');
|
||||
set([hXLabel, hYLabel] , ...
|
||||
'FontSize' , 14 );
|
||||
grid on
|
||||
|
||||
subplot(1,3,2)
|
||||
plot(W, F * 1e-4)
|
||||
text(626, F_626 * 1e-4 + 0.2, sprintf('F @ 626 = %1.f', F_626), 'FontSize' , 10)
|
||||
text(750, F_842 * 1e-4 - 0.5, sprintf('F @ 842 = %1.f', F_842), 'FontSize' , 10)
|
||||
hXLabel = xlabel('Wavelength (nm)');
|
||||
hYLabel = ylabel('Finesse (x 10^{4})');
|
||||
set([hXLabel, hYLabel] , ...
|
||||
'FontSize' , 14 );
|
||||
grid on
|
||||
|
||||
subplot(1,3,3)
|
||||
plot(W, L * 1e-3)
|
||||
text(626, L_626 * 1e-3 - 5, sprintf('L @ 626 = %1.f kHz', L_626* 1e-3), 'FontSize' , 10)
|
||||
text(750, L_842 * 1e-3 + 12, sprintf('L @ 842 = %1.f kHz', L_842* 1e-3), 'FontSize' , 10)
|
||||
hXLabel = xlabel('Wavelength (nm)');
|
||||
hYLabel = ylabel('Linewidth (kHz)');
|
||||
set([hXLabel, hYLabel] , ...
|
||||
'FontSize' , 14 );
|
||||
grid on
|
||||
|
||||
hTitle = sgtitle('SLS ULE Cavity Characterisitics');
|
||||
set( hTitle , ...
|
||||
'FontSize' , 18 );
|
||||
Loading…
Reference in New Issue
Block a user