diff --git a/@MyFit/MyFit.m b/@MyFit/MyFit.m index 961f6a2..71f9be1 100644 --- a/@MyFit/MyFit.m +++ b/@MyFit/MyFit.m @@ -1,832 +1,851 @@ classdef MyFit < dynamicprops %Note that dynamicprops classes are handle classes. properties (Access=public) Data; init_params=[]; scale_init=[]; lim_lower; lim_upper; enable_plot; plot_handle; %Calibration values supplied externally CalVals=struct(); init_color='c'; end properties (GetAccess=public, SetAccess=private) Fit; Gui; Fitdata; Gof; FitInfo; FitStruct; coeffs; fit_name='Linear' end properties (Access=private) %Structure used for initializing GUI of userpanel UserGui; Parser; enable_gui=1; hline_init; end properties (Dependent=true) fit_function; anon_fit_fun; fit_tex; fit_params; fit_param_names; valid_fit_names; n_params; scaled_params; init_param_fun; x_vec; n_user_fields; user_field_tags; user_field_names; user_field_vals; fullpath; filename; base_dir; save_path; session_name; end events NewFit; NewInitVal; end %Parser function methods (Access=private) %Creates parser for constructor function createParser(this) p=inputParser; addParameter(p,'fit_name','Linear',@ischar) addParameter(p,'Data',MyTrace()); addParameter(p,'Fit',MyTrace()); addParameter(p,'x',[]); addParameter(p,'y',[]); addParameter(p,'enable_gui',1); addParameter(p,'enable_plot',0); addParameter(p,'plot_handle',[]); addParameter(p,'base_dir',pwd); addParameter(p,'session_name','placeholder'); addParameter(p,'filename','placeholder'); this.Parser=p; end end %Public methods methods (Access=public) %Constructor function function this=MyFit(varargin) createFitStruct(this); createParser(this); parse(this.Parser,varargin{:}); parseInputs(this); initCalVals(this); if ismember('Data',this.Parser.UsingDefaults) &&... ~ismember('x',this.Parser.UsingDefaults) &&... ~ismember('y',this.Parser.UsingDefaults) this.Data.x=this.Parser.Results.x; this.Data.y=this.Parser.Results.y; end %Sets the scale_init to 1, this is used for the GUI. this.scale_init=ones(1,this.n_params); this.init_params=ones(1,this.n_params); %Creates the structure that contains variables for calibration %of fit results createUserGuiStruct(this); if this.enable_gui; createGui(this); end %If the data is appropriate, generates initial %parameters if validateData(this); genInitParams(this); end end %Deletion function of object function delete(this) if this.enable_gui %Avoids loops set(this.Gui.Window,'CloseRequestFcn',''); %Deletes the figure delete(this.Gui.Window); %Removes the figure handle to prevent memory leaks this.Gui=[]; end if ~isempty(this.hline_init); delete(this.hline_init); end if ~isempty(this.Fit.hlines); delete(this.Fit.hlines{:}); end end %Close figure callback simply calls delete function for class function closeFigure(this,~,~) delete(this); end %Saves the metadata function saveParams(this,varargin) p=inputParser; addParameter(p,'save_user_params',true); addParameter(p,'save_gof',true); parse(p,varargin{:}); save_user_params=p.Results.save_user_params; save_gof=p.Results.save_gof; assert(~isempty(this.coeffs) && ... length(this.coeffs)==this.n_params,... ['The number of calculated coefficients (%i) is not',... ' equal to the number of parameters (%i).', ... ' Perform a fit before trying to save parameters.'],... length(this.coeffs),this.n_params); %Creates combined strings of form: Linewidth (b), where %Linewidth is the parameter name and b is the parameter tag headers=cellfun(@(x,y) sprintf('%s (%s)',x,y),... this.fit_param_names, this.fit_params,'UniformOutput',0); save_data=this.coeffs; if save_user_params user_field_headers=cellfun(@(x,y) ... sprintf('%s. %s',this.UserGui.Fields.(x).parent,y),... this.user_field_tags,this.user_field_names,... 'UniformOutput',0)'; headers=[headers, user_field_headers]; save_data=[save_data,this.user_field_vals']; end if save_gof headers=[headers,fieldnames(this.Gof)']; save_data=[save_data,struct2array(this.Gof)]; end n_columns=length(headers); %Sets the column width. Pads 2 for legibility. col_width=cellfun(@(x) length(x), headers)+2; %Min column width of 24 col_width(col_width<24)=24; if ~exist(this.base_dir,'dir') mkdir(this.base_dir) end if ~exist(this.save_path,'dir') mkdir(this.save_path) end if exist(this.fullpath,'file') fileID=fopen(this.fullpath,'a'); else fileID=fopen(this.fullpath,'w'); pre_fmt_str=repmat('%%%is\\t',1,n_columns); fmt_str=sprintf([pre_fmt_str,'\r\n'],col_width); fprintf(fileID,fmt_str,headers{:}); end pre_fmt_str_nmb=repmat('%%%i.15e\\t',1,n_columns); nmb_fmt_str=sprintf([pre_fmt_str_nmb,'\r\n'],col_width); fprintf(fileID,nmb_fmt_str,save_data); fclose(fileID); end function loadFit(this,fullfilename,varargin) p=inputParser; addParameter(p,'line_no',1); parse(p,varargin{:}) n=p.Results.line_no; load_table=readtable(fullfilename); load_names=fieldnames(load_table); for i=1:this.n_params this.coeffs(i)=load_table.(load_names{i})(n); end end function setFitParams(this,coeffs) assert(length(coeffs)==this.n_params,... ['The length of the coefficient vector (currently %i) ',... 'must be equal to the number of parameters (%i)'],... length(this.coeffs),this.n_params) this.coeffs=coeffs; end %Initializes the CalVals structure. function initCalVals(this) switch this.fit_name case 'Lorentzian' %Line spacing is the spacing between all the lines, %i.e. number of lines times the spacing between each %one this.CalVals.line_spacing=1; case 'DoubleLorentzian' this.CalVals.line_spacing=1; end end %Fits the trace using currently set parameters, depending on the %model. function fitTrace(this) switch this.fit_name case 'Linear' %Fits polynomial of order 1 this.coeffs=polyfit(this.Data.x,this.Data.y,1); case 'Quadratic' %Fits polynomial of order 2 this.coeffs=polyfit(this.Data.x,this.Data.y,2); -<<<<<<< HEAD this.Fit.y=polyval(this.coeffs,this.Fit.x); - case {'Lorentzian','LorentzianGrad',... - 'DoubleLorentzian','DoubleLorentzianGrad',...... - 'Exponential','Gaussian'} + case {'Lorentzian','LorentzianGrad','Gaussian'} + offset=median(this.Data.x); + this.Data.x=this.Data.x-offset; + this.init_params(3)=this.init_params(3)-offset; + doFit(this); + this.coeffs(3)=this.coeffs(3)+offset; + case {'DoubleLorentzian','DoubleLorentzianGrad'} +% offset=median(this.Data.x); +% this.Data.x=this.Data.x-offset; +% this.init_params(3)=this.init_params(3)-offset; +% this.init_params(6)=this.init_params(6)-offset; + doFit(this); +% this.coeffs(3)=this.coeffs(3)+offset; +% this.coeffs(6)=this.coeffs(6)+offset; + case {'Exponential','Gorodetsky2000'} doFit(this); -======= - case {'Exponential','Gaussian','Lorentzian','DoubleLorentzian'} - doFit(this) ->>>>>>> 16b427aaa83c68ba5bb0ef84e3578761bbb49bc7 otherwise error('Selected fit is invalid'); end calcFit(this); calcUserParams(this); %Sets the new initial parameters to be the fitted parameters this.init_params=this.coeffs; %Resets the scale variables for the GUI this.scale_init=ones(1,this.n_params); %Updates the gui if it is enabled if this.enable_gui; updateGui(this); end %Plots the fit if the flag is on if this.enable_plot; plotFit(this); end %Triggers new fit event triggerNewFit(this); end function calcUserParams(this) switch this.fit_name case 'Lorentzian' this.mech_lw=this.coeffs(2); %#ok this.mech_freq=this.coeffs(3); %#ok this.Q=this.mech_freq/this.mech_lw; %#ok this.opt_lw=convOptFreq(this,this.coeffs(2)); %#ok this.Qf=this.mech_freq*this.Q; %#ok case 'DoubleLorentzian' this.opt_lw1=convOptFreq(this,this.coeffs(2)); %#ok this.opt_lw2=convOptFreq(this,this.coeffs(5)); %#ok splitting=abs(this.coeffs(6)-this.coeffs(3)); this.mode_split=convOptFreq(this,splitting); %#ok case 'Exponential' this.tau=1/this.coeffs(2); %#ok this.lw=this.coeffs(2)/pi; %#ok this.Q=pi*this.freq*this.tau; %#ok this.Qf=this.Q*this.freq; %#ok otherwise end end function real_freq=convOptFreq(this,freq) real_freq=freq*this.spacing*this.line_no/this.CalVals.line_spacing; end function createUserGuiStruct(this) this.UserGui=struct('Fields',struct(),'Tabs',struct()); switch this.fit_name case 'Lorentzian' %Parameters for the tab relating to mechanics this.UserGui.Tabs.Mech.tab_title='Mech.'; this.UserGui.Tabs.Mech.Children={}; addUserField(this,'Mech','mech_lw','Linewidth (Hz)',1,... 'enable_flag','off') addUserField(this,'Mech','Q',... 'Qualify Factor (x10^6)',1e6,... 'enable_flag','off','conv_factor',1e6) addUserField(this,'Mech','mech_freq','Frequency (MHz)',1e6,... 'conv_factor',1e6, 'enable_flag','off') addUserField(this,'Mech','Qf','Q\times f (10^{14} Hz)',1e14,... 'conv_factor',1e14,'enable_flag','off'); %Parameters for the tab relating to optics this.UserGui.Tabs.Opt.tab_title='Optical'; this.UserGui.Tabs.Opt.Children={}; addUserField(this,'Opt','spacing',... 'Line Spacing (MHz)',1e6,'conv_factor',1e6,... 'Callback', @(~,~) calcUserParams(this)); addUserField(this,'Opt','line_no','Number of lines',10,... 'Callback', @(~,~) calcUserParams(this)); addUserField(this,'Opt','opt_lw','Linewidth (MHz)',1e6,... 'enable_flag','off','conv_factor',1e6); case 'DoubleLorentzian' this.UserGui.Tabs.Opt.tab_title='Optical'; this.UserGui.Tabs.Opt.Children={}; addUserField(this,'Opt','spacing',... 'Line Spacing (MHz)',1e6,'conv_factor',1e6,... 'Callback', @(~,~) calcUserParams(this)); addUserField(this,'Opt','line_no','Number of lines',10,... 'Callback', @(~,~) calcUserParams(this)); addUserField(this,'Opt','opt_lw1','Linewidth 1 (MHz)',1e6,... 'enable_flag','off','conv_factor',1e6); addUserField(this,'Opt','opt_lw2','Linewidth 2 (MHz)',1e6,... 'enable_flag','off','conv_factor',1e6); addUserField(this,'Opt','mode_split',... 'Modal splitting (MHz)',1e6,... 'enable_flag','off','conv_factor',1e6); case 'Exponential' this.UserGui.Tabs.Q.tab_title='Q'; this.UserGui.Tabs.Q.Children={}; addUserField(this,'Q','tau','\tau (s)',1,... 'enable_flag','off') addUserField(this,'Q','lw','Linewidth (Hz)',1,... 'enable_flag','off') addUserField(this,'Q','Q',... 'Qualify Factor (x10^6)',1e6,... 'enable_flag','off','conv_factor',1e6) addUserField(this,'Q','freq','Frequency (MHz)',1e6,... 'conv_factor',1e6, 'enable_flag','on') addUserField(this,'Q','Qf','Q\times f (10^{14} Hz)',1e14,... 'conv_factor',1e14,'enable_flag','off'); otherwise %Do nothing if there is no defined user parameters end end %Parent is the parent tab for the userfield, tag is the tag given %to the GUI element, title is the text written next to the field, %initial value is the initial value of the property and change_flag %determines whether the gui element is enabled for writing or not. function addUserField(this, parent, tag, title, ... init_val,varargin) %Parsing inputs p=inputParser(); addRequired(p,'Parent'); addRequired(p,'Tag'); addRequired(p,'Title'); addRequired(p,'init_val'); addParameter(p,'enable_flag','on'); addParameter(p,'Callback',''); addParameter(p,'conv_factor',1); parse(p,parent,tag,title,init_val,varargin{:}); tag=p.Results.Tag; %Populates the UserGui struct this.UserGui.Fields.(tag).parent=p.Results.Parent; this.UserGui.Fields.(tag).title=p.Results.Title; this.UserGui.Fields.(tag).init_val=p.Results.init_val; this.UserGui.Fields.(tag).enable_flag=... p.Results.enable_flag; this.UserGui.Fields.(tag).conv_factor=p.Results.conv_factor; this.UserGui.Fields.(tag).Callback=... p.Results.Callback; this.UserGui.Tabs.(p.Results.Parent).Children{end+1}=tag; %Adds the new property to the class addUserProp(this, tag); end function addUserProp(this,tag) prop=addprop(this,tag); if this.enable_gui prop.GetMethod=@(this) getUserVal(this,tag); prop.SetMethod=@(this, val) setUserVal(this, val, tag); prop.Dependent=true; end end function val=getUserVal(this, tag) conv_factor=this.UserGui.Fields.(tag).conv_factor; val=str2double(this.Gui.([tag,'Edit']).String)*conv_factor; end function setUserVal(this, val, tag) conv_factor=this.UserGui.Fields.(tag).conv_factor; this.Gui.([tag,'Edit']).String=num2str(val/conv_factor); end % Callbacks %Save fit function callback function saveFitCallback(this,~,~) assert(~isempty(this.base_dir),'Save directory is not specified'); assert(ischar(this.base_dir),... ['Save directory is not specified.',... ' Should be of type char but is %s.'], ... class(this.base_dir)) try this.Fit.save('name',this.filename,... 'base_dir',this.base_dir) catch error(['Attempted to save to directory %s',... ' with file name %s, but failed'],this.base_dir,... this.filename); end end function saveParamCallback(this,~,~) saveParams(this); end %Callback functions for sliders in GUI. Uses param_ind to find out %which slider the call is coming from, this was implemented to %speed up the callback. function sliderCallback(this, param_ind, hObject, ~) %Gets the value from the slider scale=get(hObject,'Value'); %Updates the scale with a new value this.scale_init(param_ind)=10^((scale-50)/50); %Updates the edit box with the new value from the slider set(this.Gui.(sprintf('Edit_%s',this.fit_params{param_ind})),... 'String',sprintf('%3.3e',this.scaled_params(param_ind))); if this.enable_plot; plotInitFun(this); end end %Callback function for edit boxes in GUI function editCallback(this, hObject, ~) init_param=str2double(get(hObject,'String')); tag=get(hObject,'Tag'); %Finds the index where the fit_param name begins (convention is %after the underscore) fit_param=tag((strfind(tag,'_')+1):end); param_ind=strcmp(fit_param,this.fit_params); %Updates the slider to be such that the scaling is 1 set(this.Gui.(sprintf('Slider_%s',fit_param)),... 'Value',50); %Updates the correct initial parameter this.init_params(param_ind)=init_param; if this.enable_plot; plotInitFun(this); end %Triggers event for new init values triggerNewInitVal(this); end %Callback function for analyze button in GUI. Checks if the data is %ready for fitting. function analyzeCallback(this, ~, ~) assert(validateData(this),... ['The length of x is %d and the length of y is',... ' %d. The lengths must be equal and greater than ',... 'the number of fit parameters to perform a fit'],... length(this.Data.x),length(this.Data.y)) fitTrace(this); end %Callback for clearing the fits on the axis. function clearFitCallback(this,~,~) clearFit(this); end %Callback function for generate init parameters button. Updates GUI %afterwards function initParamCallback(this,~,~) genInitParams(this); updateGui(this); end %Generates model-dependent initial parameters, lower and upper %boundaries. function genInitParams(this) assert(validateData(this), ['The data must be vectors of',... ' equal length greater than the number of fit parameters.',... ' Currently the number of fit parameters is %d, the',... ' length of x is %d and the length of y is %d'],... this.n_params,length(this.Data.x),length(this.Data.y)); %Cell for putting parameters in to be interpreted in the %parser. Element 1 contains the init params, Element 2 contains %the lower limits and Element 3 contains the upper limits. params={}; switch this.fit_name case 'Exponential' [params{1},params{2},params{3}]=... initParamExponential(this.Data.x,this.Data.y); case 'Gaussian' [params{1},params{2},params{3}]=... initParamGaussian(this.Data.x,this.Data.y); case 'Lorentzian' [params{1},params{2},params{3}]=... initParamLorentzian(this.Data.x,this.Data.y); case 'LorentzianGrad' [params{1},params{2},params{3}]=... initParamLorentzianGrad(this.Data.x,this.Data.y); case 'DoubleLorentzian' [params{1},params{2},params{3}]=... initParamDblLorentzian(this.Data.x,this.Data.y); case 'DoubleLorentzianGrad' [params{1},params{2},params{3}]=... initParamDblLorentzianGrad(this.Data.x,this.Data.y); + case 'Gorodetsky2000' + [params{1},params{2},params{3}]=... + initParamGorodetsky2000(this.Data.x,this.Data.y); end %Validates the initial parameters p=createFitParser(this.n_params); parse(p,params{:}); %Loads the parsed results into the class variables this.init_params=p.Results.init_params; this.lim_lower=p.Results.lower; this.lim_upper=p.Results.upper; + %Resets scale init variables + this.scale_init=ones(1,this.n_params); %Plots the fit function with the new initial parameters if this.enable_plot; plotInitFun(this); end + updateGui(this); end function calcFit(this) this.Fit.x=this.x_vec; input_coeffs=num2cell(this.coeffs); this.Fit.y=this.anon_fit_fun(this.Fit.x,input_coeffs{:}); end %Plots the trace contained in the Fit MyTrace object after %calculating the new values function plotFit(this,varargin) calcFit(this); assert((isa(this.plot_handle,'matlab.graphics.axis.Axes')||... isa(this.plot_handle,'matlab.ui.control.UIAxes')),... 'plot_handle property must be defined to valid axis in order to plot') this.Fit.plotTrace(this.plot_handle,varargin{:}); end %Clears the plots function clearFit(this) cellfun(@(x) delete(x), this.Fit.hlines); delete(this.hline_init); this.hline_init=[]; this.Fit.hlines={}; end %Function for plotting fit model with current initial parameters. function plotInitFun(this) %Substantially faster than any alternative - generating %anonymous functions is very cpu intensive. input_cell=num2cell(this.scaled_params); y_vec=feval(this.FitStruct.(this.fit_name).anon_fit_fun,... this.x_vec,input_cell{:}); if isempty(this.hline_init) this.hline_init=plot(this.plot_handle,this.x_vec,y_vec,... 'Color',this.init_color); else set(this.hline_init,'XData',this.x_vec,'YData',y_vec); end end end %Private methods methods(Access=private) %Creates the GUI of MyFit, in separate file. createGui(this); %Creates a panel for the GUI, in separate file createTab(this, tab_tag, bg_color, button_h); %Creats two vboxes (from GUI layouts) to display values of %quantities createUnitBox(this, bg_color, h_parent, name); %Creates edit box inside a UnitDisp for showing label and value of %a quantity. Used in conjunction with createUnitBox createUnitDisp(this,varargin); %Sets the class variables to the inputs from the inputParser. function parseInputs(this) for i=1:length(this.Parser.Parameters) %Takes the value from the inputParser to the appropriate %property. if isprop(this,this.Parser.Parameters{i}) this.(this.Parser.Parameters{i})=... this.Parser.Results.(this.Parser.Parameters{i}); end end end %Does the fit with the currently set parameters function doFit(this) + ft=fittype(this.fit_function,'coefficients',this.fit_params); + opts=fitoptions('Method','NonLinearLeastSquares',... + 'Lower',this.lim_lower,... + 'Upper',this.lim_upper,... + 'StartPoint',this.init_params,... + 'MaxFunEvals',2000,... + 'MaxIter',2000,... + 'TolFun',1e-9,... + 'TolX',1e-9); %Fits with the below properties. Chosen for maximum accuracy. [this.Fitdata,this.Gof,this.FitInfo]=... - fit(this.Data.x,this.Data.y,this.fit_function,... - 'Lower',this.lim_lower,'Upper',this.lim_upper,... - 'StartPoint',this.init_params, .... - 'MaxFunEvals',2000,'MaxIter',2000,'TolFun',1e-9); + fit(this.Data.x,this.Data.y,ft,opts); %Puts the coeffs into the class variable. this.coeffs=coeffvalues(this.Fitdata); end %Triggers the NewFit event such that other objects can use this to %e.g. plot new fits function triggerNewFit(this) notify(this,'NewFit'); end function triggerNewInitVal(this) notify(this,'NewInitVal'); end %Creates the struct used to get all things relevant to the fit %model function createFitStruct(this) %Adds fits addFit(this,'Linear','a*x+b','$$ax+b$$',{'a','b'},... {'Gradient','Offset'}) addFit(this,'Quadratic','a*x^2+b*x+c','$$ax^2+bx+c$$',... {'a','b','c'},{'Quadratic coeff.','Linear coeff.','Offset'}); addFit(this,'Gaussian','a*exp(-((x-c)/b)^2/2)+d',... '$$ae^{-\frac{(x-c)^2}{2b^2}}+d$$',{'a','b','c','d'},... {'Amplitude','Width','Center','Offset'}); addFit(this,'Lorentzian','1/pi*a*b/2/((x-c)^2+(b/2)^2)+d',... '$$\frac{a}{\pi}\frac{b/2}{(x-c)^2+(b/2)^2}+d$$',{'a','b','c','d'},... {'Amplitude','Width','Center','Offset'}); - addFit(this,'LorentzianGrad','1/pi*a*b/2/((x-c)^2+(b/2)^2)+d*x+e',... + addFit(this,'LorentzianGrad','1/pi*a*b/2/((x-c)^2+(b/2)^2)+d*(x-c)+e',... '$$\frac{a}{\pi}\frac{b/2}{(x-c)^2+(b/2)^2}+d*x+e$$',... {'a','b','c','d','e'},... {'Amplitude','Width','Center','Gradient','Offset'}); addFit(this,'Exponential','a*exp(b*x)+c',... '$$ae^{bx}+c$$',{'a','b','c'},... {'Amplitude','Rate','Offset'}); addFit(this,'DoubleLorentzian',... '1/pi*b/2*a/((x-c)^2+(b/2)^2)+1/pi*e/2*d/((x-f)^2+(e/2)^2)+g',... '$$\frac{a}{\pi}\frac{b/2}{(x-c)^2+(b/2)^2}+\frac{d}{\pi}\frac{e/2}{(x-f)^2+(e/2)^2}+g$$',... {'a','b','c','d','e','f','g'},... {'Amplitude 1','Width 1','Center 1','Amplitude 2',... 'Width 2','Center 2','Offset'}); addFit(this,'DoubleLorentzianGrad',... - '1/pi*b/2*a/((x-c)^2+(b/2)^2)+1/pi*e/2*d/((x-f)^2+(e/2)^2)+g',... + '1/pi*b/2*a/((x-c)^2+(b/2)^2)+1/pi*e/2*d/((x-f)^2+(e/2)^2)+g*(x-c)+h',... '$$\frac{a}{\pi}\frac{b/2}{(x-c)^2+(b/2)^2}+\frac{d}{\pi}\frac{e/2}{(x-f)^2+(e/2)^2}+g*x+h$$',... {'a','b','c','d','e','f','g','h'},... {'Amplitude 1','Width 1','Center 1','Amplitude 2',... 'Width 2','Center 2','Gradient','Offset'}); addFit(this,'Gorodetsky2000',... - 'a*abs( (k0^2 - kex^2 + gamma^2 - (x-d).^2 + 2i*k0.*(x-d))./( (k0 + kex)^2 + gamma^2 - (x-d).^2 + 2i.*(x-d)*(k0 + kex) )).^2+e*x+f',... - '$$a\left|\frac{\kappa_0^2-\kappa_{ex}^2+\gamma^2-(x-d)^2+2i\kappa_0(x-d)}{(\kappa_0+\kappa_{ex})^2+\gamma^2-(x-d)^2+2i(x-d)(\kappa_0+\kappa_{ex})}\right|^2+ex+f$$',... - {'k_0', 'k_ex', 'gamma', 'a','d','e','f'},... - {'Intrinsic Lw','Extrinsic Lw','Loss','Amplitude','Center','Gradient','Offset'}); + 'a*abs( (k0^2 - kex^2 + gamma^2 - (x-b).^2 + 2i*k0.*(x-b))./( (k0 + kex)^2 + gamma^2 - (x-b).^2 + 2i.*(x-b)*(k0 + kex) )).^2',... + '$$a\left|\frac{\kappa_0^2-\kappa_{ex}^2+\gamma^2-(x-b)^2+2i\kappa_0(x-b)}{(\kappa_0+\kappa_{ex})^2+\gamma^2-(x-b)^2+2i(x-b)(\kappa_0+\kappa_{ex})}\right|^2$$',... + { 'a','b', 'gamma','k0', 'kex'},... + {'Background','Center','Mode splitting','Intrinsic','Extrinsic'}); end %Updates the GUI if the edit or slider boxes are changed from %elsewhere. function updateGui(this) %Converts the scale variable to the value between 0 and 100 %necessary for the slider slider_vals=50*log10(this.scale_init)+50; for i=1:this.n_params set(this.Gui.(sprintf('Edit_%s',this.fit_params{i})),... 'String',sprintf('%3.3e',this.scaled_params(i))); set(this.Gui.(sprintf('Slider_%s',this.fit_params{i})),... 'Value',slider_vals(i)); end end %Adds a fit to the list of fits function addFit(this,fit_name,fit_function,fit_tex,fit_params,... fit_param_names) this.FitStruct.(fit_name).fit_function=fit_function; this.FitStruct.(fit_name).fit_tex=fit_tex; this.FitStruct.(fit_name).fit_params=fit_params; this.FitStruct.(fit_name).fit_param_names=fit_param_names; %Generates the anonymous fit function from the above args=['@(x,', strjoin(fit_params,','),')']; - anon_fit_fun=str2func(vectorize([args,fit_function])); - this.FitStruct.(fit_name).anon_fit_fun=anon_fit_fun; + this.FitStruct.(fit_name).anon_fit_fun=... + str2func(vectorize([args,fit_function])); end %Checks if the class is ready to perform a fit function bool=validateData(this) bool=~isempty(this.Data.x) && ~isempty(this.Data.y) && ... length(this.Data.x)==length(this.Data.y) && ... length(this.Data.x)>=this.n_params; end end % Get and set functions methods % Set functions %Set function for fit_name. function set.fit_name(this,fit_name) assert(ischar(fit_name),'The fit name must be a string'); %Capitalizes the first letter fit_name=[upper(fit_name(1)),lower(fit_name(2:end))]; %Checks it is a valid fit name ind=strcmpi(fit_name,this.valid_fit_names);%#ok assert(any(ind),'%s is not a supported fit name',fit_name); this.fit_name=this.valid_fit_names{ind}; %#ok end % Get functions for dependent variables %Generates the valid fit names function valid_fit_names=get.valid_fit_names(this) valid_fit_names=fieldnames(this.FitStruct); end %Grabs the correct fit function from FitStruct function fit_function=get.fit_function(this) fit_function=this.FitStruct.(this.fit_name).fit_function; end %Grabs the correct tex string from FitStruct function fit_tex=get.fit_tex(this) fit_tex=this.FitStruct.(this.fit_name).fit_tex; end %Grabs the correct anon fit function from FitStruct function anon_fit_fun=get.anon_fit_fun(this) anon_fit_fun=this.FitStruct.(this.fit_name).anon_fit_fun; end %Grabs the correct fit parameters from FitStruct function fit_params=get.fit_params(this) fit_params=this.FitStruct.(this.fit_name).fit_params; end %Grabs the correct fit parameter names from FitStruct function fit_param_names=get.fit_param_names(this) fit_param_names=this.FitStruct.(this.fit_name).fit_param_names; end %Calculates the scaled initial parameters function scaled_params=get.scaled_params(this) scaled_params=this.scale_init.*this.init_params; end %Calculates the number of parameters in the fit function function n_params=get.n_params(this) n_params=length(this.fit_params); end %Generates a vector of x values for plotting function x_vec=get.x_vec(this) x_vec=linspace(min(this.Data.x),max(this.Data.x),1000); end function n_user_fields=get.n_user_fields(this) n_user_fields=length(this.user_field_tags); end function user_field_tags=get.user_field_tags(this) user_field_tags=fieldnames(this.UserGui.Fields); end function user_field_names=get.user_field_names(this) user_field_names=cellfun(@(x) this.UserGui.Fields.(x).title,... this.user_field_tags,'UniformOutput',0); end function user_field_vals=get.user_field_vals(this) user_field_vals=cellfun(@(x) this.(x), this.user_field_tags); end function fullpath=get.fullpath(this) fullpath=[this.save_path,this.filename,'.txt']; end function save_path=get.save_path(this) save_path=createSessionPath(this.base_dir,this.session_name); end function base_dir=get.base_dir(this) try base_dir=this.Gui.BaseDir.String; catch base_dir=pwd; end end function set.base_dir(this,base_dir) this.Gui.BaseDir.String=base_dir; end function session_name=get.session_name(this) try session_name=this.Gui.SessionName.String; catch session_name=''; end end function set.session_name(this,session_name) this.Gui.SessionName.String=session_name; end function filename=get.filename(this) try filename=this.Gui.FileName.String; catch filename='placeholder'; end end function set.filename(this,filename) this.Gui.FileName.String=filename; end end end \ No newline at end of file diff --git a/@MyPeakFinder/MyPeakFinder.m b/@MyPeakFinder/MyPeakFinder.m index 8d6070c..f9be676 100644 --- a/@MyPeakFinder/MyPeakFinder.m +++ b/@MyPeakFinder/MyPeakFinder.m @@ -1,271 +1,284 @@ classdef MyPeakFinder < handle properties Trace; Peaks; end methods function this=MyPeakFinder(varargin) p=inputParser; addParameter(p,'Trace',MyTrace()); parse(p,varargin{:}) this.Trace=p.Results.Trace; + this.Peaks=struct('Location',[],'Width',[],'Prominence',[],... + 'Value',[]); + end + + %Checks if a peak exists within the given limits + function bool=checkPeakExist(this,x_min,x_max) + assert(isnumeric(x_min) & isnumeric(x_max),... + 'x_min and x_max must be numbers'); + assert(x_max>x_min,['x_max must be greater than x_min,',... + ' currently x_min is %e while x_max is %e'],x_min,x_max); + bool=any([this.Peaks.Location]>x_min & [this.Peaks.Location] 0); isnumber=@(x) isnumeric(x) && isscalar(x); valid_sorts={'none','ascend','descend'}; sortstrcheck=@(x) any(validatestring(x,valid_sorts)); addParameter(p,'FindMinima',false); addParameter(p,'MinPeakDistance',0,ispositive); addParameter(p,'MinPeakHeight',-Inf,isnumber); addParameter(p,'MinPeakWidth',0,ispositive); addParameter(p,'MaxPeakWidth',Inf,ispositive); addParameter(p,'MinPeakProminence',0.1,isnumber); addParameter(p,'SortStr','none',sortstrcheck); addParameter(p,'Threshold',0); addParameter(p,'ClearPeaks',true); addParameter(p,'Limits',[min(this.Trace.x),max(this.Trace.x)]) addParameter(p,'NPeaks',0); parse(p,varargin{:}); %Sets the indices to be searched x_lim=p.Results.Limits; ind=(this.Trace.xx_lim(1)); %Sets the minimum peak prominence, which is always normalized min_peak_prominence=p.Results.MinPeakProminence*... peak2peak(this.Trace.y(ind)); %We must condition the Trace such that x is increasing and has %no NaNs nan_ind=isnan(this.Trace.x) | isnan(this.Trace.y); this.Trace.x(nan_ind)=[]; this.Trace.y(nan_ind)=[]; ind(nan_ind)=[]; if ~issorted(this.Trace.x) this.Trace.x=flipud(this.Trace.x); this.Trace.y=flipud(this.Trace.y); %If it is still not sorted, we sort it if ~issorted(this.Trace.x) 'Sorting'; [this.Trace.x,sort_ind]=sort(this.Trace.x); this.Trace.y=this.Trace.y(sort_ind); end end %If we are looking for minima, we invert the trace if p.Results.FindMinima y=-this.Trace.y; else y=this.Trace.y; end %As there is no way to tell it to look for infinite numbers of %peaks when you specify NPeaks, we only specify this parameter %if we need to if p.Results.NPeaks extra_args={'NPeaks',p.Results.NPeaks}; else extra_args={}; end %We now search for the peaks using the specified parameters [pks,locs,wdth,prom]=findpeaks(y(ind),... this.Trace.x(ind),... 'MinPeakDistance',p.Results.MinPeakDistance,... 'MinPeakheight',p.Results.MinPeakHeight,... 'MinPeakWidth',p.Results.MinPeakWidth,... 'MaxPeakWidth',p.Results.MaxPeakWidth,... 'SortStr',p.Results.SortStr,... 'MinPeakProminence',min_peak_prominence,... 'Threshold',p.Results.Threshold,extra_args{:}); %We invert the results back if we are looking for minima. if p.Results.FindMinima pks=-pks; prom=-prom; end PeakStruct=struct('Value',num2cell(pks),... 'Location',num2cell(locs),... 'Width',num2cell(wdth),... 'Prominence',num2cell(prom)); %If the clearpeaks flag is set, we delete the previous peaks. if p.Results.ClearPeaks clearPeaks(this); this.Peaks=PeakStruct; else this.Peaks=[this.Peaks;PeakStruct]; end end function loadTrace(this,fullfilename) %Finds type of file [~,~,ext]=fileparts(fullfilename); switch ext case '.txt' loadTrace(this.Trace,fullfilename); case '.mat' DataStruct=load(fullfilename); fields=fieldnames(DataStruct); %We try to find which vectors are x and y for the data, %first we find the two longest vectors in the .mat file vec_length=structfun(@(x) length(x), DataStruct); max_val=max(vec_length); ind=(max_val==vec_length); assert(sum(ind)==2,['There are %i vectors with length %i, ',... 'cannot determine which are the data vectors'],... sum(ind),max_val); vec_names=fields(ind); %Now we do some basic conditioning of these vectors: %Make column vectors and remove NaNs. vec{1}=DataStruct.(vec_names{1})(:); vec{2}=DataStruct.(vec_names{2})(:); nan_ind=isnan(vec{1}) | isnan(vec{2}); vec{1}(nan_ind)=[]; vec{2}(nan_ind)=[]; %We find what x is by looking for a sorted vector ind_x=cellfun(@(x) issorted(x,'monotonic'),vec); this.Trace.x=vec{ind_x}; this.Trace.y=vec{~ind_x}; otherwise error('File type %s is not supported',ext) end end function fitAllPeaks(this,varargin) valid_fits={'Lorentzian','DoubleLorentzian','Gaussian'}; validate_fit=@(x) all(ismember(x,valid_fits)); p=inputParser; - addParameter(p,'FitNames',{'LorentzianGrad'},validate_fit); + addParameter(p,'FitNames',{'DoubleLorentzianGrad'},validate_fit); addParameter(p,'base_dir',pwd); addParameter(p,'session_name','placeholder'); addParameter(p,'filename','placeholder'); parse(p,varargin{:}); fit_names=p.Results.FitNames; %We instantiate the MyFit objects used for the fitting Fits=struct(); for i=1:length(fit_names) Fits.(fit_names{i})=MyFit('fit_name',fit_names{i},... 'enable_gui',0); Fits.(fit_names{i}).base_dir=p.Results.base_dir; Fits.(fit_names{i}).session_name=p.Results.session_name; Fits.(fit_names{i}).filename=p.Results.filename; end %We fit the peaks for i=1:length(this.Peaks) %First extract the data around the peak [x_fit,y_fit]=extractPeak(this,i); for j=1:length(fit_names) Fits.(fit_names{j}).Data.x=x_fit; Fits.(fit_names{j}).Data.y=y_fit; genInitParams(Fits.(fit_names{j})); fitTrace(Fits.(fit_names{j})); saveParams(Fits.(fit_names{j}),... 'save_user_params',false,... 'save_gof',true) + Fits.(fit_names{j}).FitInfo end end end function [x_peak,y_peak]=extractPeak(this,peak_no) loc=this.Peaks(peak_no).Location; w=this.Peaks(peak_no).Width; - ind=(loc-10*wthis.Trace.x); + ind=(loc-5*wthis.Trace.x); x_peak=this.Trace.x(ind); y_peak=this.Trace.y(ind); end function save(this,varargin) %Parse inputs for saving p=inputParser; addParameter(p,'filename','placeholder',@ischar); addParameter(p,'save_dir',pwd,@ischar); addParameter(p,'overwrite_flag',false); addParameter(p,'save_prec',15); parse(p,varargin{:}); %Assign shorter names filename=p.Results.filename; save_dir=p.Results.save_dir; overwrite_flag=p.Results.overwrite_flag; save_prec=p.Results.save_prec; %Puts together the full file name fullfilename=fullfile([save_dir,filename,'.txt']); %Creates the file in the given folder write_flag=createFile(save_dir,fullfilename,overwrite_flag); %Returns if the file is not created for some reason if ~write_flag; return; end col_names={'Value','Location','Width','Prominence'}; n_columns=length(col_names); %Finds appropriate column width cw=max([cellfun(@(x) length(x),col_names), save_prec+7]); cw_vec=repmat(cw,1,4); pre_fmt_str=repmat('%%%is\\t',1,n_columns); fmt_str=sprintf([pre_fmt_str,'\r\n'],cw_vec); fileID=fopen(fullfilename,'w'); fprintf(fileID,fmt_str,col_names{:}); pre_fmt_str_nmb=repmat('%%%i.15e\\t',1,n_columns); %Removes the tab at the end pre_fmt_str_nmb((end-2):end)=[]; nmb_fmt_str=sprintf([pre_fmt_str_nmb,'\r\n'],cw_vec); fprintf(fileID,nmb_fmt_str,... [[this.Peaks.Value];[this.Peaks.Location];... [this.Peaks.Width];[this.Peaks.Prominence]]); fclose(fileID); end function loadPeaks(this,fullfilename) assert(ischar(fullfilename),... 'File name must be a char, currently it is a %s',... class(fullfilename)); if ~exist(fullfilename,'file') error('File named ''%s'' does not exist, choose a different file',... fullfilename); end LoadStruct=importdata(fullfilename); headers=regexp(LoadStruct.textdata{1},'\s*(\w*)\s*','Tokens'); this.Peaks=struct(headers{1}{1},num2cell(LoadStruct.data(:,1)),... headers{2}{1},num2cell(LoadStruct.data(:,2)),... headers{3}{1},num2cell(LoadStruct.data(:,3)),... headers{4}{1},num2cell(LoadStruct.data(:,4))); end function clearPeaks(this) - this.Peaks=struct(); + this.Peaks=struct('Location',[],'Width',[],'Prominence',[],... + 'Value',[]); end end end \ No newline at end of file diff --git a/@MyPeakFinderGui/MyPeakFinderGui.m b/@MyPeakFinderGui/MyPeakFinderGui.m index 3647b18..5f0ded8 100644 --- a/@MyPeakFinderGui/MyPeakFinderGui.m +++ b/@MyPeakFinderGui/MyPeakFinderGui.m @@ -1,219 +1,226 @@ classdef MyPeakFinderGui < handle properties PeakFinder Gui; axis_handle; end properties (Access=private) peak_color='r'; data_color='b'; peak_handle; end properties (Dependent=true) trace_handle; filename; session_name; base_dir; save_dir; end methods function this=MyPeakFinderGui() this.PeakFinder=MyPeakFinder(); createGui(this); end function delete(this) %Avoids loops set(this.Gui.Window,'CloseRequestFcn',''); %Deletes the figure delete(this.Gui.Window); %Removes the figure handle to prevent memory leaks this.Gui=[]; end function closeFigure(this, ~, ~) delete(this); end function analyzeCallback(this, src, ~) src.BackgroundColor=[1,0,0]; src.String='Analyzing..'; searchPeaks(this.PeakFinder,... 'MinPeakProminence',str2double(this.Gui.PromEdit.String),... 'MinPeakDistance',str2double(this.Gui.SepEdit.String),... 'FindMinima',this.Gui.MinimaCheck.Value); plotPeaks(this); src.BackgroundColor=[0.94,0.94,0.94]; src.String='Analyze'; end function plotPeaks(this) delete(this.peak_handle); this.peak_handle=plot(this.axis_handle,... [this.PeakFinder.Peaks.Location],... [this.PeakFinder.Peaks.Value],... 'Marker','o',... 'LineStyle','none',... 'Color',this.peak_color); end function fitPeakCallback(this,~,~) fitAllPeaks(this.PeakFinder,... 'base_dir',this.base_dir,... 'session_name',this.session_name,... 'filename',this.filename); end function clickCallback(this,~,~) switch this.Gui.Window.SelectionType case 'normal' %Left click addPeak(this); case 'alt' %Right click axis(this.axis_handle,'tight') case 'extend' %Shift click coords=this.axis_handle.CurrentPoint; removePeak(this, coords(1)); otherwise end end function windowScrollCallback(this, ~, event) coords=this.axis_handle.CurrentPoint; if event.VerticalScrollCount>0 %Zoom out zoomAxis(this,coords(1),0.1) else %Zoom in zoomAxis(this,coords(1),10); end end function removePeak(this, coord) [~,ind]=min(abs([this.PeakFinder.Peaks.Location]-coord)); this.PeakFinder.Peaks(ind)=[]; plotPeaks(this); end function addPeak(this) x_lim=[this.axis_handle.XLim(1),... this.axis_handle.XLim(2)]; + if checkPeakExist(this.PeakFinder,x_lim(1),x_lim(2)) + opts=struct('WindowStyle','modal','Interpreter','tex'); + errordlg(['A peak already exists in the window. ',... + 'To add a new peak, zoom such that there is ',... + 'no peak in the window'],'Error',opts); + return + end searchPeaks(this.PeakFinder,... 'ClearPeaks',false,... 'Limits',x_lim,... 'MinPeakProminence',str2double(this.Gui.PromEdit.String),... 'FindMinima',this.Gui.MinimaCheck.Value,... 'NPeaks',1,... 'SortStr','descend'); plotPeaks(this); end function zoomAxis(this,coords,zoom_factor) curr_width=this.axis_handle.XLim(2)-this.axis_handle.XLim(1); new_width=curr_width/zoom_factor; this.axis_handle.XLim=... [coords(1)-new_width/2,coords(1)+new_width/2]; end function clearCallback(this, ~, ~) delete(getLineHandle(this.PeakFinder.Trace,this.axis_handle)); clearData(this.PeakFinder.Trace); cla(this.axis_handle); end function loadTraceCallback(this, src, ~) %Window can find all files [fname,path]=uigetfile('*.*'); if fname==0 warning('No file was selected'); return end src.BackgroundColor=[1,0,0]; src.String='Loading..'; loadTrace(this.PeakFinder,[path,fname]); plotTrace(this.PeakFinder.Trace,this.axis_handle); this.trace_handle.ButtonDownFcn=... @(src, event) clickCallback(this, src, event); exitLoad(this); end function exitLoad(this) this.Gui.LoadTraceButton.BackgroundColor=[0.94,0.94,0.94]; this.Gui.LoadTraceButton.String='Load trace'; end function savePeaksCallback(this,~,~) save(this.PeakFinder,... 'save_dir',this.save_dir,... 'filename',this.filename); end function loadPeaksCallback(this,~,~) [fname,path]=uigetfile('*.*'); if fname==0 warning('No file was selected'); return end loadPeaks(this.PeakFinder,[path,fname]); plotPeaks(this); end function clearPeaksCallback(this,~,~) clearPeaks(this.PeakFinder); delete(this.peak_handle); end end methods function trace_handle=get.trace_handle(this) trace_handle=getLineHandle(this.PeakFinder.Trace,this.axis_handle); end function base_dir=get.base_dir(this) try base_dir=this.Gui.BaseEdit.String; catch base_dir=pwd; end end function set.base_dir(this,base_dir) this.Gui.BaseEdit.String=base_dir; end function session_name=get.session_name(this) try session_name=this.Gui.SessionEdit.String; catch session_name=''; end end function set.session_name(this,session_name) this.Gui.SessionEdit.String=session_name; end function filename=get.filename(this) try filename=this.Gui.FileNameEdit.String; catch filename='placeholder'; end end function set.filename(this,filename) this.Gui.FileNameEdit.String=filename; end %Get function from save directory function save_dir=get.save_dir(this) save_dir=createSessionPath(this.base_dir,this.session_name); end end end \ No newline at end of file diff --git a/Test functions/GorodetskyFit test/testgorodetsky.m b/Test functions/GorodetskyFit test/testgorodetsky.m new file mode 100644 index 0000000..e0fefab --- /dev/null +++ b/Test functions/GorodetskyFit test/testgorodetsky.m @@ -0,0 +1,15 @@ +load('testmat') +figure(124) +ax=gca; +testfit=MyFit('fit_name','Gorodetsky2000','x',xf,'y',yf,'plot_handle',ax,... + 'enable_gui',true,'enable_plot',true); +testfit.Data.plotTrace(ax); +hold on +testfit.genInitParams; +testfit.init_params(1) +testfit.init_params(3:end) +testfit.plotInitFun; +testfit.fitTrace; +testfit.plotFit('Color','r'); +testfit.Gof +testfit.FitInfo \ No newline at end of file diff --git a/Test functions/GorodetskyFit test/testmat.mat b/Test functions/GorodetskyFit test/testmat.mat new file mode 100644 index 0000000..426662e Binary files /dev/null and b/Test functions/GorodetskyFit test/testmat.mat differ diff --git a/Utility functions/Initial parameter estimation/initParamGorodetsky2000.m b/Utility functions/Initial parameter estimation/initParamGorodetsky2000.m new file mode 100644 index 0000000..c578c38 --- /dev/null +++ b/Utility functions/Initial parameter estimation/initParamGorodetsky2000.m @@ -0,0 +1,27 @@ +function [p_in,lim_lower,lim_upper]=initParamGorodetsky2000(x,y) +% { 'a','b','c','d', 'gamma','k0', 'kex'},... +lim_upper=[Inf,Inf,Inf,Inf,Inf]; +lim_lower=[-Inf,-Inf,-Inf,-Inf,-Inf]; + + +%Finds peaks on the negative signal (max 2 peaks) +[~,locs,widths,proms]=findpeaks(-y,x,... + 'MinPeakDistance',0.001*range(x),'SortStr','descend','NPeaks',2); + + +p_in(1)=max(y); + +%position +p_in(2)=mean(locs); + +if length(locs)==2 + p_in(3)=abs(diff(locs))/2; +else + p_in(3)=0; +end + +p_in(4)=mean(widths)/4; +%Assume critical coupling +p_in(5)=p_in(4); + +end \ No newline at end of file diff --git a/Utility functions/initParamDblLorentzianGrad.m b/Utility functions/initParamDblLorentzianGrad.m index 2c64f3b..4507077 100644 --- a/Utility functions/initParamDblLorentzianGrad.m +++ b/Utility functions/initParamDblLorentzianGrad.m @@ -1,63 +1,63 @@ function [p_in,lim_lower,lim_upper]=initParamDblLorentzianGrad(x,y) -%Assumes form a/pi*b/2/((x-c)^2+(b/2)^2)+d/pi*e/2/((x-f)^2+(e/2)^2))+g +%Assumes form a/pi*b/2/((x-c)^2+(b/2)^2)+d/pi*e/2/((x-f)^2+(e/2)^2))+g*x+h lim_upper=[Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf]; lim_lower=[-Inf,0,-Inf,-Inf,0,-Inf,-Inf,-Inf]; %Finds peaks on the positive signal (max 2 peaks) [~,locs{1},widths{1},proms{1}]=findpeaks(y,x,... 'MinPeakDistance',0.001*range(x),'SortStr','descend','NPeaks',2); %Finds peaks on the negative signal (max 2 peaks) [~,locs{2},widths{2},proms{2}]=findpeaks(-y,x,... 'MinPeakDistance',0.001*range(x),'SortStr','descend','NPeaks',2); %If the prominence of the peak in the positive signal is greater, we adapt %our limits and parameters accordingly, if negative signal has a greater %prominence, we use this for fitting. if isempty(proms{2}) || proms{1}(1)>proms{2}(1) ind=1; lim_lower(1)=0; lim_lower(4)=0; p_in(8)=min(y); else lim_upper(1)=0; lim_upper(4)=0; ind=2; p_in(8)=max(y); proms{2}=-proms{2}; end p_in(2)=widths{ind}(1); %Calculates the amplitude, as when x=c, the amplitude is 2a/(pi*b) p_in(1)=proms{ind}(1)*pi*p_in(2)/2; p_in(3)=locs{ind}(1); if length(locs{ind})==2 p_in(5)=widths{ind}(2); p_in(4)=proms{ind}(2)*pi*p_in(5)/2; p_in(6)=locs{ind}(2); else p_in(5)=widths{ind}(1); p_in(4)=proms{ind}(1)*pi*p_in(5)/2; p_in(6)=locs{ind}(1); end %If one of the lorentzians is found to be much smaller than the other, we %instead fit using only the greater lorentzian's parameters. This is an %adaption for very closely spaced lorentzians. if abs(p_in(1))>abs(10*p_in(4)) p_in(1)=p_in(1)/2; p_in(5)=p_in(2); p_in(6)=p_in(3); p_in(4)=p_in(1); end %Find gradient offset -p_in(7)=(y(end)-y())/(x(end)-x(1)); +p_in(7)=(y(end)-y(1))/(x(end)-x(1)); lim_lower(2)=0.01*p_in(2); lim_upper(2)=100*p_in(2); lim_lower(5)=0.01*p_in(5); lim_upper(5)=100*p_in(5); end \ No newline at end of file