%% Input
filename = 'XXX';

% Define Rectangular Box for Line Profile Extraction
x_Pos = 0.54;
x_size = 0.075;
x_size_background = 0.02;
y_Pos = 0;
y_size = 86;
regridding_x = 20;
regridding_y = 300;

%% Import Peaklist to Help Identify Useful Line Profiles
opts = delimitedTextImportOptions("NumVariables", 14);
opts.DataLines = [2, Inf];
opts.Delimiter = "\t";
opts.VariableNames = ["h", "k", "l", "q_x", "q_y", "q_z", "q_xy", "q",...
    "realSF", "imagSF", "absSF", "absSF2", "absSF2Rotation", "d"];
opts.VariableTypes = ["string", "string", "string", "double",...
    "double", "double", "double", "double", "double", "double",...
    "double", "double", "double", "double"];
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
tbl = readtable(".\peaklist.txt", opts);
% Transform table into column vectors
q_xy = tbl.q_xy;                                                            % q_z & q_xy for all possible peaks
q_z = tbl.q_z;
q_ref = tbl.q;
Psi_ref = 180./pi.*atan(q_xy./q_z);                                         % calculate Psi from q_xy and q_z
clear tbl opts q_xy q_z

%% Import Measured Data Transformed to Reciprocal Spherical Coordinates
opts = delimitedTextImportOptions("NumVariables", 3);
opts.DataLines = [2, Inf];
opts.Delimiter = "\t";
opts.VariableNames = ["q", "Psi", "I"];
opts.VariableTypes = ["double", "double", "double"];
opts.ExtraColumnsRule = "ignore";
opts.EmptyLineRule = "read";
tbl_data = readtable(".\"+filename+'.txt', opts);
% Transform table into column vectors
q = tbl_data.q;
Psi = tbl_data.Psi;
I = tbl_data.I;
clear opts tbl_data

q_z = q.*(tan(Psi.*pi./180).^2+1).^(-0.5);                                  % transform to cylindrical coordinates to remove data at too small alpha_f

% Apply q-dependent Part of Lorentz Correction for Spherical Reciprocal Space Coordinates
I = I.*q.^2;

% Remove shaded areas
I(I<=0) = NaN;

% Loop over all line Profiles
for k = 1:numel(x_Pos)
    %% Plot Data with Expected Peak Positions and Rectangle for Line Profile Extraction
    x_lim = [0,2];                                                          % x-axis limits
    y_lim = [-90,90];                                                       % y-axis limits
    c_lim = [40,150];                                                       % colormap limits
    
    f = figure(1);
    clf(f,'reset')
    set(f,'DefaultAxesFontSize',15)
    scatter3(q,Psi,sqrt(I),40,sqrt(I),'.')                                  % plot scattered data
    hold on
    [black_x,black_y] = meshgrid([0,max(x_lim)],[-180,180]);                % add a black underground to the plot
    surf(black_x,black_y,zeros(size(black_x)),EdgeColor="none",...
        FaceColor='k')
    
    scatter3(q_ref,-Psi_ref,max(sqrt(I))*ones(size(q_ref)),60,'wo')         % plot circles for theoretical peak positions
    scatter3(q_ref,Psi_ref,max(sqrt(I))*ones(size(q_ref)),60,'wo')
    
    % Display Box for the Line Profile Extraction
    x_rec = [x_Pos(k)-x_size(k),x_Pos(k)+x_size(k)];
    y_rec = [y_Pos(k)-y_size(k);y_Pos(k)+y_size(k)];
    xx_rec = [x_rec(1);x_rec(1);x_rec(2);x_rec(2);x_rec(1)];
    yy_rec = [y_rec(:);flipud(y_rec(:));y_rec(1)];
    plot3(xx_rec(:),yy_rec(:),max(sqrt(I))*ones(size(xx_rec(:))),'-w')      % plot a white rectangle for the box for evaluation
    
    % Display Boxes for the Undergound Fit
    x_rec_1 = [x_Pos(k)-x_size(k),x_Pos(k)-x_size(k)+x_size_background(k)];
    x_rec_2 = [x_Pos(k)-x_size_background(k)+x_size(k),x_Pos(k)+x_size(k)];
    xx_rec_1 = [x_rec_1(1);x_rec_1(1);x_rec_1(2);x_rec_1(2);x_rec_1(1)];
    xx_rec_2 = [x_rec_2(1);x_rec_2(1);x_rec_2(2);x_rec_2(2);x_rec_2(1)];
    plot3(xx_rec_1(:),yy_rec(:),max(sqrt(I))*ones(size(xx_rec_1(:))),'-w')  % plot white rectangles for the boxes for the underground fit
    plot3(xx_rec_2(:),yy_rec(:),max(sqrt(I))*ones(size(xx_rec_2(:))),'-w')
    
    % Plot Settings
    colormap(turbo)
    clim(c_lim)
    view(0,90)
    xlim(x_lim)
    ylim(y_lim)
    xlabel(['q / ',char(197)])
    ylabel('\Psi / °')
    zlabel('sqrt(I) / a.u.')
    set(f,'Position',[1925 45 980 948])
    clear black_x black_y xx_rec yy_rec xx_rec_1 xx_rec_2 x_lim y_lim c_lim
    
    %% Calculate Underground
    % Regridding
    q_size_underground = 1;
    q_size = regridding_x(k);
    Psi_size = regridding_y(k);
    index_new = x_rec_1(2)<q & q<x_rec_2(1) & y_rec(1)<Psi & Psi<y_rec(2);
    
    index_underground_l = x_rec_1(1)<q & q<x_rec_1(2) & y_rec(1)<Psi & Psi<y_rec(2);
    I_regrid_l = regrid([Psi(index_underground_l);y_rec(1);y_rec(2);y_rec(1);y_rec(2)],...
        [q(index_underground_l);x_rec_1(2);x_rec_1(2);x_rec_1(1);x_rec_1(1)],...
        [I(index_underground_l);NaN;NaN;NaN;NaN],[Psi_size,q_size_underground]);
    index_underground_r = x_rec_2(1)<q & q<x_rec_2(2) & y_rec(1)<Psi & Psi<y_rec(2);
    I_regrid_r = regrid([Psi(index_underground_r);y_rec(1);y_rec(2);y_rec(1);y_rec(2)],...
        [q(index_underground_r);x_rec_2(2);x_rec_2(2);x_rec_2(1);x_rec_2(1)],...
        [I(index_underground_r);NaN;NaN;NaN;NaN],[Psi_size,q_size_underground]);
    q_regrid_l = min(q(index_underground_l))+(max(q(index_underground_l))-min(q(index_underground_l)))./2;
    q_regrid_r = min(q(index_underground_r))+(max(q(index_underground_r))-min(q(index_underground_r)))./2;
    q_regrid = linspace(x_rec_1(2),x_rec_2(1),q_size);
    Psi_regrid = linspace(y_rec(1),y_rec(2),Psi_size);
      
    % Calculate line connecting the left and right underground
    incl_underground = (I_regrid_r-I_regrid_l)./(q_regrid_r-q_regrid_l);
    d_underground = I_regrid_r-incl_underground.*q_regrid_r;
    
    I_underground = incl_underground*q_regrid+d_underground;

    % Plot Data with UndergroundFit
    g = figure(2);
    clf(g,'reset')
    % set(g,'DefaultAxesFontSize',15)
    plot3(q(index_new),Psi(index_new),I(index_new),'k.')                    % plot data points in black
    hold on
    plot3(q(index_underground_l),Psi(index_underground_l),...
        I(index_underground_l),'r.')                                        % plot underground data points in red
    plot3(q(index_underground_r),Psi(index_underground_r),...
        I(index_underground_r),'r.')                                        % plot underground data points in red
    q_underground_plot = [q_regrid_l,q_regrid,q_regrid_r];
    I_underground_plot = incl_underground*q_underground_plot+d_underground;
    [qq_regrid,PsiPsi_regrid] = meshgrid(q_underground_plot,Psi_regrid);
    s_background = surf(qq_regrid,PsiPsi_regrid,I_underground_plot,EdgeColor="none",FaceColor="interp");
    alpha(s_background,0.8)
    clear I_underground_plot q_underground_plot qq_regrid PsiPsi_regrid
    clim([0,max(I(index_new))])
    view(60,10)
    xlim([x_rec_1(1),x_rec_2(2)])
    ylim(y_rec)
    xlabel(['q / ',char(197)])
    ylabel('\Psi / °')
    zlabel('I / a.u.')
    set(g,'Position',[1925+985 45+340 925 608])
    clear index_underground_r index_underground_l I_regrid_l I_regrid_r
   
    %% Perform Regridding, Extract Line Profile  
    I_regrid = regrid([Psi(index_new);y_rec(1);y_rec(2);y_rec(1);y_rec(2)],...
        [q(index_new);x_rec_1(2);x_rec_1(2);x_rec_2(1);x_rec_2(1)],...
        [I(index_new);NaN;NaN;NaN;NaN],[Psi_size,q_size]);                  % NaN is added in every corner to compensate for missing data at high q

    % Remove Underground from Intensity Data
    I_peak = I_regrid-I_underground;
    
    %% Plot Regridded Data
    [qq_regrid,PsiPsi_regrid] = meshgrid(q_regrid,Psi_regrid);
    s_regrid = surf(qq_regrid,PsiPsi_regrid,I_regrid,EdgeColor="none",FaceColor="interp");
    alpha(s_regrid,0.5)
    colormap(turbo)
    clear index_new qq_regrid PsiPsi_regrid s_regrid

    % Extract Line Profile
    remove = sum(~isnan(I_peak),2);                                         % determine the number of entries for each Psi value
    I_peak(isnan(I_peak)) = 0;                                              % set NaN values to 0 to perform integration
    I_line = trapz(q_regrid,I_peak,2);                                      % independent 1D integration for each Psi value, a simple summation does not work when comparing peaks at different q
    remove(remove < 0.80*q_size) = NaN;                                     % remove data where we maily sum over background i.e. less than 80 % of data is available for a specific Psi value
    I_line(isnan(remove)) = NaN;
    
    normalization = ~isnan(I_line);                                         % normalization necessary because of missing data in positive or negative Psi regime
    I_line(isnan(I_line)) = 0;                                              % set NaN values to 0 to perform summation
    I_line_sum = (I_line+flipud(I_line))./...
        (normalization+flipud(normalization));                              % average the positive and negative Psi regime to improve statistics and reduce missing data
    
    I_line(normalization == 0) = NaN;                                       % bring back NaN values where no data was available in both positive and negative Psi regime
    I_line_sum(normalization+flipud(normalization) == 0) = NaN;

    %% Plot Extracted Line Profile
    h = figure(3);
    clf(h,'reset')
    plot(Psi_regrid,I_line,'r',LineWidth=2)
    hold on
    plot(Psi_regrid,I_line_sum,'k')
    xlim([-90,90])
    set(h,'Position',[1925+985 45 925 253])
    
    clear remove normalization I_line
    
    % Save .txt file with Extracted Radial Line Profile
    LineProfile = table(Psi_regrid(Psi_size/2+1:end).',I_line_sum(Psi_size/2+1:end));
    writetable(LineProfile,"./LineProfile_q"+num2str(x_Pos(k))+"_"+filename+".txt",'Delimiter','\t')
end
% Save .txt with Extraction Settings
LineProfileSpecs = table(x_Pos,x_size,x_size_background,y_Pos,y_size,regridding_x,regridding_y);
writetable(LineProfileSpecs,"./LineProfileSpecs_"+filename+".txt",'Delimiter','\t')

