%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Solving the stripline issue
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear, close all;
global c Q num_ch gap m volt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c = 3e8;            %speed of light
Q = -1.6e-19;       %electron charge
m = 9.1093e-31;     %elctron mass
eps = 8.85e-12;     %vacuum permittivity

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Geometry %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
length = 5e-2;      %stripline length
width = 1e-3;       %stripline width
height = 1e-3;      %stripline height
gap = 1e-3;         %gap distance

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Simulation data %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
num_ch = 1;                          %number of charges considered
t_step = 20;                        %simulation time step number
t_step = 2^nextpow2(t_step);
x_step = 5000;                       %number of point on the strip
volt = 100;                          %top electrode voltage
tf = gap*sqrt(2*m/(abs(Q)*num_ch*volt));    %time of crossing for the electrons
t0 = tf/1e4;                                %starting time for the non-linear vect
length_fft = 2^(nextpow2(t_step));          %number of point for the fft

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Processing %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
strip= linspace(-length/2,length/2,x_step);                         %position on the strip

time = 1./-exp(linspace(-log(tf),-log(t0),t_step))+tf;              %time vector non-linear (+ pts w e- closer)
time_linz = (log(-1./(time-tf))+log(tf))*(tf-t0)/log(tf/t0)+t0;     %time vector linearized
time = linspace(0,tf,t_step);
time_linz = linspace(0,tf,t_step);
time_lin = linspace(0,tf,t_step);                                   %time vector linear for fft purpose

freq_fft = t_step/(2*tf)*linspace(0,1,length_fft/2+1);              %frequency vector for fft purpose

[t,x] = meshgrid(time,strip);                                       %grid for the charge calculation
[t_plot,x_plot] = meshgrid(time_linz,strip);                        %grid for the charge display
[t_fft,x_fft] = meshgrid(time_lin, strip);                          %grid for the frequency calculation
[f_freq,x_freq] = meshgrid(freq_fft,strip);                         %grid for the frequency display

[d,v,char] = charge(t,x);                                 %compute the charge on the strip for each time and location
curr = current(t,x);                                %compute the current on the strip for each time and location
curr_t = time*num_ch^2*Q^2*volt^2/(m*gap^2);        %compute the theoretical current on the whole strip

i_curr = sum(abs(curr),1);                          %compute the integral value of the current
i_volt = char/(eps*length*width/(x_step*height));   
i_volt = sum(i_volt,1)/x_step;                      %compute the integral value of the voltage

char_freq = charge(t_fft,x_fft);                    %compute the charge on the strip for frequency analysis
mini = char_freq(:,1)';                              %compute the charge value for t=0
cont = ones(t_step,1)*mini;
char_freq = subplus(char_freq-cont');               %substract to the charge its initial value
    
char_freq = [char_freq,-char_freq(:,t_step:-1:1)];  %symmetrized charge vector
v_fft = fft(char_freq',length_fft);                 %zero-padding and fft
char_fft = 2*abs(v_fft(1:length_fft/2+1,:)')/length_fft;    %real part and normalized

pulse = real(ifft(v_fft));


char = subplus(char);
mchar  = max(max(char));
char = min(char, mchar/3);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Figures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
figure(1)
surf(t_plot,x_plot,char,'EdgeColor','none');
xlabel('time')
ylabel('position')
zlabel('charge')
colorbar

figure(2)
surf(t,x,curr/max(max(curr))+(max(max(curr))/1e6),'EdgeColor','none');
xlabel('time')
ylabel('position')
zlabel('current')
colorbar

% figure(3)
% surf(t,x,log(i_volt),'EdgeColor','none');
% xlabel('time')
% % ylabel('position')
% % zlabel('volt')
% % colorbar
% 
figure(4)
hold on;
xlabel('time')
ylabel('integral of current')
plot(time,i_curr);
%plot(time,curr_t,'r');
hold off;

figure(5)
xlabel('time')
ylabel('average voltage on the strip')
plot(time_linz,i_volt);

figure(6)
surf(f_freq,x_freq,char_fft,'EdgeColor','none');
xlabel('frequencies')
ylabel('position')
zlabel('amplitude')
colorbar

figure(7)
surf(t_fft,x_fft,pulse,'EdgeColor', 'none');

figure(8)
semilogx(freq_fft,char_fft(50,:));
grid on
xlabel('Frequencies (Hz)')
ylabel('Amplitude')

figure(11)
subplot(2,1,1)
plot (t(1,:),d(1,:));
xlabel('time (s)');
ylabel('distance (m)');
subplot(2,1,2);
plot(t(1,:),v(1,:));
xlabel('time (s)');
ylabel('velocity (m/s)');