Skip to content

GitLab

Commit ba54ecbb authored by Michael Rudolf's avatar Michael Rudolf Committed by Michael Rudolf
Browse files

Minor formatting changes in VSTanalysis 

Added some more files to .gitignore
parent 06bcae10
Hide whitespace changes
Inline Side-by-side
.gitignore
View file @ ba54ecbb
......@@ -2,14 +2,21 @@
Tests/
.ipynb_checkpoints/
Examples/
.vscode/
# ignore data files
*.tdms
*.h5
*.txt
*.asc
# ignore pictures and pdfs
*.pdf
*.png
*.jpg
*.eps
*.txt
*.tdms
# Byte-compiled / optimized / DLL files
__pycache__/
......
RSTanalysis/VST/Scripts/VSTanalysis.py
View file @ ba54ecbb
# -*- coding: utf-8 -*-
"""
Spyder Editor
Authors: Michael Warsitzka, Matthias Rosenau, Malte Ritter
Date (last modified): 07.11.2018
When using the data please use the citation as given in the description of data of this data publication.
Authors: Michael Warsitzka, Matthias Rosenau, Malte Ritter, Michael Rudolf
Date (last modified): 26.06.2019
When using the data please use the citation as given in the description of
data of this data publication.
"""
#%%=============IMPORT=====================================================
# %%=============IMPORT=====================================================
import os
import numpy as np
from scipy import stats
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 \
# %%==============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 f in range(0,n):
names = [splitext(files[i])[0] for i in range(0, n)]
# Process all vst files in path:
for f in range(0, n):
data = np.loadtxt(path+files[f], 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)
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)
slope, intercept, r, p, std = stats.linregress(logvel, fric)
fitfric = np.polyval((slope, intercept), logvel)
# Standard deviations:
m=len(logvel)
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 = ({:.4f}$\pm${:.5f})log10(x) +'+'\n'+
'\t'+'\t'+ r'({:.4f}$\pm${:.5f})').format(slope,std_slope,
intercept,std_intercept)
#%%=======================PLOT VST DATA=================================================
# standard error Coef. friction
std_slope = np.sqrt(m*s/(m*np.sum(logvel**2)-np.sum(logvel)**2))
# standard error cohesion
std_intercept = (np.sqrt(np.sum(logvel**2) * s /
(m*np.sum(logvel**2)-np.sum(logvel)**2)))
name_fit = (
r'y = ({:.4f}$\pm${:.5f})log10(x) +' +
'\n\t\t' +
r'({:.4f}$\pm${:.5f})').format(slope,
std_slope,
intercept,
std_intercept)
# %%==================PLOT VST DATA==================================
# set figure parameters:
plt.rcParams['savefig.format'] = 'pdf'
plt.rcParams['font.size'] = 10
plt.rcParams['font.family'] = 'Arial Unicode MS'
plt.rcParams['savefig.dpi'] = 1000
plt.rcParams['figure.figsize'] = (16, 5)
fig=plt.figure()
plt.rcParams['savefig.format'] = 'pdf'
plt.rcParams['font.size'] = 10
plt.rcParams['font.family'] = 'Arial Unicode MS'
plt.rcParams['savefig.dpi'] = 1000
plt.rcParams['figure.figsize'] = (16, 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 = 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('Friction', 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)
# 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.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)
# 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=8,
facecolor='w',
edgecolor='k',
loc= 'upper right',
framealpha=0.5)
ax3.legend(fontsize=8,
facecolor='w',
edgecolor='k',
loc='upper right',
framealpha=0.5)
fig.suptitle(names[f], fontsize=14, y=0.95)
plt.savefig(path+names[f],
bbox_inches='tight',
edgecolor='w')
plt.savefig(path+names[f],
bbox_inches='tight',
edgecolor='w')
plt.close()
print(" ")
print("======================================================================================")
print("==================================================================")
print("Files of material: "+names[f]+" were processed.")
print("======================================================================================")
print("==================================================================")
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment