% Program tracking f1 and f2 for the elements (gallium is the default
% element) as a function of photon energy and an animation of the vector
% addition of different components of f

clear; close all;

myBlue = [0.4 0.44 0.73];

scale=1;

% Prompt user for input file containing 3 columns of hv, f_1(Q = 0), and
% f_2(Q = 0). If no file is read, the data for gallium is taken
prompt = 'Do you want to input a form factor dispersion file (hv, f1, f2)? Y/N [N]: ';
str = input(prompt,'s');
if isempty(str) || (str == 'N') || (str == 'n')
    % Data for gallium
    f12 = importdata('Gaf1f2.dat');
    vid = VideoWriter('Ga_f1f2Venergy.mp4','MPEG-4');
else
    prompt = 'Name of input file? : ';
    filestr = input(prompt,'s');
    f12 = importdata(filestr);
    prompt = 'Name of output mp4 movie file (including .mp4)? : ';
    mp4str = input(prompt,'s');
    vid = VideoWriter(mp4str,'MPEG-4');
end

datSizef12 = size(f12); % Get size of data file
N = datSizef12(1); % Number of data points
energy = f12(:,1); % Photon energy
f1plot = f12(:,2); % f_1(Q = 0) = f^0(Q = 0) + f' (f' assumes negative values)
f2plot = f12(:,3); % f_2 = f"
f1max = max(f1plot);
f2max = max(f2plot);

prompt = 'Value of f^0(Q = 0)? [default = 31, Ga]: ';
f00 = input(prompt);
if isempty(f00)
    f00 = 31;
end


prompt = 'Plot in eV [e] or keV [k]: ';
str = input(prompt,'s');
if (str == 'e') || (str == 'E')
    energy = 1000*energy;
    edatamin = min(energy);
    edatamax = max(energy);
    
    % Prompt user for photon-energy range to be plotted
    str1 = 'Minimum photon energy in eV? (> ';
    str2 = num2str(edatamin);
    str3 = ' eV) [default minimum of input data] : ';
    strTot1 = [str1,str2,str3];
    
    prompt = strTot1;
    hvmin = input(prompt);
    if isempty(hvmin)
        hvmin = edatamin;
    end
    
    str1 = 'Maximum photon energy in eV? (< ';
    str2 = num2str(edatamax);
    str3 = ' eV) [default maximum of input data] : ';
    strTot1 = [str1,str2,str3];
    streVkeV = '[eV]';
    prompt = strTot1;
    hvmax = input(prompt);
    if isempty(hvmax)
        hvmax = edatamax;
    end
else % Plot in keV
    edatamin = min(energy);
    edatamax = max(energy);
    
    % Prompt user for photon-energy range to be plotted
    str1 = 'Minimum photon energy in keV? (> ';
    str2 = num2str(edatamin);
    str3 = ' keV) [default minimum of input data] : ';
    strTot1 = [str1,str2,str3];
    streVkeV = '[keV]';
    prompt = strTot1;
    hvmin = input(prompt);
    if isempty(hvmin)
        hvmin = edatamin;
    end
    
    str1 = 'Maximum photon energy in keV? (< ';
    str2 = num2str(edatamax);
    str3 = ' keV) [default maximum of input data] : ';
    strTot1 = [str1,str2,str3];
    prompt = strTot1;
    hvmax = input(prompt);
    if isempty(hvmax)
        hvmax = edatamax;
    end
end

figure('units','pixels','position',[0 0 1920 1080],'ToolBar','none');
set(0,'defaultfigurecolor',[1 1 1]);
vid.FrameRate = 60;    % Default 30
vid.Quality = 100;    % Default 75

% Determine range to plot out
minindex = 1;
maxindex = 1;
for ii = 2:N
    if (energy(ii)<=hvmin)
        minindex = ii;
    end
    if (energy(ii)>hvmax) && (energy(ii-1)<hvmax)
        maxindex = ii-1;
    elseif (energy(ii)==hvmax)
        maxindex = ii;
    end
end

open(vid);
for ii = minindex:maxindex % all data rows
    newplot
    % subplot 1
    pos1 = [0.25 0.43 0.5 0.5]; % x, y, Dx, Dy
    subplot('Position',pos1) % Plot of f1 and f2 v hv
    
    plot(hvmin,hvmax,0.0,f1max,energy(1:N),f1plot(1:N),'color',myBlue, 'LineWidth', 2.5);
    hold on
    plot(hvmin,hvmax,0.0,f1max,energy(1:N),f2plot(1:N),'color','g', 'LineWidth', 2.5);
    plot(energy(ii),f1plot(ii),'o','MarkerFaceColor',myBlue,'MarkerSize',10,...
        'MarkerEdgeColor','none')  % marking the ith data point of x and y
    plot(energy(ii),f2plot(ii),'o','MarkerFaceColor','g','MarkerSize',10,...
        'MarkerEdgeColor','none')  % marking the ith data point of x and y
    
    hold off
    set(gca,'FontName','Helvetica','fontsize',18);
    set(gca,'TickLength',[0.016, 2]);
    xlim([hvmin hvmax])
    ylim([0 1.1*f1max])
    
    str1 = 'f_1';
    str2 = 'f_2';
    strTot = [str1,str2,str3];
    text(hvmin+(hvmax-hvmin)/10, 0.989*f1max, str1, 'FontSize',18, 'color', myBlue);
    text(hvmin+(hvmax-hvmin)/10, 0.22*f1max, str2, 'FontSize',18, 'color', 'green');
    set(gca,'linewidth',2);
    
    xlabstr1 = 'Photon energy ';
    xlabstr = [xlabstr1,streVkeV];
    xlabel(xlabstr);
    ylabel('f_1, f_2');
    
    % subplot 2
    
    pos2 = [0.2 0.1 0.6 f2max/f00]; % x, y, Dx, Dy
    subplot2 = subplot('Position',pos2); % Plot of f1 and f2 v hv
    
    xlim(subplot2,[-f00/20 21*f00/20]);
    ylim(subplot2,[-f2max/20 21*f2max/20]);
    
    % Plot quiver vectors
    axis equal
    
    quiver(0,0,f00,0,scale, 'LineWidth', 2.5, 'MaxHeadSize', 0.1, 'Color', 'red');
    hold on
    quiver(0,0,0,f2max,scale, 'LineWidth', 0.01, 'MaxHeadSize', 0.01, 'Color', 'white');
    quiver(f00,0,-f00+f1plot(ii),0,scale, 'LineWidth', 2.5, 'MaxHeadSize', 0.4, 'Color', myBlue);
    quiver(f1plot(ii), 0,0,f2plot(ii),scale, 'LineWidth', 2.5, 'MaxHeadSize', 0.4, 'Color', 'g');
    quiver(0,0,f1plot(ii),f2plot(ii),scale, 'LineWidth', 3.5, 'MaxHeadSize', 0.1, 'Color', 'black');
    
    xlim(subplot2,[-f00/20 21*f00/20]);
    ylim(subplot2,[-f2max/20 21*f2max/20]);
    hold off;
    
    axis off
    
    % Store the frame
    frame = getframe(gcf);
    writeVideo(vid,frame);
end

% Output the movie as an mp4 file
close(vid);
