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

function two_pulse_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.

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

%Get pulse 1 parameters from gui
pulse1_wave=1e-9*str2double(get(findobj('Tag','pulse1_wave_box'),'string'));
pulse1_power=1e-6*str2double(get(findobj('Tag','pulse1_power_box'),'string'));
pulse1_size=1e-4*str2double(get(findobj('Tag','pulse1_size_box'),'string'));
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'));

%Pulse 2 parameters
pulse2_wave=1e-9*str2double(get(findobj('Tag','pulse2_wave_box'),'string'));
pulse2_power=1e-6*str2double(get(findobj('Tag','pulse2_power_box'),'string'));
pulse2_size=1e-4*str2double(get(findobj('Tag','pulse2_size_box'),'string'));
pulse2_shape=get(findobj('Tag','pulse2_shape_box'),'Value');
pulse2_QY=str2double(get(findobj('Tag','pulse2_QY_box'),'string'));
pulse2_cross_sec=str2double(get(findobj('Tag','pulse2_cross_sec_box'),'string'));

%Get two-pulse parameters
pulses_delay=str2num(get(findobj('Tag','pulses_delay_box'),'string'));
pulses_recovery=str2num(get(findobj('Tag','pulses_recovery_box'),'string'));

%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
if max(isnan([pulse2_wave pulse2_power pulse2_size pulse2_QY pulse2_cross_sec]))==1
    errordlg('Please fill in all boxes to define Pulse 2.')
    error('Please fill in all boxes to define Pulse 2.')
end
if length(pulses_delay)~=length(pulses_recovery)    %Consistency check
    errordlg('Number of time delays does not match with number of ground state recovery fractions!')
    error('Number of time delays does not match with number of ground state recovery fractions!')
end

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

phot_energy2=6.626068e-34*2.99792e8/(pulse2_wave);   %Planck, light and excitation wavelength
if pulse2_shape==0
    phot2=pulse2_power/phot_energy2;   %nr photons/pulse
elseif pulse2_shape==1
    phot2=log(2)*pulse2_power/phot_energy2;   %nr photons/pulse
end

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

J2=phot2/(pi*(pulse2_size/2)^2);
N2=3*J2*pulse2_cross_sec*pulse2_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_2pulses_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_2pulses_ana_output_box'),'String',num2str(double(y_pulse1_anaML)))

%--------------
%Second pulse
n2=1- exp(-N2*x^2);
subs(n2);
pulse2_Taylor_opt=str2double(get(findobj('Tag','pulse2_Taylor_opt_box'),'string'));

for i=1:size(pulses_recovery,2)
    n2_t(i)=(1- (1-pulses_recovery(i))*n1 ) * n2;
end

for i=1:size(pulses_recovery,2)
    y_pulse2_Taylor(i)=0.5*int(taylor(n2_t(i),pulse2_Taylor_opt),-1,1); %Works with power>250!; power 31 for long delays OK; short delays: 31 OK for 2 decimal, 59 for 5decimals, inf for correct sum
end
set(findobj('Tag','pulse2_2pulses_num_output_box'),'String',num2str(double(y_pulse2_Taylor)));

for i=1:size(pulses_recovery,2)
    y_pulse2_ana(i)=0.5*int(n2_t(i),-1,1);
end
set(findobj('Tag','pulse2_2pulses_ana_output_box'),'String',num2str(double(y_pulse2_ana)));
%############################


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Plotting
if pulse1_shape==0 && pulse2_shape==0
    two_pulse=figure('Name','Two pulse results with ground state recovery, square beam profiles','NumberTitle','off','WindowStyle','docked');
elseif pulse1_shape==1 && pulse2_shape==0
    two_pulse=figure('Name','Two pulse results with ground state recovery, Gaussian beam pulse 1, square pulse 2','NumberTitle','off','WindowStyle','docked');
elseif pulse1_shape==0 && pulse2_shape==1
    two_pulse=figure('Name','Two pulse results with ground state recovery, square pulse 1, Gaussian beam pulse 2','NumberTitle','off','WindowStyle','docked');
elseif pulse1_shape==1 && pulse2_shape==1
    two_pulse=figure('Name','Two pulse results with ground state recovery, Gaussian beam profiles','NumberTitle','off','WindowStyle','docked');
