# -*- coding: utf-8 -*-
"""
Spyder Editor
Authors: Michael Warsitzka, Matthias Rosenau, Malte Ritter
Date (last modified): 22.10.2018
When using the data please use the citation as given in the description of data of this data publication.
"""
#%%=============IMPORT=====================================================
import os
import numpy as np
from scipy import stats 
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from os.path import splitext
#%%==============PARAMETERS=================================================
A =     0.022619            # area of shear zone (= surface area of lid) (m^2)
li =    0.0776              # inner lever (center cell to center of shear zone)
lo =    0.1250              # outer lever (center cell to hinge point)
v=      30                  # shear velocity (mm/min)
#%%==================NAMES=================================================
path =        './../Data files/' #define path of input data
#%%==========================READ & CONVERT DATA==============================================  
files = [i for i in os.listdir(path) \
         if i.endswith(".txt") and \
         ('VST' in i or 'vst' in i)]
n = len(files)
names = [splitext(files[i])[0] for i in range(0,n)]
#Process all vst files in path:
for i in range(0,n):
    data = np.loadtxt(path+files[i], skiprows=1)
    time= data[:,0]
    vel = data[:,1]
    normalforce = data[:,2]
    shearforce = data[:,3]
    #convert:
    normalstress = normalforce/A                 #calculate normal stress
    shearstress = (shearforce*lo)/(li*A)          #calculate shear stress
    fric = shearstress/normalstress                        #friction []
    displ = np.cumsum((time[1]-time[0])*vel)              #cumulative displacement in [mm]
#%%==========FIT SHEAR LOAD CURVE============================================
    logvel=np.log10(vel)
    # linear fit and evaluation:
    slope, intercept,r,p,std = stats.linregress(logvel,fric)
    fitfric=   np.polyval((slope,intercept),logvel) 
    # Standard deviations:
    m=len(logvel)
    s = np.sum((fric-intercept-slope*logvel)**2)/(m-2)
    std_slope = np.sqrt(m*s/(m*np.sum(logvel**2)-np.sum(logvel)**2))                    #standard error Coef. friction
    std_intercept = np.sqrt(np.sum(logvel**2) * s/ (m*np.sum(logvel**2)-np.sum(logvel)**2))       #standard error cohesion
    name_fit = r'y = ({:.3f}$\pm${:.3f})log10(x) + ({:.3f}$\pm${:.3f})'.format(slope,
                          std_slope,intercept,std_intercept) 
#%%=======================PLOT VST DATA================================================= 
    # set figure parameters:
    plt.rcParams['savefig.format'] =        'pdf'
    plt.rcParams['font.size'] =             9
    plt.rcParams['font.family'] =           'Arial Unicode MS'
    plt.rcParams['savefig.dpi'] =           1000
    plt.rcParams['figure.figsize'] =        (15, 5)
    fig=plt.figure()
    gs = gridspec.GridSpec(1, 3, wspace=0.4, hspace=0.3)
    ax1=fig.add_subplot(gs[0,:2]) 
    ax2=ax1.twinx()
    ax3=fig.add_subplot(gs[0,2]) 
    
    #Displacement vs. friction:
    ax1.plot(displ,fric,'royalblue', lw=.2)  
    ax1.set_xlabel('Shear displacement $d$ [mm]')
    ax1.set_ylabel(r'Friction ($\tau$/$\sigma_N$)', color='royalblue')
    ax1.spines['left'].set_color('royalblue')
    ax1.tick_params('y', which='both', colors='royalblue')
    ax1.set_xlim(0, max(displ)+max(displ)/100)
    
    #Displacement vs. shear velocity:
    ax2.plot(displ,vel, 'r', lw=1)
    ax2.set_xlabel('Shear displacement $d$ [mm]')
    ax2.set_ylabel('Shear velocity $v$ [$\mathregular{mm s^{-1}}$]', color='r')
    ax2.set_yscale('log')
    ax2.spines['right'].set_color('r')
    ax2.spines['left'].set_color('royalblue')
    ax2.tick_params('y', which='both',colors='r')
    ax2.set_xlim(0, max(displ)+max(displ)/100)
    
    #shear velocity vs. fricton:
    ax3.plot(vel,fric, 'k.', ms=2)
    ax3.plot(vel,fitfric, 'g-', label='curve fit: '+name_fit)
    ax3.set_xscale('log')
    ax3.set_xlabel('Shear velocity $v$ [$\mathregular{mm s^{-1}}$]')
    ax3.set_ylabel('Friction')
    ax3.set_xlim(min(vel)-(min(vel)/5), max(vel)+max(vel)/5)
    ax3.legend(fontsize=7.5, 
             facecolor='w', 
             edgecolor='k', 
             loc= 'upper right',
             framealpha=0.5)
    
    fig.suptitle(names[i])
    plt.savefig(path+names[i], 
            bbox_inches='tight', 
            edgecolor='w')
    plt.close()