matlab .m; program for fluorescence lifetime fitting and Defocused imaging of single molecules and QD

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);

matlab .m; random walk

Routines for simulating paths of stochastic processes: random walk, Poisson process, Brownian motion and their multidimensional versions, as well as birth-and-death processes, branching and reproduction models. See

http://www.math.uu.se/research/telecom/software/

2002; birthdeath.m, bm3plot.m, brownian.m, galtonwatson.m, moran.m, poisson2d.m, poisson3d.m, poissonjp.m, poissonti.m, ranwalk2d.m, ranwalk3d.m, ranwalk.m, README, rw3plot.m

comment: Intuitive, but not the most efficient way of programming.

Ref.

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=2493&objectType=file

----------------

a 2-Dimensional Random Walk process program in matlab.

comment: A couple remarks for this homework solution:
- several (most!) variables are never used (T,s,t, mx,my,mr)
- explain what the input variable n is in the help, because not "any value of time n" is valid. (it can be a positive integer only!)
- add an error check for that in the code
- The help can be improved. Explain what the user may expect (a plot!).
- why not return both X and Y?
- Have you actually typed "help two_D_Rand_walk"? What you get is not very informative for the user, isn't it? Or do you want a lot of emails ;-)

Ref.

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=14864&objectType=file

----------

http://ocw.mit.edu/OcwWeb/Biological-Engineering/20-011JSpring-2004/Simulations/

Simulations

MATLAB® software is required to run the .m and .fig files in this section.Resources

Monte Carlo (PDF)

Matlab Tutorial (PDF)

PCR simulator

pcr_sim.m (M)
Replicate a certain number of starting copies of DNA for several cycles. Try different probabilities of replication to get different histograms showing how much DNA is produced over many trials.

pcr_sim.fig (FIG)
Put this in the same directory as pcr_sim.m to be able to run "pcr_sim" from the matlab prompt.

Random Walk - 2D with 2 Particles (M)
Look at how these simulations work and try making your own 3D random walk or have your particle trace out a trajectory in Matlab.

Lab

Random Walk - 2D (M)

----------

----------tutorial of Matlab (first approach for dummies) , exemple walk 1D:

web.mit.edu/14.451/www/Matlab_Tutorial.pdf

---------

matlab .m Techniques for Visualizing Scalar Volume MRI Data

Techniques for Visualizing Scalar Volume Data

What Is Scalar Volume Data?

Typical scalar volume data is composed of a 3-D array of data and three coordinate arrays of the same dimensions. The coordinate arrays specify the x-, y-, and z-coordinates for each data point.

The units of the coordinates depend on the type of data. For example, flow data might have coordinate units of inches and data units of psi.

A number of MATLAB^® functions are useful for visualizing scalar data:

• Slice planes provide a way to explore the distribution of data values within the volume by mapping values to colors. You can orient slice planes at arbitrary angles, as well as use nonplanar slices. (For illustrations of how to use slice planes, see slice, a volume slicing example, and slice planes used to show context.) You can specify the data used to color isosurfaces, enabling you to display different information in color and surface shape (see isocolors).
• Contour slices are contour plots drawn at specific coordinates within the volume. Contour plots enable you to see where in a given plane the data values are equal. See contourslice for an example
• Isosurfaces are surfaces constructed by using points of equal value as the vertices of patch graphics objects.

updateSectionId("f5-14693");addTopOfPageButtons(); Back to Top

Example — Ways to Display MRI Data

Changing the Data Format
Displaying Images of MRI Data
Displaying a 2-D Contour Slice
Displaying 3-D Contour Slices
Displaying an Isosurface
Adding an Isocap to Show a Cutaway Surface
Defining the View
Add Lighting

An example of scalar data includes Magnetic Resonance Imaging (MRI) data. This data typically contains a number of slice planes taken through a volume, such as the human body. MATLAB includes an MRI data set that contains 27 image slices of a human head. This example illustrate the following techniques applied to MRI data:

• A series of 2-D images representing slices through the head
• 2-D and 3-D contour slices taken at arbitrary locations within the data
• An isosurface with isocaps showing a cross section of the interior

Changing the Data Format

The MRI data, D, is stored as a 128-by-128-by-1-by-27 array. The third array dimension is used typically for the image color data. However, since these are indexed images (a colormap, map, is also loaded) there is no information in the third dimension, which you can remove using the squeeze command. The result is a 128-by-128-by-27 array.

The first step is to load the data and transform the data array from 4-D to 3-D.

load mri
D = squeeze(D);

Displaying Images of MRI Data

To display one of the MRI images, use the image command, indexing into the data array to obtain the eighth image. Then adjust axis scaling, and install the MRI colormap, which was loaded along with the data.

image_num = 8;
image(D(:,:,image_num))
axis image
colormap(map)

Save the x- and y-axis limits for use in the next part of the example.

