% ShekieFactor.m

% Crystallite peak fit

% Created by Jonathan Ogle July 26, 2018
% Modified by Jonathan Ogle Aug 8, 2018
% Modified by Jonathan Ogle Feb 14, 2019

% This code is designed to quickly quantify Herman's orientation (HF),
% and the mosaicity factor (MF) for I vs Phi data

% Code suffix changed to *.txt for supplementary upload to use code
% change suffix to *.m

%-------------------------------------------------------------------------

% modifiable variables

% initial fit conditions
fit_from = [0]; % phi_range for fitting start point for two sets write in the form '[X Y];'
fit_to = [90]; % phi_range for fitting end point 
fit_phi = [0]; % what angle are you trying to quantify?

num_fits = length(fit_from);

% Non oriented data subtraction, if 1, an average intensity from the
% desired range will be collected and subtracted from the overall
% intensity, and the data will be re-normalized.

subtract_check = 0; %if you want to quantify any intensity that is non-oriented put 1, otherwise 0
sub_start = 55; % phi start point where data is assumed non-oriented
sub_end = 88; %phi end point where data is assumed non-oriented

%-------------------------------------------------------------------------

% Input phi data

A = exist('path1');
if A == 1
[file1,path1] = uigetfile('.txt','Select Phi Data',path1);
else
[file1,path1] = uigetfile('.txt','Select Phi Data');
end

pathfile1 = strcat(path1,file1);

Phi_data = dlmread(pathfile1);

Phi_data(:,2) = Phi_data(:,2)/max(Phi_data(:,2));
Phi_data

tic

A = exist('fit');
if A == 1
    clear fit
end

% If the code is written from 90 to zero, flip to be zero to 90

if Phi_data(1,1) > Phi_data(length(Phi_data),1)
    Phi_data(1:length(Phi_data),:) = flipud(Phi_data);
end

%-------------------------------------------------------------------------

% Determine how many gaussian peaks that need to be fit

idx_fit_from = zeros(num_fits,1);
idx_fit_to = zeros(num_fits,1);

for i = 1:num_fits
    idx_fit_from(i) = find(Phi_data(:,1) >= fit_from(i),1,'first');
    idx_fit_to(i) = find(Phi_data(:,1) <= fit_to(i),1,'last');
end


%-------------------------------------------------------------------------

% Subtract non-oriented area from plot using user input for phi range where
% data is non-oriented
if num_fits > 1
    Percent_oriented_fit = zeros(num_fits,1);
end


if subtract_check == 1
    
    if Phi_data(1,1) < Phi_data(length(Phi_data),1)
    sub_start_idx = find(Phi_data(:,1)>=sub_start,1,'first');
    sub_end_idx = find(Phi_data(:,1)>=sub_end,1,'first');
    else
    sub_start_idx = find(Phi_data(:,1)>=sub_end,1,'last');
    sub_end_idx = find(Phi_data(:,1)>=sub_start,1,'last'); 
    end
    
    Phi_data_raw = Phi_data;
    Phi_data_sub = Phi_data;
    Phi_data_sub(:,2) = mean(Phi_data(sub_start_idx:sub_end_idx,2));
    
    Percent_non_oriented = sum(Phi_data_sub(:,2))/sum(Phi_data_raw(:,2))
    Percent_oriented = 1-Percent_non_oriented
   
    if num_fits > 1
    for i = 1:num_fits
        Percent_oriented_fit(i) = sum(Phi_data_raw(idx_fit_from(i):idx_fit_to(i),2)-Phi_data_sub(idx_fit_from(i):idx_fit_to(i),2))/sum(Phi_data_raw(:,2));
    end
    else
        Percent_oriented_fit = Percent_oriented;
   
    end
    
    Phi_data(:,2) = (Phi_data_raw(:,2)-Phi_data_sub(:,2))./max((Phi_data_raw(:,2)-Phi_data_sub(:,2)));
    figure(1)
    plot(Phi_data_raw(:,1),Phi_data_raw(:,2),Phi_data(:,1),Phi_data(:,2),Phi_data_sub(:,1),Phi_data_sub(:,2))
else
    
    Percent_oriented = 0;
    Percent_non_oriented = 1;
    figure(1)
    plot(Phi_data(:,1),Phi_data(:,2))
    
    if num_fits > 1
    for i = 1:num_fits
        Percent_oriented_fit(i) = sum(Phi_data(idx_fit_from(i):idx_fit_to(i),2))/sum(Phi_data(:,2));
    end
    else
        Percent_oriented_fit = 1;
   
    end
end

%-------------------------------------------------------------------------

% create the amplitude and phi data for each respective region

Phi_shift = Phi_data(:,1)+1; % create a shift so that I can find the first index >= 0, by searching for 1;

for i = 1:num_fits
    
    fit(i).amp = Phi_data(idx_fit_from(i):idx_fit_to(i),2);
    fit(i).phi = Phi_data(idx_fit_from(i):idx_fit_to(i),1);
    
end

%--------------------------------------------------------------------------

