%dbstop if error;

%make everything nice and tidy
clear;

%get the folder with *.dat - Data from Espresso
pathname = uigetdir('', 'Choose Directory');
if pathname == 0
  return
end

%scan the folder
Mdir = dir(pathname);
nbentries = size(Mdir, 1);
Mfiles = [];
for entry_i = 1:nbentries
  if Mdir(entry_i).isdir == false

    filename = Mdir(entry_i).name;
    if filename(1) ~= '.'
      [p, n, ext] = fileparts(filename);
      if strcmpi(ext, '.dat')
        Mfiles = strvcat(Mfiles, filename);
      end
    end
  end
end
nbfiles = size(Mfiles, 1);

%import *.dat-Files
for i=1:nbfiles;
    filename_info = Mfiles(i,:);
    filename = strrep(filename_info, ' ', ''); %oder deblank(filename_info)
    data(i).file = filename;
    [data(i).Time,data(i).Strain,data(i).Shearstress] = importfile(filename);
end

%analyze peaks and cycles in data
for i=1:nbfiles;
    signum = sign(data(i).Strain)
    halfcycle = find(diff(signum) < 0)
    endofcycle = find(diff(signum) > 0)
    [data(i).strainpeaks,data(i).timepeaksloc] = findpeaks(data(i).Strain)
    [data(i).stresspeaks,data(i).stresspeaksloc] = findpeaks(data(i).Shearstress)
end

%determine strain, modulus, periodlength, frequency, number of cycles and
%samplingrate
for i=1:nbfiles;
    data(i).totalstrain=data(i).strainpeaks(1)
    data(i).modulus=data(i).Shearstress./data(i).Strain
    data(i).periodlength=data(i).Time(end)/numel(data(i).strainpeaks)
    data(i).frequency=1/data(i).periodlength
    data(i).numberofcycles=numel(data(i).strainpeaks)
    data(i).samplingrate=numel(data(i).Time)/data(i).numberofcycles
end

%Fourier-Analysis of Strain
m = 1
for i=1:nbfiles;

    N = numel(data(i).Time) % FFT-Length ~2^x otherwise DFT
    fn = data(i).samplingrate/2 % Nyquist-Frequency
    df = data(i).samplingrate/N %Resolution
    data(i).DFT_Strain=fft(data(i).Strain,N)
    data(i).DFT_Strain_Amplitudes=abs(data(i).DFT_Strain)
    data(i).amplitude_output_strain=fftshift(data(i).DFT_Strain_Amplitudes/N);
    x_fn = 0 : df : fn-df;
    x_fa = 0 : df : data(i).samplingrate-df;

    figure(m)
    stem(x_fa-fn, data(i).amplitude_output_strain, 'b.-')
    axis([-fn fn 0 max(data(i).Strain)/2*1.1])
    title('Amplitude-Output-Strain')
    ylabel('Amplitude')
    xlabel(['Resolution ',num2str(df),' Frequency in Hz'])
    grid

    m = m+1

    data(i).DFT_Strain_Amplitudes_2=data(i).DFT_Strain_Amplitudes'
    data(i).amplitude_output_strain_2 = [data(i).DFT_Strain_Amplitudes_2(1)/N data(i).DFT_Strain_Amplitudes_2(2:N/2)/(N/2)]; % DC-Bin auf N normieren!
    figure(m)
    stem(x_fn, data(i).amplitude_output_strain_2, 'b.-')
    axis([0 fn 0 max(data(i).Strain)*1.1])
    title('Amplitude-Output-Strain')
    ylabel('Amplitude')
    xlabel(['Resolution ',num2str(df),' Frequency in Hz'])
    grid

    m = m+1

end

%Fourier-Analysis of Shearstress
% should i add x = detrend(x) =???
n = m+1
for i=1:nbfiles;

    N = numel(data(i).Time) % FFT-Length ~2^x otherwise DFT
    fn = data(i).samplingrate/2 % Nyquist-Frequency
    df = data(i).samplingrate/N %Resolution
    data(i).DFT_Shearstress=fft(data(i).Shearstress,N)
    data(i).DFT_Shearstress_Amplitudes=abs(data(i).DFT_Shearstress)
    data(i).amplitude_output=fftshift(data(i).DFT_Shearstress_Amplitudes/N);
    x_fn = 0 : df : fn-df;
    x_fa = 0 : df : data(i).samplingrate-df;

    figure(n)
    stem(x_fa-fn, data(i).amplitude_output, 'b.-')
    axis([-fn fn 0 max(data(i).Shearstress)/2*1.1])
    title('Amplitude-Output-Shearstress')
    ylabel('Amplitude')
    xlabel(['Resolution ',num2str(df),' Frequency in Hz'])
    grid

    n = n+1

    data(i).DFT_Shearstress_Amplitudes_2=data(i).DFT_Shearstress_Amplitudes'
    data(i).amplitude_output_shearstress_2 = [data(i).DFT_Shearstress_Amplitudes_2(1)/N data(i).DFT_Shearstress_Amplitudes_2(2:N/2)/(N/2)]; % DC-Bin auf N normieren!
    figure(n)
    stem(x_fn, data(i).amplitude_output_shearstress_2, 'b.-')
    axis([0 fn 0 max(data(i).Shearstress)*1.1])
    title('Amplitude-Output-Shearstress')
    ylabel('Amplitude')
    xlabel(['Resolution ',num2str(df),' Frequency in Hz'])
    grid

    n = n+1

end

%Fourier-Analysis of Phase
o = n+1
for i=1:nbfiles;

    N = numel(data(i).Time)
    fn = data(i).samplingrate/2 % Nyquist-Frequency
    df = data(i).samplingrate/N
    x_fn = 0 : df : fn-df;
    x_fa = 0 : df : data(i).samplingrate-df;

    data(i).Phase = angle(fft(data(i).DFT_Shearstress));
    data(i).DFT_Phase = unwrap(data(i).Phase)
    data(i).DFT_Phase_2=data(i).DFT_Phase';
    data(i).phase_2 = [data(i).DFT_Phase_2(1)/N data(i).DFT_Phase_2(2:N/2)/(N/2)];

    figure(o)
    stem(x_fn, data(i).phase_2, 'b.-')
    axis([0 fn min(data(i).phase_2)*1.1 0])
    title('Phase-Shearstress')
    ylabel('Phase')
    xlabel(['Frequency in Hz'])
    grid
    o = o+1
    
end

%Determine G' and G''
for i=1:nbfiles;
    
    N = numel(data(i).Time) % FFT-Length ~2^x otherwise DFT
    fn = data(i).samplingrate/2 % Nyquist-Frequency
    df = data(i).samplingrate/N 
    data(i).frequencies = 0 : df : fn-df;
    
    data(i).Phase_2=data(i).Phase'
    data(i).phase_3 = [data(i).Phase_2(1)/N data(i).Phase_2(2:N/2)/(N/2)];
    k = find(data(i).frequencies==data(i).frequency)
    
    data(i).G_Prime = (max(data(i).Shearstress)/max(data(i).Strain))*cos(abs(data(i).phase_3(k)))
    data(i).G_DPrime = (max(data(i).Shearstress)/max(data(i).Strain))*sin(abs(data(i).phase_3(k)))
    
end