%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%% MAIN %%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc; clear all; close all;
global Box Network Fiber Particle Water_shed Character

arch_black = [0,0,0];
bulldog_red = [186, 12, 47]./255;
odyssey_grey = [200, 216, 235]/275;
olympic_blue = [0, 78, 96]/255;
creamery_beige = [214, 210, 196]/255;
lake_herrick = [0, 163, 173]/255;
glory_glory = [228, 0, 43]/255;
hedges = [183, 191, 16]/255;
mat_blue = [0, 114, 189]/255;

% Units ng, um, us 
% (E=ng*um^2/us^2 == 1 nJ)
% (mass_density=ng/um^3 == 1e9 kg/m^3)
% (P=ng/(um*us^2) == 1 GPa)

[in_num, in_text] = xlsread('input.xlsx');

% Load packages
disp("Importing dependencies"); 

ImportDependencies()

% Generate objects
disp("Generating network objects")

Box         = box();
Network     = network();
Fiber       = fiber();
Particle    = particle();
Water_shed  = water_shed();
Character   = character();

% Set up environment

Box.set_box_size(in_num(1), in_num(2), in_num(3));

Network.set_network_properties(in_num(4), in_num(5), in_num(6), in_num(7), in_num(8), in_num(9), in_num(10), in_num(11), in_text{12,2}, in_num(13), in_text{14,2}, in_num(15));

Fiber.set_fiber_properties(in_num(16), in_num(17), in_num(18));

Particle.set_particle_properties();

Character.get_initial_network_character(Box, Network, Fiber, Particle);

plot_network([89  179  173]./255, [89  179  173]./255, [89  179  173]./255, 1, [-35,15])
print(gcf,'Initial_network_configuration','-dpng','-r300');

natural_frequency = (1/2)*sqrt(max([max(Particle.k_bond)./max(Particle.particle_mass);
                                    min(Particle.k_bond)./min(Particle.particle_mass)]));                             
alpha   = 64;
beta    = 32;
delta_t = 1/(alpha*natural_frequency);
delta_E = min((Particle.k_bond.*Particle.bond_length.^2))/beta^2;
epsilon = Particle.N_particle*(delta_E/delta_t);

time    = 0;
lambda  = 2*sqrt(pi*Particle.particle_mass.*Particle.particle_diameter/6);
v       = zeros(Particle.N_particle,3);
time_step = 0;
comp_time = 0;
time_evolution = [];
time_evolution_E = [];

delta_E = epsilon*10;
E_last  = 0;
E_soft  = 2*epsilon*delta_t;

disp(['sim_time ', 'E_total ', 'E_soft' , 'E_bond ', 'E_angle ','E_kinetic ', 'delta_E' , 'comp_time'])

while delta_E > epsilon || sum(E_soft) > epsilon*delta_t
    
    tic()
    [E_soft, F_soft]    = soft_potential(Particle.N_particle, Particle.x_particle, Particle.particle_diameter, Particle.k_soft, Particle.particle_N_bond, Particle.particle_bond_neigh, Box.box_bounds);
    [E_bond, F_bond]    = harmonic_bond_potential(Particle.N_particle, Particle.N_bond, Particle.x_particle, Particle.bond_length, Particle.k_bond, Particle.bond ,Box.box_bounds);
    [E_angle, F_angle]  = harmonic_angle_potential(Particle.N_particle, Particle.N_angle, Particle.x_particle, Particle.cos_theta_0, Particle.k_angle, Particle.angle, Box.box_bounds);

    E_kinetic   = (1/2).*Particle.particle_mass.*(v(:,1).^2 + v(:,2).^2 + v(:,3).^2);
    speed       = sqrt(v(:,1).^2 + v(:,3).^2 + v(:,2).^2);
    if sum(isnan(speed)) == 0 && sum(speed==0) == 0 
        direction   = v./speed;
        F_drag = -lambda.*direction.*speed.^2;
    else
        F_drag = -lambda.*v;
    end
    
    F_drag = -lambda.*v;
    
    E_total     = E_soft + E_bond + E_angle + E_kinetic;
    F_total     = F_soft + F_bond + F_angle + F_drag;
    
    a = F_total./Particle.particle_mass;
    v = v + a.*delta_t;
    Particle.x_particle = Particle.x_particle + v*delta_t;
    Particle.x_particle_unwrapped = Particle.x_particle_unwrapped + v*delta_t;
    
    Particle.wrap_particles(Box.box_bounds);
    
    comp_time = comp_time + toc();
    time_evolution  = [time_evolution; time_step];
    time_evolution_E = [time_evolution_E; sum(E_total)];
    
    if mod(time_step,10) == 0 || time_step == 0
        disp([delta_t*time_step, sum(E_total), sum(E_soft), sum(E_bond), sum(E_angle), sum(E_kinetic), delta_E, comp_time/60])
        figure(2)
        cla
        scatter(time_evolution, time_evolution_E, 50, 'k', 'Filled') 
        if time_step > 0
            axis([0 max([max(time_evolution),1]) min(time_evolution_E) max(time_evolution_E)])
        end
        xlabel('time step')
        ylabel('Energy (nJ)')
        set(gca,'FontName','Times New Roman','FontSize',22,'fontweight','bold','Box','on')
        set(gcf,'Color','w');
        print(gcf,'Energy_minimization','-dpng','-r300');
        pause(1e-6)
    end
    
    time = time + delta_t;
    time_step = time_step + 1;
    
    delta_E = abs(sum(E_total) - E_last);
    
    E_last = sum(E_total);
    
    Particle.k_soft = Particle.k_soft*1.005;