x = xlim;
y = ylim;

Displaying a 2-D Contour Slice

You can treat this MRI data as a volume because it is a collection of slices taken progressively through the 3-D object. Use contourslice to display a contour plot of a slice of the volume. To create a contour plot with the same orientation and size as the image created in the first part of this example, adjust the y-axis direction (axis), set the limits (xlim, ylim), and set the data aspect ratio (daspect).

contourslice(D,[],[],image_num)
axis ij
xlim(x)
ylim(y)
daspect([1,1,1])
colormap('default')

This contour plot uses the figure colormap to map color to contour value.

Displaying 3-D Contour Slices

Unlike images, which are 2-D objects, contour slices are 3-D objects that you can display in any orientation. For example, you can display four contour slices in a 3-D view. To improve the visibility of the contour line, increase the LineWidth to 2 points (one point equals 1/72 of an inch).

phandles = contourslice(D,[],[],[1,12,19,27],8);
view(3); axis tight
set(phandles,'LineWidth',2)

Displaying an Isosurface

You can use isosurfaces to display the overall structure of a volume. When combined with isocaps, this technique can reveal information about data on the interior of the isosurface.

First, smooth the data with smooth3; then use isosurface to calculate the isodata. Use patch to display this data as a graphics object.

Ds = smooth3(D);
hiso = patch(isosurface(Ds,5),...
'FaceColor',[1,.75,.65],...
'EdgeColor','none');

Adding an Isocap to Show a Cutaway Surface

Use isocaps to calculate the data for another patch that is displayed at the same isovalue (5) as the surface. Use the unsmoothed data (D) to show details of the interior. You can see this as the sliced-away top of the head.

hcap = patch(isocaps(D,5),...
'FaceColor','interp',...
'EdgeColor','none');
colormap(map)

Defining the View

Define the view and set the aspect ratio (view, axis, daspect).

view(45,30) 
axis tight 
daspect([1,1,.4])

Add Lighting

Add lighting and recalculate the surface normals based on the gradient of the volume data, which produces smoother lighting (camlight, lighting, isonormals). Increase the AmbientStrength property of the isocap to brighten the coloring without affecting the isosurface. Set the SpecularColorReflectance of the isosurface to make the color of the specular reflected light closer to the color of the isosurface; then set the SpecularExponent to reduce the size of the specular spot.

lightangle(45,30); 
set(gcf,'Renderer','zbuffer'); lighting phong
isonormals(Ds,hiso)
set(hcap,'AmbientStrength',.6)
set(hiso,'SpecularColorReflectance',0,'SpecularExponent',50)

Example of an Isocap

The isocap uses interpolated face coloring, which means the figure colormap determines the coloring of the patch. This example uses the colormap supplied with the data.

To display isocaps at other data values, try changing the isosurface value or use the subvolume command. See the isocaps and subvolume reference pages for examples.

Ref.

http://www.mathworks.com/access/helpdesk/help/techdoc/index.html?/access/helpdesk/help/techdoc/visualize/f5-3291.html&http://www.mathworks.com/cgi-bin/texis/webinator/search/?db=MSS&prox=page&rorder=750&rprox=750&rdfreq=500&rwfreq=500&rlead=250&sufs=0&order=r&is_summary_on=1&ResultCount=10&query=MRI&submitButtonName=Search

Volume visualization of data sets defined on 3D grids. help matlab

http://www.mathworks.com/access/helpdesk/help/techdoc/index.html?/access/helpdesk/help/techdoc/visualize/f5-6010.html&http://www.google.fr/search?hl=fr&client=firefox-a&rls=org.mozilla%3Afr%3Aofficial&hs=PuP&q=volume+visualization+matlab&btnG=Rechercher&meta=

Volume visualization is the creation of graphical representations of data sets that are defined on three-dimensional grids. Volume data sets are characterized by multidimensional arrays of scalar or vector data. These data are typically defined on lattice structures representing values sampled in 3-D space. There are two basic types of volume data:

• Scalar volume data contains single values for each point.
• Vector volume data contains two or three values for each point, defining the components of a vector.

An example of scalar volume data is that produced by the flow M-file. The flow data represents the speed profile of a submerged jet within an infinite tank. Typing

[x,y,z,v] = flow;

produces four 3-D arrays. The x, y, and z arrays specify the coordinates of the scalar values in the array v.

The wind data set is an example of vector volume data that represents air currents over North America. You can load this data in the MATLAB^® workspace with the command

load wind

This data set comprises six 3-D arrays: x, y, and z are the coordinate data for the arrays u, v, and w, which are the vector components for each point in the volume.

updateSectionId("f5-4446");addTopOfPageButtons(); Back to Top

Selecting Visualization Techniques

The techniques you select to visualize volume data depend on what type of data you have and what you want to learn. In general,

