%% MATLAB Recipes for Earth Sciences - Chapter 10

%% Section 10.2
clear

%%
data_degrees_1 = load('directional_1.txt');

%%
data_radians_1 = pi*data_degrees_1/180;
rose(data_radians_1,12)

%%
rose(data_radians_1,12)
view(90,-90)

%%
rose_sqrt(data_radians_1,12)
view(90,-90)

%% Section 10.3
clear

%%
data_degrees_1 = load('directional_1.txt');
data_radians_1 = pi*data_degrees_1/180;

%%
x_1 = sum(sin(data_radians_1))
y_1 = sum(cos(data_radians_1))

%%
mean_radians_1 = atan(x_1/y_1)
mean_degrees_1 = 180*mean_radians_1/pi

%%
mean_degrees_1 = mean_degrees_1 + 180

%%
R_1 = sqrt(x_1^2 + y_1^2)

%%
Rm_1 = R_1 / (length(data_radians_1))

%%
sigma_1 = 1 - Rm_1

%% Section 10.4
clear

%%
mu = 0; kappa = [0 1 2 3 4]';

%%
theta = -180:1:180;

%%
mu_radians = pi*mu/180;
theta_radians = pi*theta/180;

%%
for i = 1:5 
 mises(i,:) = (1/(2*pi*besseli(0,kappa(i))))* ...
 exp(kappa(i)*cos(theta_radians-mu_radians));
 theta(i,:) = theta(1,:);
end

%%
for i = 1:5
 plot(theta(i,:),mises(i,:))
 axis([-180 180 0 max(mises(i,:))])
 hold on
end

%% Section 10.5
clear

%%
data_degrees_1 = load('directional_1.txt');

%%
counts = hist(data_degrees_1,15:30:345);

%%
expect = 40/12 * ones(1,12);

%%
chi2 = sum((counts - expect).^2 ./expect)

%%
chi2inv(0.95,12-1-1)

%% Section 10.6
clear

%%
data_degrees_1 = load('directional_1.txt');
data_radians_1 = pi*data_degrees_1/180;

%%
x_1 = sum(sin(data_radians_1));
y_1 = sum(cos(data_radians_1));

%%
mean_radians_1 = atan(x_1/y_1);
mean_degrees_1 = 180*mean_radians_1/pi;
mean_degrees_1 = mean_degrees_1 + 180;

%%
Rm_1 = 1/length(data_degrees_1) .*(x_1.^2+y_1.^2).^0.5

%%
theta_1 = 180 * atan(x_1/y_1) / pi

%% Section 10.7
clear

%%
data_degrees_1 = load('directional_1.txt');
data_degrees_2 = load('directional_2.txt');

%%
data_radians_1 = pi*data_degrees_1/180;
data_radians_2 = pi*data_degrees_2/180;

%%
x_1 = sum(sin(data_radians_1));
y_1 = sum(cos(data_radians_1));
x_2 = sum(sin(data_radians_2));
y_2 = sum(cos(data_radians_2));

%%
mean_radians_1 = atan(x_1/y_1);
mean_degrees_1 = 180*mean_radians_1/pi;
mean_radians_2 = atan(x_2/y_2);
mean_degrees_2 = 180*mean_radians_2/pi;

%%
mean_degrees_1 = mean_degrees_1 + 180
mean_degrees_2 = mean_degrees_2 + 180

%%
R_1 = sqrt(x_1^2 + y_1^2);
R_2 = sqrt(x_2^2 + y_2^2);

%%
data_radians_T = [data_radians_1;data_radians_2];

%%
x_T = sum(sin(data_radians_T));
y_T = sum(cos(data_radians_T));

%%
mean_radians_T = atan(x_T/y_T);
mean_degrees_T = 180*mean_radians_T/pi;

%%
mean_degrees_T = mean_degrees_T + 180;

%%
R_T = sqrt(x_T^2 + y_T^2)
Rm_T = R_T / (length(data_radians_T))

%%
n = length(data_radians_T);

%%
F = (1+3/(8*4.177)) * (((n-2)*(R_1+R_2-R_T))/(n-R_1-R_2))

%%
finv(0.95,1,78)



