Renamed GPE Solver to Dipolar Gas Simulator, modified ODT Calculator and Time Series Analyzer.
This commit is contained in:
parent
02bc8cb364
commit
aaf9b3ba53
Binary file not shown.
@ -1,38 +0,0 @@
|
||||
/*
|
||||
* 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]);
|
||||
}
|
||||
}
|
||||
@ -1,46 +0,0 @@
|
||||
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
|
||||
@ -1,15 +0,0 @@
|
||||
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
|
||||
@ -1,10 +0,0 @@
|
||||
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
|
||||
@ -1,6 +0,0 @@
|
||||
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
|
||||
@ -1,18 +0,0 @@
|
||||
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
|
||||
@ -1,191 +0,0 @@
|
||||
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
|
||||
|
||||
|
||||
|
||||
@ -1,92 +0,0 @@
|
||||
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.';
|
||||
@ -1,55 +0,0 @@
|
||||
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');
|
||||
@ -1,148 +0,0 @@
|
||||
% 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
|
||||
|
||||
|
||||
@ -1,820 +0,0 @@
|
||||
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
|
||||
@ -1,50 +0,0 @@
|
||||
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);
|
||||
@ -1,77 +0,0 @@
|
||||
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
|
||||
@ -1,58 +0,0 @@
|
||||
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
|
||||
@ -1,47 +0,0 @@
|
||||
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]);
|
||||
@ -1,94 +0,0 @@
|
||||
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
|
||||
@ -1,18 +0,0 @@
|
||||
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;
|
||||
@ -1,54 +0,0 @@
|
||||
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?
|
||||
@ -1,28 +0,0 @@
|
||||
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;
|
||||
@ -1,35 +0,0 @@
|
||||
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)
|
||||
@ -1,31 +0,0 @@
|
||||
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;
|
||||
@ -1,24 +0,0 @@
|
||||
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);
|
||||
@ -1,99 +0,0 @@
|
||||
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;
|
||||
@ -1,23 +0,0 @@
|
||||
%-% 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);
|
||||
@ -1,32 +0,0 @@
|
||||
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;
|
||||
@ -1,311 +0,0 @@
|
||||
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
|
||||
@ -1,98 +0,0 @@
|
||||
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');
|
||||
@ -3,6 +3,9 @@ from scipy.optimize import curve_fit
|
||||
from scipy import interpolate
|
||||
from astropy import units as u, constants as ac
|
||||
|
||||
from Potentials import *
|
||||
from Helpers import *
|
||||
|
||||
DY_POLARIZABILITY = 184.4 # in a.u, most precise measured value of Dy polarizability
|
||||
DY_MASS = 164*u.u
|
||||
DY_DIPOLE_MOMENT = 9.93 * ac.muB
|
||||
|
||||
File diff suppressed because one or more lines are too long
@ -3,6 +3,10 @@ import matplotlib.pyplot as plt
|
||||
import matplotlib.ticker as mtick
|
||||
from astropy import units as u, constants as ac
|
||||
|
||||
from Calculator import *
|
||||
from Helpers import *
|
||||
from Potentials import *
|
||||
|
||||
#####################################################################
|
||||
# PLOTTING #
|
||||
#####################################################################
|
||||
|
||||
@ -1,6 +1,8 @@
|
||||
import numpy as np
|
||||
from astropy import units as u, constants as ac
|
||||
|
||||
from Helpers import *
|
||||
|
||||
DY_POLARIZABILITY = 184.4 # in a.u, most precise measured value of Dy polarizability
|
||||
DY_MASS = 164*u.u
|
||||
DY_DIPOLE_MOMENT = 9.93 * ac.muB
|
||||
|
||||
@ -270,7 +270,7 @@ def plot_analysis(data, data_bkg, Sxx, Sxx_bkg, data_str, bkg_str, peak_find_thr
|
||||
# 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.ylabel('Power Spectral Density [dBV/Hz]', fontsize=14)
|
||||
# plt.title('SNR= %.2f dB' % (SNR_f), fontsize=14)
|
||||
|
||||
# Adjust layout
|
||||
@ -323,7 +323,7 @@ def plot_analysis(data, data_bkg, Sxx, Sxx_bkg, data_str, bkg_str, peak_find_thr
|
||||
# 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_ylabel('Power Spectral Density [dBV/Hz]', fontsize=14)
|
||||
# axs2.set_title('SNR= %.2f dB' % (SNR_f), fontsize=14)
|
||||
|
||||
# Adjust layout
|
||||
|
||||
File diff suppressed because one or more lines are too long
Loading…
Reference in New Issue
Block a user