% Program to plot out atomic form factor of an element v scattering vector
% and track the maximum accessible scattering vector for a given photon
% energy

clear; close all;

myBlue = [0.4 0.44 0.73];

% Prompt user for input file containing 9 parameters to generate form
% factor. If no file is read, the parameters for sulphur are taken
prompt = 'Do you want to input an elemental parameter file (1 x 9)? Y/N [N]: ';
str = input(prompt,'s');
if isempty(str) || (str == 'N') || (str == 'n')
    % Parameters needed to generate form factor for sulphur
    a1 = 6.9053;
    b1 = 1.4679;
    a2 = 5.2034;
    b2 = 22.2151;
    a3 = 1.4379;
    b3 = 0.2536;
    a4 = 1.5863;
    b4 = 56.172;
    c = 0.8669;
    vid = VideoWriter('fVq_sulphur.mp4','MPEG-4');
else
    prompt = 'Name of input file? : ';
    filestr = input(prompt,'s');
    fileID = fopen(filestr,'r');
    formatSpec = '%f';
    ffPars = fscanf(fileID,formatSpec);
    a1 = ffPars(1);
    b1 = ffPars(2);
    a2 = ffPars(3);
    b2 = ffPars(4);
    a3 = ffPars(5);
    b3 = ffPars(6);
    a4 = ffPars(7);
    b4 = ffPars(8);
    c = ffPars(9);
    prompt = 'Name of output mp4 movie file (including .mp4)? : ';
    mp4str = input(prompt,'s');
    vid = VideoWriter(mp4str,'MPEG-4');
end

fmax = a1+a2+a3+a4+c;

% Prompt user for input parameters
prompt = 'Minimum photon energy in eV? [default = 500 eV]: ';
hvmin = input(prompt);
if isempty(hvmin)
    hvmin = 500;
end

prompt = 'Maximum photon energy in eV? [default = 30000]: ';
hvmax = input(prompt);
if isempty(hvmax)
    hvmax = 30000;
end

prompt = 'Step size in photon energy in eV? [default = 50]: ';
hvstep = input(prompt);
if isempty(hvstep)
    hvstep = 50;
end

qmax = 4*pi*hvmax/12398.4; % Maximum accessible Q at 2theta = 180 degrees
% at maximum photon energy
q = (0:0.05:qmax);

f1 = a1*exp(-b1*(q/(4*pi)).^2);
f2 = a2*exp(-b2*(q/(4*pi)).^2);
f3 = a3*exp(-b3*(q/(4*pi)).^2);
f4 = a4*exp(-b4*(q/(4*pi)).^2);
f5 = c;
f = f1 + f2 + f3 + f4 + f5;     % Scattering form factor


N=(hvmax-hvmin)/hvstep; % Number of frames

open(vid);
figure('ToolBar','none');


for i = hvmin:hvstep:hvmax
    lambda = 12398.4/i; % in Angstrom
    qnow = (4*pi/lambda); % max accessible Q in reciprocal Angstrom
    % for present photon energy
    amin = 2*pi/qnow;
    f1 = a1*exp(-b1*(qnow/(4*pi))^2);
    f2 = a2*exp(-b2*(qnow/(4*pi))^2);
    f3 = a3*exp(-b3*(qnow/(4*pi))^2);
    f4 = a4*exp(-b4*(qnow/(4*pi))^2);
    f5 = c;
    fnow = f1 + f2 + f3 + f4 + f5;     % Max scattering vector for ith energy
    
    plot(0.0,0.0,qmax,fmax,q,f,'color',myBlue, 'LineWidth', 2.5);
    hold on
    plot(qnow,fnow,'ro','MarkerFaceColor','r','MarkerSize',10)  % marking the ith data point of x and y
    hold off
    xaxisLabel = get(gca,'XTickLabel');
    set(gca,'FontName','Helvetica','fontsize',18);
    set(gca,'TickLength',[0.016, 2]);
    xlim([0.0 qmax])
    ylim([0 fmax+1])
    
    str1 = 'Q [';
    str2 = char(197);
    str3 = '^{-1}]';
    strTot = [str1,str2,str3];
    set(gca,'linewidth',2);
    xlabel(strTot);
    ylabel('f(Q)');
    
    str1 = 'E = ';
    str2 = num2str(i,'% 4.0f');
    str3 = ' eV; ';
    str4 = '{\lambda} = ';
    str5 = num2str(lambda,'%4.3f');
    str6 = ' ';
    str7 = char(197);
    strTot = [str1,str2,str3,str4,str5,str6,str7];
    str8 = 'a_{min} = ';
    str9 = num2str(amin,'%4.3f');
    str10 = ' ';
    strTot2 = [str8,str9,str10,str7];
    hText = text(qmax/16, 0.1*fmax, strTot, 'FontSize',14);
    hText2 = text(qmax/1.4, 0.91*fmax, strTot2, 'FontSize',14);
    
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
% Output the movie as an mp4 file
close(vid);