how to make this file work "Error using unique"

1 view (last 30 days)
single_peaks=1:8; %single peaks to fit
double_peaks=[];
plot_fits=1; %set to 1 to plot raw and fitted data
azi_bin=360; %number of pre-existing eta bins to put in a single bin for fitting (1 <=azi_bin<= eta.N)
dx=8; %number of radial points on either side of guessed peak position to evaluate for fit
bx=1; %number of background points on either side to use in fitting poly background
fwhmguess=2.5; %guess for peak fwhm, in pixels
maxdrift=3; %maximum range in pixels to find the y max (keep tight as possible!)
peak_tol=0.5; %minimum peak intensity (background subtracted, evaluated over range 2*maxdrift)at which peaks will be fit - otherwise set as 1e-6
below_tol_val=999; %value to give all parameters for peaks which are below intensity tolerance
delta=0; %detector offset parameter, keep at 0 unless you have a good reason not to!
de=0; %strain amplitude (guess)
% crystallographic parameters, hkls to fit and resulting peak position guesses
a1=3.054; % refined from NiTi (austenite_B2)lattice parameters from pdf #18-0899
%a1=3.054; % refined from NiTi (austenite_B2)lattice parameters from pdf #18-0899
% % % NiTiCu B19 lattice parameters from pdf Potatov et. al.
a2=2.919;
b2=4.288;
c2=4.500;
%alpha=1.5708; %90 degrees
%beta=1.5708; %90 degrees
%gamma=1.5708; %90 degrees
% % % NiTiPd B19 lattice parameters from pdf Potatov et. al.
% a3=2.811;
% b3=4.421;
% c3=4.694;
% alpha=1.5708; % 90 degrees
% beta=1.5708; % 90 degrees
% gamma=1.5708; % 90 degrees
% a3=2.885; b3=4.622; c3=4.120; alpha=1.5708; beta=1.68948; gamma=1.5708; % NiTi (mart) B19 lattice parameters from pdf #35-1281
a3=2.898; b3=4.108; c3=4.646; %alpha=1.5708; beta=1.7066; gamma=1.5708; % NiTi (mart) B19' lattice parameters from pdf #03-065-0145
h=[1 0 1 0 1 1 1 1];
k=[0 1 1 0 1 1 2 1];
l=[0 1 0 2 1 2 0 2];
%h=[1 1 0 1 2 2];
%k=[0 1 1 1 0 1];
%l=[0 0 2 1 0 1];
for i=1:1:4; cen0(i)=round(cen0calc('cubic',h(i),k(i),l(i),a2,a2,a2,center.z,energy,pix.size,delta)); end; % NiTi (austenite) cubic peaks
%for i=5:1:11; cen0(i)=round(cen0calcmono('mono',h(i),k(i),l(i),a3,b3,c3,alpha,beta,gamma,center.z,energy,pix.size)); end; % NiTi (martensite) monoclinic peaks
%for i=1:1:6; cen0(i)=round(cen0calc('cubic',h(i),k(i),l(i),a1,a1,a1,center.z,energy,pix.size,delta)); end; % NiTi (austenite) cubic peaks
%for i=2:1:8; cen0(i)=round(cen0calcmono('mono',h(i),k(i),l(i),a3,b3,c3,alpha,beta,gamma,center.z,energy,pix.size)); end; % NiTi (martensite) monoclinic peaks
%%calculate arrays based on user defined parameters
nPeak = length(cen0); nDouble=length(double_peaks);
nAzi=size(imm,1)/azi_bin; %bin data using azi_bin
azi_deg=[0.5:nAzi-0.5]*(eta.max-eta.min)/eta.N*azi_bin; %mean values of azimuths for fit, in degrees
%%fit singlet peaks
opt = optimset('disp','off','large','on','jacobi','on','tolx',1e-4,'tolf',1e-4,'maxi',120);
fn = {'backg1' 'psv1'}; npar = [2 4]; %define peakfitting choice (PV here)
low = [-inf 0 0 0 0 0]; up = [inf inf inf inf inf 1]; %set limits for fitting parameters
pfix = [nan nan nan nan nan nan]; % fix fitting parameters (currently none are fixed)
if ~isempty(single_peaks); fprintf('\nFitting single peaks for following azimuth bins:\n'); end;
for azi = nAzi:-1:1
if ~isempty(single_peaks); fprintf('%4i',azi*azi_bin);
for peak = single_peaks;
cen0a(peak,azi)=round(cen0(peak)*(1+de/2*cos(2.*azi_deg(azi)/57.3))); %azi-dependent starting guess
x = cen0a(peak,azi) - dx: cen0a(peak,azi) + dx; Nx=length(x);
yraw = mean(imm([azi*azi_bin:-1:azi*azi_bin-azi_bin+1], x-rho.min));
xback=x([1:2 Nx-1:Nx]); yback=yraw([1:2 Nx-1:Nx]); %x,y values to use for background determination
[pback_pars,s,mu]=polyfit(xback,yback,1); % 1st order polynomial (linear) background fit
pback=polyval(pback_pars,x,[],mu); %mu are centering and scaling parameters to help polyfit
y=yraw-pback+1000;
ymax= max(max(y)); xmax = x(find(round(y)==round(ymax))); xmax = xmax(1); ymin=min(min(y));
clear x pback_pars pback y yraw xback yback;
x = xmax - dx: xmax + dx; Nx=length(x); % shift data so xmax is center point
x = cen0a(peak,azi) - dx: cen0a(peak,azi) + dx; Nx=length(x); %do not shift data
yraw = mean(imm((azi*azi_bin:-1:azi*azi_bin-azi_bin+1), x-rho.min));ymin=min(min(yraw));
xback=x([1:2 Nx-1:Nx]); yback=yraw([1:2 Nx-1:Nx]); %x,y values to use for background determination
[pback_pars,s,mu]=polyfit(xback,yback,1); % 1st order polynomial (linear) background fit
pback=polyval(pback_pars,x,[],mu); %mu are centering and scaling parameters to help polyfit
y=yraw-pback+1000;
ymax = max(max(y((length(y)-1)/2-maxdrift:(length(y)-1)/2+maxdrift)));
xmax = x(find(round(y)==round(ymax))); xmax = xmax(1); ymin=min(min(y));
%check if peak intensity is above threshold, to decide if this particular azimuth will be fitted
fom=ymax-1000;
if plot_fits
figure(6);clf;
plot(x,yraw,'b+-',x,pback,'r-');
xlabel('radial position, pixels');ylabel('intensity');
title(sprintf('Peak %i and azimuth %i of %s; blue=data, red =background',peak,azi,filename)); grid on;
if (fom>=peak_tol); text(0.05,0.95,sprintf('Will fit as FOM is %4.2f while tolerance is %4.2f',fom,peak_tol),'sc'); pause;
else text(0.05,0.95,sprintf('Will NOT fit as FOM is %4.2f while tolerance is %4.2f',fom,peak_tol),'sc'); pause; end;
end; %plotting raw data and background
if (fom>=peak_tol); % do the fit as intensity is above threshold
p0 = [0, ymin, (ymax-ymin), xmax, fwhmguess, .5];
y0 = sumfun1(p0,pfix,npar,fn,x);
[pfit,rn,rs,ex,out,lam,jac] = lsqnonlin('sumfun1',p0,low,up,opt,pfix,npar,fn,x,y,ones(size(y)));
[~, var] = confint(pfit,rs,jac); % variance of fitted params
var=full(var);
yfit = sumfun1(pfit,pfix,npar,fn,x);
[~, dum, intint]=psv1([pfit(3:end)],x);
if plot_fits
figure(7);clf; subplot('position',[0.1 0.3 0.85 0.6]);
plot(x,y0,'-r',x,y,'+',x,yfit,'-b'); ylabel('Intensity+100(offset)');
title(sprintf('Fit of peak %i and azimuth %i of %s; symbols=data, blue=fit, red =guess',peak,azi,filename)); grid on;
text(0.05,0.9,sprintf('peak int= %f +- %f',pfit(3),sqrt(var(3))),'sc'); text(0.05,0.85,sprintf('center= %f +- %f',pfit(4),sqrt(var(4))),'sc');
text(0.05,0.8,sprintf('fwhm= %f +- %f',pfit(5),sqrt(var(5))),'sc'); text(0.05,0.75,sprintf('fgauss= %f +- %f',1-pfit(6),sqrt(var(6))),'sc');
text(0.05,0.7,sprintf('integ int= %f ',intint),'sc'); text(0.05,0.65,sprintf('min int(tol)= %f ',peak_tol),'sc');
subplot('position',[0.1 0.1 0.85 0.15]); plot(x,y-yfit'); ylabel('Residual'); xlabel('Radial position'); grid on; pause;
end; %if plot_fits==1;
fit(ifile,peak).int(azi) = pfit(3);fit(ifile,peak).cen(azi) = pfit(4); fit(ifile,peak).fwhm(azi) = pfit(5); fit(ifile,peak).shape(azi)=pfit(6);
fit(ifile,peak).dint(azi) = sqrt(var(3));fit(ifile,peak).dcen(azi) = sqrt(var(4));fit(ifile,peak).dfwhm(azi) = sqrt(var(5));fit(ifile,peak).dshape(azi)=sqrt(var(6));
else %peak intensity too weak to fit
fit(ifile,peak).int(azi) = below_tol_val; fit(ifile,peak).cen(azi) = below_tol_val; fit(ifile,peak).fwhm(azi) = below_tol_val; fit(ifile,peak).shape(azi)=below_tol_val;
fit(ifile,peak).dint(azi) = below_tol_val;fit(ifile,peak).dcen(azi) = below_tol_val;fit(ifile,peak).dfwhm(azi) = below_tol_val;fit(ifile,peak).dshape(azi)=below_tol_val;
end; %if statement for thresholding
end; %looping over azimuths
end; end; %looping over single peaks and checking if there are any singlepeaks
%%fit 2 PV simultaneously
opt = optimset('disp','off','large','on','jacobi','off','tolx',1e-2,'tolf',1e-2,'maxi',80);
slguess=0; fwhmguess1 = 3; fwhmguess2=3; %some guesses for peaks
fn = {'backg1' 'psv1' 'psv1'}; npar = [2 4 4];
low = [-inf 0 0 0 0 0 0 0 0 0]; up = [inf inf inf inf inf 1 inf inf inf 1];
for igroup=1:nDouble;
fprintf('\nFitting two peaks simultaneously for following azimuth bins:\n');
peaks=double_peaks{igroup}; dpeak1=peaks(1);dpeak2=peaks(2);
for azi = nAzi:-1:1
fprintf('%4i',azi*azi_bin);
x = cen0(dpeak2) - dx2: cen0(dpeak1) + dx2; Nx=length(x);
yraw = median(imm([azi*azi_bin:-1:azi*azi_bin-azi_bin+1], x-rho.min));
xback=x([1:2 Nx-1:Nx]); yback=yraw([1:2 Nx-1:Nx]); %x,y values to use for background determination
[pback_pars,s,mu]=polyfit(xback,yback,1); % 1st order polynomial (linear) background fit
pback=polyval(pback_pars,x,[],mu); %mu are centering and scaling parameters to help polyfit
y=yraw-pback+1000;
ymax= max(max(y)); xmax = x(round(y)==round(ymax)); xmax = xmax(1);
ymin = min(y); ymin = ymin(1);
ymax1 = y(x==cen0(dpeak1)); ymax2=y(x==cen0(dpeak2));
fom=ymax-1000;
if plot_fits;
figure(6);clf;
plot(x,yraw,'b+-',x,pback,'r-');
xlabel('radial position, pixels');ylabel('intensity');
title(sprintf('Peaks %i and %i and azimuth %i of %s; blue=data, red =background',dpeak1,dpeak2,azi,filename)); grid on;
if (fom>=peak_tol); text(0.05,0.95,sprintf('Will fit as FOM is %4.2f while tolerance is %4.2f',fom,peak_tol),'sc'); pause;
else text(0.05,0.95,sprintf('Will NOT fit as FOM is %4.2f while tolerance is %4.2f',fom,peak_tol),'sc'); pause; end;
end; %plotting raw data and background
if (fom>=peak_tol); % do the fit as intensity is above threshold
p0 = [slguess, ymin, (ymax1-ymin), cen0(dpeak1), fwhmguess1, .5,(ymax2-ymin), cen0(dpeak2), fwhmguess2, .5];
y0 = sumfun1(p0,[],npar,fn,x);
[pfit,rn,rs,ex,out,lam,jac] = lsqnonlin('sumfun1',p0,low,up,opt,[],npar,fn,x,y,ones(size(y)));
[dum, var] = confint(pfit,rs,jac); % variance of fitted params
var=full(var);
yfit = sumfun1(pfit,[],npar,fn,x);
yfit1=sumfun1(pfit(1:6),[],[2 4],{'backg1' 'psv1'},x); yfit2=sumfun1(pfit([1:2 7:10]),[],[2 4],{'backg1' 'psv1'},x);
if plot_fits
figure(8);clf; subplot('position',[0.1 0.3 0.85 0.6]);
plot(x,y,'+',x,y0,'-r',x,yfit1,'-b',x,yfit2,'--b',x,yfit,'-k'); ylabel('Intensity+100(offset)');
title(sprintf('Peaks %i and %i & azi %i of %s; symb=data, red=guess, blue/black=fit',dpeak1,dpeak2,azi,filename)); grid on;
text(0.05,0.9,sprintf('peak ints= %f & %f',pfit(3),pfit(7)),'sc'); text(0.05,0.85,sprintf('centers= %f & %f',pfit(4),pfit(8)),'sc');
text(0.05,0.8,sprintf('fwhms= %f & %f',pfit(5),pfit(9)),'sc'); text(0.05,0.75,sprintf('fgauss= %f &%f',1-pfit(6),1-pfit(10)),'sc');
subplot('position',[0.1 0.1 0.85 0.15]); plot(x,y-yfit'); ylabel('Residual'); xlabel('Radial position'); grid on; pause;
end; %if plot_fits==1;
fit(ifile,dpeak1).int(azi) = pfit(3);fit(ifile,dpeak1).cen(azi) = pfit(4); fit(ifile,dpeak1).fwhm(azi) = pfit(5);
fit(ifile,dpeak1).shape(azi) = pfit(6);fit(ifile,dpeak1).dint(azi) = sqrt(var(3));
fit(ifile,dpeak1).dcen(azi) = sqrt(var(4));fit(ifile,dpeak1).dfwhm(azi) = sqrt(var(5));
fit(ifile,dpeak2).int(azi) = pfit(7); fit(ifile,dpeak2).cen(azi) = pfit(8); fit(ifile,dpeak2).fwhm(azi) = pfit(9);
fit(ifile,dpeak2).shape(azi) = pfit(10);fit(ifile,dpeak2).dint(azi) = sqrt(var(7));
fit(ifile,dpeak2).dcen(azi) = sqrt(var(8));fit(ifile,dpeak2).dfwhm(azi) = sqrt(var(9));
else %peak too weak to fit
fit(ifile,dpeak1).int(azi) = below_tol_val; fit(ifile,dpeak1).cen(azi) = below_tol_val; fit(ifile,dpeak1).fwhm(azi) = below_tol_val; fit(ifile,dpeak1).shape(azi)=below_tol_val;
fit(ifile,dpeak1).dint(azi) = below_tol_val;fit(ifile,dpeak1).dcen(azi) = below_tol_val;fit(ifile,dpeak1).dfwhm(azi) = below_tol_val;fit(ifile,dpeak1).dshape(azi)=below_tol_val;
fit(ifile,dpeak2).int(azi) = below_tol_val; fit(ifile,dpeak2).cen(azi) = below_tol_val; fit(ifile,dpeak2).fwhm(azi) = below_tol_val; fit(ifile,dpeak2).shape(azi)=below_tol_val;
fit(ifile,dpeak2).dint(azi) = below_tol_val;fit(ifile,dpeak2).dcen(azi) = below_tol_val;fit(ifile,dpeak2).dfwhm(azi) = below_tol_val;fit(ifile,dpeak2).dshape(azi)=below_tol_val;
end;
end;
end;
fprintf('\n');
  1 Comment
per isakson
per isakson on 9 Feb 2015
  • I don't think unique is called by your script.
  • The script cannot be run because of missing variables, e.g. center
  • Descibe your problem in more detail and upload the file with the "paper-clip-button"

Sign in to comment.

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!