Archive for the 'soft_interessant' Category



SPM-Haemodynamic Response Function

Generally speaking, cognitive processing is associated with increases in neuronal firing rates. The increased neural activity lead to increased metabolic requirements for the neurons. The onset of neural activity leads to a systematic series of physiological changes in the local network of blood vessels that include changes in the cerebral blood volume per unit of brain tissue (CBV), changes in the rate of cerebral blood flow (CBF), and changes in the concentration of oxyhaemoglobin and deoxyhaemoglobin.

There are different fMRI techniques that can pick up a functional signal corresponding to changes in each of the previously mentioned components of the haemodynamic response. The most common functional imaging signal is the Blood Oxygenation Level Dependent signal (BOLD), which primarily corresponds to the concentration of deoxyhaemoglobin. In simple terms, the magnetic resonance signal comes from exciting hydrogen nuclei with a radiofrequency pulse, and detecting the radio waves emitted as the nuclei return to a lower-energy configuration. Deoxyhaemoglobin has different magnetic properties than oxyhaemoglobin– it is paramagnetic, which means that it will make the local magnetic field over a microscopic domain inhomogenous. This has the effect of dephasing the signal emitted by the nuclei in this domain, causing destructive interference in the observed MR signal. Over a macroscopic domain (i.e., one functional voxel) greater amounts of deoxyhaemoglobin lead to less signal. The functional BOLD signal is seen as an increase in the MR signal that corresponding to a decrease in the concentration of deoxyhaemoglobin. The decrease of deoxy-Hb is seen because the increase in CBF following neural activity more than accounts for the effect of increased uptake of oxygen.

Image:spm_hrf.png

For the purposes of estimating the BOLD signal in an experimental paradigm, SPM makes use of a canonical haemodynamic response function (HRF). This function is assumed to be the response of the system (as reflected by the MR signal) to a brief, intense period of neural stimulation. The SPM HRF is shown above, and exhibits a rise peaking around 6 sec, followed by an undershoot that persists for a considerable period. The code for this graph is below.

>> RT = 1; hrf = spm_hrf(RT); plot(0:RT:32, hrf);

In this graph, the y-axis is in arbitrary units. A common way to plot the impulse response is in units of percent signal change from a baseline condition. A very robust stimulus (such as a contrast taken between a flickering visual stimulus and no visual stimulus) may produce changes on the order of 2%-4% in the BOLD signal. The change observed in contrasts involving higher-level cognitive processes is typically much smaller.

Ref. http://en.wikibooks.org/wiki/SPM/Haemodynamic_Response_Function

Boot camp, virtual box, VM ware fusion, parallels desktop

Dans notre série « Comment faire tourner Windows sur un Mac sans redémarrer avec BootCamp ?« ,

poursuivons notre petit tour des solutions de virtualisation !

c’est vraiment stupéfiant d’avoir les 2 mondes en basculant d’une fenêtre à une autre. On a une interopérabilité géniale.

Le mode fusion permet d’avoir l’appli dans mac os comme si c’etait natif!

Après Parallels et VMWare , voir l’analyse de janvier2009 sur VirtualBox.
Initialement développée par Innotek, une entreprise allemande rachetée par Sun en février 2008, cette solution de virtualisation gratuite se pose en concurrent sérieux de Parallels Desktop et de VMWare Fusion.

http://www.commentcamarche.net/telecharger/telecharger-3673479-virtualbox

http://www.virtualbox.org/wiki/Downloads

Petite visite en vidéo (et en français) de la VirtualBox avec Achim Hasenmueller, co-fondateur historique d’Innotek, aujourd’hui ingénieur chez Sun microsystems.

C’est bien d’avoir un tel programme gratuit, par contre il y as encore quelque progre a faire pour etre au niveau de Fusion 2.

C’est clair que niveau finition et intégration système, il ne joue pas dans la même cours que ses concurrents payant (même pas de drag & drop).
Ceci dit, il répond parfaitement à mes besoins en virtualization. Je trouve que niveau perf il s’en sort très bien.

VirtualBox permet d’émuler complètement un PC. C’est comme si vous aviez un second PC dans une simple fenêtre.
C’est utile pour tester d’autres système d’exploitation sans repartitionner et sans risque (repartitionner par exemple Linux), pour naviguer en toute sécurité ou pour tester un logiciel sans risque de rendre son système d’exploitation instable.

Vous pouvez créer autant de machines virtuelles que vous le souhaitez, et installer tous les systèmes d’exploitation que vous le voulez dedans.
Il est possible de définir – pour chaque machine virtuelle – combien de mémoire elle possède, de disque dur, si elle aura accès aux ports USB, au réseau, à la carte son, etc.

