Ref.

http://www.joerg-enderlein.de/index.html

———–

« Fast fitting of multi-exponential decay curves »

J. Enderlein, R. Erdmann

Optics Communications 134(1-6), 1997, pp.371-8

—-6 .m

First-3nd: 7 april2008

4nd: May 2007

5-6nd: feb 2004;

—first .m

function [err A z] = lsfit(param, irf, y, p)

% LSFIT(param, irf, y, p) returns the Least-Squares deviation between the data y

% and the computed values.

% LSFIT assumes a function of the form:

%

% y = yoffset + A(1)*convol(irf,exp(-t/tau(1)/(1-exp(-p/tau(1)))) + …

%

% param(1) is the color shift value between irf and y.

% param(2) is the irf offset.

% param(3:…) are the decay times.

% irf is the measured Instrumental Response Function.

% y is the measured fluorescence decay curve.

% p is the time between to laser excitations (in number of TCSPC channels).

n = length(irf);

t = 1:n;

tp = (1:p)’;

c = param(1);

tau = param(2:length(param)); tau = tau(:)’;

x = exp(-(tp-1)*(1./tau))*diag(1./(1-exp(-p./tau)));

irs = (1-c+floor(c))*irf(rem(rem(t-floor(c)-1, n)+n,n)+1) + (c-floor(c))*irf(rem(rem(t-ceil(c)-1, n)+n,n)+1);

z = convol(irs, x);

z = [ones(size(z,1),1) z];

A = z\y;

%A = lsqnonneg(z,y);

z = z*A;

if 1

semilogy(t, irs/max(irs)*max(y), t, y, ‘bo’, t, z);

drawnow

end

err = sum((z-y).^2./abs(z))/(n-length(tau))

–2nd .m

function [c, offset, A, tau, dc, dtau, irs, zz, t, chi] = fluofit(irf, y, p, dt, tau, init)

% The function FLUOFIT performs a fit of a multi-exponential decay curve.

% It is called by:

% [c, offset, A, tau, dc, doffset, dtau, irs, z, t, chi] = fluofit(irf, y, p, dt, tau, init).

% The function arguments are:

% irf = Instrumental Response Function

% y = Fluorescence decay data

% p = Time between laser exciation pulses (in nanoseconds)

% dt = Time width of one TCSPC channel (in nanoseconds)

% tau = Initial guess times

% init = Whether to use a initial guess routine or not

%

% The return parameters are:

% c = Color Shift (time shift of the IRF with respect to the fluorescence curve)

% offset = Offset

% A = Amplitudes of the different decay components

% tau = Decay times of the different decay components

% dc = Color shift error

% doffset = Offset error

% dtau = Decay times error

% irs = IRF, shifted by the value of the colorshift

% zz Fitted fluorecence component curves

% t = time axis

% chi = chi2 value

%

% The program needs the following m-files: simplex.m, lsfit.m, mlfit.m, and convol.m.

% (c) 1996 Jˆrg Enderlein

fitfun = ‘lsfit’;

close all

irf = irf(:);

offset = 0;

y = y(:);

n = length(irf);

if nargin>5

if isempty(init)

init = 1;

end

elseif nargin>4

init = 0;

else

init = 1;

end

if init>0

[cx, tau, c, c] = DistFluofit(irf, y, p, dt, [-3 3]);

cx = cx(:)’;

tmp = cx>0;

t = 1:length(tmp);

t1 = t(tmp(2:end)>tmp(1:end-1)) + 1;

t2 = t(tmp(1:end-1)>tmp(2:end));

if length(t1)==length(t2)+1

t1(end)=[];

end

if length(t2)==length(t1)+1

t2(1)=[];

end

if t1(1)>t2(1)

t1(end)=[];

t2(1)=[];

end

tmp = [];

for j=1:length(t1)

tmp = [tmp cx(t1(j):t2(j))*tau(t1(j):t2(j))/sum(cx(t1(j):t2(j)))];

end

tau = tmp;

else

c = 0;

end

p = p/dt;

tp = (1:p)’;

tau = tau(:)’/dt;

t = 1:length(y);

m = length(tau);

x = exp(-(tp-1)*(1./tau))*diag(1./(1-exp(-p./tau)));

irs = (1-c+floor(c))*irf(rem(rem(t-floor(c)-1, n)+n,n)+1) + (c-floor(c))*irf(rem(rem(t-ceil(c)-1, n)+n,n)+1);

z = convol(irs, x);

z = [ones(size(z,1),1) z];

A = z\y;

%A = lsqnonneg(z,y);

z = z*A;

close all

if init<2

disp(‘Fit = Parameters =’);

param = [c; tau’];

% Decay times and Offset are assumed to be positive.

paramin = [-Inf, zeros(1,length(tau))];

paramax = [Inf, Inf*ones(1,length(tau))];

[param, dparam] = simplex(fitfun, param, paramin, paramax, [], [], irf(:), y(:), p);

c = param(1);

dc = dparam(1);

tau = param(2:length(param))’;

dtau = dparam(2:length(param));

x = exp(-(tp-1)*(1./tau))*diag(1./(1-exp(-p./tau)));

irs = (1-c+floor(c))*irf(rem(rem(t-floor(c)-1, n)+n,n)+1) + (c-floor(c))*irf(rem(rem(t-ceil(c)-1, n)+n,n)+1);

z = convol(irs, x);

z = [ones(size(z,1),1) z];

z = z./(ones(n,1)*sum(z));

A = z\y;

%A = lsqnonneg(z,y);

zz = z.*(ones(size(z,1),1)*A’);

z = z*A;

dtau = dtau;

dc = dt*dc;

else

dtau = 0;

dc = 0;

end

chi = sum((y-z).^2./abs(z))/(n-m);

t = dt*t;

tau = dt*tau’;

c = dt*c;

offset = zz(1,1);

A(1) = [];

if 1

hold off

subplot(‘position’,[0.1 0.4 0.8 0.5])

plot(t,log10(y),t,log10(irf),t,log10(z));

v = axis;

v(1) = min(t);

v(2) = max(t);

axis(v);

xlabel(‘Time in ns’);

ylabel(‘Log Count’);

s = sprintf(‘COF = %3.3f %3.3f’, c, offset);

text(max(t)/2,v(4)-0.05*(v(4)-v(3)),s);

s = [‘AMP = ‘];

for i=1:length(A)

s = [s sprintf(‘%1.3f’,A(i)/sum(A)) ‘ ‘];

end

text(max(t)/2,v(4)-0.12*(v(4)-v(3)),s);

s = [‘TAU = ‘];

for i=1:length(tau)

s = [s sprintf(‘%3.3f’,tau(i)) ‘ ‘];

end

text(max(t)/2,v(4)-0.19*(v(4)-v(3)),s);

subplot(‘position’,[0.1 0.1 0.8 0.2])

plot(t,(y-z)./sqrt(abs(z)));

v = axis;

v(1) = min(t);

v(2) = max(t);

axis(v);

xlabel(‘Time in ns’);

ylabel(‘Residue’);

s = sprintf(‘%3.3f’, chi);

text(max(t)/2,v(4)-0.1*(v(4)-v(3)),[‘\chi^2 = ‘ s]);

set(gcf,’units’,’normalized’,’position’,[0.01 0.05 0.98 0.83])

end

—-3nd .m

function [x, dx, steps] = Simplex(fname, x, xmin, xmax, tol, steps, varargin)

% [x, dx, steps] = Simplex(‘F’, X0, XMIN, XMAX, TOL, STEPS, VARARGIN)

% attempts to return a vector x and its error dx, so that x minimzes the

% function F(x) near the starting vector X0 under the conditions that

% xmin <= x <= xmax.

% TOL is the relative termination tolerance dF/F; (default = 1e-10)

% STEPS is the maximum number of steps; (default = 200*number of parameters).

% The returned value of STEPS is the actual number of performed steps.

% Simplex allows for up to 10 additional arguments for the function F.

% Simplex uses a Nelder-Mead simplex search method.

x = x(:);

if nargin<5

tol = 1e-10;

if nargin<4

xmax = Inf*ones(length(x),1);

if nargin<3

xmin = -Inf*ones(length(x),1);

end

end

elseif isempty(tol)

tol = 1e-5;

end

if nargin<6

steps = [];

end

if isempty(xmin)

xmin = -Inf*ones(size(x));

end

if isempty(xmax)

xmax = Inf*ones(size(x));

end

xmin = xmin(:);

xmax = xmax(:);

xmax(xmax<xmin) = xmin(xmax<xmin);

x(x<xmin) = xmin(x<xmin);

x(x>xmax) = xmax(x>xmax);

xfix = zeros(size(x));

tmp = xmin==xmax;

xfix(tmp) = xmin(tmp);

mask = diag(~tmp);

mask(:, tmp) = [];

x(tmp) = [];

xmin(tmp) = [];

xmax(tmp) = [];

if isa(fname,’function_handle’)

fun = fname;

evalstr = ‘fun’;

else

evalstr = fname;

end

evalstr = [evalstr, ‘(mask*x+xfix’];

if nargin>6

evalstr = [evalstr, ‘,varargin{:}’];

end

evalstr = [evalstr, ‘)’];

n = length(x);

if n==0

x = xfix;

dx = zeros(size(xfix));

steps = 0;

return

end

if isempty(steps)

steps = 200*n;

end

xin = x(:);

%v = 0.9*xin;

v = xin;

v(v<xmin) = xmin(v<xmin);

v(v>xmax) = xmax(v>xmax);

x(:) = v; fv = eval(evalstr);

for j = 1:n

y = xin;

if y(j) ~= 0

y(j) = (1 +.2*rand)*y(j);

else

y(j) = 0.2;

end

if y(j)>=xmax(j)

y(j) = xmax(j);

end

if y(j)<=xmin(j)

y(j) = xmin(j);

end

v = [v y];

x(:) = y; f = eval(evalstr);

fv = [fv f];

end

[fv, j] = sort(fv);

v = v(:,j);

count = n+1;

% Parameter settings for Nelder-Meade

alpha = 1; beta = 1/2; gamma = 2;

% Begin of Nelder-Meade simplex algorithm

while count < steps

if 2*abs(fv(n+1)-fv(1))/(abs(fv(1))+abs(fv(n+1))) <= tol

break

end

% Reflection:

vmean = mean(v(:, 1:n),2);

vr = (1 + alpha)*vmean – alpha*v(:, n+1);

x(:) = vr;

fr = eval(evalstr);

count = count + 1;

vk = vr; fk = fr;

if fr < fv(1) && all(xmin<=vr) && all(vr<=xmax)

% Expansion:

ve = gamma*vr + (1-gamma)*vmean;

x(:) = ve;

fe = eval(evalstr);

count = count + 1;

if fe < fv(1) && all(xmin<=ve) && all(ve<=xmax)

vk = ve; fk = fe;

end

else

vtmp = v(:,n+1); ftmp = fv(n+1);

if fr < ftmp && all(xmin<=vr) && all(vr<=xmax)

vtmp = vr; ftmp = fr;

end

% Contraction:

vc = beta*vtmp + (1-beta)*vmean;

x(:) = vc;

fc = eval(evalstr);

count = count + 1;

if fc < fv(n) && all(xmin<=vc) && all(vc<=xmax)

vk = vc; fk = fc;

else

% Shrinkage:

for j = 2:n

v(:, j) = (v(:, 1) + v(:, j))/2;

x(:) = v(:, j);

fv(j) = eval(evalstr);

end

count = count + n-1;

vk = (v(:, 1) + v(:, n+1))/2;

x(:) = vk;

fk = eval(evalstr);

count = count + 1;

end

end

v(:, n+1) = vk;

fv(n+1) = fk;

[fv, j] = sort(fv);

v = v(:,j);

end

x = v(:,1);

dx = abs(v(:,n+1)-v(:,1));

x = mask*x + xfix;

dx = mask*dx;

if count>=steps

disp([‘Warning: Maximum number of iterations (‘, int2str(steps),’) has been exceeded’]);

else

steps = count;

end

——–4nd .m

function [cx, tau, offset, csh, z, t, err] = DistFluofit(irf, y, p, dt, shift, flag, bild, N)

% The function DistFluofit performs a fit of a distributed decay curve.

% It is called by:

% [cx, tau, offset, csh, z, t, err] = DistFluofit(irf, y, p, dt, shift).

% The function arguments are:

% irf = Instrumental Response Function

% y = Fluorescence decay data

% p = Time between laser exciation pulses (in nanoseconds)

% dt = Time width of one TCSPC channel (in nanoseconds)

% shift = boundaries of colorshift in channels

%

% The return parameters are:

% cx = lifetime distribution

% tau = used lifetimes

% offset = Offset

% z = Fitted fluorecence curve

% t = time axis

% err = chi2 value

%

% The program needs the following m-files: convol.m.

% (c) 2003 Jˆrg Enderlein

if nargin<6 | isempty(flag)

flag = 0;

end

if nargin<7 | isempty(bild)

bild = 1;

end

close all

if isempty(irf)

irf = zeros(size(y));

irf(1) = 1;

end

irf = irf(:);

y = y(:);

n = length(irf);

tp = dt*(1:p/dt)’;

t = (1:n)’;

if nargin<8 | isempty(N)

N = 100;

end

shifton = 1;

if nargin>4 & ~isempty(shift)

sh_min = shift(1);

sh_max = shift(2);

else

sh_min = -3;

sh_max = 3;

end

%tau = (1/dt/10)./exp((0:N)/N*log(p/dt/10)); % distribution of decay times

tau = (1/dt)./exp((0:N)/N*log(p/dt)); % distribution of decay times

M0 = [ones(size(t)) convol(irf,exp(-tp*tau))];

M0 = M0./(ones(n,1)*sum(M0));

err = [];

if sh_max-sh_min>0

for c=sh_min:sh_max

M = (1-c+floor(c))*M0(rem(rem(t-floor(c)-1, n)+n,n)+1,:) + (c-floor(c))*M0(rem(rem(t-ceil(c)-1, n)+n,n)+1,:);

ind = max([1,1+c]):min([n,n+c]);

cx = lsqnonneg(M(ind,:),y(ind));

z = M*cx;

err = [err sum((z-y).^2./abs(z))/n];

err(end)

end

shv = sh_min:0.1:sh_max;

tmp = interp1(sh_min:sh_max, err, shv);

[pos, pos] = min(tmp);

csh = shv(pos);

else

csh = sh_min;

end

M = (1-csh+floor(csh))*M0(rem(rem(t-floor(csh)-1, n)+n,n)+1,:) + (csh-floor(csh))*M0(rem(rem(t-ceil(csh)-1, n)+n,n)+1,:);

c = ceil(abs(csh))*sign(csh);

ind = max([1,1+c]):min([n,n+c]);

cx = lsqnonneg(M(ind,:),y(ind));

z = M*cx;

err = sum((z-y).^2./abs(z))/n

if bild

t = dt*t;

semilogy(t,y,’.b’,t,z,’g’,’linewidth’,2); times16

v = axis;

v(1) = min(t);

v(2) = max(t);

axis(v);

xlabel(‘time [ns]’);

ylabel(‘lg count’);

figure;

subplot(2,1,1);

plot(t,(y-z)./sqrt(z));

v = axis;

v(1) = min(t);

v(2) = max(t);

axis(v);

xlabel(‘time [ns]’);

ylabel(‘weighted residual’);

ind=1:length(cx)-2;

len = length(ind);

tau = 1./tau;

fac = sqrt(tau(1:end-1)/tau(2:end));

subplot(2,1,2)

semilogx(reshape([fac*tau(ind);fac*tau(ind);tau(ind)/fac;tau(ind)/fac],4*len,1),reshape([0*tau(ind);cx(ind+1)’;cx(ind+1)’;0*tau(ind)],4*len,1));

patch(reshape([fac*tau(ind);fac*tau(ind);tau(ind)/fac;tau(ind)/fac],4*len,1),reshape([0*tau(ind);cx(ind+1)’;cx(ind+1)’;0*tau(ind)],4*len,1),’b’);

xlabel(‘decay time [ns]’);

ylabel(‘distribution’);

end

tau = tau’;

offset = cx(1);

cx(1) = [];

if flag>0

cx = cx’;

tmp = cx>0;

t = 1:length(tmp);

t1 = t(tmp(2:end)>tmp(1:end-1)) + 1;

t2 = t(tmp(1:end-1)>tmp(2:end));

if t1(1)>t2(1)

t2(1)=[];

end

if t1(end)>t2(end)

t1(end)=[];

end

if length(t1)==length(t2)+1

t1(end)=[];

end

if length(t2)==length(t1)+1

t2(1)=[];

end

tmp = []; bla = [];

for j=1:length(t1)

tmp = [tmp cx(t1(j):t2(j))*tau(t1(j):t2(j))/sum(cx(t1(j):t2(j)))];

bla = [bla sum(cx(t1(j):t2(j)))];

end

cx = bla.*tmp;

cx = cx/sum(cx);

tau = tmp;

end

—5nd .m

function y = convol(irf, x)

% convol(irf, x) performs a convolution of the instrumental response

% function irf with the decay function x. Periodicity (=length(x)) is assumed.

mm = mean(irf(end-10:end));

if size(x,1)==1 | size(x,2)==1

irf = irf(:);

x = x(:);

end

p = size(x,1);

n = length(irf);

if p>n

irf = [irf; mm*ones(p-n,1)];

else

irf = irf(1:p);

end

y = real(ifft((fft(irf)*ones(1,size(x,2))).*fft(x)));

t = rem(rem(0:n-1,p)+p,p)+1;

y = y(t,:);

% t = (1:length(irf))’;

% x = x(:)’;

% y = cumsum((irf(:)*ones(1,length(x))).*exp(t*x));

% y = y.*exp(-t*x) + (ones(length(t),1)*((mm*(exp(p*x)-exp(t(end)*x))./(exp(x)-1) + y(end,:))./(1-exp(-p*x)))).*exp(-(p+t)*x);

—-6nd .m; the last ;o)

