%Copyright (c) 2010, Luuk J.G.W. van Wilderen
%
%single_pulse1_frac.m

function single_pulse1_frac
%This function requires the Symbolic Toolbox
%The first time this function is called, it is rather slow as Matlab needs
%to load the Toolbox. Runs after that will go much faster.
%Try these values for realistic plots (Example is Phytochrome Cph1)
% Excitation wavelength = 640 nm;
% Laser pulse power = 1.65 microJoule;
% Diameter of laser beam in micrometers = 120;
% Quantum yield of reaction = 0.1;
% Absorption cross-section of molecule = 3.3e-17;

%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%Don't change anything below
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
syms x N1 n1 n2;  %Symbolic paramaters

pulse1_wave=1e-9*str2double(get(findobj('Tag','pulse1_wave_box'),'string')); %nanometers
pulse1_power=1e-6*str2double(get(findobj('Tag','pulse1_power_box'),'string')); %microJoules
pulse1_size=1e-4*str2double(get(findobj('Tag','pulse1_size_box'),'string')); %centimeters
%beam diameter in micrometer, so converted do cm (because of
% cross-section)
pulse1_shape=get(findobj('Tag','pulse1_shape_box'),'Value');
pulse1_QY=str2double(get(findobj('Tag','pulse1_QY_box'),'string'));
pulse1_cross_sec=str2double(get(findobj('Tag','pulse1_cross_sec_box'),'string'));
% pulse1_pH=str2double(get(findobj('Tag','pulse1_pH_box'),'string'));
% %Not needed for any calculation

%Input checking
if max(isnan([pulse1_wave pulse1_power pulse1_size pulse1_QY pulse1_cross_sec]))==1
    errordlg('Please fill in all boxes to define Pulse 1.')
    error('Please fill in all boxes to define Pulse 1.')
end

%##################################################
phot_energy=6.626068e-34*2.99792e8/(pulse1_wave);   %Planck, speed of light and excitation wavelength
if pulse1_shape==0
    phot=pulse1_power/phot_energy;   %nr photons/pulse
elseif pulse1_shape==1    
    phot=log(2)*pulse1_power/phot_energy;   %nr photons/pulse
    
end

%Calculating numerical (via Taylor series), and analytical solution (via 
%symbolic toolbox for averaged photolysed fraction
J=phot/(pi*(pulse1_size/2)^2);
N1=3*J*pulse1_cross_sec*pulse1_QY;

%First pulse
n1=1- exp(-N1*x^2); %Remember that this is a symbolic formula

%Ensemble averaged photolysed fraction
subs(n1); %Replaces all the variables in the symbolic expression n1
pulse1_Taylor_opt=str2double(get(findobj('Tag','pulse1_Taylor_opt_box'),'string'));

disp('Calculating ensemble averaged photolysed fractions numerically...Please wait...');
y_pulse1_Taylor=0.5*int(taylor(n1,pulse1_Taylor_opt),-1,1); %Works with power>250!;
set(findobj('Tag','pulse1_num_output_box'),'String',num2str(double(y_pulse1_Taylor)));

disp('Calculating ensemble averaged photolysed fractions analytically...Please wait...');
y_pulse1_anaML=0.5*int(n1,-1,1);
set(findobj('Tag','pulse1_ana_output_box'),'String',num2str(double(y_pulse1_anaML)))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plotting
if pulse1_shape==0
    single_pulse1=figure('Name','Single pulse 1 excitation results for square beam profile','NumberTitle','off','WindowStyle','docked');
elseif pulse1_shape==1
    single_pulse1=figure('Name','Single pulse 1 excitation results for Gaussian beam profile','NumberTitle','off','WindowStyle','docked');
end
color_index='krbgcmkrbgcm';

%Plot of excitation power dependence of photolysed fraction
subplot(1,3,1)
plot_power_factor_orient=str2num(get(findobj('Tag','plot_power_factor_orient_box'),'string')); %Grap factors to multiply power with
orientation_x_spacing=str2num(get(findobj('Tag','orientation_x_spacing_box'),'string')); %Grab spacing value from GUI to generate x-axis coordinates
x=-1:orientation_x_spacing:1;

for i=1:size(plot_power_factor_orient,2)
    n1_power(i,:)=1- exp(-plot_power_factor_orient(i)*N1*x.^2);
    figure(single_pulse1)
    plot(x,n1_power(i,:),color_index(i))
    hold on
end
grid on
legend(num2str(plot_power_factor_orient'))
title({'Photolysed fraction as function of ';'orientation and excitation power';'(multiples of laser power in legend)'});
xlabel('{\itx} (\mu_{laser}\bullet\mu_{molecule})');
ylabel('Fraction photolysed');

%Plot of integrated population as function of actual power density
syms x;
disp('Calculating power dependence of integrated photolysed fraction...Please wait...');
plot_power_factor_integr_start=str2double(get(findobj('Tag','plot_power_factor_integr_start_box'),'string'));
plot_power_factor_integr_spacing=str2double(get(findobj('Tag','plot_power_factor_integr_spacing_box'),'string'));
plot_power_factor_integr_end=str2double(get(findobj('Tag','plot_power_factor_integr_end_box'),'string'));
power_fac=plot_power_factor_integr_start:plot_power_factor_integr_spacing:plot_power_factor_integr_end;

power_factors_waitbar=waitbar(0,'Please wait. Calculating power factors...');
for i=1:size(power_fac,2)
    n1_power_pop(i)=0.5*int(1- exp(-power_fac(i)*N1*x.^2),-1,1);
    waitbar(i/size(power_fac,2));
end
close(power_factors_waitbar) 

subplot(1,3,2)
figure(single_pulse1)
semilogx(power_fac*pulse1_power*1e6,double(n1_power_pop))
axis tight
grid on
title({'Integrated population as function';'of actual power density'})
xlabel('Excitation power (\muJ)')
ylabel('Integrated photolysed fraction')

%Plot of ensemble averaged photolysed fraction on beam diameter dependence
plot_beam_factor_start=str2double(get(findobj('Tag','plot_beam_factor_start_box'),'string'));
plot_beam_factor_spacing=str2double(get(findobj('Tag','plot_beam_factor_spacing_box'),'string'));
plot_beam_factor_end=str2double(get(findobj('Tag','plot_beam_factor_end_box'),'string'));
beam_fac=plot_beam_factor_start:plot_beam_factor_spacing:plot_beam_factor_end;

disp('Calculating beam size dependence of ensemble averaged photolysed fraction (analytical only)...Please wait...');
diam_factors_waitbar=waitbar(0,'Please wait. Calculating diameter factors...');
for j=1:size(beam_fac,2)
    J_size_dep=phot/(pi*(beam_fac(j)*pulse1_size/2)^2);
    N1_size_dep=3*J_size_dep*pulse1_cross_sec*pulse1_QY;
    n1_size_dep=1- exp(-N1_size_dep*x^2); %Remember that this is a symbolic formula
    subs(n1_size_dep); %Ensemble averaged photolysed fraction
    waitbar(j/size(beam_fac,2));
    y_pulse1_anaML_size_dep(j)=0.5*int(n1_size_dep,-1,1);
end
close(diam_factors_waitbar)

subplot(1,3,3)
figure(single_pulse1)
semilogx(beam_fac*pulse1_size*1e4,double(y_pulse1_anaML_size_dep))
axis tight
grid on
title({'Photolysed fraction as function';'of actual beam size'});
xlabel('Beam diameter (\mum)');
ylabel('Averaged photolysed fraction');

disp('Done!')