end

color_index='krbgcmkrbgcm';

%Plot of excitation power dependence of photolysed fraction
subplot(2,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;

recovery_populations_waitbar=waitbar(0,'Please wait. Calculating power dependence of orientational distrubution...');
for i=1:size(plot_power_factor_orient,2)
    waitbar(i/size(plot_power_factor_orient,2));
    n1_power(i,:)=1- exp(-plot_power_factor_orient(i)*N1*x.^2);
    figure(two_pulse)
    plot(x,n1_power(i,:),color_index(i))
    hold on
end
close (recovery_populations_waitbar)
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 n_1');


%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(2,3,2)
figure(two_pulse)
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 <n_1>')

%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=phot1/(pi*(beam_fac(j)*pulse1_size/2)^2);
    N0_size_dep=3*J_size_dep*pulse1_cross_sec*pulse1_QY;
    n1_size_dep=1- exp(-N0_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(2,3,3)
figure(two_pulse)
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('<n_1>');

%Plot of recovery of the orientation function of the ensemble averaged 
%photolysed fraction as function of user-given ground state recovery 
%fractions at given excitation power of pulse 1
subplot(2,3,4)
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; %Get x again from gui as has been redefined as symbolic before

for i=1:size(pulses_recovery,2)
    n1_rec(i,:)=subs(n1)*(pulses_recovery(i)); %n1 is calculated for excitation power
    figure(two_pulse)
    plot(x,n1_rec(i,:),color_index(i))
    hold on
end
grid on
axis tight
legend(num2str(pulses_recovery'))
title({'Recovery of the orientational distribution of different fractions';'excited by pulse 1 for different ground state recovery fractions';['(Power = ',num2str(1e6*pulse1_power),' \muJ )']});
xlabel('{\itx} (\mu_{laser}\bullet\mu_{molecule})');
ylabel('$(1-\beta) \times n_1$','Interpreter','latex');


%Plot of recovering orientation distribution of pulse 1 x  orientation
%distribution of pulse 2 at given excitation power of pulse 1
subplot(2,3,5)
for i=1:size(pulses_recovery,2)
    figure(two_pulse)
    plot(x,n1_rec(i,:).*subs(n2),color_index(i))
    hold on
end
grid on
axis tight
title({'Orientation distribution after';' interaction with pulse 2';['(Power pulse 1 = ',num2str(1e6*pulse1_power),' \muJ, Power pulse 2 = ',num2str(1e6*pulse2_power),' \muJ)']})
xlabel('{\itx} (\mu_{laser}\bullet\mu_{molecule})');
ylabel('$(1-\beta) \times n_1 \times n_2$','Interpreter','latex');

%Plot of recovering populations as function of given times (as set by user
%in Delay between pulses panel)
syms x
recovery_populations_waitbar=waitbar(0,'Please wait. Calculating recovery populations...');
for i=1:size(pulses_recovery,2)
    waitbar(i/size(pulses_recovery,2));
    n2_rec_pop(i)=0.5*int((1- (1-pulses_recovery(i)) *n1 ) * n2 ,-1,1);
end
close(recovery_populations_waitbar)
subplot(2,3,6)
semilogx(pulses_delay,double(n2_rec_pop));
hold on;
semilogx(pulses_delay,(1-pulses_recovery)*(subs(y_pulse1_anaML)),'k');
axis tight
grid on
axis tight
legend('<n_2>','\beta \times <n_1>','location','NE')
title({'Time course of ground state recovered';'photolysed fraction of pulse 1';'and of pulse 2 as function of delay time';['(Power pulse 1 = ',num2str(1e6*pulse1_power),' \muJ, Power pulse 2 = ',num2str(1e6*pulse2_power),' \muJ)']})
axis tight
xlabel ('Time')
ylabel('$<n>$','Interpreter','latex');

%############################
%Display in gui
set(findobj('Tag','n_2_time_course_output_box'),'String',num2str(double(n2_rec_pop)));
set(findobj('Tag','n_1_beta_time_course_output_box'),'String',num2str(double((1-pulses_recovery)*(subs(y_pulse1_anaML)))));
%############################

disp('Done!')