clc
clear all
close all

tic

cd 'C:\Users\odin_\OneDrive\Documents\HVL- Byggenginør\6. semester\Bacheloroppgave\Målinger kloridsensor\cvs reader'

dager = 54;
molinger_per_sekund = 2;
tid_per_moling = 720;
molinger_totalt = molinger_per_sekund * tid_per_moling + 1;
gjeldende_molinger = 240;
betong = ["M90","M60"];
salt = ["0%","1,75%","3,5%","5,25%","7%"];
kube = ["nr.1","nr.2","nr.3"];
kuber_per_type = length(salt) * length(kube);
M90_M60 = zeros(dager,kuber_per_type,length(betong));
M90_M60_std = zeros(dager,kuber_per_type,length(betong));
defekt = [10,23]; %[10,23] (5,25_1,3,5_2)

for dag = 1:dager
    data = les_fil(dag);
    if not(isempty(data))
        if rem(dag,2) == 0
            M90_M60(dag,:,1) = mean(data(molinger_totalt - gjeldende_molinger:end,:));
            M90_M60_std(dag,:,1) = std(data(molinger_totalt - gjeldende_molinger:end,:));
        else
            M90_M60(dag,:,2) = mean(data(molinger_totalt - gjeldende_molinger:end,:));
            M90_M60_std(dag,:,2) = std(data(molinger_totalt - gjeldende_molinger:end,:));
        end
    end
end

if not(defekt ==0)
    for i = 1:length(defekt)
        if defekt(i)<=15
            M90_M60(:,defekt(i),1) = 0;
        else
            M90_M60(:,defekt(i)-15,2) = 0;
        end
    end
end

minste_moling = min(M90_M60(find(M90_M60)));
storste_moling = max(M90_M60(find(M90_M60)));
storste_std = max(M90_M60_std(find(M90_M60_std)))/sqrt(gjeldende_molinger);
def = 0;


for i = 1:length(betong)
    figure('Name',betong(i),'WindowState','maximized')
    figure('Name',betong(i)+'gjennomsnitt','WindowState','maximized')
    figure(1+(i-1)*2)
    [x,y] = find(M90_M60(:,:,i));
    X = unique(x).';
    n = 0;
    for sal = 1:length(salt)
        if sal == 3
            subplot(3,2,[3])
        elseif sal > 3
            subplot(3,2,sal+1)
        else
            subplot(3,2,sal)
        end
        hold on
        title(betong(i)+'-'+salt(sal))
        if i == 1
            xlim([0.75,dager+0.25-5])
        else
            xlim([0.75,dager+0.25])
        end
        ylim([minste_moling - storste_std,storste_moling + storste_std])
        xticks(1:2:dager)
        grid on
        yticks(0:1:round(storste_moling + storste_std))
        fjern = [false,0];
        for kub = 1:length(kube)
            n = n + 1;
            def = def + 1;
            if not(ismember(def,defekt))
                if i == 1
                    errorbar(X-5,M90_M60(X,n,i),M90_M60_std(X,n,i)/sqrt(gjeldende_molinger))
                else
                    errorbar(X,M90_M60(X,n,i),M90_M60_std(X,n,i)/sqrt(gjeldende_molinger))
                end
            else
                fjern = [true,rem(n,3)];
            end
        end
        if fjern(1)
            legend(kube([1:fjern(2)-1,fjern(2)+1:3]),'Location','northwest')
        else
            legend(kube,'Location','northwest')
        end
        figure(i*2)
        if sal == 3
            subplot(3,2,[3,4])
        elseif sal > 3
            subplot(3,2,sal+1)
        else
            subplot(3,2,sal)
        end
        if i == 1
            xlim([0.75,dager+0.25-5])
        else
            xlim([0.75,dager+0.25])
        end
        ylim([minste_moling - storste_std,storste_moling + storste_std])
        xticks(1:2:dager)
        yticks(0:1:round(storste_moling + storste_std))
        grid on
        title(betong(i)+'-'+salt(sal))
        hold on
        if fjern(1)
            if i == 1
                errorbar(X-5,mean(M90_M60(X,[1:fjern(2)-1,fjern(2)+1:3]+(n-3),i),2),std(M90_M60(X,[1:fjern(2)-1,fjern(2)+1:3]+(n-3),i),0,2)/sqrt(2))
            else
                errorbar(X,mean(M90_M60(X,[1:fjern(2)-1,fjern(2)+1:3]+(n-3),i),2),std(M90_M60(X,[1:fjern(2)-1,fjern(2)+1:3]+(n-3),i),0,2)/sqrt(2))
            end
            b(sal) = linreg(X,M90_M60(:,[1:fjern(2)-1,fjern(2)+1:3]+(n-3),:),i);
        else
            if i == 1
                errorbar(X-5,mean(M90_M60(X,n-2:n,i),2),std(M90_M60(X,n-2:n,i),0,2)/sqrt(3))
            else
                errorbar(X,mean(M90_M60(X,n-2:n,i),2),std(M90_M60(X,n-2:n,i),0,2)/sqrt(3))
            end
            b(sal) = linreg(X,M90_M60(:,n-2:n,:),i);
        end
        legend(["gjennomsnitt","Linjær regresjon: Stigningstall = "+num2str(round(b(sal),3))],'Location','northwest')
        figure(1+(i-1)*2)
    end
    sgtitle(betong(i))
    figure(i*2)
    sgtitle(betong(i)+' Gjennomsnitt')
end
toc

function b = linreg(X,M90_M60,i)
    Y = mean(M90_M60(X,:,i),2);
    Y = Y.';
    XY = X.*Y;
    X2 = X.^2;
    Y2 = Y.^2;
    a = (sum(Y)*sum(X2)-sum(X)*sum(XY))/(length(X)*sum(X2)-sum(X)^2);
    b = (length(X)*sum(XY)-sum(X)*sum(Y))/(length(X)*sum(X2)-sum(X)^2);
    Y(1) = a + b*(X(1));
    Y(2) = a + b*(X(end));
    if i == 1
        plot([X(1),X(end)]-5,[Y(1),Y(2)])
    else
        plot([X(1),X(end)],[Y(1),Y(2)])
    end
end
%lager filnavn og henter ut data fra filen
function data = les_fil(dag)
    filnavn = "dag" + string(dag) + ".csv";
    if isfile(filnavn)
        data = readmatrix(filnavn);
        data = data(:,2:2:end);
    else
        data = [];
    end
end