• Scalar data is best viewed with isosurfaces, slice planes, and contour slices.
• Vector data represents both a magnitude and direction at each point, which is best displayed by stream lines (particles, ribbons, and tubes), cone plots, and arrow plots. Most visualizations, however, employ a combination of techniques to best reveal the content of the data.

The material in these sections describes how to apply a variety of techniques to typical volume data.

updateSectionId("f5-5598");addTopOfPageButtons(); Back to Top

Steps to Create a Volume Visualization

Creating an effective visualization requires a number of steps to compose the final scene. These steps fall into four basic categories:

1. Determine the characteristics of your data. Graphing volume data usually requires knowledge of the range of both the coordinates and the data values.
2. Select an appropriate plotting routine. The information in this section helps you select the right methods.
3. Define the view. The information conveyed by a complex three-dimensional graph can be greatly enhanced through careful composition of the scene. Viewing techniques include adjusting camera position, specifying aspect ratio and project type, zooming in or out, and so on.
4. Add lighting and specify coloring. Lighting is an effective means to enhance the visibility of surface shape and to provide a three-dimensional perspective to volume graphs. Color can convey data values, both constant and varying.

updateSectionId("brbl8tm");addTopOfPageButtons(); Back to Top

Volume Visualization Functions

MATLAB functions enable you to apply a variety of volume visualization techniques. The following tables group these functions into two categories based on the type of data (scalar or vector) that each is designed to work with. The reference page for each function provides examples of the intended use.

Functions for Scalar Data

Function	Purpose
contourslice	Draw contours in volume slice planes
isocaps	Compute isosurface end-cap geometry
isocolors	Compute the colors of isosurface vertices
isonormals	Compute normals of isosurface vertices
isosurface	Extract isosurface data from volume data
patch	Create a patch (multipolygon) graphics object
reducepatch	Reduce the number of patch faces
reducevolume	Reduce the number of elements in a volume data set
shrinkfaces	Reduce the size of each patch face
slice	Draw slice planes in volume
smooth3	Smooth 3-D data
surf2patch	Convert surface data to patch data
subvolume	Extract subset of volume data set

Functions for Vector Data

Function	Purpose
coneplot	Plot velocity vectors as cones in 3-D vector fields
curl	Compute the curl and angular velocity of a 3-D vector field
divergence	Compute the divergence of a 3-D vector field
interpstreamspeed	Interpolate streamline vertices from vector-field magnitudes
streamline	Draw stream lines from 2-D or 3-D vector data
streamparticles	Draw stream particles from vector volume data
streamribbon	Draw stream ribbons from vector volume data
streamslice	Draw well-spaced stream lines from vector volume data
streamtube	Draw stream tubes from vector volume data
stream2	Compute 2-D stream line data
stream3	Compute 3-D stream line data
volumebounds	Return coordinate and color limits for volume (scalar and vector)

matlab .m , voxel render a 3-D array using 2-D OpenGL texture maps.

VOL3D Volume (voxel) render a 3-D array using 2-D OpenGL texture maps.

Useful for visualizing and exploring 3-D data such as MRI images.

This function requires OpenGL hardware acceleration. See 'opengl' command for more information.

Use vol3dtool for editing the colormap and alphamap. Adjusting these maps will allow you to explore your 3-D volume data at various intensity levels.

Example:

v = flow(50);
h = vol3d('cdata',v,'texture','3D');

See vol3d.m m-help for more information.

Ref.

http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=4927&objectType=FILE

Matlab .m Unstructured Mesh Generation (Mesh2d v2.3)

Unstructured Mesh Generation (Mesh2d v2.3)

MESH2D is a toolbox for the generation and manipulation of unstructured triangular meshes in MATLAB. High quality meshes can be generated automatically for user defined geometries. These meshes are suitable for subsequent FEM or FVM analysis. MESH2D will automatically adapt the element size to ensure that the geometry is adequately resolved, allowing very complex geometries to be meshed with no additional user input. User defined size functions are also available, allowing the user to control mesh resolution if desired. MESH2D is based on an iterative continuous smoothing method, and generally results in very high quality meshes with no small angles and smooth element size variations. WHAT'S NEW IN v2.3? MESH2D v2.3 allows connected polygons to be meshed, with a compatible mesh generated along internal boundaries. MESH 2D v2.3 should also produce meshes with a higher mean element quality in most cases. - MESH2D can produce large scale meshes, with meshes incorporating over 500,000 elements produced based on satellite coastline data. - MESH2D allows fully user defined size functions, giving the user full control over mesh resolution if desired. A template for boundary layer regions is also included. - MESH2D includes functions for uniform and non-uniform mesh refinement, allowing existing meshes to be refined without costly re-triangulation. - MESH2D includes functions to assemble the connectivity information necessary for FE or FV methods. Tested on MATLAB 6.5, 7.0, R2006b and R2007a. DEMO FUNCTIONS: Please use the following functions to launch MESH2D demos: - meshdemo(); - mesh_collection(n); - facedemo(n); MESH2D is (always!) under development, so any problems or suggestions are welcome via email. MESH2D is distributed under the GNU GPL.