VirtualBox contient un gestionnaire de disques qui vous permet de créer des disques virtuel sous forme de fichiers .vdi qui apparaîtront comme de vrais disques dans les machines virtuelles.
Cela vous permet donc de « créer » à volonté des disques, et cela sans jamais avoir à repartitionner votre disque dur.

Vous pouvez également utiliser directement des images ISO de CD et DVD, ce qui permet de tester des distributions Linux sans avoir à les graver.

Cerise sur le gâteau, VirtualBox possède un serveur RDP intégré, ce qui permet de démarrer une machine virtuelle sur un ordinateur, et utiliser cette machine virtuelle à partir d’un autre ordinateur.

——-

Comparaison des soft sur un paln techno

en particulier sur un macpro intel x86  tiger  osx 10.5.x

Version du système :    Mac OS X 10.5.6 (9G55)
Version du noyau :    Darwin 9.6.0

VMWare est plus rapide dans l’ensemble, notamment au démarrage mais Parallel est plus confortable à utiliser.

Les outils VMWare n’affecte pas le système windows (au cas où c’est via bootcamp) en revanche parallel s’implique un peu plus. On va dire que VMWare reste à la surface et pas Parallel.

dans une optique « serveur ». virtualBox s’en sort très bien! Même beaucoup mieux que Fusion et Parrallels. Son mode en ligne de commande, et le fait de pouvoir démarrer les VM en arrière plan (sans fenêtre)…

—–

qq remarques:

les logiciels de virtualisation prennent beaucoup de mémoire vive, surtout si les OS invités en demandent beaucoup (surtout Vista argh). Donc le maitre-mot, c’est la ram. Il faut au strict minimum avoir 2 Go de ram.

les logiciels de virtualisation sont incapables de gérer la CG hôte, et n’émulent qu’une CG de base, par exemple les effets 3D sont très mal gérés.

——

Je viens fraîchement d’installer la dernière version de Virtual Box et je me demandais comment utiliser la partition Windows creer par bootCamp.

Auparavant j’utilisait Parrallels Desktop, lors de la configuration j’ai simplement eu à préciser ma partition windows, il me la toute suite reconnue et windows à démarré.

Virtual Desktop me demande de creer une image de cette partition qui bien sur me consommera de la place, est t’il possible d’utiliser directement la partition windows sans passer par la création d’une image ?

—–

Je n’ai jamais essayé VMWare par j’utilise paralelle et auparavant virtual box
Le mode cohérence est sympa même si je ne l’utilise que rarement, en effet j’ai installé windows pour un logiciel que j’utilise seul qui n’a pas de lien avec d’autres applications.
details marrant on a accès à la logithèque mac via le menu démarrer de win.

J’en profite pour glisser quelque mots sur virtual box .
Ce logiciel gratuit est vraiment bien construit et facile à prendre en main il n’est pas à grand chose des logiciel payants.
Cependant son interface est moins bien fini que paralelle, ne paramètre pas le réseau pour avoir accès aux dossier mac automatiquement, nécessite un plugin pour avoir le passage de souris de win à OSX ,de plus le mode cohérence ou fusion n’existe pas.
Niveau vitesse il est un peu plus lent mais ce n’est qu’une sensation.

une mise à jour de virtualbox de 2008. http://www.macgeneration.com/news/voir/131628/sun-lance-sa-virtualbox-en-version-2.0

SPM version 8b statistical parametric mapping MATLAB

les 2 liens les plus « usefull »:

http://www.fil.ion.ucl.ac.uk/spm/software/spm8b/#Introduction

http://www.fil.ion.ucl.ac.uk/spm/doc/intro/

http://www.fil.ion.ucl.ac.uk/spm/ext/

The SPM approach in brief

The Statistical Parametric Mapping approach is voxel based:

  • Images are realigned, spatially normalised into a standard space, and smoothed.
  • Parametric statistical models are assumed at each voxel, using the General Linear Model GLM to describe the data in terms of experimental and confounding effects, and residual variability.
  • For fMRI the GLM is used in combination with a temporal convolution model.
  • Classical statistical inference is used to test hypotheses that are expressed in terms of GLM parameters. This uses an image whose voxel values are statistics, a Statistic Image, or Statistical Parametric Map (SPM{t}, SPM{Z}, SPM{F})
  • For such classical inferences, the multiple comparisons problem is addressed using continuous random field theory RFT, assuming the statistic image to be a good lattice representation of an underlying continuous stationary random field. This results in inference based on corrected p-values.
  • Bayesian inference can be used in place of classical inference resulting in Posterior Probability Maps PPMs .
  • For fMRI, analyses of effective connectivity can be implemented using Dynamic Causal Modelling DCM.

