Hello All,
> On May 25, 2012, at 9:56 PM, Luis Emiro Linares Garcia wrote:
>
> > Jack,
> >
> > Would you mind to share the details about your fitting code in Mathematica?
I have matlab code that does essentially the same:
Given data file 3 columns: wavelength [m], Re(N), Im(N) where N - complex
refractive index; Fits the data using guess vector and the limits vectors. Only
one specific fit function from matlab is used: lsqnonlin. See the code below:
_________________________
WBR, Serge
_________________________
% Main function: %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ------------------------------------------------------------------
% Fitting of measured dispersion (n,k files) to get Meep dispersion parameters.
% ------------------------------------------------------------------
% Parameters:
a=1e-6; % Meep unit in [m]
NKfile = 'ITO.nk'; % File with wl, n,k
% Fit region:
wl_l=350; %lower wavelength, [nm]
wl_h=1000; %upper wavelength, [nm]
only_plot=0; % set to 1 to only plot the fit
% Initial fit vector (guess):
% Format of guess and low, upper boundaries (lb,ub):
% guess: ['constant epsilon'
% 'polarization 1 (sigma omega gamma: 3 numbers)'
% 'polarization 2'
% 'polarisation N'
% 'electric conductivity']
guess = [4.4517 ...
2.5012 0.361 0.0274 ...
1.3517 2.4131 0.4832 ...
0.5522 2.1597 0.4142 ...
0.00001]; % Initial fit vector
lb = [1 ...
0 0.1 0 ...
0]; % low boundary
ub = [20 ...
1000 10 5 ...
0.0001] ; % high --"--
% End of parameters ------------------------------------------------
% ------------------------------------------------------------------
% Calculates p - MEEP dispersion parameters, the same format as guess.
% Note: 'guess' is used when 'p' is not exists, so delete 'p' from matlab
% variables to start again with 'quess'.
% ------------------------------------------------------------------
% ------------------------------------------------------------------
aSi=dlmread(NKfile,'',0,0); % read file with wl, n, k
wl=aSi(:,1);
wl = wl/a; % change from metres to Meep units
wl_l= wl_l/1e9/a; % change from nm to Meep units
wl_h= wl_h/1e9/a;
for j=1:length(wl) % cut [wl_l:wl_h] range
if wl(j)>=wl_l
b=j;
break
end
end
for j=1:length(wl)
if wl(j)>=wl_h
wl=wl(b:j);
n=aSi(b:j,2);
k=aSi(b:j,3);
break
end
end
f = 1./wl; %freq = 1/wavelegth
N = n + 1i*k; % complex refr. index
Eps = N.^2; % complex dielectic permittivity
f1 = max(f):((min(f)-max(f))/100):min(f);
Eps = interp1(f,Eps,f1,'spline');
ImReRatio = mean(real(Eps))/mean(imag(Eps)); % Im/Re ratio in fit function
fcen=(1/wl_l+1/wl_h)/2; % pulse parameters for Meep
df=(1/wl_l-1/wl_h);
fprintf(';----------------------------------------------------\n');
fprintf('(define-param fcen %f)\n', fcen);
fprintf('(define-param df %f)\n', df);
fprintf(';----------------------------------------------------\n');
switch length(guess)
case {0,1,2,3}
disp('Wrong guess length!')
otherwise
if(not(exist('p','var')))
p = guess;
end
p = Lorenzian_fit_N_Econd(f1, Eps, p, lb, ub, only_plot, ImReRatio);
fprintf(';----------------------------------------------------\n');
fprintf('; fitted "%s" for wl=[%g:%g]*%g [m]\n', NKfile, wl_l, wl_h, a);
fprintf('(define TCO (make medium (epsilon %g) (D-conductivity
%g)\n',...
p(1), p(length(p)));
fprintf('\t(E-polarizations\n');
for n=1:((length(p)-2)/3)
fprintf('\t\t(make polarizability (sigma %g) (omega %g)
(gamma %g))\n',...
p(3*(n-1)+1+1), p(3*(n-1)+1+2), p(3*(n-1)+1+3));
end
fprintf(')))\n');
fprintf(';----------------------------------------------------\n');
end
% Fit function: -----------------------------------------------------------
function [p, fEps] = Lorenzian_fit_N_Econd(freq, epsilon, p, lb, ub, only_plot,
ImReRatio)
%#codegen
%p(1) - epsilon_infinity
%p(2),p(5) - Delta epsilon_n : S(n)
%p(3),p(6) - rez. freq. : w(n)
%p(4),p(7) - extin.gamma : y(n)
%p(last) - conductivity : d_cond
coder.extrinsic('int2str', 'strcat', 'inline', 'optimset', 'lsqnonlin',...
'fprintf', 'figure', 'plot');
lrntz = (['(1+i*p(' int2str(length(p)) ')./(2*pi*f)).*(p(1)']);
lbl = zeros(length(p),1);
ubl = zeros(length(p),1);
lbl(1) = lb(1);
ubl(1) = ub(1);
for n=1:((length(p)-2)/3)
% if (length(lb) < 3*n) % expand lb to length(p)
lbl(3*(n-1)+1+1) = lb(2);
% lb(3*(n-1)+1+2) = p(3*(n-1)+1+2);
lbl(3*(n-1)+1+2) = lb(3);
lbl(3*(n-1)+1+3) = lb(4);
% end
% if (length(ub) < 3*n) % expand ub to length(p)
ubl(3*(n-1)+1+1) = ub(2);
% ub(3*(n-1)+1+2) = p(3*(n-1)+1+2)*1.001;
ubl(3*(n-1)+1+2) = ub(3);
ubl(3*(n-1)+1+3) = ub(4);
% end
lrntz = strcat(lrntz,'+','p(',int2str(3*(n-1)+1+1),')*p(',...
int2str(3*(n-1)+1+2),')^2./(p(',int2str(3*(n-1)+1+2),...
')^2-f.^2-i*f*p(',int2str(3*(n-1)+1+3),'))');
end
lbl(end) = lb(end);
ubl(end) = ub(end);
lrntz = strcat(lrntz,')');
lorentzian=inline(lrntz, 'p','f');
if (not(only_plot))
lrntz = strcat('[real(',lrntz,')-re,',int2str(ImReRatio),'*(imag(',lrntz,')-
im)]');
% lrntz = strcat('[real(', lrntz, ')./re-1,imag(', lrntz, ')./im-1]');
fitfunc=inline(lrntz,'p','f','re','im');
options = optimset('MaxFunEvals',100000,'MaxIter',10000);
[p,resnorm]=lsqnonlin(fitfunc,p, lbl,ubl, options,
freq,real(epsilon),imag(epsilon));
err = resnorm/(2*length(freq))^2; % print the residual at x:
sum(fun(x).^2)/(2*npoints)^2.
fprintf('ResNorm/(2*Nfreqs)^2=%g\n', err);
end
f= min(freq):((max(freq)-min(freq))/100):max(freq);
Eps = interp1(freq,epsilon,f);
ReEps = real(Eps);
ImEps = imag(Eps);
%ReEps = interp1(freq,real(epsilon),f);
%ImEps = interp1(freq,imag(epsilon),f);
EpsFit = lorentzian(p,f);
maxStd = max(std(Eps-EpsFit));
createfigureEpsFitting(1, f, [ReEps;ImEps.*ImReRatio], ...
f, [real(lorentzian(p,f));imag(lorentzian(p,f)).*ImReRatio],...
ImReRatio) % Eps
createfigureEpsFitting_wl(2, 1./f, [ReEps;ImEps.*ImReRatio], ...
1./f, [real(lorentzian(p,f));imag(lorentzian(p,f)).*ImReRatio],...
ImReRatio, maxStd) % Eps
figure(3)
plot(1./freq,real(sqrt(epsilon)),'bo', 1./f,real(sqrt(lorentzian(p,f))),'b');
hold on
plot(1./freq,imag(sqrt(epsilon)),'ro', 1./f,imag(sqrt(lorentzian(p,f))),'r');
hold off
ylabel('n & k')
xlabel('wavelength')
legend('meas n','fit n','meas k','fit k',0)
fEps = [f; lorentzian(p,f)].';
%fnk = flipud([f; real(sqrt(lorentzian(p,f))); imag(sqrt(lorentzian(p,f)))].');
end
% Plotting functions: ------------------------------------------------------
function createfigureEpsFitting(nFig, X1, YMatrix1, X2, YMatrix2, ImReRatio)
% X1: vector of x data
% YMATRIX1: matrix of y data
% X1: vector of x data
% YMATRIX2: matrix of y data
% Auto-generated by MATLAB on 04-Apr-2011 11:46:52
% Create figure
figure1 = figure(nFig);
% Create axes
axes1 = axes('Parent',figure1,'YGrid','on','XGrid','on',...
'Units','centimeters','Position',[1.3 1.25 7 7]);
box(axes1,'on');
hold(axes1,'all');
% Create multiple lines using matrix input to plot
plot1 = plot(X1,YMatrix1,'Parent',axes1,'LineWidth',2,'LineStyle','--');
set(plot1(1),'MarkerSize',5,'Color',[0 0 1],'DisplayName','\Re(\epsilon)');
set(plot1(2),'Color',[1 0 0],'DisplayName',['\Im(\epsilon)*'
num2str(ImReRatio)]);
% Create multiple lines using matrix input to plot
plot2 = plot(X2,YMatrix2,'Parent',axes1,'LineWidth',2);
set(plot2(1),'Color',[0 0 1],'DisplayName','fitted \Re(\epsilon)');
set(plot2(2),'Color',[1 0 0],'DisplayName',['fitted \Im(\epsilon)*'
num2str(ImReRatio)]);
%plot3 = plot(X2,abs(YMatrix2(2,:)-YMatrix1(2,:)),'Parent',axes1,'LineWidth',2);
%set(plot3(1),'Color',[0 1 1],'DisplayName','err \Im(\epsilon)');
% Create ylabel
ylabel('\Re(\epsilon) and \Im(\epsilon)','FontSize',10);
% Create xlabel
xlabel('scaled frequency, [1/\mum]','FontSize',10);
% Create legend
legend1 = legend(axes1,'show');
set(legend1,'Units','centimeters','Position',[1.4 5.6 3.4 2.5],'FontSize',10);
% Create textbox
annotation(figure1,'textbox',[0.02 0.024 0.1 0.07],...
'String',{'b)'}, 'FontSize',10, 'FitBoxToText','off','LineStyle','none');
hold(axes1, 'off');
function createfigureEpsFitting_wl(nFig, X1, YMatrix1, X2, YMatrix2, ImReRatio,
maxStd)
% X1: vector of x data
% YMATRIX1: matrix of y data
% X1: vector of x data
% YMATRIX2: matrix of y data
% Create figure
figure1 = figure(nFig);
% Create axes
axes1 = axes('Parent',figure1,'YGrid','on','XGrid','on',...
'Units','centimeters','Position',[1.3 1.25 7 7]);
box(axes1,'on');
hold(axes1,'all');
% Create multiple lines using matrix input to plot
plot1 = plot(X1,YMatrix1,'Parent',axes1,'LineWidth',2,'LineStyle','--');
set(plot1(1),'MarkerSize',5,'Color',[0 0 1],'DisplayName','\Re(\epsilon)');
set(plot1(2),'Color',[1 0 0],'DisplayName',['\Im(\epsilon)*'
num2str(ImReRatio)]);
% Create multiple lines using matrix input to plot
plot2 = plot(X2,YMatrix2,'Parent',axes1,'LineWidth',2);
set(plot2(1),'Color',[0 0 1],'DisplayName','fitted \Re(\epsilon)');
set(plot2(2),'Color',[1 0 0],'DisplayName',['fitted \Im(\epsilon)*'
num2str(ImReRatio)]);
%plot3 = plot(X2,abs(YMatrix2(2,:)-YMatrix1(2,:)),'Parent',axes1,'LineWidth',2);
%set(plot3(1),'Color',[0 1 1],'DisplayName','err \Im(\epsilon)');
% Create ylabel
ylabel('\Re(\epsilon) and \Im(\epsilon)','FontSize',10);
% Create xlabel
xlabel('wavelength, [\mum]','FontSize',10);
% Create legend
legend1 = legend(axes1,'show');
set(legend1,'Units','centimeters','Position',[1.4 5.6 3.4 2.5],'FontSize',10);
% Create textbox
annotation(figure1,'textbox',[0.02 0.024 0.1 0.07],...
'String',{'b)'}, 'FontSize',10, 'FitBoxToText','on','LineStyle','none');
annotation(figure1,'textbox',[0.12 0.124 0.1 0.07],...
'String',{['max(std)=' num2str(maxStd,3)]}, 'FontSize',10,
'FitBoxToText','on','LineStyle','none');
hold(axes1, 'off');
__________________________
_______________________________________________
meep-discuss mailing list
meep-discuss@ab-initio.mit.edu
http://ab-initio.mit.edu/cgi-bin/mailman/listinfo/meep-discuss