Ref.

http://www.mathworks.de/matlabcentral/fileexchange/loadFile.do?objectId=10307&objectType=file

matlab .m DISTMESH A Simple Mesh Generator in MATLAB

HomePage Per-Olof Persson

http://www.mit.edu/~persson/software.html

Software

Meshing:

• DistMesh - A Simple Mesh Generator in MATLAB

Educational MATLAB codes:

• Tridiagonal Eigenvalues in MATLAB - Interface to LAPACK routines for computing eigenvalues of tridiagonal matrices and singular values of bidiagonal matrices
• Level Set Demo - Simple MATLAB scripts for illustration of explicit/implicit interface tracking, reinitialization, and the fast marching method (undocumented, but see presentations for slides and notes).
• fempoisson.m - Solves the Poisson equation on an unstructured grid (square in this example but easy to change) using linear finite elements. Good start to learn about implementation of FEM.
• poiunit.m - Fourier solution of Poisson's equation on the unit line, square, or cube. Good for verification of Poisson solvers, but slow if many Fourier terms are used (high accuracy).
• laplacefft.m - Solve the Laplace equation on a rectangular domain using the FFT. Supports Dirichlet or Dirichlet/Neumann conditions. Contains the following short functions for discrete Sine and Cosine transforms:
  
  • dst.m - Discrete Sine Transform DST-I
  • idst.m - Inverse Discrete Sine Transform IDST-I
  • dct.m - Discrete Cosine Transform DCT-I
  • idct.m - Inverse Discrete Cosine Transform DCT-I
  
  
• Implementation of Finite Element-Based Navier-Stokes Solver

--------------------------end HomePage

http://www-math.mit.edu/~persson/mesh/gallery_images.html

3-D meshes, the left plots show surface meshes and the right plots show cross sections.



-----------

http://people.scs.fsu.edu/~burkardt/m_src/distmesh/distmesh.html

DISTMESH is a MATLAB library which generates and manipulates unstructured meshes in 2D, 3D and general ND. The code is relatively simple, and the user is able to define a variety of geometric shapes, and desired mesh densities.

DISTMESH can be a very quick and flexible means of computing a set of points in a region. However, keep in mind the following flaws:

• Especially if you are have specified some fixed points which must appear in the mesh, it is possible for DISTMESH to return multiple instances of the same point. For finite element applications, in particular, this can result in catastrophe. The program TABLE_MERGE can fix this.
• The nodes produced by DISTMESH are not ordered or sorted in any way whatsoever.
• Because the nodes are not ordered in any way, the triangular elements produced by DISTMESH will typically contain nodes with widely ranging indices. For finite element applications, this can result in a system matrix with an unnecessarily outrageous bandwidth. The program TRIANGULATION_RCM can fix this, after the fact.
• The triangles produced by DISTMESH are not necessarily oriented; they are just as likely to have positive or negative orientation. Some finite element programs insist that all triangles have positive orientation. The program TRIANGULATION_ORIENT can fix this, after the fact.
• When DISTMESH is trying to approximate a boundary, particularly a long straight boundary, it is possible for several points that really belong on the boundary to be slightly out of line. This means that they will be used to form a triangle, of very small area, and terrible conditioning. This can result in perplexing problems near the boundary.
• Because DISTMESH uses tolerances, it is possible for some nodes to lie outside the boundary of the region; it is possible for some triangles on the boundary to lie partially outside the region; it is possible for triangles that lie partly on the boundary, but entirely within the region, to be deleted, leaving a triangular hole.
• Once you have the mesh, you may want to know which nodes lie on the boundary. You'll want this information, for instance, if you need to impose boundary conditions on such nodes. You can get a list of the boundary nodes using the program TRIANGULATION_BOUNDARY_NODES.

Usage:

[ p, t ] = distmesh_2d ( fd, fh, h, box, iteration_max, fixed );

takes:

• fd, the name of a distance function defining the region;
• fh, the name of a mesh density function;
• h, the nominal mesh spacing;
• box, defining a box that contains the region;
• iteration_max, which limits the number of iterations;
• fixed, a list of points which must be included in the mesh.

and returns a triangulation defined by:

• p, a list of node coordinates;
• t, a list of node indices forming triangles;

Related Data and Programs:

DIST_PLOT is a MATLAB program which creates a color contour plot of the distance functions that are used by DISTMESH.

DISTMESH_3D is a MATLAB program which is a subset of the DISTMESH routines, exclusively for 3D problems.

