%% Problem 1
fprintf('\nProblem 1 Results\n\n')

% load data
load('HW10_Data')

% The function that I need to find its parameters
y = @(A,tau,t)  A*exp(-t./tau);

% Guess two values for the parameters of given function
A=max(Counts);
tau=2.5;
Theta0=[A ; tau];

% Theta_guess=HW10_Solution_Exponential(Theta0,T,Counts);

% Using fminsearch to find the parameters
options= optimset();
Theta_found=fminsearch(@HW10_Solution_Exponential,Theta0,options,T,Counts);

% Found model
A_found=Theta_found(1);
tau_found=Theta_found(2);

% plot data from HW10_Data
figure(1)
plot(T,Counts,'.');  % plot the real Data
hold on
% plot the guess model
plot(T,y(A,tau,T),'g-')
% plot the found model
plot(T,y(A_found,tau_found,T),'r-')
hold off
% set the plot
legend('Data','Guess','Model')
xlabel('Time(second)')
ylabel('Counts')
title ('Decay of Muons')

% print the results
fprintf('Look at the figure to compare the results\n\n')
fprintf('The estimated parameters for exponential function are A=%d and tau=%d\n',A_found,tau_found)

%% Problem 2
fprintf('\nProblem 2 Results\n\n')
% check if the parameters exist
Prob1 = exist('A_found') && exist('tau_found');
if (~Prob1)
    fprintf('Please run problem 1 to get the prameters then do this problem\n')
end
fprintf('Wait around 1 min till this problem will be finished for Ntrail=1000\n\n')

% the given function to find the uncertainity of its parameters
y = @(A,tau,t)  A*exp(-t./tau);

% give the parameters to find their uncertainity
A=A_found;
tau=tau_found;

Options=optimset();
% Use 'waitbar' to show where the program will be finished
w=waitbar(0,'Please wait...');

% Loop stores calculated parameters each iteration in Avec and Bvec
Ntrials=1000;
for nn=1:Ntrials
    % Create random counts for exponential decay process
    y_rand=poissrnd(y(A,tau,T));
    
    % Using previously found parameters to calculate new parameters
    Theta0=[A,tau];
    [Theta_Found]=fminsearch(@HW10_Solution_Exponential,Theta0,Options,T,y_rand);
    
    % Store the parameters in appriopate vectors
    A_vec(nn)=Theta_Found(1);
    tau_vec(nn)=Theta_Found(2);
    waitbar(nn/Ntrials,w)
    
    
end
close (w)
% show the results

% Look at the hist of A_vec and tau_vec. They should be gaussian distribution
figure(2)
hist(A_vec,100)
title('Histogram for A')
figure(3)
hist(tau_vec,100)
title('Histogram for \tau')

sprintf('A = %d+/- %.7g\',mean(A_vec),std(A_vec))
sprintf('tau = %d +/- %.5g\n',mean(tau_vec),std(tau_vec))


%% Problem 3
fprintf('\nProblem 3 Results\n\n')


% You can find the wavefuncions of a quantum harmonic oscillators in 1D in
% any Quantumn textbook

% Let a=m*(omega)/(2*h-bar)
% Let y=sqrt(a)x
% Psi_0 (y) = (a/pi)^.25 * exp(-y^2/2)
% Psi_1 (y) = (a/pi)^.25 * sqrt(2)y        * exp(-y^2/2)
% Psi_2 (y) = (a/pi)^.25 * (2y^2-1)/sqrt(2)* exp(-y^2/2)

% Initilize x values
x=linspace(-5,5,500);

% Initilize parameters
m=1; w=1; h=1;
a=(m*w)/(h);

% Initilize new variable
y=sqrt(a)*x;

% Wave Functions and their probability
%  Psi_n    = ...
% (Psi_n)^2 = ...

Psi_0  = ((a/pi)^0.25)*exp(-y.^2./2);
Psi2_0 =  Psi_0.^2;

Psi_1  = ((a/pi)^0.25)*sqrt(2)*y.*exp(-y.^2./2);
Psi2_1 = Psi_1.^2;

Psi_2  = ((a/pi)^0.25)*(1/sqrt(2))*(2*y.^2-1).*exp(-y.^2./2);
Psi2_2 =  Psi_2.^2;


% Create subplots
figure
subplot(3,2,1);
plot (x,Psi_0,'r.');
title('\Psi_0');
xlabel('Position');
ylabel('Wavefunction');

subplot(3,2,2);
plot (x,Psi2_0,'r.');
title('\Psi_0^2');
xlabel('Position');
ylabel('Probability');

subplot(3,2,3);
plot (x,Psi_1,'b.');
title('\Psi_1');
xlabel('Position');
ylabel('Wavefunction');

subplot(3,2,4);
plot (x,Psi2_1,'b.');
title('\Psi_1^2');
xlabel('Position');
ylabel('Probability');

subplot(3,2,5);
plot (x,Psi_2,'g.');
title('\Psi_2');
xlabel('Position');
ylabel('Wavefunction');

subplot(3,2,6);
plot (x,Psi2_2,'g.');
title('\Psi_2^2');
xlabel('Position');
ylabel('Probability');

%% Problem 4
fprintf('\nProblem 4 Results\n\n')

mean_x=0;  % Mean of x
mean_y=0;  % Mean of y
sigma_x=2; % standard deviation in x-direction
sigma_y=5; % standard deviation in y-direction

% Create scaling factor to create a maxmium of 10
scale=10/(normpdf(0,mean_x,sigma_x).*normpdf(0,mean_y,sigma_y));

% Create meshgrid
xmax=15;
ymax=15;
[X,Y]=meshgrid(-xmax:0.1:xmax,-ymax:0.1:ymax);

% Calculate Z values
Z=scale*(normpdf(X,mean_x,sigma_x).*normpdf(Y,mean_y,sigma_y));

% Plot using surf
figure
surf(X,Y,Z)
shading interp;

% Create labels and title
xlabel('X')
ylabel('Y')
zlabel('f(X,Y)')
title (' 2D Gaussian shape')