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);
it's very help me... THANKS
ReplyDelete