TABLE_IO is a MATLAB library which reads and writes files using the TABLE format; these routines are used by DISTMESH when creating some output files.

TABLE_MERGE is a FORTRAN90 program which removes duplicate points from a TABLE file; it can also remove points that are "close" to each other;

TEST_TRIANGULATION is a MATLAB library which defines some test regions for triangulation.

TRIANGLE is a C program which triangulates a region.

TRIANGULATION_BOUNDARY_NODES is a MATLAB program which reads data defining a triangulation and determines which nodes lie on the boundary.

TRIANGULATION_DISPLAY_OPEN_GL is a C++ program which reads files defining a triangulation and displays an image using Open GL.

TRIANGULATION_L2Q is a MATLAB program which reads data defining a linear triangulation and adds midpoint nodes to create a quadratic triangulation.

TRIANGULATION_MASK is a MATLAB program which is compiled with a user routine that defines a region; it then reads data defining a triangulation and removes all triangles that are outside the region. This is one way to clean up an unconstrained Delaunay triangulation of a nonconvex region.

TRIANGULATION_ORDER3 is a data directory which discusses order 3 triangulations; The node and triangle files output by DISTMESH are an example of such triangulations.

TRIANGULATION_ORIENT is a MATLAB program which reads data defining a triangulation, makes sure that every triangle has positive orientation, and if not, writes a corrected triangle file.

TRIANGULATION_PLOT is a MATLAB program which plots a triangulation.

TRIANGULATION_RCM is a MATLAB program which reads data defining a triangulation, and uses the Reverse Cuthill McKee algorithm to reorder the nodes so as the reduce the bandwidth of the corresponding adjacency matrix. This can be very helpful for cases where the data is to be handled by a frontal technique, or treated as a banded linear system.

TRIANGULATION_REFINE is a MATLAB program which reads data defining a triangulation and creates a refinement of the triangulation by subdividing each triangle.

Author:

DISTMESH is Copyright (C) 2004 Per-Olof Persson.

This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

If you use DISTMESH in any program or publication, please acknowledge its authors by citing the reference.

Reference:

1. http://math.mit.edu/~persson/
   Per-Olof Persson's web site.
2. Per-Olof Persson and Gilbert Strang,
   A Simple Mesh Generator in MATLAB,
   SIAM Review,
   Volume 46, Number 2, June 2004, pages 329-345,
   Available online at ../../pdf/persson_distmesh.pdf

Tar File:

A GZIP'ed TAR file of the contents of this directory is available. This is only done as a convenience for users who want ALL the files, and don't want to download them individually. This is not a convenience for me, so don't be surprised if the tar file is somewhat out of date.

• distmesh_m_src.tar.csh the commands used to create the GZIP'ed TAR file.
• distmesh_m_src.tar.gz the GZIP'ed TAR file.

Source Code:

• boundedges.m finds the surface edges in a triangular mesh.
• circumcenter.m computes the circumcenters of the triangles that form a triangulation of a set of nodes.
• dcircle.m returns the signed distance of one or more points to a circle.
• ddiff.m returns the signed distance of one or more points to a region defined as the set difference of two regions.
• dexpr.m returns the signed distance of one or more points to a region defined by a general symbolic expression.
• dintersect.m returns the signed distance to a region that is the intersection of two regions.
• distmesh_2d.m computes a mesh of a given 2D region.
• distmesh_nd.m computes a mesh of a given ND region.
• dmatrix.m returns the signed distance to a region by interpolation of known distance values on a Cartesian grid.
• dpoly.m returns the signed distance of one or more points to a polygon.
• drectangle.m returns the signed distance of one or more points to a rectangle.
• drectangle0.m, returns the signed distance of one or more points to a rectangle.
• dsegment.cpp (C++ file), a version of the algorithm for the signed distance of one or more points to a set of line segments, for use with a non-Windows version of MATLAB. This file must be compiled, and the corresponding MEX file must be available to define the MATLAB/C++ interface.
• dsegment.dll (binary file), a Windows DLL file, returns the signed distance of one or more points to a set of line segments.
• dsegment.m a pure MATLAB version of DSEGMENT, which returns the signed distance of one or more points to a set of line segments. This routine will be significantly slower than the MEXGLX or DLL versions.
• dsegment.mexglx (binary file), used on non-Windows machines to allow MATLAB to calculate the DSEGMENT algorithm (distance to a line segment) by calling a compiled C++ routine (dsegment.cpp).
• dsphere.m returns the signed distance of one or more points to a sphere.
• dunion.m returns the signed distance to a region that is the union of two regions.
• hmatrix.m computes the mesh size function by interpolation from values specified on a Cartesian grid.
• huniform.m computes a uniform mesh size function.
• meshdemo_nd.m demonstrates the use of the program for higher dimensional problems.
• post_2d.m performs postprocessing for output from DISTMESH_2D.
• protate.m rotates a set of points by a given angle.
• pshift.m shifts a set of points by a given increment.
• r8_epsilon.m returns the R8 arithmetic precision.
• simp_plot_2d.m displays a plot of the triangles that form a mesh in 2D.
• simp_plot_2d_demo.m reads and displays the node and triangle data for each problem.
• simp_plot_3d.m displays a plot of the tetrahedrons that form a mesh in 3D.
• simpqual.m computes the simplex quality of the mesh.
• simpvol.m computes the volume of a simplex.
• surftri.m finds the surface triangles in a tetrahedral mesh.
• timestamp.m prints the current YMDHMS time as a timestamp.
• timestring.m returns the current YMDHMS time as a string.
• triangulation_order3_plot.m writes a PostScript file containing an image of the mesh.
• uniformity.m computes the uniformity of the mesh.

