%--------------------------------------------------------------------------
% function [F,fsr,fwhm] = FT_cavity_finesse_from_scan(xdata,ydata)
%
% A function for Matalb which measures the free spectral range (FSR) and
% full-width-half-maxima (FWHM) of data representing a cavity scan.  These
% values are used to calculate the finesse of the cavity.  This function
% scan data over at least 1 FSR, i.e. at least two peaks are
% required.  It also requires at least one full peak within the data, with
% atleast 1/4*fsr data either side of one of the main peaks, in order to 
% measure the fwhm (i.e. the minimum Finesse measureable by this function 
% is 4. 
%
% N.B. This function will not work with very 'dirty' data, i.e. when there
% are large peaks (greater than 90% of the main peak power) for lots of
% different order modes.
%
% xdata:    Vector of tuning values for the cavity scan.  Generally either
%           in degrees (i.e. tuning of mirror position) or hertz (i.e.
%           tuning of the laser frequency).
% ydata:    Power data comprising the cavity scan
%
% F:        Finesse
% fsr:      Free spectral range 
% fwhm:     Full width half maxima
%
% Part of the Simtools package, http://www.gwoptics.org/simtools
% Charlotte Bond    17.08.2012
%--------------------------------------------------------------------------
%

function [F,fsr,fwhm] = FT_cavity_finesse_from_scan(xdata,ydata)
    
    baseid = 'cavity_finesse_from_scan';

    % Maximum power in the scan
    PMAX = max(ydata);

    % Locate main peaks.  Use findpeaks as it is unlikley the subsequent
    % peaks are exactly equal to the maximum.
    [peaks,locs] = findpeaks(ydata,'minpeakheight',0.9*PMAX);
    
    % Error message if less than two peaks are detected
    if length(peaks)<2
        result='Invalid data length: scan data must include at least two main peaks';
        msgid=[baseid,':checkarguments'];
        error(msgid,result);
    end
    
    % Calculate fsr
    fsr = xdata(locs(2))-xdata(locs(1));
    
    % Cutting out data around main peak
    i_fsr = locs(2)-locs(1);
    i_end = length(xdata);
    
    i = 0;
    peak_no = 1;
    
    % Find peak in data with clear 1/4 fsr either side
    while i==0
    
        % Select peak indices
        peak_start = locs(peak_no)-round(i_fsr/4);
        peak_end = locs(peak_no)+round(i_fsr/4);
        
        % Check indices are valid
        if peak_start>=1 && peak_end<=i_end
            i=1;
        else
            peak_no = peak_no+1;
            % Error message if clear peak is not included in data
            if peak_no > length(peaks)
                result='Invalid data range: require one peak with 1/4 fsr either side';
                msgid=[baseid,':checkarguments'];
                error(msgid,result);
            end
        end
    end
         
    % Select complete peak data 
    peak_y = ydata(peak_start:peak_end);
    peak_x = xdata(peak_start:peak_end);

    R = 1;
    F_ = 0;
    j = 0;
    
    while j==0
        % Interpolate
        peak_y = interp(peak_y,R);
        peak_x = interp(peak_x,R);
    
        % Find scan data greater than half the maximum power
        [idx] = find(peak_y >= 0.5*PMAX(1));

        % Calculate fwhm
        fwhm = peak_x(idx(length(idx)))-peak_x(idx(1));

        % Calculate finesse
        F = fsr/fwhm;
        
        % Check if previous finesse matches the current value (to 2dp)
        if round(F*100)==round(F_*100)
            j=1;
        else
            R = R+1;
        end
        
        % Reset previous value of finesse
        F_ = F;
        
    end
        


end

