% Program to plot out atomic form factor of an element v scattering angle
% in a polar plot as a function of photon energy

clear; close all;

% Coefficients to determine dependence of f(Q) on Q/4pi (sin theta/lambda)
% for sulphur. Change according to element: look up in the International
% Tables of Crystallography

% 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('PolarDependence_qV2theta_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

% 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

figure('ToolBar','none');

open(vid);
set(0,'defaultfigurecolor',[1 1 1]);

theta = (-pi/2:pi/1000:pi/2)';
twoTheta = 2*theta;
fmax = a1+a2+a3+a4+c;

% Create rgb progression through a rainbow
bot = 8:8:128; % 0 to 7f in 16 steps
top = 128:8:248; % 80 to ff
all = 0:16:240; % 0 to ff
botc = 0*all; % vector of zeros
midc = botc + 127; % vector of 127s
topc = botc + 255; % vector of 255s
rcomp = [top topc topc flip(all) botc botc bot flip(bot)]/256; % R in rgb
gcomp = [botc bot top flip(top) midc flip(bot) botc botc]/256; % G in rgb
bcomp = [botc botc botc botc bot midc midc bot]/256; % B in rgb

%N = (hvmax - hvmin)/hvstep;
for i = hvmin:hvstep:hvmax
    lambda = 12398.4/i;
    q = (4*pi/lambda).*sin(theta);
    
    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 vector q
    
    clii = 1+floor(127*(i-hvmin)/(hvmax - hvmin)); % Index to look up rgb components
    
    polarplot(twoTheta,f,'color',[rcomp(clii), gcomp(clii), bcomp(clii)], 'LineWidth', 2.5);
    rlimits = [0 fmax*1.025];
    rlim(rlimits);
    %axis off;
    
    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];
    hText = text(-pi/2, 1.25*fmax, strTot, 'FontSize',14);
    frame = getframe(gcf);
    writeVideo(vid,frame);
end
% Output the movie as an mp4 file
close(vid);