Routines to read and write data to files (borrowed from TABLE_IO) include:

• dtable_data_read.m reads the data of a double precision table file.
• dtable_data_write.m writes the data of a double precision table file.
• dtable_header_read.m reads the header of a double precision table file.
• dtable_header_write.m writes the header of a double precision table file.
• dtable_read.m reads a double precision table file.
• dtable_write.m writes a double precision table file.
• file_column_count.m counts the number of columns in a table file.
• file_row_count.m counts the number of rows in a table file.
• itable_data_read.m reads the data of an integer table file.
• itable_data_write.m writes the data of an integer table file.
• itable_header_read.m reads the header of an integer table file.
• itable_header_write.m writes the header of an integer table file.
• itable_read.m reads an integer table file.
• itable_write.m writes an integer table file.
• s_len_trim.m returns the length of a string to the last nonblank.
• s_word_count.m counts the number of "words" in a string.

----------

DistMesh Function Reference

Back

---

boundedges

Syntax:	e=boundedges(p,t)
Description:	Find all the boundary edges e in triangular mesh p,t.
Comments:	Useful for implementation of boundary conditions for PDE solvers. See surftri for 3-D version.

circumcenter

Syntax:	[pc,r]=circumcenter(p,t)
Description:	Compute the circumcenters pc and the circumradii r for all triangles in the mesh p,t.
Comments:	Not vectorized.

dcircle

Syntax:	d=dcircle(p,xc,yc,r)
Description:	Compute signed distance function for circle centered at xc,yc with radius r.
Comments:

ddiff

Syntax:	d=ddiff(d1,d2)
Description:	Compute signed distance function for set difference of two regions described by signed distance functions d1,d2.
Comments:	Not exactly the true signed distance function for the difference, for example around corners.

dellipse

Syntax:	d=dellipse(p,axes)
Description:	Compute distance from points p to the ellipse centered at the origin with axes=[a,b].
Comments:	C++ code, uses LAPACK for eigenvalue problem.

dellipsoid

Syntax:	d=dellipsoid(p,axes)
Description:	Compute distance from points p to the ellipsoid centered at the origin with axes=[a,b,c].
Comments:	C++ code, uses LAPACK for eigenvalue problem.

dexpr

Syntax:	d=dexpr(p,fin,nit,alpha)
Description:	Compute signed distance function for general implicit expression fin. The parameters nit and alpha have the default values 20 and 0.1.
Comments:	Requires the Symbolic Toolbox, although easy to rewrite to accept derivatives of fin as inputs. The performance is poor, a simple C implementation makes a big difference.

dintersect

Syntax:	d=dintersect(d1,d2)
Description:	Compute signed distance function for set intersection of two regions described by signed distance functions d1,d2.
Comments:	Not exactly the true signed distance function for the intersection, for example around corners.

distmesh2d

Syntax:	[p,t]=distmesh2d(fd,fh,h0,bbox,pfix,fparams)
Description:	2-D Mesh Generator. See other documentation for details on usage.
Comments:

distmeshnd

Syntax:	[p,t]=distmeshnd(fd,fh,h0,bbox,pfix,fparams)
Description:	3-D Mesh Generator. See other documentation for details on usage.
Comments:

dmatrix

Syntax:	d=dmatrix(p,xx,yy,dd)
Description:	Compute signed distance function by interpolation of the values dd on the Cartesian grid xx,yy.
Comments:	xx,yy can be created with meshgrid.

dmatrix3d

Syntax:	d=dmatrix3d(p,xx,yy,zz,dd)
Description:	Compute signed distance function by interpolation of the values dd on the Cartesian grid xx,yy,zz.
Comments:	xx,yy,zz can be created with ndgrid.

dpoly

Syntax:	d=dpoly(p,pv)
Description:	Compute signed distance function for polygon with vertices pv.
Comments:	Uses dsegment and inpolygon. It is usually good to provide pv as fix points to distmesh2d.

drectangle

