%Script Testing Case Study for CIP, 17/08/2017, Matthias Sinnesael

clear all
close all
clc
tic
% 1) Load Insolation Curve

data = load('Laskar_Insolation_8_6Ma_55N_21July_meandaily.txt');
time=data(:,1);
Insolation=data(:,2);

figure();
plot(time,Insolation);
title('Insolation 55°N, 21 July');
xlabel('ka');ylabel('W/m2');

% 2) Introduce threshold & nonlinearity, e.g. (X-440)² W m-2 for this record
% 2a) Threshold

Insolation_Threshold = Insolation-440;
figure();
plot(time,Insolation_Threshold);
title('Insolation 55°N, 21 July - Threshold');
xlabel('ka');ylabel('W/m2');

% 2b) Non-lin
length = length(Insolation_Threshold);
for i = 1 : length
Insolation_nonlin(i) = Insolation_Threshold(i)^2;
end
Insolation_nonlin=Insolation_nonlin'
figure();
plot(time,Insolation_nonlin);
title('Insolation 55°N, 21 July - Threshold + Nonlinear (squared)');
xlabel('ka');ylabel('W/m2');

% 3) "Create Lithology using thresholds"
%length2 = length(Insolation_nonlin);
for j = 1 : length
    if Insolation_nonlin(j)<500
        Lithology(j)=1;
    elseif Insolation_nonlin(j)<1000
        Lithology(j)=2;
    elseif Insolation_nonlin(j)<1500
        Lithology(j)=3;
    elseif Insolation_nonlin(j)<2100
        Lithology(j)=4;
    else
        Lithology(j)=5;
    end
end
Lithology=Lithology';

figure();
plot(time,Insolation/150);
hold on
plot(time,Lithology,'Color','gr');
title('Lithology & Insolation 55°N, 21 July');
xlabel('ka');% ylabel('W/m2');
ylim([0 6])

% 4) Transition from the time to the distance domain
%starting at 0m (-8 Ma) and havind a constant sed. of 10 m/Myr

distance = time + 8000; 

figure();
plot(distance,Lithology)
xlabel('cm');% ylabel('W/m2');
ylim([0 6])

%%%---------------------------------------------------------------------%%%

% 5) Thickness distortion

%Create new matrix (table) for rescalling, number of rows
rows=50000;

%make rescaled empty
data_rescaled=zeros(rows,2);

%determine sampling_interval
sampling_interval=data(2,1)-data(1,1);

%fill up new array with depths
for j=1:rows
    data_rescaled(j,1)=j*sampling_interval;
end

%save new distance scale
new_distance_scale=data_rescaled(:,1);

%make number of rows an input variable to declare
[m,n]=size(data); n_rows=m; n_col=n;

%scaling thicknesses
scaling_L1 = 1; 
scaling_L2 = 2;
scaling_L3 = 3;
scaling_L4 = 4; 
scaling_L5 = 5;


%CALCULATIONS%
%__________________________________________________________________________


for k = 1 : n_rows
    if Lithology(k) == 1
        next=min(find(not(data_rescaled)))-rows;
        data_rescaled(next:next+(scaling_L1-1),2)=1;
    end
        if Lithology(k) == 2
        next=min(find(not(data_rescaled)))-rows;
        data_rescaled(next:next+(round(rand(1))*scaling_L2-1),2)=2;
        end
    if Lithology(k) == 3
        next=min(find(not(data_rescaled)))-rows;
        data_rescaled(next:next+(round(rand(1))*scaling_L3-1),2)=3;
    end
        if Lithology(k) == 4
        next=min(find(not(data_rescaled)))-rows;
        data_rescaled(next:next+(round(rand(1))*scaling_L4-1),2)=4;
        end
    if Lithology(k) == 5
        next=min(find(not(data_rescaled)))-rows;
        data_rescaled(next:next+(round(rand(1))*scaling_L5-1),2)=5;
    end
end

%OUTPUT%
%__________________________________________________________________________

figure()
plot(new_distance_scale,data_rescaled(:,2))
title('Distorted Lithology');
xlim([0 next]);ylim([0 6]);
xlabel('cm'); ylabel('Lithology');

%%%---------------------------------------------------------------------%%%

% 6) Filling up with color according to lithology

for k = 1 : next
    if data_rescaled(k,2) == 1 
        colorplot_L1(k)=1;
    else
        colorplot_L1(k)=0;
    end
    if data_rescaled(k,2) == 2 
        colorplot_L2(k)=1;
    else
        colorplot_L2(k)=0;
    end
    if data_rescaled(k,2) == 3 
        colorplot_L3(k)=1;
    else
        colorplot_L3(k)=0;
    end
    if data_rescaled(k,2) == 4 
        colorplot_L4(k)=1;
    else
        colorplot_L4(k)=0;
    end
    if data_rescaled(k,2) == 5 
        colorplot_L5(k)=1;
    else
        colorplot_L5(k)=0;
    end
end

% figure()
% area(colorplot_L1,'FaceColor',[1 0 0],'LineStyle','none');%red
% hold on
% area(colorplot_L2,'FaceColor',[0 1 0],'LineStyle','none');%green
% hold on
% area(colorplot_L3,'FaceColor',[0 0 1],'LineStyle','none');%blue

figure()
area(colorplot_L1,'FaceColor',[0.6 0.6 0.6],'LineStyle','none');%red
hold on
area(colorplot_L2,'FaceColor',[0.5 0.5 0.5],'LineStyle','none');%green
hold on
area(colorplot_L3,'FaceColor',[0.4 0.4 0.4],'LineStyle','none');%blue
hold on
area(colorplot_L4,'FaceColor',[0.3 0.3 0.3],'LineStyle','none');%green
hold on
area(colorplot_L5,'FaceColor',[0.2 0.2 0.2],'LineStyle','none');%blue


%xlswrite('L1',colorplot_L1')

toc