% PROGRAMMA PER L'ELABORAZIONE DEI DATI IN OUTPUT
% DEL SOFTWARE DI CONTROLLO DELLA CELLA EDOMETRICA
% NB:
% Questo programma cerca i file di testo 
% contenenti i dati delle prove edometriche
% a partire dalla directory dove si trova

clc
clear all
close all

% nomefile = input('Inserisci il nome del file: ','s');
% nomefile = strcat(nomefile,'.txt');
% eval(['fid = fopen(',nomefile,');']);
% fid = fopen('edom0084.txt');
% a = fscanf(fid,'%g %g',[2 inf]) % It has two rows now.
% a = a';
% fclose(fid)


[filename,pathname]=uigetfile({'*.txt',  'All Files (*.*)'}, 'Seleziona il file della prova edometrica')
% Carico il file in una matrice (l'istruzione %q toglie le "")
[a,b,c,d,e,f,g] = textread([pathname, filename],'%q %q %q %q %q %q %q', 'delimiter',',','emptyvalue',NaN);

% Sistemo il dato di profondità
prof = b{7,1}
[Profmin,Profmax] = strread(prof,'%f %f','delimiter','-')
Profmed = (Profmin+Profmax)/2;

% Leggo gli altri dati
SezProv = str2double(b(13));
Hiniz = str2double(b(14));
Hfin = str2double(b(15));
PesoTara1 = str2double(b(17));
TaraPesoUmidoIniz = str2double(b(18));
PesoTara2 = str2double(b(20));
TaraPesoUmidoFin = str2double(b(21));
TaraPesoProvSecco = str2double(b(22));
PesoSpecGran = str2double(b(23));

% Calcolo e0
e0 = (Hiniz*SezProv/10-((TaraPesoProvSecco-PesoTara2)/PesoSpecGran))/((TaraPesoProvSecco-PesoTara2)/PesoSpecGran)

% Salto le righe di intestazione
inum = 29
numpassi = str2double(b(24));

% Memorizzo i dati in vettori
for i =1:numpassi
    numdati(i)= str2double(b(inum));
    Sigmav(i) = str2double(b(inum-1));
    in = inum+3;
    fin = inum+3+numdati(i)-1;
    for j=in:fin
        if (str2double(c(j))~=0);
            Tempo(j-in+1,i) = str2double(b(j));
            Dati(j-in+1,i) = str2double(c(j));
            EpsAssSupp(j-in+1,i) = Dati(j-in+1,i)*100./Hiniz;
        else
            Dati(j-in+1,i) = NaN;
        end
    end
    inum = fin+4;
end

nummax = max(numdati);
for i=1:numpassi
    % Assegno NaN ai valori vuoti finali di ciascun vettore
    Dati(numdati(i)+1:nummax,i) = NaN;
    Tempo(numdati(i)+1:nummax,i) = NaN;
    EpsAssSupp(numdati(i)+1:nummax,i) = NaN;
    % Calcolo le deformazioni e gli indici dei vuoti
    EpsAss(i) = Dati(numdati(i),i)*100./Hiniz;
    e = e0-(Dati./Hiniz.*(1+e0));
    efin(i) = e0-((Dati(numdati(i),i)./Hiniz).*(1+e0));
end

% Raccolgo solo i dati dei passi di carico
j= 1;
k =0;
Carico(:,1) = Dati(:,1);
Tempoc(:,1) = Tempo(:,1);
EpsAssCar(:,1) = EpsAssSupp(:,1);
Sigmavc(1) = Sigmav(1);
EpsAssc(1) = EpsAss(1);

for i = 1:numpassi-1
    % Separo i carichi dagli scarichi
    if Sigmav(i+1)>Sigmav(i)
        j=j+1;
        Carico(:,j) = Dati(:,i+1);
        Tempoc(:,j) = Tempo(:,i+1);
        EpsAssCar(:,j) = EpsAssSupp(:,i+1);
        % Registro le sole tensioni di carico
        Sigmavc(j) = Sigmav(i+1);
        EpsAssc(j) = EpsAss(i+1);
    else
        k=k+1;
        Scarico(:,k) = Dati(:,i+1);
        Tempos(:,k) = Tempo(:,i+1);
    end
end

% Decido se calcolare i cv o caricarli da file
k = menu('Scegli','Determina i c_v manualmente','Caricali da un file di testo preparato')
switch k
    case 1
        %Determino i cv
        for i =1:j
            figure(5)
            subplot(121)
            plot(sqrt(Tempoc(:,i)),Carico(:,i),'.-r')
            set(gca,'YDir','reverse')
            title('Metodo di Taylor')
            subplot(122)
            plot(log10(Tempoc(:,i)),Carico(:,i),'.-r')
            set(gca,'YDir','reverse')
            title('Metodo di Casagrande')
            z = menu(['Scegli per lo step \sigma''_v = ',num2str(Sigmavc(i)),' kPa'],'Metodo di Taylor','Metodo di Casagrande')
            switch z
                case 1
                    % Metodo di Terzaghi
                    figure(4)
                    clf
                    plot(sqrt(Tempoc(:,i)),Carico(:,i),'.-r')
                    set(gca,'YDir','reverse')
                    title('Traccia la linea di compressione e premi INVIO')
                    xlabel('Tempo, \surd{t} [s]');
                    ylabel('Spostamento, dH [mm]');
                    grid on
                    q=gline(4)
                    pause
                    xq = get(q,'XData');
                    yq = get(q,'YData');
                    xq = 1.1*xq;
                    cv(i) = -diff(yq)/diff(xq)
                    text((xq(1)+xq(2))/2,(yq(1)+yq(2))/2,[' c_v = ',num2str(cv(i),'%1.4f')],'Color',[0 0 1],'VerticalAlignment','bottom');
                    pause

                case 2
                    % Metodo grafico di Casagrande per valutare il cv
                    figure(4)
                    clf
                    semilogx(Tempoc(:,i),Carico(:,i),'.-r')
                    set(gca,'YDir','reverse')
                    title({'Traccia la linea di compressione secondaria','passante per (t_{100}, DH_{100}) e premi INVIO'})
                    xlabel('Tempo, t [s]');
                    ylabel('Spostamento, dH [mm]');
                    grid on
                    q=gline(4)
                    pause
                    xq = get(q,'XData');
                    yq = get(q,'YData');
                    cq = diff(yq)/log10(xq(2)/xq(1));
                    title({'Traccia una linea tangente e passante per',' il punto di flesso e (t_{100}, DH_{100}) e premi INVIO'})
                    r=gline(4)
                    pause
                    xr = get(r,'XData');
                    yr = get(r,'YData');
                    cr = diff(yr)/log10(xr(2)/xr(1));
                    t100(i) = 10^((cq*log10(xq(1))-cr*log10(xr(1))+yr(1)-yq(1))/(cq-cr));
                    DH100(i) = cq*(log10(t100(i))-log10(xq(1)))+yq(1);
                    hold on
                    semilogx(t100(i),DH100(i),'+r');
                    t50(i) = t100(i)/2;
                    semilogx([t50(i),t50(i)],[min(Carico(:,i)),max(Carico(:,i))],'-r')
                    text(t50(i),max(Carico(:,i)),['t_{50}= ',num2str(t50(i),'%4.2f'), ' s'],'Color',[1 0 0],'VerticalAlignment','bottom','HorizontalAlignment','left','BackgroundColor',[1 1 1]);
                    title({'Individua un punto a sinistra della retta rossa','e premi INVIO'})
                    [xa,ya] = ginput(1);
                    pause
                    for m=1:length(Carico(:,i))
                        if Tempoc(m,i)>xa
                            break, end
                    end
                    for n=1:length(Carico(:,i))
                        if Tempoc(n,i)>xa/4
                            break, end
                    end
                    ycurv=((Carico(m,i)-Carico(m-1,i))*(log10(xa)-log10(Tempoc(m-1,i)))/(log10(Tempoc(m,i)/Tempoc(m-1,i))))+Carico(m-1,i)
                    ycurv2=((Carico(n,i)-Carico(n-1,i))*(log10(xa/4)-log10(Tempoc(n-1,i)))/(log10(Tempoc(n,i)/Tempoc(n-1,i))))+Carico(n-1,i)
                    % Determino DH0 e H0
                    DH0(i) = 2*ycurv2-ycurv;
                    H0(i) = Hiniz-DH0(i);
                    % Determino H100
                    H100(i) = Hiniz-DH100(i);
                    % Determino H50 ipotetico
                    H50(i) = H100(i)-H0(i);
                    semilogx([xa,xa],[max(Carico(:,i)),ycurv],'--k')
                    text(xa,max(Carico(:,i)),['t_{1}= ',num2str(xa,'%4.2f'), ' s'],'Color',[0 0 0],'VerticalAlignment','bottom','HorizontalAlignment','right','BackgroundColor',[1 1 1]);
                    semilogx([xa/4,xa/4],[max(Carico(:,i)),ycurv2],'--k')
                    semilogx([xa/4,xa],[ycurv2,ycurv2],'--k')
                    semilogx([min(Tempoc(:,i)),xa],[DH0(i),DH0(i)],'--k')
                    text(xa/4,max(Carico(:,i)),['t_{1}/4= ',num2str(xa/4,'%4.2f'), ' s'],'Color',[0 0 0],'VerticalAlignment','bottom','HorizontalAlignment','right','BackgroundColor',[1 1 1]);
                    semilogx([xa,xa],[ycurv,ycurv2],'.-b','LineWidth',3);
                    semilogx([xa,xa],[ycurv2,DH0(i)],'.-r','LineWidth',3);

                    % Dovrei verificare ora che U50 = metà tra U0 e U100!!!
                    % Calcolo c_v
                    cv(i) = (0.197*(H50(i)^2))/(t50(i)*4);          
                    title(['Risulta c_v = ',num2str(cv(i),'%4.4f'),'cm^2/s'])
                    %calfa(i) = -diff(yq)/diff(xq)
                    % Determino il t50
                    %text((xr(1)+xr(2))/2,(yr(1)+yr(2))/2,[' c_{alpha} = ',num2str(calpha(i),'%1.4f')],'Color',[0 0 1],'VerticalAlignment','bottom');
                    pause
            end
            clear z
        end
    case 2
        % Carico i cv dal file *.txt
        cd (pathname)
        cv = load (uigetfile({'*.txt',  'All Files (*.*)'}, 'Seleziona il file contenente i cv'));
end

% Gamma naturale (peso/volume) [g/cm3]
GammaNat= (TaraPesoUmidoIniz-PesoTara1)/(SezProv*Hiniz/10)

%Trasformo Gamma nat in [KN/m3 = kPa]
GammaNat = GammaNat*9.806

% Determino il carico effettivo in sito del campione a quella profondità
ProfFalda = 1;
sigmavin = (GammaNat*Profmed) - 9.806*1*(Profmed-ProfFalda);
% L'incremento di carico è:
DeltaSigma = 0;

% La tensione finale effettiva è [kPa]
sigmavfin = sigmavin+DeltaSigma;

% Plotto il grafico dell'edometrica con tutti i cv e l'intervallo
% tensionale di riferimento del carico in sito

h = figure
subplot(211)
semilogx(Sigmav,EpsAss,'.-r')
set(gca,'YDir','rev');
hold on
for i=1:length(cv) 
    hold on
    text(Sigmavc(i),EpsAssc(i),[' c_v = ',num2str(cv(i))]);
end
hold on
semilogx([sigmavin sigmavin],[0 max(EpsAssc)],'-b');
hold on
semilogx([sigmavfin sigmavfin],[0 max(EpsAssc)],'-b');
xlabel('Tensione verticale, \sigma''_v [kPa]');
ylabel('Deformazione assiale, \epsilon_a [%]');
grid on

subplot(212)
semilogx(Sigmavc,cv,'.-r')
xlabel('Tensione verticale, \sigma''_v [kPa]');
ylabel('c_v delle sole fasi di carico, [cm^2/s]');
grid on



% Determino il cc, il cr e il cs con una interfaccia grafica

    
figure(2)
semilogx(Sigmav,efin,'.-r');
title('Traccia la linea di compressione vergine e premi INVIO')
xlabel('Tensione verticale, \sigma''_v [kPa]');
ylabel('Indice dei vuoti, e');
grid on

k = menu('Scegli come determinare la retta di compressione vergine','Tracciala manualmente','Tracciala in automatico')
switch k
    case 1
        m=gline(2)
        pause
        xc = get(m,'XData');
        yc = get(m,'YData');
        cc = -diff(yc)/log10(xc(2)/xc(1))
    case 2
        [lrt,i] = max(Sigmav);
        xc = [Sigmav(i) Sigmav(i-1)];
        yc = [efin(i) efin(i-1)];
        cc = -diff(yc)/log10(xc(2)/xc(1))
        xe0 = 10^(log10(xc(1))+((e0-yc(1))*log10(xc(2)/xc(1))/(yc(2)-yc(1))));
        hold on
        plot([xe0 xc(1)],[e0 yc(1)],'-k','Linewidth',2.3)
end

hold on
text(10^((log10(xc(1))+log10(xc(2)))/2),(yc(1)+yc(2))/2,['c_c = ',num2str(cc,'%1.4f')],'Color',[0 0 1],'VerticalAlignment','bottom');

title('Traccia la linea di ricompressione e premi INVIO')
n=gline(2)
pause
xr = get(n,'XData');
yr = get(n,'YData');
cr = -diff(yr)/log10(xr(2)/xr(1))
hold on
text(10^((log10(xr(1))+log10(xr(2)))/2),(yr(1)+yr(2))/2,['c_r = ',num2str(cr,'%1.4f')],'Color',[0 0 1],'VerticalAlignment','bottom','BackgroundColor',[1 1 1]);

title('Traccia la linea di scarico e premi INVIO')
l=gline(2)
pause
xs = get(l,'XData');
ys = get(l,'YData');
cs = -diff(ys)/log10(xs(2)/xs(1))
hold on
text(10^((log10(xs(1))+log10(xs(2)))/2),(ys(1)+ys(2))/2,['c_s = ',num2str(cs,'%1.4f')],'Color',[0 0 1],'VerticalAlignment','bottom');

% Plotto la costruzione grafica di Pacheco-Silva (alternativa a Casagrande) per la determinazione della
% tensione di preconsolidazione
% Determino l'intersezione con l'asse di e0
xe0 = 10^(log10(xc(1))+((e0-yc(1))*log10(xc(2)/xc(1))/(yc(2)-yc(1))));
% Determino l'intersezione con la curva edometrica
for i=1:length(Sigmav)
    if Sigmav(i)>xe0
        break,end
end
ycurv=((efin(i)-efin(i-1))*(log10(xe0)-log10(Sigmav(i-1)))/(log10(Sigmav(i)/Sigmav(i-1))))+efin(i-1)
hold on
plot([xe0 xe0],[e0 ycurv],'--k');
% Infine determino la sigma di preconsolidazione
xprecons=10^(log10(xc(1))+((ycurv-yc(1))*log10(xc(2)/xc(1))/(yc(2)-yc(1))));
plot([xe0 xprecons],[ycurv ycurv],'--k');
plot([xprecons xprecons],[ycurv min(efin)],'--k');
plot([Sigmav(1) max(Sigmav)],[e0 e0],'--k','Linewidth',2.5);
text(xprecons,min(efin),['\sigma_p = ',num2str(xprecons,'%1.2f'), ' kPa'],'Color',[1 0 0],'VerticalAlignment','top','HorizontalAlignment','right','BackgroundColor',[1 1 1]);
text(Sigmav(1),e0,['e_0 = ',num2str(e0,'%1.3f')],'Color',[0 0 0],'VerticalAlignment','bottom')

title(['Prova edometrica. Sond. ',char(b(8)), ' Camp. ',char(b(9)),' ',char(b(12)) ' prof. ',char(b(7))]);


figure(3)
semilogx(Tempoc,EpsAssCar,'.-r');
set(gca,'YDir','rev');
for i = 1:length(Sigmavc)
    hold on
    text(max(Tempoc(:,i)),max(EpsAssCar(:,i)),[' \sigma''_v = ',num2str(Sigmavc(i)),' kPa'],'BackgroundColor',[1 1 1]);
end
for i = 1:length(Sigmavc)
title('Traccia la linea di compressione secondaria e premi INVIO')
xlabel('Tempo, t [s]');
ylabel('Deformazione assiale, \epsilon_a [%]');
grid on
set(gca,'XLim',[min(Tempoc(:,i)),max(Tempoc(:,i))],'YLim',[min(EpsAssCar(:,i)),max(EpsAssCar(:,i))])
g=gline(3)
pause
xcalfa = get(g,'XData');
ycalfa = get(g,'YData');
calfa(i) = diff(ycalfa)/100/log10(xcalfa(2)/xcalfa(1))
hold on
text(10^((log10(xcalfa(1))+log10(xcalfa(2)))/2),(ycalfa(1)+ycalfa(2))/2,[' c_{\alpha,\epsilon} = ',num2str(calfa(i),'%1.4f'),' [-]'],'Color',[0 0 1],'VerticalAlignment','bottom','BackgroundColor',[1 1 1]);

title('Posiziona il punto e_{100} e premi INVIO')
[xe100,ye100]=ginput(1);
pause
e100(i) = e0-((ye100/100)*(1+e0));
Eps100(i) = ye100;
plot (xe100,ye100,'+r');
text(xe100,ye100,[' e_{100} = ',num2str(e100(i),'%1.4f'),' [-]'],'Color',[0 0 1],'VerticalAlignment','bottom','BackgroundColor',[1 1 1]);
title(['Determinazione c_{\alpha,\epsilon} Sond. ',char(b(8)), ' Camp. ',char(b(9)),' ',char(b(12)), ' prof. ',char(b(7))]);
end
xlim('auto')
ylim('auto')


% Preparo il cell array dei dati da portare in Excel

A(10,3) = b(8)     % Sondaggio
A(10,7) = b(9)     % Campione
A(10,8) = b(12)    % Note
A(10,11) = b(7)    % Profondità
A(16,4) = b(2)     % Data

Diametro = 2*sqrt(SezProv/pi)
A(18,4) = num2cell(Diametro)       % Diametro provino [cm]
A(19,4) = b(14)                    % Altezza iniziale provino [mm]
A(20,4) = b(15)                    % Altezza finale provino [mm]
A(21,4) = b(13)                    % Sezione provino [cm2]

A(22,4) = b(16)                    % N. tara1
A(23,4) = b(17)                    % Peso anello + peso prov. umido iniz. [g]
A(24,4) = b(18)                    % Peso anello [g]
A(25,4) = b(19)                    % N. tara2
A(26,4) = b(20)                    % Peso tara [g]
A(27,4) = b(21)                    % Peso tara + peso prov. umido fin. [g]
A(28,4) = b(22)                    % Peso tara + peso provino secco [g]
A(30,4) = b(23)                    % Peso specifico grani [g/cm3]

PProvUmidIniz = TaraPesoUmidoIniz-PesoTara1;
A(18,10) = num2cell(PProvUmidIniz)   % Peso provino umido iniz. [g]
PPSecco = TaraPesoProvSecco-PesoTara2;
A(19,10) = num2cell(PPSecco)         % Peso provino secco [g]
PPUmidoFin = TaraPesoUmidoFin-PesoTara2;
A(20,10) = num2cell(PPUmidoFin)      % Peso provino umido finale [g]

WIniz = ((TaraPesoUmidoIniz - PesoTara1) -(TaraPesoProvSecco-PesoTara2))*100/(TaraPesoProvSecco-PesoTara2)
WFin = ((TaraPesoUmidoFin - PesoTara2) - (TaraPesoProvSecco-PesoTara2))*100/(TaraPesoProvSecco-PesoTara2)
A(23,10) = num2cell(WIniz)           % Contenuto d'acqua iniziale [%]
A(24,10) = num2cell(WFin)            % Contenuto d'acqua finale [%]

A(25,10) = num2cell(GammaNat)        % Peso di volume iniziale [kN/m3]
GammaFin = (TaraPesoUmidoFin - PesoTara2)*9.806/(SezProv*Hfin/10)
A(26,10) = num2cell(GammaFin)        % Peso di volume finale [kN/m3]
GammaDryIniz = (TaraPesoProvSecco - PesoTara2)*9.806/(SezProv*Hiniz/10)
A(27,10) = num2cell(GammaDryIniz)    % Peso di volume secco iniziale [kN/m3]
GammaDryFin = (TaraPesoProvSecco - PesoTara2)*9.806/(SezProv*Hfin/10)
A(28,10) = num2cell(GammaDryFin)     % Peso di volume secco finale [kN/m3]

A(33,10) = num2cell(e0)              % Indice dei vuoti iniziale
eFin = (Hfin*SezProv/10-((TaraPesoProvSecco-PesoTara2)/PesoSpecGran))/((TaraPesoProvSecco-PesoTara2)/PesoSpecGran)
A(34,10) = num2cell(eFin)            % Indice dei vuoti finale

% Gli indici di compressione
A(39,4) = num2cell(cr)
A(39,7) = num2cell(cc)
A(39,10) = num2cell(cs)

% La tabella riassuntiva dei dati di ciascun passo di carico

A(47:length(Sigmav)+46,3) = num2cell(Sigmav')       % Tensione [kPa]
A(47:length(EpsAss)+46,4) = num2cell(EpsAss')       % Deformazione [%]
A(47:length(efin)+46,5) = num2cell(efin')           % Indice dei vuoti [-]

% M viene calcolato a partire dalle differenze di deformazione al 100%
% della consolidazione primaria (e100 ed Eps100)
M = (diff(Sigmavc)./(diff(Eps100)/100))/1000;
A(48:length(M)+47,6) = num2cell(M')                 % Modulo Edometrico [MPa]

           
% k = cv.*mv.*1;
% A(48:length(k)+47,7) = num2cell(k)                % Coefficiente di filtrazione
A(48:length(cv)+47,8) = num2cell(cv)                % Coeff. di consolidazione
A(48:length(calfa)+47,10) = num2cell(calfa)         % Coeff. di compressione secondaria


% Scrivo il cell array in un foglio Excel

xlswrite(['Alleg ',char(b(9)),' ',char(b(12)),'2.xls'],A,filename, 'a1')


% Salvo le figure in eps
% saveas(1,[char(b(8)),' ',char(b(9)),' ',char(b(12)),'_cv2.eps'],'psc2')
% saveas(2,[char(b(8)),' ',char(b(9)),' ',char(b(12)),'_edo2.eps'],'psc2')
% saveas(3,[char(b(8)),' ',char(b(9)),' ',char(b(12)),'_calfa2.eps'],'psc2')