Syntax:	d=drectangle(p,x1,x2,y1,y2)
Description:	Compute signed distance function for rectangle with corners (x1,y1), (x2,y1), (x1,y2), (x2,y2).
Comments:	Incorrect distance to the four corners, see drectangle0 for a true distance function.

drectangle0

Syntax:	d=drectangle0(p,x1,x2,y1,y2)
Description:	Compute signed distance function for rectangle with corners (x1,y1), (x2,y1), (x1,y2), (x2,y2).
Comments:	See drectangle for simpler version ignoring corners.

dsegment

Syntax:	ds=dsegment(p,pv)
Description:	Compute distance from points p to the line segments in pv.
Comments:	C++ code, used by dpoly.

dsphere

Syntax:	d=dsphere(p,xc,yc,zc,r)
Description:	Compute signed distance function for sphere centered at xc,yc,zc with radius r.
Comments:

dunion

Syntax:	d=dunion(d1,d2)
Description:	Compute signed distance function for set union of two regions described by signed distance functions d1,d2.
Comments:	Not exactly the true signed distance function for the union, for example around corners.

fixmesh

Syntax:	[p,t]=fixmesh(p,t)
Description:	Remove duplicated and unused nodes from p and update t correspondingly. Also make all elements orientations equal.
Comments:

hmatrix

Syntax:	h=hmatrix(p,xx,yy,dd,hh)
Description:	Compute mesh size function by interpolation of the values hh on the Cartesian grid xx,yy.
Comments:	xx,yy can be created with meshgrid. The parameter dd is not used, but included to get a syntax consistent with dmatrix.

hmatrix3d

Syntax:	h=hmatrix3d(p,xx,yy,zz,dd,hh)
Description:	Compute mesh size function by interpolation of the values hh on the Cartesian grid xx,yy,zz.
Comments:	xx,yy,zz can be created with ndgrid. The parameter dd is not used, but included to get a syntax consistent with dmatrix.

huniform

Syntax:	h=huniform(p)
Description:	Implements the trivial uniform mesh size function h=1.
Comments:

meshdemo2d

Syntax:	meshdemo2d
Description:	Demonstration of distmesh2d.
Comments:

meshdemond

Syntax:	meshdemond
Description:	Demonstration of distmeshnd.
Comments:

mkt2t

Syntax:	[t2t,t2n]=mkt2t(t)
Description:	Compute element connectivities from element indices.
Comments:

protate

Syntax:	p=protate(p,phi)
Description:	Rotate points p the angle phi around origin.
Comments:

pshift

Syntax:	p=pshift(p,x0,y0)
Description:	Move points p by (x0,y0).
Comments:

simpplot

Syntax:	simpplot(p,t,expr,bcol,icol)
Description:	Plot 2-D or 3-D mesh p,t. The parameters expr, bcol, icol are only used in 3-D and they have default values.
Comments:

simpqual

Syntax:	q=simpqual(p,t,type)
Description:	Compute qualities of triangular or tetrahedral elements in the mesh p,t. If type==1 (default) the inradius/outradius expression is used. If type==2 a slightly different expression is used.
Comments:

simpvol

Syntax:	v=simpvol(p,t)
Description:	Compute the signed volumes of the simplex elements in the mesh p,t.
Comments:

surftri

Syntax:	tri=surftri(p,t)
Description:	Find all the surface triangles tri in tetrahedral mesh p,t.
Comments:	Used by simpplot. Also useful for implementation of boundary conditions for PDE solvers. See boundedges for 2-D version.

uniformity

Syntax:	u=uniformity(p,t,fh,fparams)
Description:	Computes "uniformity measure", that is, how close the element sizes in the mesh p,t are to the desired mesh size function fh.
Comments:

------------

Examples and Tests:

MESHDEMO_2D runs all the 2D tests:

• meshdemo_2d.m demonstrates the use of the program for 2D problems.
• meshdemo_2d_output.txt the output from a run of MESHDEMO_2D.

P01 is the circle:

• p01_demo.m runs the problem.
• p01_fd.m the signed distance function.
• p01_fh.m a uniform mesh density function.
• p01_mesh.png a PNG image of the mesh.
• p01_nodes.txt the mesh nodes of the problem.
• p01_triangles.txt the mesh triangles of the problem.

P02 is the circle with a hole:

• p02_demo.m runs the problem.
• p02_fd.m the signed distance function.
• p02_fh.m a uniform mesh density function.
• p02_mesh.png a PNG image of the mesh.
• p02_nodes.txt the mesh nodes of the problem.
• p02_triangles.txt the mesh triangles of the problem.

P03 is the square with a hole:

• p03_demo.m runs the problem.
• p03_fd.m the signed distance function.
• p03_fh.m a uniform mesh density function.
• p03_mesh.png a PNG image of the mesh.
• p03_nodes.txt the mesh nodes of the problem.
• p03_triangles.txt the mesh triangles of the problem.