end

%Character.get_final_network_character(Box, Network, Fiber, Particle);

%plot_network([89  179  173]./255, [89  179  173]./255, [89  179  173]./255, 3, [-35,15])
%print(gcf,'Final_network_configuration','-dpng','-r300');

length_experimental = csvread('xiaole_chen_fiber_length_distribution.csv');
LData = csvread('xiaole_chen_fiber_length_distribution.csv');
LData(:,2) = LData(:,2)*(Network.N_fiber/simps(LData(:,1),LData(:,2)));

a = 4.97;
b = 0.212;

length = LData(:,1);
length_frequency = LData(:,2);

x = [0:0.01:100];
y = b^a*x.^(a-1).*exp(-b*x);
y_sum = simps(x,y);
y = y*(Network.N_fiber/y_sum);

[y_cdf_l1,x_cdf_l1] = ecdf(Character.init_fiber_length);
[y_cdf_l2,x_cdf_l2] = ecdf(Character.final_fiber_length);

figure(4)
hold on
yyaxis left
scatter(LData(:,1),LData(:,2)*7,100,odyssey_grey,'Filled')
plot(x,y*6.85,'LineWidth',2,'Color',bulldog_red)
histogram(Character.init_fiber_length, 'BinWidth', 2, 'FaceColor', 'none', 'EdgeColor', arch_black, 'Normalization','countdensity')
histogram(Character.final_fiber_length, 'BinWidth', 2, 'FaceColor', 'none', 'EdgeColor', glory_glory, 'Normalization','countdensity')
xlabel('Fiber Length (\mum)')
ylabel('Count Density')
set(gca,'FontName','Times New Roman','FontSize',28,'fontweight','bold','Box','on','ycolor','k')
axis([0 90 0 28])
yyaxis right
plot(x_cdf_l1,y_cdf_l1,'LineWidth', 1.5, 'Color', 'k')
plot(x_cdf_l2,y_cdf_l2,'LineWidth', 1.5, 'Color', 'k')
ylabel('CDF')
set(gcf,'Color','w');
legend('Exp','Fit','Initial','Final','Initial','Final')
set(gca,'FontName','Times New Roman','FontSize',28,'fontweight','bold','Box','on','ycolor','k')
print(gcf,'Length_distribution','-dpng','-r300');

ml1 = mean(Character.init_fiber_length);
sl1 = std(Character.init_fiber_length);

ml2 = mean(Character.final_fiber_length);
sl2 = std(Character.final_fiber_length);

[hl,pl] = kstest2(Character.init_fiber_length, Character.final_fiber_length);

tortuosity_experimental = xlsread('LogNormalDist3.xlsx');
experimental_tau = tortuosity_experimental(:,1);
experimental_frequency = tortuosity_experimental(:,2);
experimental_frequency = experimental_frequency*(Network.N_fiber/simps(experimental_tau,experimental_frequency));
experimental_frequency = experimental_frequency + 0*randn(size(experimental_frequency,1),size(experimental_frequency,2));

mu = 1.18;
sigma = 0.058;

m  = log(mu/sqrt(1 + sigma^2/mu^2));
s  = sqrt(log(1 + sigma^2/mu^2));

x = [1:1/1000:3.5];
y = (1./x).*(1/(s*sqrt(2*pi))).*exp(-((log(x)-m).^2)/(2*s^2));
y_sum = simps(x,y);
y = y*(Network.N_fiber/y_sum);

[y_cdf_tau1,x_cdf_tau1] = ecdf(Character.init_fiber_tortuosity);
[y_cdf_tau2,x_cdf_tau2] = ecdf(Character.final_fiber_tortuosity);
[y_cdf_texp,x_cdf_texp] = ecdf(experimental_tau);