function err = mlfit(param, irf, y, p)

% MLFIT(param, irf, y, p) returns the ML error between the data y

% and the computed values.

% MLFIT assumes a function of the form:

%

% y = yoffset + A(1)*convol(irf,exp(-t/tau(1)/(1-exp(-p/tau(1)))) + …

%

% param(1) is the color shift value between irf and y.

% param(2) is the irf offset.

% param(3:…) are the decay times.

% irf is the measured Instrumental Response Function.

% y is the measured fluorescence decay curve.

% p is the time between to laser excitations (in number of TCSPC channels).

global bild Plothandle

n = length(irf);

t = 1:n;

tp = (1:p)’;

c = param(1);

tau = param(2:length(param)); tau = tau(:)’;

x = exp(-(tp-1)*(1./tau))*diag(1./(1-exp(-p./tau)));

irs = (1-c+floor(c))*irf(rem(rem(t-floor(c)-1, n)+n,n)+1) + (c-floor(c))*irf(rem(rem(t-ceil(c)-1, n)+n,n)+1);

z = convol(irs, x);

z = [ones(size(z,1),1) z];

A = lsqnonneg(z,y);

z = z*A;

if bild

set(Plothandle,’ydata’,z)

drawnow

end

ind = y>0;

err = sum(y(ind).*log(y(ind)./z(ind))-y(ind)+z(ind))/(n-length(tau))

————————————————–

————one .m defocusing

and one .fig

http://www.joerg-enderlein.de/qdots/QDControl.fig

—-one .m

function varargout = QDControl(varargin)
% QDCONTROL M-file for QDControl.fig
% QDCONTROL, by itself, creates a new QDCONTROL or raises the existing
% singleton*.
%
% H = QDCONTROL returns the handle to a new QDCONTROL or the handle to
% the existing singleton*.
%
% QDCONTROL('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in QDCONTROL.M with the given input arguments.
%
% QDCONTROL('Property','Value',...) creates a new QDCONTROL or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before QDControl_OpeningFunction gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to QDControl_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES
% Edit the above text to modify the response to help QDControl
% Last Modified by GUIDE v2.5 22-Feb-2005 12:49:48
% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @QDControl_OpeningFcn, ...
'gui_OutputFcn', @QDControl_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin & isstr(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
% End initialization code - DO NOT EDIT
% --- Executes just before QDControl is made visible.
function QDControl_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to QDControl (see VARARGIN)
% Choose default command line output for QDControl
handles.output = hObject;
% UIWAIT makes QDControl wait for user response (see UIRESUME)
% uiwait(handles.figure1);
handles.n0 = 1.52;
handles.n1 = 1.0;
handles.lambda = 0.57;
nn = 20;
handles.pix = 6.45;
[x,y] = meshgrid(-nn:nn,-nn:nn);
handles.p = angle(x+i*y);
handles.r = sqrt(x.^2+y.^2);
handles.kappa = str2double(get(handles.kappa_edit,'String'));
handles.ratio = str2double(get(handles.ratio_edit,'String'));
handles.omega = str2double(get(handles.omega_edit,'String'));
handles.theta = str2double(get(handles.theta_edit,'String'));
handles.phi = 0;
if nargin
pos = find(strcmpi('kappa', varargin));
if ~isempty(pos)
handles.kappa = varargin{pos+1};
set(handles.kappa_edit,'String',num2str(handles.kappa));
end
pos = find(strcmpi('ratio', varargin));
if ~isempty(pos)
handles.ratio = varargin{pos+1};
set(handles.ratio_edit,'String',num2str(handles.ratio));
end
pos = find(strcmpi('omega', varargin));
if ~isempty(pos)
handles.omega = varargin{pos+1};
set(handles.omega_edit,'String',num2str(handles.omega));
end
pos = find(strcmpi('theta', varargin));
if ~isempty(pos)
handles.theta = varargin{pos+1};
set(handles.theta_edit,'String',num2str(handles.theta));
end
pos = find(strcmpi('phi', varargin));
if ~isempty(pos)
handles.phi = varargin{pos+1};
end
end
set(handles.theta_slider,'sliderstep',1/90*[1 10],'max',90,'min',0,'Value',handles.theta);
set(handles.omega_slider,'sliderstep',1/180*[1 10],'max',180,'min',0,'Value',handles.omega);
set(handles.kappa_slider,'sliderstep',0.01*[1 5],'max',1,'min',-1,'Value',handles.kappa);
set(handles.ratio_slider,'sliderstep',0.02*[1 5],'max',1,'min',0,'Value',handles.ratio);
set(handles.phi_slider,'sliderstep',1/360*[1 10],'max',360,'min',0,'Value',handles.phi);
handles.mag = str2double(get(handles.mag_edit,'String'));
handles.focus = str2double(get(handles.foc_edit,'String'));
handles.na = str2double(get(handles.na_edit,'String'));
if nargin
pos = find(strcmpi('mag', varargin));
if ~isempty(pos)
handles.mag = varargin{pos+1};
set(handles.mag_edit,'String',num2str(handles.mag));
end
pos = find(strcmpi('focus', varargin));
if ~isempty(pos)
handles.focus = varargin{pos+1};
set(handles.foc_edit,'String',num2str(handles.focus));
end
pos = find(strcmpi('na', varargin));
if ~isempty(pos)
handles.na = varargin{pos+1};
set(handles.na_edit,'String',num2str(handles.na));
end
end
[tmp, tmp, tmp, rho, tmp, fxx0, fxx2, fxz, byx0, byx2, byz] = ...
SEPDipole([0 2],0,handles.na,handles.n0,handles.n1,handles.n1,[],0,[],handles.lambda,handles.mag,handles.focus);
f00 = fxx0.*conj(byx0);
f01 = -fxx0.*conj(byz);
f02 = -fxx0.*conj(byx2);
f10 = -fxz.*conj(byx0);
f11 = fxz.*conj(byz);
f12 = fxz.*conj(byx2);
f20 = -fxx2.*conj(byx0);
f21 = fxx2.*conj(byz);
f22 = fxx2.*conj(byx2);
handles.rho = rho/handles.pix;
handles.f00 = f00;
handles.f01 = f01;
handles.f02 = f02;
handles.f10 = f10;
handles.f11 = f11;
handles.f12 = f12;
handles.f20 = f20;
handles.f21 = f21;
handles.f22 = f22;
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
if nargin
pos = find(strcmpi('image',varargin));
if ~isempty(pos)
axes(handles.imm);
mim(varargin{pos+1});
end
end
axes(handles.pic);
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
% Update handles structure
guidata(hObject, handles);
% --- Outputs from this function are returned to the command line.
function varargout = QDControl_OutputFcn(hObject, eventdata, handles)
% varargout cell array for returning output args (see VARARGOUT);
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Get default command line output from handles structure
varargout{1} = handles.output;
% --- Executes during object creation, after setting all properties.
function na_edit_CreateFcn(hObject, eventdata, handles)
% hObject handle to na_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '1.2');
function na_edit_Callback(hObject, eventdata, handles)
% hObject handle to na_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of na_edit as text
% str2double(get(hObject,'String')) returns contents of na_edit as a double
handles.na = str2double(get(hObject,'String'));
[tmp, tmp, tmp, rho, tmp, fxx0, fxx2, fxz, byx0, byx2, byz] = ...
SEPDipole([0 2],0,handles.na,handles.n0,handles.n1,handles.n1,[],0,[],handles.lambda,handles.mag,handles.focus);
f00 = fxx0.*conj(byx0);
f01 = -fxx0.*conj(byz);
f02 = -fxx0.*conj(byx2);
f10 = -fxz.*conj(byx0);
f11 = fxz.*conj(byz);
f12 = fxz.*conj(byx2);
f20 = -fxx2.*conj(byx0);
f21 = fxx2.*conj(byz);
f22 = fxx2.*conj(byx2);
handles.f00 = f00;
handles.f01 = f01;
handles.f02 = f02;
handles.f10 = f10;
handles.f11 = f11;
handles.f12 = f12;
handles.f20 = f20;
handles.f21 = f21;
handles.f22 = f22;
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function mag_edit_CreateFcn(hObject, eventdata, handles)
% hObject handle to mag_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '110');
function mag_edit_Callback(hObject, eventdata, handles)
% hObject handle to mag_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of mag_edit as text
% str2double(get(hObject,'String')) returns contents of mag_edit as a double
handles.mag = str2double(get(hObject,'String'));
[tmp, tmp, tmp, rho, tmp, fxx0, fxx2, fxz, byx0, byx2, byz] = ...
SEPDipole([0 2],0,handles.na,handles.n0,handles.n1,handles.n1,[],0,[],handles.lambda,handles.mag,handles.focus);
f00 = fxx0.*conj(byx0);
f01 = -fxx0.*conj(byz);
f02 = -fxx0.*conj(byx2);
f10 = -fxz.*conj(byx0);
f11 = fxz.*conj(byz);
f12 = fxz.*conj(byx2);
f20 = -fxx2.*conj(byx0);
f21 = fxx2.*conj(byz);
f22 = fxx2.*conj(byx2);
handles.rho = rho/handles.pix;
handles.f00 = f00;
handles.f01 = f01;
handles.f02 = f02;
handles.f10 = f10;
handles.f11 = f11;
handles.f12 = f12;
handles.f20 = f20;
handles.f21 = f21;
handles.f22 = f22;
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function foc_edit_CreateFcn(hObject, eventdata, handles)
% hObject handle to foc_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '1.2');
function foc_edit_Callback(hObject, eventdata, handles)
% hObject handle to foc_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of foc_edit as text
% str2double(get(hObject,'String')) returns contents of foc_edit as a double
handles.focus = str2double(get(hObject,'String'));
[tmp, tmp, tmp, rho, tmp, fxx0, fxx2, fxz, byx0, byx2, byz] = ...
SEPDipole([0 2],0,handles.na,handles.n0,handles.n1,handles.n1,[],0,[],handles.lambda,handles.mag,handles.focus);
f00 = fxx0.*conj(byx0);
f01 = -fxx0.*conj(byz);
f02 = -fxx0.*conj(byx2);
f10 = -fxz.*conj(byx0);
f11 = fxz.*conj(byz);
f12 = fxz.*conj(byx2);
f20 = -fxx2.*conj(byx0);
f21 = fxx2.*conj(byz);
f22 = fxx2.*conj(byx2);
handles.f00 = f00;
handles.f01 = f01;
handles.f02 = f02;
handles.f10 = f10;
handles.f11 = f11;
handles.f12 = f12;
handles.f20 = f20;
handles.f21 = f21;
handles.f22 = f22;
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function omega_edit_CreateFcn(hObject, eventdata, handles)
% hObject handle to omega_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '0');
function omega_edit_Callback(hObject, eventdata, handles)
% hObject handle to omega_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of omega_edit as text
% str2double(get(hObject,'String')) returns contents of omega_edit as a double
handles.omega = str2double(get(hObject,'String'));
set(handles.omega_slider,'Value',handles.omega);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function theta_edit_CreateFcn(hObject, eventdata, handles)
% hObject handle to theta_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '0');
function theta_edit_Callback(hObject, eventdata, handles)
% hObject handle to theta_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of theta_edit as text
% str2double(get(hObject,'String')) returns contents of theta_edit as a double
handles.theta = str2double(get(hObject,'String'));
set(handles.theta_slider,'Value',handles.theta);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function kappa_edit_CreateFcn(hObject, eventdata, handles)
% hObject handle to kappa_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '0');
function kappa_edit_Callback(hObject, eventdata, handles)
% hObject handle to kappa_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of kappa_edit as text
% str2double(get(hObject,'String')) returns contents of kappa_edit as a double
handles.kappa = str2double(get(hObject,'String'));
set(handles.kappa_slider,'Value',handles.kappa);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function ratio_edit_CreateFcn(hObject, eventdata, handles)
% hObject handle to ratio_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: edit controls usually have a white background on Windows.
% See ISPC and COMPUTER.
if ispc
set(hObject,'BackgroundColor','white');
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
set(hObject, 'String', '0');
function ratio_edit_Callback(hObject, eventdata, handles)
% hObject handle to ratio_edit (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'String') returns contents of ratio_edit as text
% str2double(get(hObject,'String')) returns contents of ratio_edit as a double
handles.ratio = str2double(get(hObject,'String'));
set(handles.ratio_slider,'Value',handles.ratio);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function pic_CreateFcn(hObject, eventdata, handles)
% hObject handle to pic (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: place code in OpeningFcn to populate pic
% --- Executes during object creation, after setting all properties.
function theta_slider_CreateFcn(hObject, eventdata, handles)
% hObject handle to theta_slider (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background, change
% 'usewhitebg' to 0 to use default. See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
set(hObject,'BackgroundColor',[.9 .9 .9]);
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on slider movement.
function theta_slider_Callback(hObject, eventdata, handles)
% hObject handle to theta_slider (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'Value') returns position of slider
% get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.theta = get(hObject,'Value');
set(handles.theta_edit,'String',num2str(handles.theta));
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function omega_slider_CreateFcn(hObject, eventdata, handles)
% hObject handle to omega_slider (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background, change
% 'usewhitebg' to 0 to use default. See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
set(hObject,'BackgroundColor',[.9 .9 .9]);
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on slider movement.
function omega_slider_Callback(hObject, eventdata, handles)
% hObject handle to omega_slider (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'Value') returns position of slider
% get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.omega = get(hObject,'Value');
set(handles.omega_edit,'String',num2str(handles.omega));
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function kappa_slider_CreateFcn(hObject, eventdata, handles)
% hObject handle to kappa_slider (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background, change
% 'usewhitebg' to 0 to use default. See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
set(hObject,'BackgroundColor',[.9 .9 .9]);
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on slider movement.
function kappa_slider_Callback(hObject, eventdata, handles)
% hObject handle to kappa_slider (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'Value') returns position of slider
% get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.kappa = get(hObject,'Value');
set(handles.kappa_edit,'String',num2str(handles.kappa));
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function ratio_slider_CreateFcn(hObject, eventdata, handles)
% hObject handle to ratio_slider (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background, change
% 'usewhitebg' to 0 to use default. See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
set(hObject,'BackgroundColor',[.9 .9 .9]);
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on slider movement.
function ratio_slider_Callback(hObject, eventdata, handles)
% hObject handle to ratio_slider (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'Value') returns position of slider
% get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.ratio = get(hObject,'Value');
set(handles.ratio_edit,'String',num2str(handles.ratio));
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function phi_slider_CreateFcn(hObject, eventdata, handles)
% hObject handle to phi_slider (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: slider controls usually have a light gray background, change
% 'usewhitebg' to 0 to use default. See ISPC and COMPUTER.
usewhitebg = 1;
if usewhitebg
set(hObject,'BackgroundColor',[.9 .9 .9]);
else
set(hObject,'BackgroundColor',get(0,'defaultUicontrolBackgroundColor'));
end
% --- Executes on slider movement.
function phi_slider_Callback(hObject, eventdata, handles)
% hObject handle to phi_slider (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hints: get(hObject,'Value') returns position of slider
% get(hObject,'Min') and get(hObject,'Max') to determine range of slider
handles.phi = get(hObject,'Value');
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, handles.phi, get(handles.radiobutton1,'Value'));
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function radiobutton1_CreateFcn(hObject, eventdata, handles)
% hObject handle to radiobutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
set(hObject,'Value',1);
% --- Executes during object creation, after setting all properties.
function radiobutton2_CreateFcn(hObject, eventdata, handles)
% hObject handle to radiobutton2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
set(hObject,'Value',0);
% --- Executes on button press in radiobutton1.
function radiobutton1_Callback(hObject, eventdata, handles)
% hObject handle to radiobutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hint: get(hObject,'Value') returns toggle state of radiobutton1
set(handles.radiobutton2,'Value',0);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, ...
handles.phi, 1);
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes on button press in radiobutton2.
function radiobutton2_Callback(hObject, eventdata, handles)
% hObject handle to radiobutton2 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% Hint: get(hObject,'Value') returns toggle state of radiobutton2
set(handles.radiobutton1,'Value',0);
handles.int = int_compute(handles.theta, handles.omega, handles.kappa, handles.ratio, handles.rho, handles.r, handles.p, ...
handles.f00, handles.f01, handles.f02, handles.f10, handles.f11, handles.f12, handles.f20, handles.f21, handles.f22, ...
handles.phi, 0);
mim(handles.int); title(int2str(round(handles.phi)),'fontsize',12)
guidata(hObject, handles);
% --- Executes during object creation, after setting all properties.
function imm_CreateFcn(hObject, eventdata, handles)
% hObject handle to imm (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called
% Hint: place code in OpeningFcn to populate imm
return
function int = int_compute(theta, omega, kappa, ratio, rho, r, p, f00, f01, f02, f10, f11, f12, f20, f21, f22, phi, flag);
theta = theta/180*pi;
omega = omega/180*pi;
p = p - phi/180*pi;
if flag==1
psi = 0;
intc1(1,:) = 3*f00+3*f22+f11 - f00*cos(2*psi)*(1-cos(2*theta)+kappa*(3+cos(2*theta))*cos(2*omega)) + ...
(f00-f11+f22)*(cos(2*theta)+2*kappa*cos(2*omega)*sin(theta)^2) + 4*kappa*f00*cos(theta)*sin(2*psi)*sin(2*omega);
intc1(2,:) = -2*(2*f01-f21+2*f10-f12)*sin(theta)*(cos(psi)*cos(theta)*(1-kappa*cos(2*omega))+kappa*sin(psi)*sin(2*omega));
intc1(3,:) = -3*f20+f11-3*f02 + (f20+f02)*cos(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2) - ...
(f20+f11+f02)*(cos(2*theta)+2*kappa*cos(2*omega)*sin(theta)^2) - 4*(f20+f02)*kappa*cos(theta)*sin(2*psi)*sin(2*omega);
intc1(4,:) = 2*(f21+f12)*sin(theta)*(cos(psi)*cos(theta)*(1-kappa*cos(2*omega))+kappa*sin(psi)*sin(2*omega));
intc1(5,:) = -f22*(cos(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2)-4*kappa*cos(theta)*sin(2*psi)*sin(2*omega));
ints1(1,:) = 2*(f21+f12)*sin(theta)*(cos(theta)*(1-kappa*cos(2*omega))*sin(psi)-kappa*cos(psi)*sin(2*omega));
ints1(2,:) = (f20+f02)*(sin(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2)+4*kappa*cos(2*psi)*cos(theta)*sin(2*omega));
ints1(3,:) = ints1(1,:);
ints1(4,:) = -f22*(sin(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)-2*sin(theta)^2)+4*kappa*cos(2*psi)*cos(theta)*sin(2*omega));
psi = pi/2;
intc2(1,:) = 3*f00+3*f22+f11 - f00*cos(2*psi)*(1-cos(2*theta)+kappa*(3+cos(2*theta))*cos(2*omega)) + ...
(f00-f11+f22)*(cos(2*theta)+2*kappa*cos(2*omega)*sin(theta)^2) + 4*kappa*f00*cos(theta)*sin(2*psi)*sin(2*omega);
intc2(2,:) = -2*(2*f01-f21+2*f10-f12)*sin(theta)*(cos(psi)*cos(theta)*(1-kappa*cos(2*omega))+kappa*sin(psi)*sin(2*omega));
intc2(3,:) = -3*f20+f11-3*f02 + (f20+f02)*cos(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2) - ...
(f20+f11+f02)*(cos(2*theta)+2*kappa*cos(2*omega)*sin(theta)^2) - 4*(f20+f02)*kappa*cos(theta)*sin(2*psi)*sin(2*omega);
intc2(4,:) = 2*(f21+f12)*sin(theta)*(cos(psi)*cos(theta)*(1-kappa*cos(2*omega))+kappa*sin(psi)*sin(2*omega));
intc2(5,:) = -f22*(cos(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2)-4*kappa*cos(theta)*sin(2*psi)*sin(2*omega));
ints2(1,:) = 2*(f21+f12)*sin(theta)*(cos(theta)*(1-kappa*cos(2*omega))*sin(psi)-kappa*cos(psi)*sin(2*omega));
ints2(2,:) = (f20+f02)*(sin(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)+2*sin(theta)^2)+4*kappa*cos(2*psi)*cos(theta)*sin(2*omega));
ints2(3,:) = ints2(1,:);
ints2(4,:) = -f22*(sin(2*psi)*(kappa*(3+cos(2*theta))*cos(2*omega)-2*sin(theta)^2)+4*kappa*cos(2*psi)*cos(theta)*sin(2*omega));
% int = interp1(rho,real(intc1(1,:)+intc2(1,:)),r,'cubic');
int = interp1(rho,real(intc1(1,:)),r,'cubic') + interp1(rho,real(intc2(1,:)),r,'cubic');
for jj=1:4
int = int + cos(jj*p).*interp1(rho,real(intc1(jj+1,:)),r,'cubic') + ...
sin(jj*p).*interp1(rho,real(ints1(jj,:)),r,'cubic');
int = int + cos(jj*(p+pi/2)).*interp1(rho,real(intc2(jj+1,:)),r,'cubic') + ...
sin(jj*(p+pi/2)).*interp1(rho,real(ints2(jj,:)),r,'cubic');
end
intc1(1,:) = 4*(f00+f22)*sin(theta)^2+4*f11*cos(theta)^2;
intc1(2,:) = -2*(f01-f21+f10-f12)*sin(2*theta);
intc1(3,:) = -4*(f20+f02)*sin(theta)^2;
tmp = interp1(rho,real(intc1(1,:)),r,'cubic');
for jj=1:2
tmp = tmp + cos(jj*p).*interp1(rho,real(intc1(jj+1,:)),r,'cubic');
end
else
intc1(1,:) = (f00+f22)*(cos(omega)^2*cos(theta)^2 + sin(omega)^2) + f11*cos(omega)^2*sin(theta)^2;
intc1(2,:) = -0.5*(f01-f21+f10-f12)*sin(2*theta)*cos(omega)^2;
intc1(3,:) = 0.25*(f20+f02)*(1-3*cos(2*omega)-2*cos(omega)^2*cos(2*theta));
ints1(1,:) = -0.5*(f01-f21+f10-f12)*sin(2*omega)*sin(theta);
ints1(2,:) = -(f20+f02)*sin(2*omega)*cos(theta);
omega = omega + pi/2;
intc2(1,:) = (f00+f22)*(cos(omega)^2*cos(theta)^2 + sin(omega)^2) + f11*cos(omega)^2*sin(theta)^2;
intc2(2,:) = -0.5*(f01-f21+f10-f12)*sin(2*theta)*cos(omega)^2;
intc2(3,:) = 0.25*(f20+f02)*(1-3*cos(2*omega)-2*cos(omega)^2*cos(2*theta));
ints2(1,:) = -0.5*(f01-f21+f10-f12)*sin(2*omega)*sin(theta);
ints2(2,:) = -(f20+f02)*sin(2*omega)*cos(theta);
int = (1-kappa)/2*interp1(rho,real(intc1(1,:)),r,'cubic') + (1+kappa)/2*interp1(rho,real(intc2(1,:)),r,'cubic');
for jj=1:2
int = int + (1-kappa)/2*cos(jj*p).*interp1(rho,real(intc1(jj+1,:)),r,'cubic') + ...
(1-kappa)/2*sin(jj*p).*interp1(rho,real(ints1(jj,:)),r,'cubic');
int = int + (1+kappa)/2*cos(jj*p).*interp1(rho,real(intc2(jj+1,:)),r,'cubic') + ...
(1+kappa)/2*sin(jj*p).*interp1(rho,real(ints2(jj,:)),r,'cubic');
end
intc1(1,:) = (f00+f22)*sin(theta)^2+f11*cos(theta)^2;
intc1(2,:) = -0.5*(f01-f21+f10-f12)*sin(2*theta);
intc1(3,:) = -(f20+f02)*sin(theta)^2;
tmp = interp1(rho,real(intc1(1,:)),r,'cubic');
for jj=1:2
tmp = tmp + cos(jj*p).*interp1(rho,real(intc1(jj+1,:)),r,'cubic');
end
end
int = (1-ratio)*int + ratio*tmp;
return
% --- Executes during object deletion, before destroying properties.
function pic_DeleteFcn(hObject, eventdata, handles)
% hObject handle to pic (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% --- Executes during object deletion, before destroying properties.
function imm_DeleteFcn(hObject, eventdata, handles)
% hObject handle to imm (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% --- Executes during object deletion, before destroying properties.
function figure1_DeleteFcn(hObject, eventdata, handles)
% hObject handle to figure1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
function [intx, inty, intz, rho, phi, fxx0, fxx2, fxz, byx0, byx2, byz, s0, dphase, v, pc, ps, psi] = SEPDipole(rho, z, NA, n1, n, n2, d1, d, d2, lambda, mag, focpos, atf, alpha)
% SEPDipole calculates the emission pattern of dipoles while imaged with a camera
%
% [intx, inty, intz, rho, phi, fxx0, fxx2, fxz, byx0, byx2, byzs0, dphase, v, pc, ps, psi] = SEPDipole(rho, z, NA, n1, n, n2, d1, d, d2, lambda, mag, focpos)
%
% input arguments:
% r - two-element vector containing the minimum and maximum value of the rho coordinate
% z - distance of molecule from surface
% NA - numerical aperture
% n1 - refractive index vector (bottom to top) below molecule
% n - refractive index of molecule's layer
% n2 - refractive index vector (bottom to top) above molecule
% d1 - thickness vector for layers below molecule
% d - thickness of embedding layer
% d2 - thickness vector for layers above molecule
% lambda - wavelength
% mag - magnification
% focpos - focus position along optical axis
%
% output arguments:
% int* - intensity distributions
% f** - electric field component distributions
% b** - magnetic field component distributions
% rho, phi - 2d coordinate fields
% s0 - aperture variable
% dphase - wave front deviation
% v - angular distribution of radiation for vertical dipole
% pc, ps - cosine and sine part of angular distribution of radiation for parallel dipole
% psi - emission angle for angular distribution of radiation
maxnum = 1e3;
ni = 1.0; % it is assumed that imaging medium is air
ki = 2*pi/lambda*ni;
if length(rho)==2
drho = lambda/50;
rho = rho(1) + (0:(rho(2)-rho(1))/drho)*drho;
elseif length(rho)>2
drho = diff(rho(1:2));
end
rho = mag*rho; rho = rho(:)';
row = ones(size(rho));
chimin = 0;
chimax = asin(NA/mag/ni);
dchi = (chimax-chimin)/maxnum;
chi = chimin + (0.5:maxnum)*dchi;
ci = cos(chi);
si = sin(chi);
s0 = ni/n1(1)*mag*si;
c0 = sqrt(1-s0.^2);
psi = asin(s0);
[v,pc,ps] = DipoleL(psi,2*pi*z/lambda,n1,n,n2,2*pi*d1/lambda,2*pi*d/lambda,2*pi*d2/lambda);
v = v.'; ps = ps.'; pc = pc.';
if nargin>12 & ~isempty(atf)
if length(atf)==1 % accounting for reflection losses when using water immersion; atf = ref index of cover slide
[tmp, tms, tmp, tms] = Fresnel(n1(1)*c0,n1(1),atf);
[mp, ms, mp, ms] = Fresnel(sqrt(atf^2-n1(1)^2*s0.^2),atf,n1(1));
else % + aberration for water immersion; atf(2) = thickness mismatch
[tmp, tms, tmp, tms] = Fresnel(n1(1)*c0,[n1(1) atf(1) atf(1)], 2*pi*atf(2)/lambda);
[mp, ms, mp, ms] = Fresnel(sqrt(atf(1)^2-n1(1)^2*s0.^2),[atf(1) n1(1) n1(1)], -2*pi*atf(2)/lambda);
end
v = tmp.*mp.*v;
pc = tmp.*mp.*pc;
ps = tms.*ms.*ps;
end
phase = -n1(1).*c0*focpos;
if nargin>13 & ~isempty(alpha)
dphase = 0;
for j=1:length(alpha)
dphase = dphase + alpha(j)*(s0/max(s0)).^(2*j);
end
phase = phase + dphase;
plot(s0/max(s0),dphase);
end
%alternative formulation but same result:
%phase = ni*mag*focpos*s0./c0.*si-n1(1)./c0*focpos;
ez = exp(i*2*pi/lambda*phase);
fac = si.*sqrt(ci./c0); % aplanatic objective
barg = ki*si'*rho;
j0 = besselj(0,barg); j1 = besselj(1,barg); j2 = besselj(2,barg);
ezi = fac.*ci.*ez;
ezr = fac.*ez;
fxx0 = (ezi.*pc+ezr.*ps)*j0; % cos(0*phi)-component
fxx2 = -(ezi.*pc-ezr.*ps)*j2; % cos(2*phi)-component
fxz = -2*i*(ezi.*v)*j1; % cos(1*phi)-component
byx0 = (ezr.*pc+ezi.*ps)*j0; % cos(0*phi)-component
byx2 = -(ezr.*pc-ezi.*ps)*j2; % cos(2*phi)-component
byz = -2*i*(ezr.*v)*j1; % cos(1*phi)-component
phi = 0:pi/100:2*pi; phi = phi'; col = ones(size(phi));
intx = real((col*fxx0+cos(2*phi)*fxx2).*conj(col*byx0+cos(2*phi)*byx2)); % int for x-dipole along x-pol
inty = real((sin(2*phi)*fxx2).*conj(sin(2*phi)*byx2)); % int for x-dipole along y-pol
intz = col*real(fxz.*conj(byz)); % int for z-dipole along x-pol
s0 = s0/max(s0);
function [v,pc,ps,tp,ts,tp1,ts1,fac] = DipoleL(theta,z,n1,n,n2,d1,d,d2)
% [v,pc,ps] = DipoleL(theta,z,n,d) calculates the electric field amplitudes of dipole radiation along emission angle theta
% of a dipole at distance z from an interface within a layer
% theta - direction of radiation downwards
% z - molecule's distance from the bottom of its layer
% n1 - vector of refractive indices of the stack below the molecule's layer
% n - refracive index of the molecule's layer
% n2 - vector of refractive indices of the stack above the the molecule's layer
% d1 - vector of layer thickness values of the stack below the molecule's layer ( length(d1)=length(n1)-1 )
% d - thickness of molecule's layer
% d2 - vector of layer thickness values of the stack above the molecule's layer ( length(d2)=length(n2)-1 )
z = z(:)';
col = ones(size(z));
theta = abs(theta(:));
ind = theta<pi/2;
v = zeros(length(theta),length(z));
pc = v; ps = v;
if sum(ind)>0
tmp = theta(ind);
w = sqrt(n^2-n1(1)^2*sin(tmp).^2);
[rpu, rsu, tpu, tsu] = Fresnel(w,[n n2],d2);
[rpd, rsd, tpd, tsd] = Fresnel(w,[n n1(end:-1:1)],d1(end:-1:1));
tp = (tpd./(1-rpu.*rpd.*exp(2*i*w.'*d))).';
ts = (tsd./(1-rsu.*rsd.*exp(2*i*w.'*d))).';
tp1 = tp.*(rpu.*exp(2*i*w.'*d)).';
ts1 = ts.*(rsu.*exp(2*i*w.'*d)).';
fac = sqrt(n1(1))*n1(1)*cos(tmp)./w;
ez = exp(i*w*z);
v(ind,:) = ((fac.*n1(1)/n.*sin(tmp).*tp)*col).*ez + ((fac.*n1(1)/n.*sin(tmp).*tp1)*col)./ez;
pc(ind,:) = ((fac.*w/n.*tp)*col).*ez - ((fac.*w/n.*tp1)*col)./ez;
ps(ind,:) = ((fac.*ts)*col).*ez + ((fac.*ts1)*col)./ez;
end
function [rp,rs,tp,ts] = Fresnel(w1,n1,n2);
maxval = exp(100);
minval = exp(-100);
w1 = w1(:).'; n1 = n1(:).'; n2 = n2(:).';
if length(n1)==1 & length(n2)==1
w2 = sqrt(n2^2-n1^2+w1.^2);
w2(imag(w2)<0) = conj(w2(imag(w2)<0));
rp = (w1*n2^2-w2*n1^2)./(w1*n2^2+w2*n1^2);
rs = (w1-w2)./(w1+w2);
tp = 2*n1*n2*w1./(w1*n2^2+w2*n1^2);
ts = 2*w1./(w1+w2);
elseif length(n1)==length(n2)+2
if length(n1)==2
[rp,rs,tp,ts] = Fresnel(w1,n1(1),n1(2));
else
n = n1; d = [0 n2 0];
w(1,:) = w1;
for j = 2:length(n)
w(j,:) = sqrt(n(j)^2-n(1)^2+w1.^2);
w(j,imag(w(j,:))<0) = conj(w(j,imag(w(j,:))<0));
end
j = length(n);
M11 = (w(j,:)./w(j-1,:)*n(j-1)/n(j)+n(j)/n(j-1))/2;
M12 = (-w(j,:)./w(j-1,:)*n(j-1)/n(j)+n(j)/n(j-1))/2;
M21 = M12;
M22 = M11;
for j=length(n)-1:-1:2
M11 = exp(-i*w(j,:)*d(j)).*M11;
M12 = exp(-i*w(j,:)*d(j)).*M12;
M21 = exp(i*w(j,:)*d(j)).*M21;
M22 = exp(i*w(j,:)*d(j)).*M22;
N11 = (w(j,:)./w(j-1,:)*n(j-1)/n(j)+n(j)/n(j-1))/2;
N12 = (-w(j,:)./w(j-1,:)*n(j-1)/n(j)+n(j)/n(j-1))/2;
N21 = N12;
N22 = N11;
tmp11 = N11.*M11+N12.*M21;
tmp12 = N11.*M12+N12.*M22;
tmp21 = N21.*M11+N22.*M21;
tmp22 = N21.*M12+N22.*M22;
M11 = tmp11;
M12 = tmp12;
M21 = tmp21;
M22 = tmp22;
end
rp = M21./M11;
tp = 1./M11;
j = length(n);
M11 = (w(j,:)./w(j-1,:)+1)/2;
M12 = (-w(j,:)./w(j-1,:)+1)/2;
M21 = M12;
M22 = M11;
for j=length(n)-1:-1:2
M11 = exp(-i*w(j,:)*d(j)).*M11;
M12 = exp(-i*w(j,:)*d(j)).*M12;
M21 = exp(i*w(j,:)*d(j)).*M21;
M22 = exp(i*w(j,:)*d(j)).*M22;
N11 = (w(j,:)./w(j-1,:)+1)/2;
N12 = (-w(j,:)./w(j-1,:)+1)/2;
N21 = N12;
N22 = N11;
tmp11 = N11.*M11+N12.*M21;
tmp12 = N11.*M12+N12.*M22;
tmp21 = N21.*M11+N22.*M21;
tmp22 = N21.*M12+N22.*M22;
M11 = tmp11;
M12 = tmp12;
M21 = tmp21;
M22 = tmp22;
end
rs = M21./M11;
ts = 1./M11;
end
end
function mim(x,p1,p2)
if nargin==1
nd = ndims(x);
switch nd
case 2
imagesc(x);
axis image
colormap hot
axis off
case 3
nj = size(x,3);
for j=1:nj
subplot(1,nj,j);
mim(x(:,:,j))
end
case 4
nj = size(x,3);
nk = size(x,4);
for j=1:nj
for k=1:nk
subplot(nk,nj,j+(k-1)*nj);
mim(x(:,:,j,k))
end
end
end
end
if nargin==2
if prod(size(p1))==2
imagesc(x,p1);
axis image
colormap hot
axis off
elseif p1=='h'
mim(x);
colorbar('h');
elseif p1=='v'
mim(x);
colorbar;
elseif isstr(p1)
mim(x);
[a,b]=size(x);
text(b*(1-0.025*length(p1)),0.06*a,p1,'FontName','Times','FontSize',16,'Color','w');
else
c = zeros(size(x,1),size(x,2),3);
mm = min(p1(isfinite(p1)));
tmp = max(p1(isfinite(p1)))-mm;
if tmp>0
p1 = (p1-mm)/tmp;
else
p1 = zeros(size(p1));
end
mm = min(x(isfinite(p1)));
tmp = max(x(isfinite(p1)))-mm;
if tmp>0
x = 1 + (x-mm)/tmp*63;
else
x = 64*ones(size(x));
end
x(isnan(x)) = 1;
x(x<1) = 1;
x(x>64) = 64;
col = double(jet);
for j=1:3
c(:,:,j) = ((floor(x)+1-x).*reshape(col(floor(x),j),size(x)) + (x-floor(x)).*reshape(col(ceil(x),j),size(x))).*p1;
end
if tmp>0
subplot('position',[0.1 0.1 0.7 0.8]);
imagesc(c);
axis image
axis off
subplot('position',[0.8 0.1 0.1 0.8]);
pcolor([0 0.1]*tmp,mm+(1:size(x,1))/size(x,1)*tmp,(mm+(1:size(x,1))'/size(x,1)*tmp)*ones(1,2))
shading flat
axis image
set(gca,'xtick',[],'yaxislocation','right')
colormap(jet)
else
imagesc(c);
axis image
axis off
end
end
end
if nargin==3
if size(x,1)==size(p1,1) & size(x,2)==size(p1,2)
if prod(size(p2))==2
c = zeros(size(x,1),size(x,2),3);
mm = min(p1(isfinite(p1)));
tmp = max(p1(isfinite(p1)))-mm;
if tmp>0
p1 = (p1-mm)/tmp;
else
p1 = zeros(size(p1));
end
mm = p2(1);
tmp = p2(2) - mm;
if tmp>0
x = 1 + (x-mm)/tmp*63;
else
x = 64*ones(size(x));
end
x(isnan(x)) = 1;
x(x<1) = 1;
x(x>64) = 64;
col = double(jet);
for j=1:3
c(:,:,j) = ((floor(x)+1-x).*reshape(col(floor(x),j),size(x)) + (x-floor(x)).*reshape(col(ceil(x),j),size(x))).*p1;
end
if tmp>0
subplot('position',[0.1 0.1 0.7 0.8]);
imagesc(c);
axis image
axis off
subplot('position',[0.8 0.1 0.1 0.8]);
pcolor([0 0.1]*tmp,mm+(1:size(x,1))/size(x,1)*tmp,(mm+(1:size(x,1))'/size(x,1)*tmp)*ones(1,2))
shading flat
axis image
set(gca,'xtick',[],'yaxislocation','right')
colormap(jet)
else
imagesc(c);
axis image
axis off
end
end
elseif prod(size(p1))==2
mim(x,p1);
if p2=='h'
colorbar('h');
elseif p2=='v'
colorbar;
end
elseif prod(size(p2))==2
mim(x,p2);
if p1=='h'
colorbar('h');
elseif p1=='v'
colorbar;
end
end
end
figure(gcf);

## Commentaires récents