P04 is the hexagon with hexagonal hole:

• p04_demo.m runs the problem.
• p04_fd.m the signed distance function.
• p04_fh.m a uniform mesh density function.
• p04_mesh.png a PNG image of the mesh.
• p04_nodes.txt the mesh nodes of the problem.
• p04_triangles.txt the mesh triangles of the problem.

P05 is the horn:

• p05_demo.m runs the problem.
• p05_fd.m the signed distance function.
• p05_fh.m a nonuniform mesh density function.
• p05_mesh.png a PNG image of the mesh.
• p05_nodes.txt the mesh nodes of the problem.
• p05_triangles.txt the mesh triangles of the problem.

P06 is the superellipse:

• p06_demo.m runs the problem.
• p06_fd.m the signed distance function.
• p06_fh.m a nonuniform mesh density function.
• p06_mesh.png a PNG image of the mesh.
• p06_nodes.txt the mesh nodes of the problem.
• p06_triangles.txt the mesh triangles of the problem.

P07 is the bicycle seat:

• p07_demo.m runs the problem.
• p07_fd.m the signed distance function.
• p07_fh.m a uniform mesh density function.
• p07_mesh.png a PNG image of the mesh.
• p07_nodes.txt the mesh nodes of the problem.
• p07_triangles.txt the mesh triangles of the problem.

P08 is the holey pie slice:

• p08_demo.m runs the problem.
• p08_fd.m the signed distance function.
• p08_fh.m a nonuniform mesh density function.
• p08_mesh.png, a PNG image of the mesh.
• p08_nodes.txt the mesh nodes of the problem.
• p08_triangles.txt the mesh triangles of the problem.

P09 is Jeff Borggaard's square with two hexagonal holes:

• p09_demo.m runs the problem.
• p09_fd.m the signed distance function.
• p09_fh.m a uniform mesh density function.
• p09_mesh.png a PNG image of the mesh.
• p09_nodes.txt the mesh nodes of the problem.
• p09_triangles.txt the mesh triangles of the problem.

P10 is the unit square:

• p10_demo.m runs the problem.
• p10_fd.m the signed distance function.
• p10_fh.m a uniform mesh density function.
• p10_mesh.png a PNG image of the mesh.
• p10_nodes.txt the mesh nodes of the problem.
• p10_triangles.txt the mesh triangles of the problem.

P11 is the L-shaped region:

• p11_demo.m runs the problem.
• p11_fd.m the signed distance function.
• p11_fh.m a uniform mesh density function.
• p11_mesh.png a PNG image of the mesh.
• p11_nodes.txt the mesh nodes of the problem.
• p11_triangles.txt the mesh triangles of the problem.

P12 is the John Shadid's H-shaped region:

• p12_demo.m runs the problem.
• p12_fd.m the signed distance function.
• p12_fh.m a uniform mesh density function.
• p12_mesh.png a PNG image of the mesh.
• p12_nodes.txt the mesh nodes of the problem.
• p12_triangles.txt the mesh triangles of the problem.

P13 is the Sandia fork:

• p13_demo.m runs the problem.
• p13_fd.m the signed distance function.
• p13_fh.m a uniform mesh density function.
• p13_mesh.png a PNG image of the mesh.
• p13_nodes.txt the mesh nodes of the problem.
• p13_triangles.txt the mesh triangles of the problem.

P14 is Marcus Garvie's Lake Alpha, with Beta Island:

• p14_demo.m runs the problem.
• p14_fd.m the signed distance function.
• p14_fh.m a nonuniform mesh density function.
• p14_mesh.png a PNG image of the mesh.
• p14_nodes.txt the mesh nodes of the problem.
• p14_triangles.txt the mesh triangles of the problem.

P15 is Sangbum Kim's forward step region:

• p15_demo.m runs the problem.
• p15_fd.m the signed distance function.
• p15_fh.m a nonuniform mesh density function.
• p15_mesh.png a PNG image of the mesh.
• p15_nodes.txt the mesh nodes of the problem.
• p15_triangles.txt the mesh triangles of the problem.

P16 is Kevin Pond's elbow, a quarter of a circular annulus:

• p16_demo.m runs the problem.
• p16_fd.m the signed distance function.
• p16_fh.m the mesh density function.
• p16_mesh.png a PNG image of the mesh.
• p16_nodes.txt the mesh nodes of the problem.
• p16_triangles.txt the mesh triangles of the problem.

P17 is a rectangular region with a Reuleaux triangle obstacle.

• p17_demo.m runs the problem.
• p17_fd.m the signed distance function.
• p17_fh.m the mesh density function.
• p17_mesh.png a PNG image of the mesh.
• p17_nodes.txt the mesh nodes of the problem.
• p17_triangles.txt the mesh triangles of the problem.