figure(5)
hold on
yyaxis left
scatter(experimental_tau,experimental_frequency*4.5,100,odyssey_grey,'Filled')
plot(x,y*4.5,'LineWidth',2,'Color',bulldog_red)
histogram(Character.init_fiber_tortuosity, 'BinWidth', 0.02, 'FaceColor', 'none', 'EdgeColor', arch_black, 'Normalization','countdensity')
histogram(Character.final_fiber_tortuosity, 'BinWidth', 0.02, 'FaceColor', 'none', 'EdgeColor', glory_glory, 'Normalization','countdensity')
axis([1 1.6 0 3000])
xlabel('Fiber Tortuosity')
ylabel('Count Density')
set(gca,'FontName','Times New Roman','FontSize',28,'fontweight','bold','Box','on','ycolor','k')
yyaxis right
plot(x_cdf_tau1,y_cdf_tau1,'LineWidth', 1.5, 'Color', 'k')
plot(x_cdf_tau2,y_cdf_tau2,'LineWidth', 1.5, 'Color', 'k')
set(gca,'FontName','Times New Roman','FontSize',28,'fontweight','bold','Box','on','ycolor','k')
set(gcf,'Color','w');
ylabel('CDF')
axis([1 1.6 0 1])
print(gcf,'Tortuosity_distribution','-dpng','-r300');

mtau1 = mean(Character.init_fiber_tortuosity);
stau1 = std(Character.init_fiber_tortuosity);

mtau2 = mean(Character.final_fiber_tortuosity);
stau2 = std(Character.final_fiber_tortuosity);

[htau,ptau] = kstest2(Character.init_fiber_tortuosity, Character.final_fiber_tortuosity);

[y_cdf_cl1,x_cdf_cl1] = ecdf(Character.init_fiber_crosslink_distance);
[y_cdf_cl2,x_cdf_cl2] = ecdf(Character.final_fiber_crosslink_distance/2);

figure(6)
hold on
yyaxis left
histogram(Character.init_fiber_crosslink_distance, 'BinWidth', 0.12, 'FaceColor', 'none', 'EdgeColor', arch_black, 'Normalization','countdensity')
histogram(Character.final_fiber_crosslink_distance/2, 'BinWidth', 0.12, 'FaceColor', 'none', 'EdgeColor', glory_glory, 'Normalization','countdensity')
xlabel('Crosslink Distance (\mum)')
ylabel('Count Density')
set(gca,'FontName','Times New Roman','FontSize',28,'fontweight','bold','Box','on','ycolor','k')
yyaxis right
plot(x_cdf_cl1,y_cdf_cl1,'LineWidth', 1.5, 'Color', 'k')
plot(x_cdf_cl2,y_cdf_cl2,'LineWidth', 1.5, 'Color', 'k')
set(gca,'FontName','Times New Roman','FontSize',28,'fontweight','bold','Box','on','ycolor','k')
set(gcf,'Color','w');
ylabel('CDF')
axis tight
set(gcf,'Color','w');
print(gcf,'Crosslink_distance_distribution','-dpng','-r300');

mcl1 = mean(Character.init_fiber_crosslink_distance);
scl1 = std(Character.init_fiber_crosslink_distance);

mcl2 = mean(Character.final_fiber_crosslink_distance/2);
scl2 = std(Character.final_fiber_crosslink_distance/2);

orientation_experimental = xlsread('O1.xlsx');
experimental_theta = orientation_experimental(:,1);
experimental_frequency = orientation_experimental(:,2);
experimental_frequency = experimental_frequency*(Network.N_fiber/simps(experimental_theta,experimental_frequency));

k = 2;
x = [0:1/100:180];
y = ((2*k+1)/(8*pi))*(cosd(x).^(2*k) + sind(x).^(2*k));
y_sum = simps(x,y);
y = y*(Network.N_fiber/y_sum);

[y_cdf_o1,x_cdf_o1] = ecdf(Character.init_fiber_orientation_theta);
[y_cdf_o2,x_cdf_o2] = ecdf(Character.final_fiber_orientation_theta);

figure(7)
hold on
yyaxis left
scatter(orientation_experimental(:,1),orientation_experimental(:,2)/50,100,odyssey_grey,'Filled')
plot(x,y*6,'LineWidth',2,'Color',bulldog_red)
histogram(Character.init_fiber_orientation_theta, 'BinWidth', 5, 'FaceColor', 'none', 'EdgeColor', arch_black, 'Normalization','countdensity')
histogram(Character.final_fiber_orientation_theta, 'BinWidth', 5, 'FaceColor', 'none', 'EdgeColor', glory_glory, 'Normalization','countdensity')
xlabel('Theta (degrees)')
ylabel('Count Density')
set(gca,'FontName','Times New Roman','FontSize',28,'fontweight','bold','Box','on','ycolor','k')
yyaxis right
plot(x_cdf_o1,y_cdf_o1,'LineWidth', 1.5, 'Color', 'k')
plot(x_cdf_o2,y_cdf_o2,'LineWidth', 1.5, 'Color', 'k')
set(gca,'FontName','Times New Roman','FontSize',28,'fontweight','bold','Box','on','ycolor','k')
set(gcf,'Color','w');
ylabel('CDF')
axis tight
print(gcf,'Orientation_distribution','-dpng','-r300');

[htau,ptau] = kstest2(Character.init_fiber_orientation_theta, Character.final_fiber_orientation_theta);