% quantify the MF

MF = zeros(num_fits,1);

for i = 1:num_fits
    
    % correct for amplitude at the edges (signal overlap, only 50% of the
    % signal contributes)
    if fit(i).phi(1) == 0
        fit(i).amp(1) = fit(i).amp(1)/2;
    end
    if fit(i).phi(length(fit(i).phi)) == 90
        fit(i).amp(length(fit(i).phi)) = fit(i).amp(length(fit(i).phi))/2;
    end
    
    MF(i) = sum(fit(i).amp./sum(fit(i).amp).*(45-abs(fit_phi(i)-fit(i).phi))./45);
    
    % correct data back to original form
    
    if fit(i).phi(1) == 0
        fit(i).amp(1) = fit(i).amp(1)*2;
    end
    if fit(i).phi(length(fit(i).phi)) == 90
        fit(i).amp(length(fit(i).phi)) = fit(i).amp(length(fit(i).phi))*2;
    end
    
end    

%--------------------------------------------------------------------------

 figure(2)
 plot(Phi_data(:,1),Phi_data(:,2),'k')
 
%--------------------------------------------------------------------------

% quantify the HF on each section, and on the total

if subtract_check == 1
    HOF_data = Phi_data_raw;
else
    HOF_data = Phi_data;
end

if min(HOF_data(:,2))<0
    HOF_data(:,2) = (HOF_data(:,2)-min(HOF_data(:,2)))./max(HOF_data(:,2)+min(HOF_data(:,2)));
end

HF = (3*sum(HOF_data(:,2).*sind(HOF_data(:,1)).*cosd(HOF_data(:,1)).^2)/sum(HOF_data(:,2).*sind(HOF_data(:,1)))-1)/2;

HOF_fit_data = Phi_data;

if min(HOF_fit_data(:,2))<0
    HOF_fit_data(:,2) = (HOF_fit_data(:,2)-min(HOF_fit_data(:,2)))./max(HOF_fit_data(:,2)+min(HOF_fit_data(:,2)));
end

HF_fit = zeros(num_fits,1);

for i = 1:num_fits
    
    HF_fit(i) = (3*sum(HOF_fit_data(idx_fit_from(i):idx_fit_to(i),2).*sind(HOF_fit_data(idx_fit_from(i):idx_fit_to(i),1)).*cosd(HOF_fit_data(idx_fit_from(i):idx_fit_to(i),1)).^2)/sum(HOF_fit_data(idx_fit_from(i):idx_fit_to(i),2).*sind(HOF_fit_data(idx_fit_from(i):idx_fit_to(i),1)))-1)/2;

end

Phi_data
Percent_oriented_fit
MF
HF
HF_fit

 
%--------------------------------------------------------------------------

% Combine data and export

A = exist('path3');
if A == 1
[file3,path3] = uiputfile('*.txt','Save data',path3);
else
[file3,path3] = uiputfile('*.txt','Save data',path1);
end

% Create header information

if subtract_check ==1
header_text = [strcat("Data file 1: ",file1);...
                strcat("Fit from: ",num2str(fit_from));...
                strcat("Fit to: ",num2str(fit_to));...
                strcat("Fit phi: ",num2str(fit_phi));...
                strcat("Subtract start: ",num2str(sub_start));...
                strcat("Subtract end: ",num2str(sub_end));...
                strcat("Percent non-oriented: ",num2str(Percent_non_oriented));...
                strcat("Percent oriented: ",num2str(Percent_oriented));...
                strcat("Hermans orientation: ",num2str(HF));...
                "-------------------------";...
                "Fit	OP	MF	HF (fit)"];
else
header_text = [strcat("Data file 1: ",file1);...
                strcat("Fit from: ",num2str(fit_from));...
                strcat("Fit to: ",num2str(fit_to));...
                strcat("Fit phi: ",num2str(fit_phi));...
                strcat("Subtract start: 0");...
                strcat("Subtract end: 0");...
                strcat("Percent non-oriented: ",num2str(Percent_non_oriented));...
                strcat("Percent oriented: ",num2str(Percent_oriented));...
                strcat("Hermans orientation: ",num2str(HF));...
                "-------------------------";...
                "Fit#	OP	MF	HF(fit)"];

end

% Create orientation information

orientation_data = zeros(num_fits,4);

for i = 1:num_fits
    orientation_data(i,:) = [i,Percent_oriented_fit(i),MF(i),HF_fit(i)];
end

Phi_Data_Fit = [Phi_data(:,1) Phi_data(:,2)]';

pathfile3 = strcat(path3,file3);

fileID = fopen(strcat(pathfile3),'w+');
fprintf(fileID,'%s\n',header_text);
fprintf(fileID,'%d\t%6.3f\t%6.3f\t%6.3f\n',orientation_data');
fprintf(fileID,'-------------------------\n');
fprintf(fileID,'Phi\tData\n');
fprintf(fileID,'%6.4f\t%6.4f\n',Phi_Data_Fit);
fclose(fileID);