LiberKey = 200 logiciels disponibles pour votre clé USB

200 logiciels disponibles pour votre clé USB

La LiberKey regroupe de nombreux logiciels gratuits. Véritable « couteau suisse », elle s’utilise aussi bien en local qu’en nomade en téléchargeant les logiciels sur une clé USB.

Parmi les logiciels libres ou gratuits, les meilleurs ont été sélectionnés pour répondre pratiquement à tous les besoins imaginables: OpenOffice, Firefox, Filezilla, Gimp, VLC, …

Vous retrouvez alors, instantanément, l’intégralité de vos logiciels sans aucune installation sur votre PC.

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

freeware 2D 3D Data graphics, code, mac OSX

http://www.vvidget.org/builder/description/screenshots.html

Vvidget Builder is a powerful point-and-click graph building tool which enables you to fine tune, add graphics, buildup specialized graphs and layout entire pages with multiple figures. Palettes provide immediate drag-and-drop access to the charts you use most often and because you can define your own palettes you can drag and drop your new graphs and graphics in addition to the premade ones.

With Vvidget Builder your graphs are only limited by your imagination because you are given complete freedom and are empowered with the capability to layout results the way you want them.

Features available with Vvidget Builder include:

Feature Description
General Graphics Line, Rectangle, Oval, Parallelogram, Circle, Ellipse, Pie Wedge, Curve, Polygon, Cubic Bezier Sections, Path, Label and Image graphics.
Data Graphics Function, Scatter, Trajectory, Point Map, 3D Surface and 3D Scatter graphics.
Graphs Linear, 2-Y Linear, Log, Semi-Log (X and Y directions), Log-Log, Full-Cycle and Sub-Cycle Log, Polar, R-Log Polar, Date, 2-Y Date, Date-Log and 3D Rectilinear Perspective.
Chart Wizards Chart wizards take data from tables and make graphs. You can make a simple line graph to complex polar-point-map graphs. Use the optional VVidget Code programmable system to automate and animate the charts found on the wizards.
General Attributes Shadow, diffusion, gradients and standard fill and stroke type effects.
Real-Time Attributes edit and change in real time. 3D graphics rotate and update in realtime.
Data Attributes Dashes, markers and labels. Unlimited capabilities to make new markers and end caps.
Graph Attributes There are a large amount of tunable graph attributes including grid and subgrid color, width, dash, limit values, tick and subtick placements, color, length, width, label format, frame and background color, title values and placement, font type and size, label angles, offsets and gaps.
Editing Point-wise vertex and spline knot editing, smoothing, a large quantity of direct attribute inspector editing, rotate, skew, copy, paste, delete, etc.
Palettes Includes standard palettes for drag and drop creation and the ability to create your own palettes.
Export and Import Export figures to PDF, print them, or copy the data directly. Importing data is as simple as pasting a space delimited sequence of numeric text values. Export high resolution images into your favorite applications like Power Point.
Documents Stores results in a Vvidget Builder document.
High-Level
Design
You are never trapped into a limited design. For example, you can have any marker because markers are graphics and the editing and graphic design is recursive. When you use Vvidget Builder you are exercising the tools of VVidget Code, VVidget Pro and the Peer Visual server so your knowledge is applicable to more sophisticated applications.

open source and scientific 3D and utilities cmake

VTK The Visualization ToolKit (VTK) is an open source, freely available software system for 3D computer graphics, image processing, and visualization used by thousands of researchers and developers around the world.
ITK The Insight Segmentation and Registration Toolkit (ITK) is an open-source software system to support the VisibleHuman Project. Currently under active development, ITK employs leading-edge segmentation and registration algorithms in two, three, and more dimensions.
CMake CMake, a cross-platform, open-source make system is used to control the software compilation process using simple platform and compiler independent configuration files. CMake generates native makefiles and workspaces that can be used in the compiler environment of your choice.
CDash CDash, an open source, web-based software testing server. CDash aggregates, analyzes and displays the results of software testing processes submitted from clients located around the world.
ParaView ParaView is an open-source, multi-platform, parallel scientific visualization application. Although it can be run on a single processor, ParaView is specifically designed to run on distributed parallel computers.

Ref.

http://www.kitware.com/opensource/opensource.html


I MOVED THIS BLOG FROM WORDPRESS TO BLOGGER. Ce blog est à
ex-ample.blogspot.com

Blog Stats

  • 221 100 hits

localization

Flickr Photos

juin 2019
L M M J V S D
« Oct    
 12
3456789
10111213141516
17181920212223
24252627282930