From 82c5f0e5d5040d0f741ffff3e238c7358878df51 Mon Sep 17 00:00:00 2001 From: eghummel <127259592+eghummel@users.noreply.github.com> Date: Tue, 7 Mar 2023 18:11:50 -0500 Subject: [PATCH] Initial Upload Initial upload of BlueDrop code to github. --- AccPickTest.py | 211 ++++++ BD code v3.py | 223 ++++++ BD code v4.2.py | 258 +++++++ BD code v4.py | 411 ++++++++++ BD code v5.py | 574 ++++++++++++++ BD code v6.py | 979 ++++++++++++++++++++++++ BD code v7.py | 1241 +++++++++++++++++++++++++++++++ BD code_FTF.py | 1004 +++++++++++++++++++++++++ BDcodev4-jupyter.ipynb | 613 +++++++++++++++ Bootleg Dropview.py | 173 +++++ Partially Saturated Sands BC.py | 79 ++ 11 files changed, 5766 insertions(+) create mode 100644 AccPickTest.py create mode 100644 BD code v3.py create mode 100644 BD code v4.2.py create mode 100644 BD code v4.py create mode 100644 BD code v5.py create mode 100644 BD code v6.py create mode 100644 BD code v7.py create mode 100644 BD code_FTF.py create mode 100644 BDcodev4-jupyter.ipynb create mode 100644 Bootleg Dropview.py create mode 100644 Partially Saturated Sands BC.py diff --git a/AccPickTest.py b/AccPickTest.py new file mode 100644 index 0000000..7c372fb --- /dev/null +++ b/AccPickTest.py @@ -0,0 +1,211 @@ +import numpy +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from pathlib import Path +from scipy.signal import find_peaks + +# SETUP VARIABLES - USER INPUTS +BD = 3 +atype = 'p' # m = mantle area, p = projected area +tiptype = 'c' # c = cone, p = parabolic, b = blunt +# paste the filepath to the folder where the bd data is stored +binFilepath = Path("H:\My Drive\CEE 5904 - Project & Report\FRF Data/test data") +# write the bin file number you want to analyze (do not include 'bLog' or '.bin') +fileNum = '02F4' +outputFile = 'data.xlsx' # this currently doesn't do anything, but eventually all data will be printed out into an excel sheet +outputPath = Path("H:\My Drive\CEE 5904 - Project & Report\FRF Data/test data") # Path for new files + +def accPick(d): + # each of the following are the same; if an accelerometer breaks on a BD, can edit that section + # the 200g accelerometer is ignored in all cases + if BD == 3: + maxAcc = d["250g (g)"].max() + global acc + if maxAcc < 5: + if d["2g (g)"].max() < 1.8: # does an extra check for to 200g because of noise + acc = g2g #d["2g (g)"] + else: + acc = d["18g (g)"] + elif maxAcc < 18: + acc = g18g #d["18g (g)"] + elif maxAcc < 50: + acc = g50g #d["50g (g)"] + else: + acc = d["250g (g)"] + + if BD == 2: + maxAcc = d["250g (g)"].max() + if maxAcc < 5: + if d["2g (g)"].max() < 1.8: # does an extra check for to 200g because of noise + acc = d["2g (g)"] + else: + acc = d["18g (g)"] + elif maxAcc < 18: + acc = d["18g (g)"] + elif maxAcc < 50: + acc = d["50g (g)"] + else: + acc = d["250g (g)"] + + if BD == 1: + maxAcc = d["250g (g)"].max() + if maxAcc < 5: + if d["2g (g)"].max() < 1.8: # does an extra check for to 200g because of noise + acc = d["2g (g)"] + else: + acc = d["18g (g)"] + elif maxAcc < 18: + acc = d["18g (g)"] + elif maxAcc < 50: + acc = d["50g (g)"] + else: + acc = d["250g (g)"] + + +# READ BD DATA IN +data_array = [] # creates an empty array for us to fill with bd data +fileName = 'bLog'+fileNum+".bin" +# print(fileName) +newPath = binFilepath / fileName +print(newPath) +file = open(newPath, 'rb') # read file +element = file.read(3) # create a byte list with each element having 3 bytes + +while element: + # Convert to signed integer before adding to data array + iVAl = int.from_bytes(element, byteorder='big', signed=True) + data_array.append(iVAl) # adds the reshaped data from the bd file to the data frame + element = file.read(3) + +np_array = np.array(data_array) # create numpy array from the list +np_array = np.reshape(np_array, (-1, 10)) # convert the 1d array to 2d array with 10 cols + +#print(np_array.shape) +# print(np_array) + +df = pd.DataFrame(np_array) # Creates a Dataframe in pandas from the bd data +df.columns = ['Count', 'no clue', 'g2g', 'g18g', 'g50g', 'ppm', 'g200g', 'gX55g', 'gY55g', 'g250g'] # names columns +# print(dfCal) + +# APPLY CALIBRATION FACTORS +if BD == 3: # calibration factors from July 2019 + g2g = (df['g2g']-34426.5)/1615925.8 # accelerometers are in g + g18g = (df['g18g']+12322.1)/163530.7 + g50g = (df['g50g']-237384.9)/63651 - 0.1120 + ppm = ((df['ppm']+62496.7)/20583.0) * 6.89475729 # converts to kPa + g200g = ((df['g200g'] -248943.7)/39009.4)+0.5518 + gX55g = (df['gX55g']-59093.7)/66674.3 + gY55g = (df['gY55g']-140224.6)/66674.3 + g250g = (df['g250g']-40536.1)/13631.6 + +if BD == 2: # calibration factors from Aug 26, 2021 + g2g = (df['g2g']+31384.7)/1624987.2-0.035 # accelerometers are in g + g18g = (df['g18g']-26631.0)/159945.4 + g50g = (df['g50g']+92987.0)/63783.5 + ppm = ((df['ppm']-35170.6)/12922.9) * 6.89475729 # converts to kPa + g200g = (df['g200g']-16264.8)/26042.8 -0.277 + gX55g = (df['gX55g']+89890.3)/63897.1 + gY55g = (df['gY55g']+14993.0)/64118.0 + g250g = (df['g250g']+17362.1)/13533.5+0.0656 + +if BD == 1: # calibration factors from July 2020 + g2g = (df['g2g']+277743.2)/1637299.6 # accelerometers are in g + g18g = (df['g18g']-3755.9)/159932.2 + g50g = (df['g50g']+92817.6)/63237.1 + ppm = ((df['ppm']-33154.0)/14763.5) * 6.89475729 # this is kPa + g200g = (df['g200g'] -1155309.9)/28368.5 - 1.464 + gX55g = (df['gX55g'] +97138.4)/62023.7 + gY55g = (df['gY55g']-9921.7)/62669.2 + g250g = (df['g250g']+59211.3)/13276.9 + +time = (df['Count']-df['Count'].iloc[0]+1)/2000 # gives time in s + +dfCal = pd.DataFrame([time, g2g, g18g, g50g, g200g, g250g, gX55g, gY55g, ppm]) # creates a calibrated data frame +dfCal = dfCal.T +dfCal.columns = ['Time (s)', '2g (g)', '18g (g)', '50g (g)', '200g (g)', '250g (g)', 'X55g (g)', 'Y55g (g)', 'Pore Pressure (kPa)'] # names columns + +#print(dfCal) + +#AUTOMATIC PEAK FINDING +x = np.array(g200g) # what accelerometer to get the peaks from +peaks, _ = find_peaks(x, height = 5, distance=10000) # finds the largest peaks more than 5g spaced at least 10000 counts apart +plt.plot(x) +plt.plot(peaks, x[peaks], "x") +plt.show() + +# CREATION OF INDIVIDUAL DROP FILES +peaksArray = np.array(peaks) # prints a list of the count where the peaks occur +print(peaksArray) +q = (peaksArray.shape) #gives number of peaks +nDrops = int(q[0]) #number of drops in the file +print(nDrops) + +a = 0 +n = 1 + + # MAKE INDIVIDUAL DATAFRAMES FOR EACH DROP AND PLOT DECELERATION VS TIME +while n <= nDrops : + b = int(peaksArray[a]) # count at the ath drop + dropstart = b - 100 # offset in counts before impact to include in sliced df + dropend = b + 100 # offset in counts after impact to include in sliced df + + if n == 1 : + drop1 = dfCal[dropstart:dropend] + drop1 = pd.DataFrame(drop1) # makes dataframe including all data within the start and end points of the drop + d = drop1 + accPick(d) + print(acc) + drop1.plot(x="Time (s)", y=acc, ylabel="Deceleration (g)") + plt.show() + #print(drop1) + + if n == 2 : + drop2 = dfCal[dropstart:dropend] + drop2 = pd.DataFrame(drop2) # makes dataframe including all data within the start and end points of the drop + d = drop2 + accPick(d) + print(acc) + drop2.plot(x="Time (s)", y="50g (g)", ylabel="Deceleration (g)") + plt.show() + #print(drop2) + + if n == 3 : + drop3 = dfCal[dropstart:dropend] + drop3 = pd.DataFrame(drop3) # makes dataframe including all data within the start and end points of the drop + d = drop3 + accPick(d) + print(acc) + drop3.plot(x="Time (s)", y="50g (g)", ylabel="Deceleration (g)") + plt.show() + #print(drop3) + + if n == 4 : + drop4 = dfCal[dropstart:dropend] + drop4 = pd.DataFrame(drop4) # makes dataframe including all data within the start and end points of the drop + d = drop4 + drop4.plot(x="Time (s)", y="50g (g)", ylabel="Deceleration (g)") + plt.show() + #print(drop4) + + if n == 5 : + drop5 = dfCal[dropstart:dropend] + drop5 = pd.DataFrame(drop5) # makes dataframe including all data within the start and end points of the drop + d = drop5 + drop5.plot(x="Time (s)", y="50g (g)", ylabel="Deceleration (g)") + plt.show() + #print(drop5) + + if n == 6 : + drop6 = dfCal[dropstart:dropend] + drop6 = pd.DataFrame(drop6) # makes dataframe including all data within the start and end points of the drop + d = drop6 + drop6.plot(x="Time (s)", y="50g (g)", ylabel="Deceleration (g)") + plt.show() + #print(drop6) + + #outputName = fileNum + " Drop " + str(n) + #print(outputName) + n = n + 1 + a = a + 1 + diff --git a/BD code v3.py b/BD code v3.py new file mode 100644 index 0000000..324c942 --- /dev/null +++ b/BD code v3.py @@ -0,0 +1,223 @@ +import numpy +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from pathlib import Path + +def tellme(s): + #print(s) + ax1.set(title=s) + plt.draw() + +# SETUP VARIABLES - USER INPUTS +BD = 3 +atype = 'p' # m = mantle area, p = projected area +tiptype = 'c' # c = cone, p = parabolic, b = blunt +# paste the filepath to the folder where the bd data is stored +binFilepath = Path("H:\My Drive\CEE 5904 - Project & Report\FRF Data/test data") +# write the bin file number you want to analyze (do not include 'bLog' or '.bin') +fileNum = '02F4' +outputFile = 'data.xlsx' # this currently doesn't do anything, but eventually all data will be printed out into an excel sheet +outputPath = Path("H:\My Drive\CEE 5904 - Project & Report\FRF Data/test data") # Path for new files + +# READ BD DATA IN +data_array = [] # creates an empty array for us to fill with bd data +fileName = 'bLog'+fileNum+".bin" +# print(fileName) +newPath = binFilepath / fileName +print(newPath) +file = open(newPath, 'rb') # read file +element = file.read(3) # create a byte list with each element having 3 bytes + +while element: + # Convert to signed integer before adding to data array + iVAl = int.from_bytes(element, byteorder='big', signed=True) + data_array.append(iVAl) # adds the reshaped data from the bd file to the data frame + element = file.read(3) + +np_array = np.array(data_array) # create numpy array from the list +np_array = np.reshape(np_array, (-1, 10)) # convert the 1d array to 2d array with 10 cols + +print(np_array.shape) +# print(np_array) + +df = pd.DataFrame(np_array) # Creates a Dataframe in pandas from the bd data +df.columns = ['Count', 'no clue', 'g2g', 'g18g', 'g50g', 'ppm', 'g200g', 'gX55g', 'gY55g', 'g250g'] # names columns +# print(dfCal) + +# APPLY CALIBRATION FACTORS +if BD == 3: # calibration factors from July 2019 + g2g = (df['g2g']-34426.5)/1615925.8 # accelerometers are in g + g18g = (df['g18g']+12322.1)/163530.7 + g50g = (df['g50g']-237384.9)/63651 - 0.1120 + ppm = ((df['ppm']+62496.7)/20583.0) * 6.89475729 # converts to kPa + g200g = ((df['g200g'] -248943.7)/39009.4)+0.5518 + gX55g = (df['gX55g']-59093.7)/66674.3 + gY55g = (df['gY55g']-140224.6)/66674.3 + g250g = (df['g250g']-40536.1)/13631.6 + +if BD == 2: # calibration factors from Aug 26, 2021 + g2g = (df['g2g']+31384.7)/1624987.2-0.035 # accelerometers are in g + g18g = (df['g18g']-26631.0)/159945.4 + g50g = (df['g50g']+92987.0)/63783.5 + ppm = ((df['ppm']-35170.6)/12922.9) * 6.89475729 # converts to kPa + g200g = (df['g200g']-16264.8)/26042.8 -0.277 + gX55g = (df['gX55g']+89890.3)/63897.1 + gY55g = (df['gY55g']+14993.0)/64118.0 + g250g = (df['g250g']+17362.1)/13533.5+0.0656 + +if BD == 1: # calibration factors from July 2020 + g2g = (df['g2g']+277743.2)/1637299.6 # accelerometers are in g + g18g = (df['g18g']-3755.9)/159932.2 + g50g = (df['g50g']+92817.6)/63237.1 + ppm = ((df['ppm']-33154.0)/14763.5) * 6.89475729 # this is kPa + g200g = (df['g200g'] -1155309.9)/28368.5 - 1.464 + gX55g = (df['gX55g'] +97138.4)/62023.7 + gY55g = (df['gY55g']-9921.7)/62669.2 + g250g = (df['g250g']+59211.3)/13276.9 + +time = (df['Count']-df['Count'].iloc[0]+1)/2000 # gives time in s + +dfCal = pd.DataFrame([time, g2g, g18g, g50g, g200g, g250g, gX55g, gY55g, ppm]) # copies the data frame; this version will be edited with calibration factors +dfCal= dfCal.T +dfCal.columns = ['Time (s)', '2g (g)', '18g (g)', '50g (g)', '200g (g)', '250g (g)', 'X55g (g)', 'Y55g (g)', 'Pore Pressure (kPa)'] # names columns +#print(dfCal) + +# GENERATE PLOTS + +fig, (ax1, ax2) = plt.subplots(2) + +# plot of all deceleration data +ax1.plot(time, g2g, label="2g") +ax1.plot(time, g18g, label="18g") +ax1.plot(time, g50g, label="50g") +#plt.plot(time, ppm) +ax1.plot(time, g200g, label="200g") +#plt.plot(time, gX55g, label="2g") +#plt.plot(time, gY55g, label="2g") +ax1.plot(time, g250g, label="250g") +ax1.legend() +ax1.set(ylabel="Deceleration (g)") +ax1.set(xlabel="Time (s)") + +# Plot pore pressure data +ax2.plot(time, ppm, label="Pore Pressure") +ax2.set(ylabel="Pore Pressure (kPa)") +ax2.set(xlabel="Time (s)") +#plt.show() + +# SELECT START AND END POINTS +tellme("How many drops?") +#plt.show() + +nDrops = int(input("Input the number of Drops in file:")) +#print(nDrops) + +n = 1 +b = 0 +a = 0 + +dropTimes = np.array([[0]*2]*nDrops) + +while n <= nDrops : + tellme('Select start and end of drop #' + str(n)) + pts = plt.ginput(2) + (x0, y0), (x1, y1) = pts + xmin, xmax = sorted([x0, x1]) + dropTimes[b,a] = xmin + a = a + 1 + dropTimes[b,a] = xmax + b = b + 1 + a = a - 1 + n = n + 1 + +print(dropTimes) +''' +#Make a new dataframe for each drop +n = 1 +x1r = 0 +x1c = 0 +x2r = 0 +x2c = 1 + +def slicedDf(dStart,dEnd): + dfCal.iloc[(dfCal["Time (s)"] == dStart : dfCal["Time (s)"] == dEnd),:] + +while n 119500: + dropstart = peak - 1500 + dropend = 120000 + else: + dropstart = peak - 1500 + dropend = peak + 500 + +def accPick(dg, d): #this function picks the smallest accelerometer that's not maxed out to perform the integration on + maxAcc = dg["250g (g)"].max() + global acc + global accName + global accNameg + global accg + if maxAcc < 5 - offset: + if dg["2g (g)"].max() < 1.8 - offset: # does an extra check for the 2g because of noise + acc = d["2g (m/s^2)"] + accg = dg["2g (g)"] + accName = "2g (m/s^2)" + accNameg = "2g (g)" + else: + acc = d["18g (m/s^2)"] + accg = dg["18g (g)"] + accName = "18g (m/s^2)" + accNameg = "18g (g)" + elif maxAcc < 18 - offset: + acc = d["18g (m/s^2)"] + accg = dg["18g (g)"] + accName = "18g (m/s^2)" + accNameg = "18g (g)" + elif maxAcc < 50 - offset: + acc = d["50g (m/s^2)"] + accg = dg["50g (g)"] + accName = "50g (m/s^2)" + accNameg = "50g (g)" + else: + acc = d["250g (m/s^2)"] + accg = dg["50g (g)"] + accName = "250g (m/s^2)" + accNameg = "250g (g)" + +def findchangepts(): #This function picks the moment that the BD impacts the ground + global drop + jlist = list() + global jindex + print("finding start of drop...") + for i in range(4,len(accg)-4): + p1 = 1 + #print(p1) + p2 = i + #print(p2) + p3 = len(accg) + #print(p3) + sample1 = list(accg[p1:p2-1]) + #print(sample1) + sample2 = list(accg[p2:p3]) + #print(sample2) + stat1 = math.log(statistics.variance(sample1)) + stat2 = math.log(statistics.variance(sample2)) + #print(stat1) + + j1 = (i-1)*stat1 + j2 = ((len(accg)-1)-i+1)*stat2 + j = j1+j2 + #print(j) + jlist.append(j) + + drop = min(jlist) + #print("drop is") + #print(drop) + jlist = np.array(jlist) + #print(jlist) + #print(jlist.size) + jlist = np.insert(jlist, 0, (0,0,0,0)) #reshape to match up with dataframe d + jlist = np.append(jlist, (0,0,0,0)) #reshape to match up with dataframe d + #print(jlist.size) #should be 2000 + jindex = np.where(jlist==drop) #finds the index of the drop start + jindex = int(jindex[0]) #converts the index into a number from a tuple + +def finddropend(): #finds the location where the deceleration is 1-offset after the peak + global num1 + global num2 + below0list = list() + for i in range(dropstart+jindex, dropend, 1): + if accg[i] < 1 - offset: + num1 = i - dropstart + #num2 = i-jindex-1 + below0list = np.append(below0list, num1) + num1=int(min(below0list)) + +def integration(d): #integrates the deceleration data to solve for velocity and penetration depth + global vel + global maxvel + global dep + global maxdep + accint = acc[jindex:num1] + vel = scipy.integrate.cumtrapz(accint, x=d["Time (s)"]) # solves for velocity + vel = np.array(vel) + vel = numpy.insert(vel, 0, 0) #not sure what to insert here, but it makes it the right size + vel = np.flip(vel) + maxvel = vel.max() + dep = scipy.integrate.cumtrapz(vel, x=d["Time (s)"]) # solves for penetration depth + dep = numpy.insert(dep, 0, 0) # not sure what to insert here, but it makes it the right size + maxdep = dep.max() + d.insert(9, "Velocity (m/s)", vel) + d.insert(10, "Penetration Depth (m)", dep) + +def areafind(): #finds the embedded area based on the penetration depth, area type, and the tip + global area + a1 = list() #placeholder for the penetrated area at that specific depth + d1 = dep*100 #penetration depth array, in cm + print(len(d1)) + r = list() #placeholder for the radius at that specific depth + if tiptype == 'c': + if atype == 'm': + for k in range(0,len(d1)): + if d1[k]=length: + r.append(4.375) + a1.append(pi*r[k]*(sqrt((r[k]^2)+(length^2)))) + a1[k] = a1[k]/10000 + area = a1 + elif atype == 'p': + for k in range(0,len(d1)): + if d1[k]=length: + r.append(4.375) + a1.append(pi*(r[k])**2) + a1[k] = a1[k]/10000 + area = a1 + + elif tiptype == 'b': + if atype =='m': + for k in range(0,len(d1)): + if d1[k]=length: + r.append(4.375) + a1.append(pi*r[k]**2 + 2*pi*r[k]*length) + a1[k]=a1[k]/10000 + area = a1 + elif atype == 'p': + for k in range(0,len(d1)): + a1.append(pi*4.375^2) + a1[k]=a1[k]/10000 + area = a1 + + '''elif tiptype == "p": + if atype == 'm': + for k in range(1,len(d)): + if d1[k] 119500: + dropstart = peak - 1500 + dropend = 120000 + else: + dropstart = peak - 1500 + dropend = peak + 500 + +def accPick(dg, d): #this function picks the smallest accelerometer that's not maxed out to perform the integration on + maxAcc = dg["250g (g)"].max() + global acc + global accName + global accNameg + global accg + if maxAcc < 5 - offset: + if dg["2g (g)"].max() < 1.8 - offset: # does an extra check for the 2g because of noise + acc = d["2g (m/s^2)"] + accg = dg["2g (g)"] + accName = "2g (m/s^2)" + accNameg = "2g (g)" + else: + acc = d["18g (m/s^2)"] + accg = dg["18g (g)"] + accName = "18g (m/s^2)" + accNameg = "18g (g)" + elif maxAcc < 18 - offset: + acc = d["18g (m/s^2)"] + accg = dg["18g (g)"] + accName = "18g (m/s^2)" + accNameg = "18g (g)" + elif maxAcc < 50 - offset: + acc = d["50g (m/s^2)"] + accg = dg["50g (g)"] + accName = "50g (m/s^2)" + accNameg = "50g (g)" + else: + acc = d["250g (m/s^2)"] + accg = dg["50g (g)"] + accName = "250g (m/s^2)" + accNameg = "250g (g)" + +def findchangepts(): #This function picks the moment that the BD impacts the ground + global drop + jlist = list() + global jindex + print("finding start of drop...") + for i in range(4,len(accg)-4): + p1 = 1 + #print(p1) + p2 = i + #print(p2) + p3 = len(accg) + #print(p3) + sample1 = list(accg[p1:p2-1]) + #print(sample1) + sample2 = list(accg[p2:p3]) + #print(sample2) + stat1 = math.log(statistics.variance(sample1)) + stat2 = math.log(statistics.variance(sample2)) + #print(stat1) + + j1 = (i-1)*stat1 + j2 = ((len(accg)-1)-i+1)*stat2 + j = j1+j2 + #print(j) + jlist.append(j) + + drop = min(jlist) + #print("drop is") + #print(drop) + jlist = np.array(jlist) + #print(jlist) + #print(jlist.size) + jlist = np.insert(jlist, 0, (0,0,0,0)) #reshape to match up with dataframe d + jlist = np.append(jlist, (0,0,0,0)) #reshape to match up with dataframe d + #print(jlist.size) #should be 2000 + jindex = np.where(jlist==drop) #finds the index of the drop start + jindex = int(jindex[0]) #converts the index into a number from a tuple + +def finddropend(n): #finds the location where the deceleration is 1-offset after the peak + global num1 + global num2 + below0list = list() + #for i in range(dropstart+jindex, dropend, 1): + for i in range(peaksArray[n], dropend, 1): + if accg[i] < 1 - offset: + num1 = i - dropstart + #num2 = i-jindex-1 + below0list = np.append(below0list, num1) + num1=int(min(below0list)) + +def integration(d): #integrates the deceleration data to solve for velocity and penetration depth + global vel + global maxvel + global dep + global maxdep + accint = acc[jindex:num1] + print(len(d)) + print(len(accint)) + vel = scipy.integrate.cumtrapz(accint, x=d["Time (s)"]) # solves for velocity + vel = np.array(vel) + vel = numpy.insert(vel, 0, 0) #not sure what to insert here, but it makes it the right size + vel = np.flip(vel) + maxvel = vel.max() + maxvel = round(maxvel,1) + dep = scipy.integrate.cumtrapz(vel, x=d["Time (s)"]) # solves for penetration depth + dep = numpy.insert(dep, 0, 0) # not sure what to insert here, but it makes it the right size + maxdep = dep.max() + maxdep = round(maxdep,4) + d.insert(9, "Velocity (m/s)", vel) + d.insert(10, "Penetration Depth (m)", dep) + +def areafind(): #finds the embedded area based on the penetration depth, area type, and the tip + global area + a1 = list() #placeholder for the penetrated area at that specific depth + d1 = dep*100 #penetration depth array, in cm + #print(len(d1)) + r = list() #placeholder for the radius at that specific depth + if tiptype == 'c': + if atype == 'm': + for k in range(0,len(d1)): + if d1[k]=length: + r.append(4.375) + a1.append(pi*r[k]*(sqrt((r[k]**2)+(length**2)))) + a1[k] = a1[k]/10000 + area = a1 + #print(r) + print(np.column_stack((d1,r, area))) + elif atype == 'p': + for k in range(0,len(d1)): + if d1[k]=length: + r.append(4.375) + a1.append(pi*(r[k])**2) + a1[k] = a1[k]/10000 + area = a1 + + elif tiptype == 'b': + if atype =='m': + for k in range(0,len(d1)): + if d1[k]=length: + r.append(4.375) + a1.append(pi*r[k]**2 + 2*pi*r[k]*length) + a1[k]=a1[k]/10000 + area = a1 + elif atype == 'p': + for k in range(0,len(d1)): + a1.append(pi*4.375**2) + a1[k]=a1[k]/10000 + area = a1 + + '''elif tiptype == "p": + if atype == 'm': + for k in range(1,len(d)): + if d1[k] 119500: + dropstart = peak - 1500 + dropend = 120000 + else: + dropstart = peak - 1500 + dropend = peak + 500 + +def accPick(dg, d): #this function picks the smallest accelerometer that's not maxed out to perform the integration on + maxAcc = dg["250g (g)"].max() + global acc + global accName + global accNameg + global accg + if maxAcc < 5 - offset: + if dg["2g (g)"].max() < 1.8 - offset: # does an extra check for the 2g because of noise + acc = d["2g (m/s^2)"] + accg = dg["2g (g)"] + accName = "2g (m/s^2)" + accNameg = "2g (g)" + else: + acc = d["18g (m/s^2)"] + accg = dg["18g (g)"] + accName = "18g (m/s^2)" + accNameg = "18g (g)" + elif maxAcc < 18 - offset: + acc = d["18g (m/s^2)"] + accg = dg["18g (g)"] + accName = "18g (m/s^2)" + accNameg = "18g (g)" + elif maxAcc < 50 - offset: + acc = d["50g (m/s^2)"] + accg = dg["50g (g)"] + accName = "50g (m/s^2)" + accNameg = "50g (g)" + else: + acc = d["250g (m/s^2)"] + accg = dg["50g (g)"] + accName = "250g (m/s^2)" + accNameg = "250g (g)" + print("acc: ", acc) + +def findchangepts(): #This function picks the moment that the BD impacts the ground + global drop + jlist = list() + global jindex + print("finding start of drop...") + for i in range(4,len(accg)-4): + p1 = 1 + #print(p1) + p2 = i + #print(p2) + p3 = len(accg) + #print(p3) + sample1 = list(accg[p1:p2-1]) + #print(sample1) + sample2 = list(accg[p2:p3]) + #print(sample2) + stat1 = math.log(statistics.variance(sample1)) + stat2 = math.log(statistics.variance(sample2)) + #print(stat1) + + j1 = (i-1)*stat1 + j2 = ((len(accg)-1)-i+1)*stat2 + j = j1+j2 + #print(j) + jlist.append(j) + + drop = min(jlist) + #print("drop is") + #print(drop) + jlist = np.array(jlist) + #print(jlist) + #print(jlist.size) + jlist = np.insert(jlist, 0, (0,0,0,0)) #reshape to match up with dataframe d + jlist = np.append(jlist, (0,0,0,0)) #reshape to match up with dataframe d + #print(jlist.size) #should be 2000 + jindex = np.where(jlist==drop) #finds the index of the drop start + jindex = int(jindex[0]) #converts the index into a number from a tuple + print("Jlist: ", jlist) + +def finddropend(n): #finds the location where the deceleration is 1-offset after the peak + global num1 + global num2 + below0list = list() + #for i in range(dropstart+jindex, dropend, 1): + for i in range(peaksArray[n], dropend, 1): + if accg[i] < 1 - offset: + num1 = i - dropstart + #num2 = i-jindex-1 + below0list = np.append(below0list, num1) + num1=int(min(below0list)) + +def integration(d): #integrates the deceleration data to solve for velocity and penetration depth + global vel + global maxvel + global dep + global maxdep + accint = acc[jindex:num1] + print(len(d)) + print(len(accint)) + vel = scipy.integrate.cumtrapz(accint, x=d["Time (s)"]) # solves for velocity + vel = np.array(vel) + vel = numpy.insert(vel, 0, 0) #not sure what to insert here, but it makes it the right size + vel = np.flip(vel) + maxvel = vel.max() + maxvel = round(maxvel,1) + dep = scipy.integrate.cumtrapz(vel, x=d["Time (s)"]) # solves for penetration depth + dep = numpy.insert(dep, 0, 0) # not sure what to insert here, but it makes it the right size + maxdep = dep.max() + maxdep = round(maxdep,4) + d.insert(9, "Velocity (m/s)", vel) + d.insert(10, "Penetration Depth (m)", dep) + +def peakpick(): + global penstart + global penend + fig, (ax1) = plt.subplots(1) + plt.plot(dfCalg['Count'], dfCalg[accNameg], label=accNameg) + ax1.set_xlim(left=dropstart, right=dropend) + ax1.legend() + ax1.set(ylabel="Deceleration (g)") + ax1.set(xlabel="Time") + #ax1.set_title('Zoom into drop start') + #time.sleep(5) + ax1.set_title('Select start and end of drop #' + str(n)) + startendpt = plt.ginput(2, 0) + pentimes = [] + for t in startendpt: + for x in t: + pentimes.append(x) + + penstart = int(pentimes[0]) + penend = int(pentimes[len(pentimes)-2]) + print("start of penetration: ", penstart) + + print("end of penetration ", penend) + +def CheckingFunction(): + #This function is added to check the values of the start and end of the drop if the autofinding function doesn't work + + global jindex + global num1 + + xvalues=np.linspace(0,len(dg)-1,len(dg)) + plt.figure(num=1) + temp_dg=dg['250g (g)'] + plt.plot(xvalues,temp_dg) + plt.grid(visible=True,which='both') + plt.show() + print ('The start of the drop is', jindex) + print ('The end of the drop is', num1) + print('Are the start and end of the drop correct?(y/n)') + drop_input=input() + if drop_input=='n': + peakpick() + jindex = penstart + num1 = penend + +def areafind(): #finds the embedded area based on the penetration depth, area type, and the tip + global area + global trunc_index + a1 = list() #placeholder for the penetrated area at that specific depth + d1 = dep*100 #penetration depth array, in cm + #print(len(d1)) + r = list() #placeholder for the radius at that specific depth + if tiptype == 'c': + if atype == 'm': + for k in range(0,len(d1)): + if d1[k]=length: + r.append(4.375) + a1.append(pi*r[k]*(sqrt((r[k]**2)+(length**2)))) + a1[k] = a1[k]/10000 + area = a1 + #print(r) + print(np.column_stack((d1,r, area))) + elif atype == 'p': + for k in range(0,len(d1)): + if d1[k]=length: + r.append(4.375) + if r[k-1]=length: + r.append(4.375) + a1.append(pi*r[k]**2 + 2*pi*r[k]*length) + a1[k]=a1[k]/10000 + area = a1 + elif atype == 'p': + for k in range(0,len(d1)): + a1.append(pi*4.375**2) + a1[k]=a1[k]/10000 + area = a1 + + '''elif tiptype == "p": + if atype == 'm': + for k in range(1,len(d)): + if d1[k]=length: + r[i]=0.04375 + A[i]=np.pi*r[i]**2 + #A[i]=A[i]/10000 + return A + + +def dropstartend(peak): #after locating the peaks, this function chops the minute long file into a smaller segment immediately before and after the peak + global dropstart + global dropend + if peak <= 1500: + dropstart = 1 + dropend = peak + 500 + elif peak > 119500: + dropstart = peak - 1500 + dropend = 120000 + else: + dropstart = peak - 1500 + dropend = peak + 500 + return dropstart,dropend + + +def accPick(dg, d): #this function picks the smallest accelerometer that's not maxed out to perform the integration on + maxAcc = dg["250g (g)"].max() + global acc + global accName + global accNameg + global accg + if maxAcc < 5 - offset: + if dg["2g (g)"].max() < 1.8 - offset: # does an extra check for the 2g because of noise + acc = d["2g (m/s^2)"] + accg = dg["2g (g)"] + accName = "2g (m/s^2)" + accNameg = "2g (g)" + else: + acc = d["18g (m/s^2)"] + accg = dg["18g (g)"] + accName = "18g (m/s^2)" + accNameg = "18g (g)" + elif maxAcc < 18 - offset: + acc = d["18g (m/s^2)"] + accg = dg["18g (g)"] + accName = "18g (m/s^2)" + accNameg = "18g (g)" + elif maxAcc < 50 - offset: + acc = d["50g (m/s^2)"] + accg = dg["50g (g)"] + accName = "50g (m/s^2)" + accNameg = "50g (g)" + else: + acc = d["250g (m/s^2)"] + accg = dg["250g (g)"] #this was set to 50g (g), I assume this should be 250g to match the rest?_FTF + accName = "250g (m/s^2)" + accNameg = "250g (g)" + +def findchangepts(): #This function picks the moment that the BD impacts the ground + + global drop + jlist = list() + global jindex + print("finding start of drop...") + for i in range(4,len(accg)-4): + p1 = 1 + #print(p1) + p2 = i + #print(p2) + p3 = len(accg) + #print(p3) + sample1 = list(accg[p1:p2-1]) + #print(sample1) + sample2 = list(accg[p2:p3]) + #print(sample2) + stat1 = math.log(statistics.variance(sample1)) + stat2 = math.log(statistics.variance(sample2)) + #print(stat1) + + j1 = (i-1)*stat1 + j2 = ((len(accg)-1)-i+1)*stat2 + j = j1+j2 + #print(j) + jlist.append(j)#adds the value of j to the end of the list + + drop = min(jlist) + #print("drop is") + #print(drop) + jlist = np.array(jlist) + #print(jlist) + #print(jlist.size) + jlist = np.insert(jlist, 0, (0,0,0,0)) #reshape to match up with dataframe d + jlist = np.append(jlist, (0,0,0,0)) #reshape to match up with dataframe d + #print(jlist.size) #should be 2000 + jindex = np.where(jlist==drop) #finds the index of the drop start + jindex = int(jindex[0]) #converts the index into a number from a tuple + #jindex = int(jindex) + #print("jindex is") + #print(jindex) + +def finddropend(): #finds the location where the deceleration is 1-offset after the peak + global num1 + global num2 + below0list = list() + for i in range(dropstart+jindex, dropend): + if accg[i] < 1 - offset: + num1 = i - dropstart + #num2 = i-jindex-11 + below0list = np.append(below0list, num1) + num1=int(min(below0list)) + #print("num 1 is ") + #print(num1) + + + +def integration(d): #integrates the deceleration data to solve for velocity and penetration depth + global vel + global maxvel + global dep + global maxdep + accint = acc[jindex:num1] + vel = scipy.integrate.cumtrapz(accint, x=d["Time (s)"]) # solves for velocity + vel = np.array(vel) + vel = np.insert(vel, 0, 0) #not sure what to insert here, but it makes it the right size + vel = np.flip(vel) + maxvel = vel.max() + dep = scipy.integrate.cumtrapz(vel, x=d["Time (s)"]) # solves for penetration depth + dep = np.insert(dep, 0, 0) # not sure what to insert here, but it makes it the right size + maxdep = dep.max() + d.insert(9, "Velocity (m/s)", vel) + d.insert(10, "Penetration Depth (m)", dep) + +#This function is added to check the values of the start and end of the drop for softer drops where the autofinding +#fucntion is not currently working FTF 12/19/2022 +#This function can be commented out so that it does not do anything if the autofinding function is working for the series of drops being worked on +#This function is placed within the while loop for each drop +def CheckingFunction(): + global jindex + global num1 + + xvalues=np.linspace(0,len(dg)-1,len(dg)) + plt.figure(num=1) + temp_dg=dg['250g (g)'] + plt.plot(xvalues,temp_dg) + plt.grid(visible=True,which='both') + plt.show() + print ('The start of the drop is', jindex) + print ('The end of the drop is', num1) + print('Are the start and end of the drop correct?(y/n)') + drop_input=input() + if drop_input=='n': +#these values should be slightly before the start and slightly after the end + print('What is the approximate drop start and drop end?') + approx_start=input() + approx_end=input() + approx_end=int(approx_end) + approx_start=int(approx_start) + test=pd.DataFrame.copy(temp_dg[approx_start:approx_end]) + #this is a little weird, im sure there is a better way to do it but this will pick the start and end of the penetration + #after the correct jindex and num1 values are input_FTF. + for i in range(0,len(test)): + if test[i+test.index[0]]<0: + test[i+test.index[0]]=np.nan + else: test[i+test.index[0]]=test[i+test.index[0]] + jindex_add=pd.Series.first_valid_index(test) + num1_add=pd.Series.last_valid_index(test) + jindex=jindex_add-temp_dg.index[0] + num1=num1_add-temp_dg.index[0] + + plt.plot(xvalues[jindex:num1],temp_dg[jindex:num1]) + plt.grid(visible=True,which='both') + plt.show() + #This is what I originally had for the corrections, this can be removed, it is not accurate enough. I am going to keep it + #in my file for now_FTF + # plt.figure(num=2) + # plt.plot(xvalues[approx_start:approx_end],temp_dg[approx_start:approx_end]) + # plt.grid(visible=True) + # plt.show() + # print('Select the correct drop start') + # jindex=input() + # jindex=int(jindex) + # plt.figure(num=1) + # print ('Select the correct drop end') + # num1=input() + # num1=int(num1) + return jindex,num1 +#The following variable and function are to calculate the QSBC for the drops. The QSBC values do not always look correct +#I need to look into this more_FTF + +#define the cone factors to use +srfk =[0.2, 0.4, 1.0, 1.5] +#computed the quasistartic bearing capacity for the drops for each cone factor selected +def QSBC(): + global qsbc_1 + global qsbc_2 + global qsbc_3 + global qsbc_4 + + force=(mass*d[accName]) + qdyn=force/A #dynamic bearing capacity [Pa] + srcv=np.array(np.log10(d['Velocity (m/s)']/0.02))#velocity portion of the strain correction, 0.02 is the 2cm/s push speed of a CPT + # srfK = [0.2, 0.4 ,1 ,1.5] + frs_1=1+srfk[0]*srcv + qsbc_1=(qdyn/frs_1)/1000 + frs_2=1+srfk[1]*srcv + qsbc_2=(qdyn/frs_2)/1000 + frs_3=1+srfk[2]*srcv + qsbc_3=(qdyn/frs_3)/1000 + frs_4=1+srfk[3]*srcv + qsbc_4=(qdyn/frs_4)/1000 + #This section removes negative values of bearing capacity. This is probably not the best way to do it, I need to check the + #matlab code a littler closer to figure out how this function should be changed_FTF + for i in range (qsbc_1.index[0],qsbc_1.index[0]+len(qsbc_1)): + if qsbc_1[i]<0 or qsbc_1[i]==-np.inf: + qsbc_1[i]=0 + for i in range (qsbc_2.index[0],qsbc_2.index[0]+len(qsbc_2)): + if qsbc_2[i]<0 or qsbc_2[i]==-np.inf: + qsbc_2[i]=0 + for i in range (qsbc_3.index[0],qsbc_3.index[0]+len(qsbc_3)): + if qsbc_3[i]<0 or qsbc_3[i]==-np.inf: + qsbc_3[i]=0 + for i in range (qsbc_4.index[0],qsbc_4.index[0]+len(qsbc_4)): + if qsbc_4[i]<0 or qsbc_4[i]==-np.inf: + qsbc_4[i]=0 + + return qsbc_1, qsbc_2, qsbc_3, qsbc_4 + +# READ BD DATA IN +data_array = [] # creates an empty array for us to fill with bd data +fileName = 'Copy of bLog'+fileNum+".bin" +# print(fileName) +newPath = binFilepath / fileName +print(newPath) +file = open(newPath, 'rb') # read file +element = file.read(3) # create a byte list with each element having 3 bytes + +while element: + # Convert to signed integer before adding to data array + iVAl = int.from_bytes(element, byteorder='big', signed=True) + data_array.append(iVAl) # adds the reshaped data from the bd file to the data frame + element = file.read(3) + +np_array = np.array(data_array) # create numpy array from the list +np_array = np.reshape(np_array, (-1, 10)) # convert the 1d array to 2d array with 10 cols + +print(np_array.shape) +# print(np_array) + +df = pd.DataFrame(np_array) # Creates a Dataframe in pandas from the bd data +df.columns = ['Count', 'no clue', 'g2g', 'g18g', 'g50g', 'ppm', 'g200g', 'gX55g', 'gY55g', 'g250g'] # names columns +# print(dfCal) + +# APPLY CALIBRATION FACTORS +if BD == 3: # calibration factors from July 2019 + g2g = (df['g2g']-34426.5)/1615925.8 - offset# accelerometers are in g + g18g = (df['g18g']+12322.1)/163530.7 - offset + g50g = (df['g50g']-237384.9)/63651 - 0.1120 - offset + ppm = ((df['ppm']+62496.7)/20583.0) * 6.89475729 # converts to kPa + g200g = ((df['g200g'] -248943.7)/39009.4)+0.5518 - offset + gX55g = (df['gX55g']-59093.7)/66674.3 - offset #check if lateral accelerometers also need to be offset + gY55g = (df['gY55g']-140224.6)/66674.3- offset + g250g = (df['g250g']-40536.1)/13631.6 - offset + +if BD == 2: # calibration factors from Aug 26, 2021 + g2g = (df['g2g']+31384.7)/1624987.2-0.035 - offset# accelerometers are in g + g18g = (df['g18g']-26631.0)/159945.4 - offset + g50g = (df['g50g']+92987.0)/63783.5 - offset + ppm = ((df['ppm']-35170.6)/12922.9) * 6.89475729 # converts to kPa + g200g = (df['g200g']-16264.8)/26042.8 -0.277 - offset + gX55g = (df['gX55g']+89890.3)/63897.1 - offset + gY55g = (df['gY55g']+14993.0)/64118.0 - offset + g250g = (df['g250g']+17362.1)/13533.5+0.0656 - offset +#These calibration factors do not match what julie has in her bluedrop code +if BD == 1: # calibration factors from July 2020 + # trying with values from Julies File + g2g = (df['g2g']-42590.9)/1626361.1 - offset # accelerometers are in g + g18g = (df['g18g']-44492.9)/161125.5 - offset + g50g = (df['g50g']-171656.1)/64020.3 - offset + ppm = ((df['ppm']-33154.0)/14763.5) * 6.89475729 # this is kPa + g200g = (df['g200g'] -723404.8)/32209.7 - offset# - 1.464 + gX55g = (df['gX55g'] -54881.1)/64858.6 - offset + gY55g = (df['gY55g']-28735.5)/63839.9 - offset + g250g = (df['g250g']+13299.7)/13697.1 - offset + + # g2g = (df['g2g']+277743.2)/1637299.6 - offset # accelerometers are in g + # g18g = (df['g18g']-3755.9)/159932.2 - offset + # g50g = (df['g50g']+92817.6)/63237.1 - offset + # ppm = ((df['ppm']-33154.0)/14763.5) * 6.89475729 # this is kPa + # g200g = (df['g200g'] -1155309.9)/28368.5 - 1.464 - offset + # gX55g = (df['gX55g'] +97138.4)/62023.7 - offset + # gY55g = (df['gY55g']-9921.7)/62669.2 - offset + # g250g = (df['g250g']+59211.3)/13276.9 - offset + +time = (df['Count']-df['Count'].iloc[0]+1)/2000 # gives time in s +count = df["Count"] + +# make a new dataframe of the calibrated values in units of g +dfCalg = pd.DataFrame([time, g2g, g18g, g50g, g200g, g250g, gX55g, gY55g, ppm]) +dfCalg = dfCalg.T +dfCalg.columns = ['Time (s)', '2g (g)', '18g (g)', '50g (g)', '200g (g)', '250g (g)', 'X55g (g)', 'Y55g (g)', 'Pore Pressure (kPa)'] # names columns +#print(dfCalg) + +#make a new dataframe of the calibrated values in units of m/s^2 +dfCal = pd.DataFrame([time, g2g, g18g, g50g, g200g, g250g, gX55g, gY55g, ppm]) +dfCal = dfCal.T +dfCal.columns = ['Time (s)', '2g (m/s^2)', '18g (m/s^2)', '50g (m/s^2)', '200g (m/s^2)', '250g (m/s^2)', 'X55g (m/s^2)', 'Y55g (m/s^2)', 'Pore Pressure (kPa)'] # names columns +dfCal['2g (m/s^2)'] = dfCal['2g (m/s^2)'] * 9.80665 +dfCal['18g (m/s^2)'] = dfCal['18g (m/s^2)'] * 9.80665 +dfCal['50g (m/s^2)'] = dfCal['50g (m/s^2)'] * 9.80665 +dfCal['200g (m/s^2)'] = dfCal['200g (m/s^2)'] * 9.80665 +dfCal['250g (m/s^2)'] = dfCal['250g (m/s^2)'] * 9.80665 +dfCal['X55g (m/s^2)'] = dfCal['X55g (m/s^2)'] * 9.80665 +dfCal['Y55g (m/s^2)'] = dfCal['Y55g (m/s^2)'] * 9.80665 +#print(dfCal) + +#Locate the drops +x = np.array(g250g) # what accelerometer to get the peaks from + +peaks, _ = find_peaks(x, height = 5, distance=10000) # finds the largest peaks more than 5g spaced at least 10000 counts apart + +peaksArray = np.array(peaks) # prints a list of the count where the peaks occur +#print(peaksArray) +q = (peaksArray.shape) #gives number of peaks +nDrops = int(q[0]) #number of drops in the file +#print(nDrops) + +# For each drop, find the start and end points and integrate to solve for velocity and acceleration + + +# #%% + +# a = 0 +# n = 1 +# for i in range(nDrops): +# peak=peaksArray[i] +# start=dropstartend(peak) +# print('dropstart=',dropstart,'dropend=',dropend) +# print('peak=',peak) +#%% +masslength(tiptype) +a=0 +n=1 +while n <= nDrops : + peak = int(peaksArray[a]) # count at the ath drop + dropstartend(peak) #zooms in the drop file to only consider 500 counts before and 1500 counts after the peak deceleration + print('drop start =',dropstart,'drop end=', dropend) + + if n == 1 : + print('a=',a) + drop1 = dfCal[dropstart:dropend] # start and end points of the drop in m/s^2 + drop1g = dfCalg[dropstart:dropend] # start and end points of the drop in g + drop1 = pd.DataFrame(drop1) # makes dataframe including all data within the start and end points of the drop + drop1g = pd.DataFrame(drop1g) + dg = drop1g + d = drop1 + accPick(dg, d) # chooses what accelerometer to use + acc1 = acc + acc1Name = accName + acc1Nameg = accNameg + findchangepts() #finds the start of the drop + finddropend() #finds the end of the drop + #print(drop) + + CheckingFunction() + + d = d[jindex:num1] #shortens the dataframe to only include the data during penetration (jindex = start, num1 = end) + dg = dg[jindex:num1] + #print(d) + #print(np.size(d)) + drop1 = d + drop1g = dg + integration(d) #solves for velocity and displacement + drop1 = d #this dataframe now includes velocity and acceleration data + A=areafind() + drop1['Cone Area (m^2)']=A #add cone area to the dataframe + + QSBC()#run the QSBC function + #add the QSBC values to the dataframe + drop1['QSBC (kPa) K=0.2']=qsbc_1 + drop1['QSBC (kPa) K=0.4']=qsbc_2 + drop1['QSBC (kPa) K=1.0']=qsbc_3 + drop1['QSBC (kPa) K=1.5']=qsbc_4 + #print(drop1) + + + + if n == 2 : + print('a=',a) + drop2 = dfCal[dropstart:dropend] # start and end points of the drop in m/s^2 + drop2g = dfCalg[dropstart:dropend] # start and end points of the drop in g + drop2 = pd.DataFrame(drop2) # makes dataframe including all data within the start and end points of the drop + drop2g = pd.DataFrame(drop2g) + dg = drop2g # chooses what accelerometer to use based on the max g + d = drop2 + accPick(dg, d) # chooses what accelerometer to use + acc2 = acc + acc2Name = accName + acc2Nameg = accNameg + #print(num1, num2) + #print(acc) + #print(acc.iloc[1]) + findchangepts() + finddropend() + #print(drop) + + CheckingFunction() + + d = d[jindex:num1] + dg = dg[jindex:num1] + drop2 = d + drop2g = dg + # drop1plot = drop1.plot(y=accName, ylabel="Deceleration (g)", title="drop 1") + # drop1plot = plt.plot(acc1Name, acc1Name[num1], "x") + integration(d) + drop2 = d + A=areafind() + drop2['Cone Area (m^2)']=A #add cone area to the dataframe + QSBC()#run the QSBC function + #add the QSBC values to the dataframe + drop2['QSBC (kPa) K=0.2']=qsbc_1 + drop2['QSBC (kPa) K=0.4']=qsbc_2 + drop2['QSBC (kPa) K=1.0']=qsbc_3 + drop2['QSBC (kPa) K=1.5']=qsbc_4 + + + if n == 3 : + print ('a=', a) + drop3 = dfCal[dropstart:dropend] # start and end points of the drop in m/s^2 + drop3g = dfCalg[dropstart:dropend] # start and end points of the drop in g + drop3 = pd.DataFrame(drop3) # makes dataframe including all data within the start and end points of the drop + drop3g = pd.DataFrame(drop3g) + dg = drop3g + d = drop3 + accPick(dg, d) # chooses what accelerometer to use + acc3 = acc + acc3Name = accName + acc3Nameg = accNameg + findchangepts() #finds the start of the drop + finddropend() #finds the end of the drop + #this section of code goes with the ChangeFunction to confirm that the points are selected correctly. + + CheckingFunction() + + d = d[jindex:num1] #shortens the dataframe to only include the data during penetration (jindex = start, num1 = end) + dg = dg[jindex:num1] + drop3 = d + drop3g = dg + integration(d) #solves for velocity and acceleration + drop3 = d #this dataframe now includes velocity and acceleration data + A=areafind() + drop3['Cone Area (m^2)']=A #add cone area to the dataframe + QSBC()#run the QSBC function + #add the QSBC values to the dataframe + drop3['QSBC (kPa) K=0.2']=qsbc_1 + drop3['QSBC (kPa) K=0.4']=qsbc_2 + drop3['QSBC (kPa) K=1.0']=qsbc_3 + drop3['QSBC (kPa) K=1.5']=qsbc_4 + + if n == 4 : + print('a=',a) + drop4 = dfCal[dropstart:dropend] # start and end points of the drop in m/s^2 + drop4g = dfCalg[dropstart:dropend] # start and end points of the drop in g + drop4 = pd.DataFrame(drop4) # makes dataframe including all data within the start and end points of the drop + drop4g = pd.DataFrame(drop4g) + dg = drop4g # chooses what accelerometer to use based on the max g + d = drop4 + accPick(dg, d) # chooses what accelerometer to use + acc4 = acc + acc4Name = accName + acc4Nameg = accNameg + #print(num1, num2) + #print(acc) + #print(acc.iloc[1]) + findchangepts() + finddropend() + #print(drop) + #this section of code goes with the ChangeFunction to confirm that the points are selected correctly. + + CheckingFunction() + + d = d[jindex:num1] + dg = dg[jindex:num1] + #print(np.size(d)) + drop4 = d + drop4g = dg + # drop1plot = drop1.plot(y=accName, ylabel="Deceleration (g)", title="drop 1") + #drop1plot = plt.plot(acc1Name, acc1Name[num1], "x") + integration(d) + drop4 = d + A=areafind() + drop4['Cone Area (m^2)']=A #add cone area to the dataframe + QSBC()#run the QSBC function + #add the QSBC values to the dataframe + drop4['QSBC (kPa) K=0.2']=qsbc_1 + drop4['QSBC (kPa) K=0.4']=qsbc_2 + drop4['QSBC (kPa) K=1.0']=qsbc_3 + drop4['QSBC (kPa) K=1.5']=qsbc_4 + + if n == 5 : + print('a=',a) + drop5 = dfCal[dropstart:dropend] # start and end points of the drop in m/s^2 + drop5g = dfCalg[dropstart:dropend] # start and end points of the drop in g + drop5 = pd.DataFrame(drop5) # makes dataframe including all data within the start and end points of the drop + drop5g = pd.DataFrame(drop5g) + dg = drop5g # chooses what accelerometer to use based on the max g + d = drop5 + accPick(dg, d) # chooses what accelerometer to use + acc5 = acc + acc5Name = accName + acc5Nameg = accNameg + finddropend() + #print(num1, num2) + #print(acc) + #print(acc.iloc[1]) + findchangepts() + + CheckingFunction() + + #print(drop) + d = d[jindex:num1] + dg = dg[jindex:num1] + drop5 = d + drop5g = dg + # drop1plot = drop1.plot(y=accName, ylabel="Deceleration (g)", title="drop 1") + #drop1plot = plt.plot(acc1Name, acc1Name[num1], "x") + integration(d) + drop5 = d + A=areafind() + drop5['Cone Area (m^2)']=A #add cone area to the dataframe + QSBC()#run the QSBC function + #add the QSBC values to the dataframe + drop5['QSBC (kPa) K=0.2']=qsbc_1 + drop5['QSBC (kPa) K=0.4']=qsbc_2 + drop5['QSBC (kPa) K=1.0']=qsbc_3 + drop5['QSBC (kPa) K=1.5']=qsbc_4 + + + if n == 6 : + print('a=',a) + drop6 = dfCal[dropstart:dropend] # start and end points of the drop in m/s^2 + drop6g = dfCalg[dropstart:dropend] # start and end points of the drop in g + drop6 = pd.DataFrame(drop6) # makes dataframe including all data within the start and end points of the drop + drop6g = pd.DataFrame(drop6g) + dg = drop6g # chooses what accelerometer to use based on the max g + d = drop6 + accPick(dg, d) # chooses what accelerometer to use + acc6 = acc + acc6Name = accName + acc6Nameg = accNameg + finddropend() + #print(num1, num2) + #print(acc) + #print(acc.iloc[1]) + findchangepts() + #print(drop) + + CheckingFunction() + + d = d[jindex:num1] + dg = dg[jindex:num1] + drop6 = d + drop6g = dg + # drop1plot = drop1.plot(y=accName, ylabel="Deceleration (g)", title="drop 1") + #drop1plot = plt.plot(acc1Name, acc1Name[num1], "x") + integration(d) + drop6 = d + A=areafind() + drop6['Cone Area (m^2)']=A #add cone area to the dataframe + QSBC()#run the QSBC function + #add the QSBC values to the dataframe + drop6['QSBC (kPa) K=0.2']=qsbc_1 + drop6['QSBC (kPa) K=0.4']=qsbc_2 + drop6['QSBC (kPa) K=1.0']=qsbc_3 + drop6['QSBC (kPa) K=1.5']=qsbc_4 + + + if n == 7 : + print('a=',a) + drop7 = dfCal[dropstart:dropend] # start and end points of the drop in m/s^2 + drop7g = dfCalg[dropstart:dropend] # start and end points of the drop in g + drop7 = pd.DataFrame(drop7) # makes dataframe including all data within the start and end points of the drop + drop7g = pd.DataFrame(drop7g) + dg =drop7g # chooses what accelerometer to use based on the max g + d = drop7 + accPick(dg, d) # chooses what accelerometer to use + acc7 = acc + acc7Name = accName + acc7Nameg = accNameg + finddropend() + #print(num1, num2) + #print(acc) + #print(acc.iloc[1]) + findchangepts() + #print(drop) + + CheckingFunction() + + d = d[jindex:num1] + drop7 = d + drop7g = dg + # drop1plot = drop1.plot(y=accName, ylabel="Deceleration (g)", title="drop 1") + #drop1plot = plt.plot(acc1Name, acc1Name[num1], "x") + integration(d) + drop7 = d + A=areafind() + drop7['Cone Area (m^2)']=A #add cone area to the dataframe + QSBC()#run the QSBC function + #add the QSBC values to the dataframe + drop7['QSBC (kPa) K=0.2']=qsbc_1 + drop7['QSBC (kPa) K=0.4']=qsbc_2 + drop7['QSBC (kPa) K=1.0']=qsbc_3 + drop7['QSBC (kPa) K=1.5']=qsbc_4 + + + n = n + 1 + a = a + 1 +#%% +#combine the dataframes and export to excel csv +if nDrops==1: + frames=[drop1] +elif nDrops==2: + frames=[drop1,drop2] +elif nDrops==3: + frames=[drop1, drop2, drop3] +elif nDrops==4: + frames=[drop1, drop2, drop3, drop4] +elif nDrops==5: + frames=[drop1, drop2, drop3, drop4, drop5] +elif nDrops==6: + frames=[drop1, drop2, drop3, drop4, drop5, drop6] +elif nDrops==7: + frames=[drop1, drop2, drop3, drop4, drop5, drop6, drop7] +data_file=pd.concat(frames) +data_file.to_csv(fileNum+'.csv') +#show plots +#%% +# GENERATE PLOTS +# If there are less than 6 drops higher drop plots will plot blank +# Add more plot sections if there are more than 6 drops + + +#PLot showing all accellerometers and pore pressure readings +fig, (ax1, ax2) = plt.subplots(2, sharex=True) + +ax1.plot(time, g2g, label="2g", ) +ax1.plot(time, g18g, label="18g") +ax1.plot(time, g50g, label="50g") +plt.plot(time, ppm) +ax1.plot(time, g200g, label="200g") +#plt.plot(time, gX55g, label="2g") +#plt.plot(time, gY55g, label="2g") +#ax1.plot(time, g250g, label="250g") +ax1.legend() +ax1.set(ylabel="Deceleration (g)") +# ax1.set(xlabel="Time (s)") + +ax2.plot(time, ppm, label="Pore Pressure") +ax2.set(ylabel="Pore Pressure (kPa)") +ax2.set(xlabel="Time (s)") +plt.show() +#%% +# Plot showing peak deceleration +peakplot = plt.plot(x) +peakplot = plt.plot(peaks, x[peaks], "x") +plt.show() + +#I edited the following plot outputs from what was originally in the file for what I wanted the plots to output as_FTF +#%% +#Deceleration,Velocity,and penetration depth vs time plots for Drop 1 + +fig, (ax1, ax2, ax3) = plt.subplots(3, sharex=True) +ax1.plot(drop1["Time (s)"], drop1[acc1Name], color = "k", marker = 11) +ax1.set(ylabel="Deceleration (m/s^2)") +ax2.plot(drop1["Time (s)"], drop1['Velocity (m/s)'] , color = "k", marker = 11) +ax2.set(ylabel="Velocity (m/s)") +ax3.plot(drop1["Time (s)"], drop1["Penetration Depth (m)"] , color = "k", marker = 11) +ax3.set(ylabel="Penetration Depth (m)", xlabel="Time(s)") +ax1.plot(drop1["Time (s)"], drop1["250g (m/s^2)"], color = "k", marker = 11) +plt.suptitle(fileNum+' Drop 1') + +plt.show() +#%% +#Deceleration,Velocity,and penetration depth vs time plots for Drop 2 +fig, (ax1, ax2, ax3) = plt.subplots(3,sharex=True) +ax1.set(ylabel="Deceleration (m/s^2)") +ax2.plot(drop2["Time (s)"], drop2['Velocity (m/s)'], color = "k", marker = 11) +ax2.set(ylabel="Velocity (m/s)",) +ax3.plot(drop2["Time (s)"], drop2["Penetration Depth (m)"], color = "k", marker = 11) +ax3.set(ylabel="Penetration Depth (m)", xlabel="Time(s)") +ax1.plot(drop2["Time (s)"], drop2["250g (m/s^2)"], color = "k", marker = 11) +plt.suptitle(fileNum+' Drop 2') + +plt.show() +#%% +#Deceleration,Velocity,and penetration depth vs time plots for Drop 3 +fig, (ax1, ax2, ax3) = plt.subplots(3,sharex=True) +ax1.plot(drop3["Time (s)"], drop3["250g (m/s^2)"], marker = 11, color = "k") +ax1.set(ylabel="Deceleration (m/s^2)") + +ax2.plot(drop3["Time (s)"], drop3['Velocity (m/s)'], marker = 11, color = "k") +ax2.set(ylabel="Velocity (m/s)") +ax3.plot(drop3["Time (s)"], drop3["Penetration Depth (m)"], marker = 11, color = "k") +ax3.set(ylabel="Penetration Depth (m)", xlabel="Time(s)") +ax1.plot(drop3["Time (s)"], drop3["250g (m/s^2)"], color = "k", marker = 11) +plt.suptitle(fileNum+' Drop 3') +plt.show() + +#%% +#Deceleration,Velocity,and penetration depth vs time plots for Drop 4 +fig, (ax1, ax2, ax3) = plt.subplots(3,sharex=True) +ax1.plot(drop4["Time (s)"], drop4["250g (m/s^2)"], marker = 11, color = "k") +ax1.set(ylabel="Deceleration (m/s^2)") + +ax2.plot(drop4["Time (s)"], drop4['Velocity (m/s)'], marker = 11, color = "k") +ax2.set(ylabel="Velocity (m/s)", xlabel="Time(s)") +ax3.plot(drop4["Time (s)"], drop4["Penetration Depth (m)"], marker = 11, color = "k") +ax3.set(ylabel="Penetration Depth (m)", xlabel="Time(s)") +ax1.plot(drop4["Time (s)"], drop4["250g (m/s^2)"], color = "k", marker = 11) +plt.suptitle(fileNum+' Drop 4') +plt.show() + +#%% +#Deceleration,Velocity,and penetration depth vs time plots for Drop 5 +fig, (ax1, ax2, ax3) = plt.subplots(3,sharex=True) +ax1.plot(drop5["Time (s)"], drop5["250g (m/s^2)"], marker = 11, color = "k") +ax1.set(ylabel="Deceleration (m/s^2)") + +ax2.plot(drop5["Time (s)"], drop5['Velocity (m/s)'], marker = 11, color = "k") +ax2.set(ylabel="Velocity (m/s)") +ax3.plot(drop5["Time (s)"], drop5["Penetration Depth (m)"], marker = 11, color = "k") +ax3.set(ylabel="Penetration Depth (m)", xlabel="Time(s)") +ax1.plot(drop5["Time (s)"], drop5["250g (m/s^2)"], color = "k", marker = 11) +plt.suptitle(fileNum+' Drop 5') +plt.show() + +#%% +#Deceleration,Velocity,and penetration depth vs time plots for Drop 6 +fig, (ax1, ax2, ax3) = plt.subplots(3,sharex=True) +ax1.plot(drop6["Time (s)"], drop6["250g (m/s^2)"], marker = 11, color = "k") +ax1.set(ylabel="Deceleration (m/s^2)") + +ax2.plot(drop6["Time (s)"], drop6['Velocity (m/s)'], marker = 11, color = "k") +ax2.set(ylabel="Velocity (m/s)") +ax3.plot(drop6["Time (s)"], drop6["Penetration Depth (m)"], marker = 11, color = "k") +ax3.set(ylabel="Penetration Depth (m)", xlabel="Time(s)") +ax1.plot(drop6["Time (s)"], drop6["250g (m/s^2)"], color = "k", marker = 11) +plt.suptitle(fileNum+' Drop 6') +plt.show() + + +#%% +#Deceleration and Velocity profile vs penetration depth for Drop 1 + + +fig, (ay1, ay2) = plt.subplots(1, 2, sharey=True) +ay1.plot(drop1g[acc1Nameg], drop1["Penetration Depth (m)"]*100, color = "k", linestyle = "solid") #marker = 11 +ay1.plot(drop1["Velocity (m/s)"], drop1["Penetration Depth (m)"]*100, color = "k", linestyle = "dashed") +ay1.set(xlabel="Deceleration (g) and Velocity (m/s)", ylabel="Penetration Depth (cm)") +ay1.invert_yaxis() +ay1.legend(["Deceleration (g)", "Velocity (m/s)"],prop={'size': 8}) +ay1.title.set_text('Deceleration/Velocity') +ay1.grid(visible=True) + +ay2.plot(drop1['QSBC (kPa) K=0.2'],drop1["Penetration Depth (m)"]*100) +ay2.plot(drop1['QSBC (kPa) K=0.4'],drop1["Penetration Depth (m)"]*100) +ay2.plot(drop1['QSBC (kPa) K=1.0'],drop1["Penetration Depth (m)"]*100) +ay2.plot(drop1['QSBC (kPa) K=1.5'],drop1["Penetration Depth (m)"]*100) +ay2.set(xlabel="QSBC (kPa)") +ay2.grid(visible=True) +ay2.title.set_text('QSBC') +ay2.legend(["K=0.2", "K=0.4", "K=1.0", "K=1.5"],prop={'size': 8}) + +plt.suptitle(fileNum+' Drop 1', fontsize=10) +plt.show() + +title=fileNum+' Drop 1 Deceleration+Velocity Profile' +fig.savefig(title) + +#%% +#Deceleration and Velocity profile vs penetration depth for Drop 2 +fig, (ay1, ay2) = plt.subplots(1, 2, sharey=True) +ay1.plot(drop2g[acc1Nameg], drop2["Penetration Depth (m)"]*100, color = "k", linestyle = "solid") #marker = 11 +ay1.plot(drop2["Velocity (m/s)"], drop2["Penetration Depth (m)"]*100, color = "k", linestyle = "dashed") +ay1.set(xlabel="Deceleration (g) and Velocity (m/s)", ylabel="Penetration Depth (cm)") +ay1.invert_yaxis() +ay1.legend(["Deceleration (g)", "Velocity (m/s)"],prop={'size': 8}) +ay1.title.set_text('Deceleration/Velocity') +ay1.grid(visible=True) + +ay2.plot(drop2['QSBC (kPa) K=0.2'],drop2["Penetration Depth (m)"]*100) +ay2.plot(drop2['QSBC (kPa) K=0.4'],drop2["Penetration Depth (m)"]*100) +ay2.plot(drop2['QSBC (kPa) K=1.0'],drop2["Penetration Depth (m)"]*100) +ay2.plot(drop2['QSBC (kPa) K=1.5'],drop2["Penetration Depth (m)"]*100) +ay2.set(xlabel="QSBC (kPa)") +ay2.grid(visible=True) +ay2.title.set_text('QSBC') +ay2.legend(["K=0.2", "K=0.4", "K=1.0", "K=1.5"],prop={'size': 8}) + +plt.suptitle(fileNum+' Drop 2', fontsize=10) +plt.show() + +title=fileNum+' Drop 2 Deceleration+Velocity Profile' +fig.savefig(title) + +#%% +#Deceleration and Velocity profile vs penetration depth for Drop 3 +fig, (ay1, ay2) = plt.subplots(1, 2, sharey=True) +ay1.plot(drop3g[acc1Nameg], drop3["Penetration Depth (m)"]*100, color = "k", linestyle = "solid") #marker = 11 +ay1.plot(drop3["Velocity (m/s)"], drop3["Penetration Depth (m)"]*100, color = "k", linestyle = "dashed") +ay1.set(xlabel="Deceleration (g) and Velocity (m/s)", ylabel="Penetration Depth (cm)") +ay1.invert_yaxis() +ay1.legend(["Deceleration (g)", "Velocity (m/s)"],prop={'size': 8}) +ay1.title.set_text('Deceleration/Velocity') +ay1.grid(visible=True) + +ay2.plot(drop3['QSBC (kPa) K=0.2'],drop3["Penetration Depth (m)"]*100) +ay2.plot(drop3['QSBC (kPa) K=0.4'],drop3["Penetration Depth (m)"]*100) +ay2.plot(drop3['QSBC (kPa) K=1.0'],drop3["Penetration Depth (m)"]*100) +ay2.plot(drop3['QSBC (kPa) K=1.5'],drop3["Penetration Depth (m)"]*100) +ay2.set(xlabel="QSBC (kPa)") +ay2.grid(visible=True) +ay2.title.set_text('QSBC') +ay2.legend(["K=0.2", "K=0.4", "K=1.0", "K=1.5"],prop={'size': 8}) + +plt.suptitle(fileNum+' Drop 3', fontsize=10) +plt.show() + +title=fileNum+' Drop 3 Deceleration+Velocity Profile' +fig.savefig(title) + + +#%% +#Deceleration and Velocity profile vs penetration depth for Drop 4 +fig, (ay1, ay2) = plt.subplots(1, 2, sharey=True) +ay1.plot(drop4g[acc1Nameg], drop4["Penetration Depth (m)"]*100, color = "k", linestyle = "solid") #marker = 11 +ay1.plot(drop4["Velocity (m/s)"], drop4["Penetration Depth (m)"]*100, color = "k", linestyle = "dashed") +ay1.set(xlabel="Deceleration (g) and Velocity (m/s)", ylabel="Penetration Depth (cm)") +ay1.invert_yaxis() +ay1.legend(["Deceleration (g)", "Velocity (m/s)"],prop={'size': 8}) +ay1.title.set_text('Deceleration/Velocity') +ay1.grid(visible=True) + +ay2.plot(drop4['QSBC (kPa) K=0.2'],drop4["Penetration Depth (m)"]*100) +ay2.plot(drop4['QSBC (kPa) K=0.4'],drop4["Penetration Depth (m)"]*100) +ay2.plot(drop4['QSBC (kPa) K=1.0'],drop4["Penetration Depth (m)"]*100) +ay2.plot(drop4['QSBC (kPa) K=1.5'],drop4["Penetration Depth (m)"]*100) +ay2.set(xlabel="QSBC (kPa)") +ay2.grid(visible=True) +ay2.title.set_text('QSBC') +ay2.legend(["K=0.2", "K=0.4", "K=1.0", "K=1.5"],prop={'size': 8}) + +plt.suptitle(fileNum+' Drop 4', fontsize=10) +plt.show() + +title=fileNum+' Drop 4 Deceleration+Velocity Profile' +fig.savefig(title) + +#%% +#Deceleration and Velocity profile vs penetration depth for Drop 5 +fig, (ay1, ay2) = plt.subplots(1, 2, sharey=True) +ay1.plot(drop5g[acc1Nameg], drop5["Penetration Depth (m)"]*100, color = "k", linestyle = "solid") #marker = 11 +ay1.plot(drop5["Velocity (m/s)"], drop5["Penetration Depth (m)"]*100, color = "k", linestyle = "dashed") +ay1.set(xlabel="Deceleration (g) and Velocity (m/s)", ylabel="Penetration Depth (cm)") +ay1.invert_yaxis() +ay1.legend(["Deceleration (g)", "Velocity (m/s)"],prop={'size': 8}) +ay1.title.set_text('Deceleration/Velocity') +ay1.grid(visible=True) + +ay2.plot(drop5['QSBC (kPa) K=0.2'],drop5["Penetration Depth (m)"]*100) +ay2.plot(drop5['QSBC (kPa) K=0.4'],drop5["Penetration Depth (m)"]*100) +ay2.plot(drop5['QSBC (kPa) K=1.0'],drop5["Penetration Depth (m)"]*100) +ay2.plot(drop5['QSBC (kPa) K=1.5'],drop5["Penetration Depth (m)"]*100) +ay2.set(xlabel="QSBC (kPa)") +ay2.grid(visible=True) +ay2.title.set_text('QSBC') +ay2.legend(["K=0.2", "K=0.4", "K=1.0", "K=1.5"],prop={'size': 8}) + +plt.suptitle(fileNum+' Drop 5', fontsize=10) +plt.show() + +title=fileNum+' Drop 5 Deceleration+Velocity Profile' +fig.savefig(title) + +#%% +#Deceleration and Velocity profile vs penetration depth for Drop 6 +fig, (ay1, ay2) = plt.subplots(1, 2, sharey=True) +ay1.plot(drop6g[acc1Nameg], drop6["Penetration Depth (m)"]*100, color = "k", linestyle = "solid") #marker = 11 +ay1.plot(drop6["Velocity (m/s)"], drop6["Penetration Depth (m)"]*100, color = "k", linestyle = "dashed") +ay1.set(xlabel="Deceleration (g) and Velocity (m/s)", ylabel="Penetration Depth (cm)") +ay1.invert_yaxis() +ay1.legend(["Deceleration (g)", "Velocity (m/s)"],prop={'size': 8}) +ay1.title.set_text('Deceleration/Velocity') +ay1.grid(visible=True) + +ay2.plot(drop6['QSBC (kPa) K=0.2'],drop6["Penetration Depth (m)"]*100) +ay2.plot(drop6['QSBC (kPa) K=0.4'],drop6["Penetration Depth (m)"]*100) +ay2.plot(drop6['QSBC (kPa) K=1.0'],drop6["Penetration Depth (m)"]*100) +ay2.plot(drop6['QSBC (kPa) K=1.5'],drop6["Penetration Depth (m)"]*100) +ay2.set(xlabel="QSBC (kPa)") +ay2.grid(visible=True) +ay2.title.set_text('QSBC') +ay2.legend(["K=0.2", "K=0.4", "K=1.0", "K=1.5"],prop={'size': 8}) + +plt.suptitle(fileNum+' Drop 6', fontsize=10) +plt.show() + +title=fileNum+' Drop 6 Deceleration+Velocity Profile' +fig.savefig(title) + + +#%% +#Deceleration and Velocity profile vs penetration depth for Drop 7 +fig, (ay1, ay2) = plt.subplots(1, 2, sharey=True) +ay1.plot(drop7g[acc1Nameg], drop7["Penetration Depth (m)"]*100, color = "k", linestyle = "solid") #marker = 11 +ay1.plot(drop7["Velocity (m/s)"], drop7["Penetration Depth (m)"]*100, color = "k", linestyle = "dashed") +ay1.set(xlabel="Deceleration (g) and Velocity (m/s)", ylabel="Penetration Depth (cm)") +ay1.invert_yaxis() +ay1.legend(["Deceleration (g)", "Velocity (m/s)"],prop={'size': 8}) +ay1.title.set_text('Deceleration/Velocity') +ay1.grid(visible=True) + +ay2.plot(drop7['QSBC (kPa) K=0.2'],drop7["Penetration Depth (m)"]*100) +ay2.plot(drop7['QSBC (kPa) K=0.4'],drop7["Penetration Depth (m)"]*100) +ay2.plot(drop7['QSBC (kPa) K=1.0'],drop7["Penetration Depth (m)"]*100) +ay2.plot(drop7['QSBC (kPa) K=1.5'],drop7["Penetration Depth (m)"]*100) +ay2.set(xlabel="QSBC (kPa)") +ay2.grid(visible=True) +ay2.title.set_text('QSBC') +ay2.legend(["K=0.2", "K=0.4", "K=1.0", "K=1.5"],prop={'size': 8}) + +plt.suptitle(fileNum+' Drop 7', fontsize=10) +plt.show() + +title=fileNum+' Drop 7 Deceleration+Velocity Profile' +fig.savefig(title) + + +xvalues=np.linspace(0,len(dg)-1,len(dg)) + +temp_dg=dg['250g (g)'] +test=np.zeros(len(temp_dg)) diff --git a/BDcodev4-jupyter.ipynb b/BDcodev4-jupyter.ipynb new file mode 100644 index 0000000..2ebb23f --- /dev/null +++ b/BDcodev4-jupyter.ipynb @@ -0,0 +1,613 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 2, + "id": "34b268be", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "import numpy\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib\n", + "matplotlib.use('TkAgg')\n", + "import matplotlib.pyplot as plt\n", + "from pathlib import Path\n", + "\n", + "import scipy.integrate\n", + "from scipy.signal import find_peaks\n", + "from numpy import trapz\n", + "from scipy.integrate import cumtrapz" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "456bff57", + "metadata": {}, + "outputs": [], + "source": [ + "# SETUP VARIABLES - USER INPUTS\n", + "BD = 3\n", + "atype = 'p' # m = mantle area, p = projected area\n", + "tiptype = 'c' # c = cone, p = parabolic, b = blunt\n", + "# paste the filepath to the folder where the bd data is stored\n", + "binFilepath = Path(\"H:\\My Drive\\CEE 5904 - Project & Report\\FRF Data/test data\")\n", + "# write the bin file number you want to analyze (do not include 'bLog' or '.bin')\n", + "fileNum = '02F4'\n", + "outputFile = 'data.xlsx' #this currently doesn't do anything, but eventually all data will be printed out into an excel sheet\n", + "outputPath = Path(\"H:\\My Drive\\CEE 5904 - Project & Report\\FRF Data/test data\") # Path for new files\n", + "offset = 1 # this value is subtracted from the accelerometer readings" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "3909e0dc", + "metadata": {}, + "outputs": [], + "source": [ + "# SETUP FUNCTIONS\n", + "\n", + "#def tellme(s):\n", + " # #print(s)\n", + " # ax1.set(title=s)\n", + " # plt.draw()\n", + "\n", + "def accPick(dg, d):\n", + " # each of the following are the same; if an accelerometer breaks on a BD, can edit that section\n", + " # the 200g accelerometer is ignored in all cases\n", + " maxAcc = dg[\"250g (g)\"].max()\n", + " global acc\n", + " global accName\n", + " global accNameg\n", + " if maxAcc < 5 - offset:\n", + " if dg[\"2g (g)\"].max() < 1.8 - offset: # does an extra check for the 2g because of noise\n", + " acc = d[\"2g (m/s^2)\"]\n", + " accName = \"2g (m/s^2)\"\n", + " accNameg = \"2g (g)\"\n", + " else:\n", + " acc = d[\"18g (m/s^2)\"]\n", + " accName = \"18g (m/s^2)\"\n", + " accNameg = \"18g (g)\"\n", + " elif maxAcc < 18 - offset:\n", + " acc = d[\"18g (m/s^2)\"]\n", + " accName = \"18g (m/s^2)\"\n", + " accNameg = \"18g (g)\"\n", + " elif maxAcc < 50 - offset:\n", + " acc = d[\"50g (m/s^2)\"]\n", + " accName = \"50g (m/s^2)\"\n", + " accNameg = \"50g (g)\"\n", + " else:\n", + " acc = d[\"250g (m/s^2)\"]\n", + " accName = \"250g (m/s^2)\"\n", + " accNameg = \"250g (g)\"\n", + "\n", + "'''def peakpicker(dnum, dse):\n", + " global acc\n", + " #global dnum\n", + " #global dse\n", + " pts = plt.ginput(2)\n", + " (x0, y0), (x1, y1) = pts\n", + " xmin, xmax = sorted([x0, x1])\n", + " dropTimes[dnum, dse] = xmin\n", + " xmin = round(xmin)\n", + " print(xmin)\n", + " dse = dse + 1\n", + " dropTimes[dnum, dse] = xmax\n", + " xmax = round(xmax)\n", + " print(xmax)\n", + " print(dropTimes)\n", + " d = dfCal[xmin:xmax]\n", + " acc = d[accName]\n", + " print(d)'''\n", + "\n", + "def integration(d):\n", + " global vel\n", + " global maxvel\n", + " global dep\n", + " global maxdep\n", + " vel = scipy.integrate.cumtrapz(acc, x=d[\"Time (s)\"]) # solves for velocity\n", + " vel = np.array(vel)\n", + " vel = numpy.insert(vel, 0, 0) #not sure what to insert here, but it makes it the right size\n", + " vel = np.flip(vel)\n", + " maxvel = vel.max()\n", + " dep = scipy.integrate.cumtrapz(vel, x=d[\"Time (s)\"]) # solves for penetration depth\n", + " dep = numpy.insert(dep, 0, 0) # not sure what to insert here, but it makes it the right size\n", + " maxdep = dep.max()\n", + " d.insert(9, \"Velocity (m/s)\", vel)\n", + " d.insert(10, \"Penetration Depth (m)\", dep)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "32343298", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "H:\\My Drive\\CEE 5904 - Project & Report\\FRF Data\\test data\\bLog02F4.bin\n", + "(120000, 10)\n" + ] + } + ], + "source": [ + "# READ BD DATA IN\n", + "\n", + "data_array = [] # creates an empty array for us to fill with bd data\n", + "fileName = 'bLog'+fileNum+\".bin\"\n", + "# print(fileName)\n", + "newPath = binFilepath / fileName\n", + "print(newPath)\n", + "file = open(newPath, 'rb') # read file\n", + "element = file.read(3) # create a byte list with each element having 3 bytes\n", + "\n", + "while element:\n", + " # Convert to signed integer before adding to data array\n", + " iVAl = int.from_bytes(element, byteorder='big', signed=True)\n", + " data_array.append(iVAl) # adds the reshaped data from the bd file to the data frame\n", + " element = file.read(3)\n", + "\n", + "np_array = np.array(data_array) # create numpy array from the list\n", + "np_array = np.reshape(np_array, (-1, 10)) # convert the 1d array to 2d array with 10 cols\n", + "\n", + "print(np_array.shape)\n", + "# print(np_array)\n", + "\n", + "df = pd.DataFrame(np_array) # Creates a Dataframe in pandas from the bd data\n", + "df.columns = ['Count', 'no clue', 'g2g', 'g18g', 'g50g', 'ppm', 'g200g', 'gX55g', 'gY55g', 'g250g'] # names columns\n", + "# print(dfCal)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "c9050842", + "metadata": {}, + "outputs": [], + "source": [ + "# APPLY CALIBRATION FACTORS\n", + "\n", + "if BD == 3: # calibration factors from July 2019\n", + " g2g = (df['g2g']-34426.5)/1615925.8 - offset# accelerometers are in g\n", + " g18g = (df['g18g']+12322.1)/163530.7 - offset\n", + " g50g = (df['g50g']-237384.9)/63651 - 0.1120 - offset\n", + " ppm = ((df['ppm']+62496.7)/20583.0) * 6.89475729 # converts to kPa\n", + " g200g = ((df['g200g'] -248943.7)/39009.4)+0.5518 - offset\n", + " gX55g = (df['gX55g']-59093.7)/66674.3 - offset #check if lateral accelerometers also need to be offset\n", + " gY55g = (df['gY55g']-140224.6)/66674.3 - offset\n", + " g250g = (df['g250g']-40536.1)/13631.6 - offset\n", + "\n", + "if BD == 2: # calibration factors from Aug 26, 2021\n", + " g2g = (df['g2g']+31384.7)/1624987.2-0.035 - offset# accelerometers are in g\n", + " g18g = (df['g18g']-26631.0)/159945.4 - offset\n", + " g50g = (df['g50g']+92987.0)/63783.5 - offset\n", + " ppm = ((df['ppm']-35170.6)/12922.9) * 6.89475729 # converts to kPa\n", + " g200g = (df['g200g']-16264.8)/26042.8 -0.277 - offset\n", + " gX55g = (df['gX55g']+89890.3)/63897.1 - offset\n", + " gY55g = (df['gY55g']+14993.0)/64118.0 - offset\n", + " g250g = (df['g250g']+17362.1)/13533.5+0.0656 - offset\n", + "\n", + "if BD == 1: # calibration factors from July 2020\n", + " g2g = (df['g2g']+277743.2)/1637299.6 - offset # accelerometers are in g\n", + " g18g = (df['g18g']-3755.9)/159932.2 - offset\n", + " g50g = (df['g50g']+92817.6)/63237.1 - offset\n", + " ppm = ((df['ppm']-33154.0)/14763.5) * 6.89475729 # this is kPa\n", + " g200g = (df['g200g'] -1155309.9)/28368.5 - 1.464 - offset\n", + " gX55g = (df['gX55g'] +97138.4)/62023.7 - offset\n", + " gY55g = (df['gY55g']-9921.7)/62669.2 - offset\n", + " g250g = (df['g250g']+59211.3)/13276.9 - offset\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "2b9d0290", + "metadata": {}, + "outputs": [], + "source": [ + "#CREATE CALIBRATED DATAFRAMES\n", + "\n", + "time = (df['Count']-df['Count'].iloc[0]+1)/2000 # gives time in s\n", + "count = df[\"Count\"]\n", + "\n", + "# make a new dataframe of the calibrated values in units of g\n", + "dfCalg = pd.DataFrame([time, g2g, g18g, g50g, g200g, g250g, gX55g, gY55g, ppm])\n", + "dfCalg = dfCalg.T\n", + "dfCalg.columns = ['Time (s)', '2g (g)', '18g (g)', '50g (g)', '200g (g)', '250g (g)', 'X55g (g)', 'Y55g (g)', 'Pore Pressure (kPa)'] # names columns\n", + "#print(dfCalg)\n", + "\n", + "#make a new dataframe of the calibrated values in units of m/s^2\n", + "dfCal = pd.DataFrame([time, g2g, g18g, g50g, g200g, g250g, gX55g, gY55g, ppm])\n", + "dfCal = dfCal.T\n", + "dfCal.columns = ['Time (s)', '2g (m/s^2)', '18g (m/s^2)', '50g (m/s^2)', '200g (m/s^2)', '250g (m/s^2)', 'X55g (m/s^2)', 'Y55g (m/s^2)', 'Pore Pressure (kPa)'] # names columns\n", + "dfCal['2g (m/s^2)'] = dfCal['2g (m/s^2)'] * 9.80665\n", + "dfCal['18g (m/s^2)'] = dfCal['18g (m/s^2)'] * 9.80665\n", + "dfCal['50g (m/s^2)'] = dfCal['50g (m/s^2)'] * 9.80665\n", + "dfCal['200g (m/s^2)'] = dfCal['200g (m/s^2)'] * 9.80665\n", + "dfCal['250g (m/s^2)'] = dfCal['250g (m/s^2)'] * 9.80665\n", + "dfCal['X55g (m/s^2)'] = dfCal['X55g (m/s^2)'] * 9.80665\n", + "dfCal['Y55g (m/s^2)'] = dfCal['Y55g (m/s^2)'] * 9.80665\n", + "#print(dfCal)" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "6072d29f", + "metadata": {}, + "outputs": [ + { + "ename": "NameError", + "evalue": "name 'np' is not defined", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mNameError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn [1], line 2\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[38;5;66;03m#AUTOMATIC PEAK FINDING\u001b[39;00m\n\u001b[1;32m----> 2\u001b[0m x \u001b[38;5;241m=\u001b[39m \u001b[43mnp\u001b[49m\u001b[38;5;241m.\u001b[39marray(g250g) \u001b[38;5;66;03m# what accelerometer to get the peaks from\u001b[39;00m\n\u001b[0;32m 3\u001b[0m peaks, _ \u001b[38;5;241m=\u001b[39m find_peaks(x, height \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m2\u001b[39m, distance\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m10000\u001b[39m) \u001b[38;5;66;03m# finds the largest peaks more than 5g spaced at least 10000 counts apart\u001b[39;00m\n\u001b[0;32m 5\u001b[0m \u001b[38;5;66;03m# CREATION OF INDIVIDUAL DROP FILES\u001b[39;00m\n", + "\u001b[1;31mNameError\u001b[0m: name 'np' is not defined" + ] + } + ], + "source": [ + "#AUTOMATIC PEAK FINDING\n", + "x = np.array(g250g) # what accelerometer to get the peaks from\n", + "peaks, _ = find_peaks(x, height = 2, distance=10000) # finds the largest peaks more than 5g spaced at least 10000 counts apart\n", + "\n", + "# CREATION OF INDIVIDUAL DROP FILES\n", + "peaksArray = np.array(peaks) # prints a list of the count where the peaks occur\n", + "print(peaksArray)\n", + "q = (peaksArray.shape) #gives number of peaks\n", + "nDrops = int(q[0]) #number of drops in the file\n", + "print(nDrops)\n", + "\n", + "a = 0\n", + "n = 1\n", + "\n", + "dropTimes = np.array([[0]*2]*nDrops)\n", + "dnum = 0 #row index for drop; 0 = drop 1, etc\n", + "dse = 0 #column index for drop 0 = start, 1 = end\n", + "\n", + " # MAKE INDIVIDUAL DATAFRAMES FOR EACH DROP AND PLOT DECELERATION VS TIME\n", + "while n <= nDrops :\n", + " b = int(peaksArray[a]) # count at the ath drop\n", + " dropstart = b - 100 # offset in counts before impact to include in sliced df\n", + " dropend = b + 100 # offset in counts after impact to include in sliced df\n", + "\n", + " if n == 1 :\n", + " drop1 = dfCal[dropstart:dropend] #start and end points of the drop in m/s^2\n", + " drop1g = dfCalg[dropstart:dropend] #start and end points of the drop in g\n", + " drop1 = pd.DataFrame(drop1) # makes dataframe including all data within the start and end points of the drop\n", + " drop1g = pd.DataFrame(drop1g)\n", + " dg = drop1g # chooses what accelerometer to use based on the max g\n", + " d = drop1\n", + " accPick(dg, d) # chooses what accelerometer to use\n", + " acc1 = acc\n", + " acc1Name = accName\n", + " acc1Nameg = accNameg\n", + " drop1.plot(y=accName, ylabel=\"Deceleration (g)\", title='Select start and end of drop #' + str(n))\n", + " #plt.show()\n", + " #plt.set(title = 'Select start and end of drop #' + str(n))\n", + " #peakpicker(dnum, dse)\n", + " pts = plt.ginput(2)\n", + " (x0, y0), (x1, y1) = pts\n", + " xmin, xmax = sorted([x0, x1])\n", + " dropTimes[dnum, dse] = xmin\n", + " xmin = round(xmin)\n", + " print(xmin)\n", + " dse = dse + 1\n", + " dropTimes[dnum, dse] = xmax\n", + " xmax = round(xmax)\n", + " print(xmax)\n", + " print(dropTimes)\n", + " d = dfCal[xmin:xmax]\n", + " acc = d[accName]\n", + " print(d)\n", + " integration(d)\n", + " drop1 = d\n", + " dnum = dnum + 1\n", + " dse = dse - 1\n", + "\n", + " if n == 2 :\n", + " drop2 = dfCal[dropstart:dropend]\n", + " drop2g = dfCalg[dropstart:dropend]\n", + " drop2 = pd.DataFrame(drop2)\n", + " drop2g = pd.DataFrame(drop2g)\n", + " dg = drop2g\n", + " d = drop2\n", + " accPick(dg, d) # chooses what accelerometer to use\n", + " acc2 = acc\n", + " acc2Name = accName\n", + " acc2Nameg = accNameg\n", + " drop2.plot(y=accName, ylabel=\"Deceleration (g)\", title='Select start and end of drop #' + str(n))\n", + " #plt.show()\n", + " pts = plt.ginput(2)\n", + " (x0, y0), (x1, y1) = pts\n", + " xmin, xmax = sorted([x0, x1])\n", + " dropTimes[dnum, dse] = xmin\n", + " xmin = round(xmin)\n", + " print(xmin)\n", + " dse = dse + 1\n", + " dropTimes[dnum, dse] = xmax\n", + " xmax = round(xmax)\n", + " print(xmax)\n", + " print(dropTimes)\n", + " d = dfCal[xmin:xmax]\n", + " acc = d[accName]\n", + " print(d)\n", + " drop2 = d\n", + " integration(d)\n", + " dnum = dnum + 1\n", + " dse = dse - 1\n", + "\n", + " if n == 3 :\n", + " drop3 = dfCal[dropstart:dropend] #start and end points of the drop in m/s^2\n", + " drop3g = dfCalg[dropstart:dropend] #start and end points of the drop in g\n", + " drop3 = pd.DataFrame(drop3) # makes dataframe including all data within the start and end points of the drop\n", + " drop3g = pd.DataFrame(drop3g)\n", + " dg = drop3g # chooses what accelerometer to use based on the max g\n", + " d = drop3\n", + " accPick(dg, d) # chooses what accelerometer to use\n", + " acc3 = acc\n", + " acc3Name = accName\n", + " acc3Nameg = accNameg\n", + " drop3.plot(y=accName, ylabel=\"Deceleration (g)\", title='Select start and end of drop #' + str(n))\n", + " pts = plt.ginput(2)\n", + " (x0, y0), (x1, y1) = pts\n", + " xmin, xmax = sorted([x0, x1])\n", + " dropTimes[dnum, dse] = xmin\n", + " xmin = round(xmin)\n", + " print(xmin)\n", + " dse = dse + 1\n", + " dropTimes[dnum, dse] = xmax\n", + " xmax = round(xmax)\n", + " print(xmax)\n", + " print(dropTimes)\n", + " d = dfCal[xmin:xmax]\n", + " acc = d[accName]\n", + " print(d)\n", + " #drop3.plot(x=\"Time (s)\", y=accName, ylabel=\"Deceleration (g)\")\n", + " #plt.show()\n", + " drop3 = d\n", + " integration(d)\n", + " dnum = dnum + 1\n", + " dse = dse - 1\n", + "\n", + " if n == 4 :\n", + " drop4 = dfCal[dropstart:dropend] #start and end points of the drop in m/s^2\n", + " drop4g = dfCalg[dropstart:dropend] #start and end points of the drop in g\n", + " drop4 = pd.DataFrame(drop4) # makes dataframe including all data within the start and end points of the drop\n", + " drop4g = pd.DataFrame(drop4g)\n", + " dg = drop4g # chooses what accelerometer to use based on the max g\n", + " d = drop4\n", + " accPick(dg, d) # chooses what accelerometer to use\n", + " acc4 = acc\n", + " #drop4.plot(x=\"Time (s)\", y=accName, ylabel=\"Deceleration (g)\")\n", + " #plt.show()\n", + " integration(d)\n", + "\n", + " if n == 5 :\n", + " drop5 = dfCal[dropstart:dropend] #start and end points of the drop in m/s^2\n", + " drop5g = dfCalg[dropstart:dropend] #start and end points of the drop in g\n", + " drop5 = pd.DataFrame(drop5) # makes dataframe including all data within the start and end points of the drop\n", + " drop5g = pd.DataFrame(drop5g)\n", + " dg = drop5g # chooses what accelerometer to use based on the max g\n", + " d = drop5\n", + " accPick(dg, d) # chooses what accelerometer to use\n", + " acc5 = acc\n", + " #drop5.plot(x=\"Time (s)\", y=accName, ylabel=\"Deceleration (g)\")\n", + " #plt.show()\n", + " integration(d)\n", + "\n", + " if n == 6 :\n", + " drop6 = dfCal[dropstart:dropend] #start and end points of the drop in m/s^2\n", + " drop6g = dfCalg[dropstart:dropend] #start and end points of the drop in g\n", + " drop6 = pd.DataFrame(drop6) # makes dataframe including all data within the start and end points of the drop\n", + " drop6g = pd.DataFrame(drop6g)\n", + " dg = drop6g # chooses what accelerometer to use based on the max g\n", + " d = drop6\n", + " accPick(dg, d) # chooses what accelerometer to use\n", + " acc6 = acc\n", + " #drop6.plot(x=\"Time (s)\", y=accName, ylabel=\"Deceleration (g)\")\n", + " #plt.show()\n", + " integration(d)\n", + "\n", + " n = n + 1\n", + " a = a + 1\n", + "\n", + "\n", + "\n", + "#Deceleration,Velocity,and penetration depth vs time plots\n", + "\n", + "\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "5c6e4e7f", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot showing peak deceleration\n", + "peakplot = plt.plot(x)\n", + "peakplot = plt.plot(peaks, x[peaks], \"x\")\n", + "plt.show()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "12f9524b", + "metadata": {}, + "outputs": [], + "source": [ + "#Plot showing all accelerometers and pore pressure readings\n", + "fig, (ax1, ax2) = plt.subplots(2)\n", + "\n", + "ax1.plot(time, g2g, label=\"2g\" )\n", + "ax1.plot(time, g18g, label=\"18g\")\n", + "ax1.plot(time, g50g, label=\"50g\")\n", + "#plt.plot(time, ppm)\n", + "#ax1.plot(time, g200g, label=\"200g\")\n", + "#plt.plot(time, gX55g, label=\"2g\")\n", + "#plt.plot(time, gY55g, label=\"2g\")\n", + "ax1.plot(time, g250g, label=\"250g\")\n", + "ax1.legend()\n", + "ax1.set(ylabel=\"Deceleration (g)\")\n", + "ax1.set(xlabel=\"Time (s)\")\n", + "\n", + "ax2.plot(time, ppm, label=\"Pore Pressure\")\n", + "ax2.set(ylabel=\"Pore Pressure (kPa)\")\n", + "ax2.set(xlabel=\"Time (s)\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "d668bd2f", + "metadata": {}, + "outputs": [], + "source": [ + "#Deceleration,Velocity,and penetration depth vs time plots - Drop #1\n", + "\n", + "fig, (ax1, ax2, ax3) = plt.subplots(3)\n", + "ax1.plot(drop1[\"Time (s)\"], drop1[acc1Name], marker = 11, color = \"k\")\n", + "ax1.set(ylabel=\"Deceleration (m/s^2)\", xlabel=\"Time(s)\",title=\"Drop 1\")\n", + "ax2.plot(drop1[\"Time (s)\"], drop1['Velocity (m/s)'], marker = 11, color = \"k\")\n", + "ax2.set(ylabel=\"Velocity (m/s)\", xlabel=\"Time(s)\")\n", + "ax3.plot(drop1[\"Time (s)\"], drop1[\"Penetration Depth (m)\"], marker = 11, color = \"k\")\n", + "ax3.set(ylabel=\"Penetration Depth (m)\", xlabel=\"Time(s)\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "6d4de205", + "metadata": {}, + "outputs": [], + "source": [ + "#Deceleration,Velocity,and penetration depth vs time plots - Drop #2\n", + "fig, (ax1, ax2, ax3) = plt.subplots(3)\n", + "ax1.plot(drop2[\"Time (s)\"], drop2[acc2Name], marker = 11, color = \"k\")\n", + "ax1.set(ylabel=\"Deceleration (m/s^2)\", xlabel=\"Time(s)\",title=\"Drop 2\")\n", + "ax2.plot(drop2[\"Time (s)\"], drop2['Velocity (m/s)'], marker = 11, color = \"k\")\n", + "ax2.set(ylabel=\"Velocity (m/s)\", xlabel=\"Time(s)\")\n", + "ax3.plot(drop2[\"Time (s)\"], drop2[\"Penetration Depth (m)\"], marker = 11, color = \"k\")\n", + "ax3.set(ylabel=\"Penetration Depth (m)\", xlabel=\"Time(s)\")\n", + "plt.show()\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "1aed80ae", + "metadata": {}, + "outputs": [], + "source": [ + "#Deceleration,Velocity,and penetration depth vs time plots - Drop #3\n", + "fig, (ax1, ax2, ax3) = plt.subplots(3)\n", + "ax1.plot(drop3[\"Time (s)\"], drop3[acc3Name], marker = 11, color = \"k\")\n", + "ax1.set(ylabel=\"Deceleration (m/s^2)\", xlabel=\"Time(s)\", title=\"Drop 3\")\n", + "ax2.plot(drop3[\"Time (s)\"], drop3['Velocity (m/s)'], marker = 11, color = \"k\")\n", + "ax2.set(ylabel=\"Velocity (m/s)\", xlabel=\"Time(s)\")\n", + "ax3.plot(drop3[\"Time (s)\"], drop3[\"Penetration Depth (m)\"], marker = 11, color = \"k\")\n", + "ax3.set(ylabel=\"Penetration Depth (m)\", xlabel=\"Time(s)\")\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "a318b707", + "metadata": {}, + "outputs": [], + "source": [ + "#Deceleration,Velocity,and penetration depth vs time plots - Drop #4\n" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "55e7ab96", + "metadata": {}, + "outputs": [], + "source": [ + "#Deceleration,Velocity,and penetration depth vs time plots - Drop #5\n" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "9243da4a", + "metadata": {}, + "outputs": [], + "source": [ + "#Deceleration,Velocity,and penetration depth vs time plots - Drop #6\n" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "857522b5", + "metadata": {}, + "outputs": [ + { + "ename": "ValueError", + "evalue": "x and y must have same first dimension, but have shapes (200,) and (46,)", + "output_type": "error", + "traceback": [ + "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m", + "\u001b[1;31mValueError\u001b[0m Traceback (most recent call last)", + "Cell \u001b[1;32mIn [17], line 4\u001b[0m\n\u001b[0;32m 1\u001b[0m \u001b[38;5;66;03m#Deceleration vs. Penetration Depth, Velocity vs. Penetration Depth\u001b[39;00m\n\u001b[0;32m 3\u001b[0m fig, (ax1) \u001b[38;5;241m=\u001b[39m plt\u001b[38;5;241m.\u001b[39msubplots(\u001b[38;5;241m1\u001b[39m)\n\u001b[1;32m----> 4\u001b[0m \u001b[43max1\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mplot\u001b[49m\u001b[43m(\u001b[49m\u001b[43mdrop3g\u001b[49m\u001b[43m[\u001b[49m\u001b[43macc3Nameg\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mdrop3\u001b[49m\u001b[43m[\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mPenetration Depth (m)\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m]\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m100\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmarker\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;241;43m11\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mcolor\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mk\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[0;32m 5\u001b[0m ax1\u001b[38;5;241m.\u001b[39mplot(drop3[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mVelocity (m/s)\u001b[39m\u001b[38;5;124m\"\u001b[39m], drop3[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mPenetration Depth (m)\u001b[39m\u001b[38;5;124m\"\u001b[39m]\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m100\u001b[39m, marker \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m11\u001b[39m, color \u001b[38;5;241m=\u001b[39m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mk\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m 6\u001b[0m ax1\u001b[38;5;241m.\u001b[39mset(xlabel\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mDeceleration (g) and Velocity (m/s)\u001b[39m\u001b[38;5;124m\"\u001b[39m, ylabel\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mPenetration Depth (cm)\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n", + "File \u001b[1;32m~\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\matplotlib\\axes\\_axes.py:1635\u001b[0m, in \u001b[0;36mAxes.plot\u001b[1;34m(self, scalex, scaley, data, *args, **kwargs)\u001b[0m\n\u001b[0;32m 1393\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m 1394\u001b[0m \u001b[38;5;124;03mPlot y versus x as lines and/or markers.\u001b[39;00m\n\u001b[0;32m 1395\u001b[0m \n\u001b[1;32m (...)\u001b[0m\n\u001b[0;32m 1632\u001b[0m \u001b[38;5;124;03m(``'green'``) or hex strings (``'#008000'``).\u001b[39;00m\n\u001b[0;32m 1633\u001b[0m \u001b[38;5;124;03m\"\"\"\u001b[39;00m\n\u001b[0;32m 1634\u001b[0m kwargs \u001b[38;5;241m=\u001b[39m cbook\u001b[38;5;241m.\u001b[39mnormalize_kwargs(kwargs, mlines\u001b[38;5;241m.\u001b[39mLine2D)\n\u001b[1;32m-> 1635\u001b[0m lines \u001b[38;5;241m=\u001b[39m [\u001b[38;5;241m*\u001b[39m\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_get_lines(\u001b[38;5;241m*\u001b[39margs, data\u001b[38;5;241m=\u001b[39mdata, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs)]\n\u001b[0;32m 1636\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m line \u001b[38;5;129;01min\u001b[39;00m lines:\n\u001b[0;32m 1637\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39madd_line(line)\n", + "File \u001b[1;32m~\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\matplotlib\\axes\\_base.py:312\u001b[0m, in \u001b[0;36m_process_plot_var_args.__call__\u001b[1;34m(self, data, *args, **kwargs)\u001b[0m\n\u001b[0;32m 310\u001b[0m this \u001b[38;5;241m+\u001b[39m\u001b[38;5;241m=\u001b[39m args[\u001b[38;5;241m0\u001b[39m],\n\u001b[0;32m 311\u001b[0m args \u001b[38;5;241m=\u001b[39m args[\u001b[38;5;241m1\u001b[39m:]\n\u001b[1;32m--> 312\u001b[0m \u001b[38;5;28;01myield from\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_plot_args\u001b[49m\u001b[43m(\u001b[49m\u001b[43mthis\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n", + "File \u001b[1;32m~\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\matplotlib\\axes\\_base.py:498\u001b[0m, in \u001b[0;36m_process_plot_var_args._plot_args\u001b[1;34m(self, tup, kwargs, return_kwargs)\u001b[0m\n\u001b[0;32m 495\u001b[0m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39maxes\u001b[38;5;241m.\u001b[39myaxis\u001b[38;5;241m.\u001b[39mupdate_units(y)\n\u001b[0;32m 497\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m x\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m0\u001b[39m] \u001b[38;5;241m!=\u001b[39m y\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m0\u001b[39m]:\n\u001b[1;32m--> 498\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mx and y must have same first dimension, but \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m 499\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mhave shapes \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mx\u001b[38;5;241m.\u001b[39mshape\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m and \u001b[39m\u001b[38;5;132;01m{\u001b[39;00my\u001b[38;5;241m.\u001b[39mshape\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n\u001b[0;32m 500\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m x\u001b[38;5;241m.\u001b[39mndim \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m2\u001b[39m \u001b[38;5;129;01mor\u001b[39;00m y\u001b[38;5;241m.\u001b[39mndim \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m2\u001b[39m:\n\u001b[0;32m 501\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mx and y can be no greater than 2D, but have \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[0;32m 502\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mshapes \u001b[39m\u001b[38;5;132;01m{\u001b[39;00mx\u001b[38;5;241m.\u001b[39mshape\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m and \u001b[39m\u001b[38;5;132;01m{\u001b[39;00my\u001b[38;5;241m.\u001b[39mshape\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m)\n", + "\u001b[1;31mValueError\u001b[0m: x and y must have same first dimension, but have shapes (200,) and (46,)" + ] + } + ], + "source": [ + "#Deceleration vs. Penetration Depth, Velocity vs. Penetration Depth\n", + "\n", + "fig, (ax1) = plt.subplots(1)\n", + "ax1.plot(drop3g[acc3Nameg], drop3[\"Penetration Depth (m)\"]*100, marker = 11, color = \"k\")\n", + "ax1.plot(drop3[\"Velocity (m/s)\"], drop3[\"Penetration Depth (m)\"]*100, marker = 11, color = \"k\")\n", + "ax1.set(xlabel=\"Deceleration (g) and Velocity (m/s)\", ylabel=\"Penetration Depth (cm)\")\n", + "ax1.invert_yaxis()\n", + "plt.show()\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.6" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Bootleg Dropview.py b/Bootleg Dropview.py new file mode 100644 index 0000000..613180d --- /dev/null +++ b/Bootleg Dropview.py @@ -0,0 +1,173 @@ +import numpy +import numpy as np +import pandas as pd +import matplotlib.pyplot as plt +from pathlib import Path +import math +from math import pi +from math import sqrt +import statistics +import scipy.integrate +from scipy.signal import find_peaks +from scipy.signal import argrelmin +from numpy import trapz +from scipy.integrate import cumtrapz + +# SETUP VARIABLES - USER INPUTS +BD = 2 #Bluedrop file is from +fileNum = '02B7' # write the bin file number you want to analyze (do not include 'bLog' or '.bin') +soiltype = "s" #s = sand, c=clay, m=mixed, u=unknown +atype = 'p' # m = mantle area (best for sands), p = projected area (best for clays) +tiptype = 'c' # c = cone, p = parabolic, b = blunt +offset = 1 # this value is subtracted from the accelerometer readings +droptype = 'w' #w = water, #a = air +sign = "uk" #enter an effective unit weight value in kg/m^3 or "uk" if unknown +# paste the filepath to the folder where the BD data is stored +binFilepath = Path("H:\My Drive\CEE 5904 - Project & Report/2021 FRF Data\Bluedrop\October/14 October 2021 AM\Pier - BD2") +#paste the filepath that you want the files containing drops copied to +outputPath = Path("H:\My Drive\CEE 5904 - Project & Report/2021 FRF Data\Bluedrop\October/14 October 2021 AM\Analysis Results-14 October 2021 AM.xlsx") # Path to pre-existing Excel File +plotPath = Path("H:\My Drive\CEE 5904 - Project & Report/2021 FRF Data\Bluedrop\October/14 October 2021 AM\Analysis Figures- 14 October 2021 AM") +#if applicable, paste the filepath to an excel file that troubleshooting data will be printed in +troubleshootingPath = Path("H:\My Drive\CEE 5904 - Project & Report/2021 FRF Data\Bluedrop\October/14 October 2021 AM\Troubleshooting.xlsx") + +def overviewplot(): #Plot showing all accellerometers and pore pressure readings + fig, (ax1, ax2, ax3) = plt.subplots(3) + fig.set_size_inches(14, 7) + ax1.plot(time, g250g, label="250g") + ax1.plot(time, g50g, label="50g") + ax1.plot(time, g18g, label="18g") + ax1.plot(time, g2g, label="2g") + ax1.legend(loc = "upper right") + ax1.set(ylabel="Deceleration (g)") + ax1.set_title("BD file "+fileNum) + ax1.grid() + + ax2.plot(time, ppm, label="Pore Pressure") + ax2.set(ylabel="Pore Pressure (kPa)") + ax2.grid() + + ax3.plot(time, gX55g, label="X 55g") + ax3.plot(time, gY55g, label="Y 55g") + ax3.legend(loc = "upper right") + ax3.set(ylabel="Deceleration (g)") + ax3.set(xlabel="Time (s)") + ax3.grid() + + fig.subplots_adjust(bottom=.1, left = .1) + #plotName = fileNum+" Overview.png" + #plt.savefig(plotPath / plotName) + plt.show() + +# READ BD DATA IN +data_array = [] # creates an empty array for us to fill with bd data +fileName = 'bLog'+fileNum+".bin" +# print(fileName) +newPath = binFilepath / fileName +#print(newPath) +file = open(newPath, 'rb') # read file +element = file.read(3) # create a byte list with each element having 3 bytes + +while element: + # Convert to signed integer before adding to data array + iVAl = int.from_bytes(element, byteorder='big', signed=True) + data_array.append(iVAl) # adds the reshaped data from the bd file to the data frame + element = file.read(3) + +np_array = np.array(data_array) # create numpy array from the list +np_array = np.reshape(np_array, (-1, 10)) # convert the 1d array to 2d array with 10 cols + +#print(np_array.shape) +# print(np_array) +''' +df = pd.DataFrame(np_array) # Creates a Dataframe in pandas from the bd data +df.columns = ['Count', 'no clue', 'g2g', 'g18g', 'g50g', 'ppm', 'g200g', 'gX55g', 'gY55g', 'g250g'] # names columns +# print(dfCal) + +# APPLY CALIBRATION FACTORS +if BD == 3: # calibration factors from July 2019 + g2g = (df['g2g']-38285.6)/1615800.9 - offset# accelerometers are in g + g18g = (df['g18g']+13738)/163516.8 - offset + g50g = (df['g50g']-238520.6)/63666 - offset + ppm = ((df['ppm']-139040.1)/20705) * 6.89475729 # converts to kPa + g200g = ((df['g200g'] +12142.6)/27751.9) - offset + gX55g = (df['gX55g']-90237)/65351.5 + gY55g = (df['gY55g']-57464.2)/65545. + g250g = (df['g250g']-40420.3)/13636.9 - offset + +if BD == 2: # calibration factors from Aug 26, 2021 + g2g = (df['g2g']+37242.2)/1639250.2 - offset# accelerometers are in g + g18g = (df['g18g']-26867.0)/160460.5 - offset + g50g = (df['g50g']-213923.3)/64080.7- offset + ppm = ((df['ppm']+55518.9)/18981.7) * 6.89475729 # converts to kPa + g200g = (df['g200g']-171448.6)/30334.2 - offset + gX55g = (df['gX55g']-54242.6)/64767.7 + gY55g = (df['gY55g']-40574.2)/66343.1 + g250g = (df['g250g']-40614.9)/13654.6 - offset + +if BD == 1: # calibration factors from July 2020 + g2g = (df['g2g']-42590.9)/1626361.1 - offset # accelerometers are in g + g18g = (df['g18g']-44492.9)/161125.5 - offset + g50g = (df['g50g']-171656.1)/64020.3 - offset + ppm = ((df['ppm']+31776.1)/20679.7) * 6.89475729 # this is kPa + g200g = (df['g200g'] -723404.8)/32209.7 - offset + gX55g = (df['gX55g'] -54881.1)/64858.6 + gY55g = (df['gY55g']-28735.5)/63839.9 + g250g = (df['g250g']+13299.7)/13697.1 - offset + +time = (df['Count']-df['Count'].iloc[0]+1)/2000 # gives time in s +count = df["Count"]-df['Count'].iloc[0] + +# make a new dataframe of the calibrated values in units of g +dfCalg = pd.DataFrame([count, time, g2g, g18g, g50g, g200g, g250g, gX55g, gY55g, ppm]) +dfCalg = dfCalg.T +dfCalg.columns = ['Count', 'Time (s)', '2g (g)', '18g (g)', '50g (g)', '200g (g)', '250g (g)', 'X55g (g)', 'Y55g (g)', 'Pore Pressure (kPa)'] # names columns +#print(dfCalg) + +#make a new dataframe of the calibrated values in units of m/s^2 +dfCal = pd.DataFrame([count, time, g2g, g18g, g50g, g200g, g250g, gX55g, gY55g, ppm]) +dfCal = dfCal.T +dfCal.columns = ['Count','Time (s)', '2g (m/s^2)', '18g (m/s^2)', '50g (m/s^2)', '200g (m/s^2)', '250g (m/s^2)', 'X55g (m/s^2)', 'Y55g (m/s^2)', 'Pore Pressure (kPa)'] # names columns +dfCal['2g (m/s^2)'] = dfCal['2g (m/s^2)'] * 9.80665 +dfCal['18g (m/s^2)'] = dfCal['18g (m/s^2)'] * 9.80665 +dfCal['50g (m/s^2)'] = dfCal['50g (m/s^2)'] * 9.80665 +dfCal['200g (m/s^2)'] = dfCal['200g (m/s^2)'] * 9.80665 +dfCal['250g (m/s^2)'] = dfCal['250g (m/s^2)'] * 9.80665 +dfCal['X55g (m/s^2)'] = dfCal['X55g (m/s^2)'] * 9.80665 +dfCal['Y55g (m/s^2)'] = dfCal['Y55g (m/s^2)'] * 9.80665 +#print(dfCal) + +overviewplot() + +print("Does the file contain a drop? (y/n)") +drop_input=input() + if drop_input=='y': + files = Path(directory).glob('*') + for file in files: + print(file) +''' +name = [*fileNum] +print(name) +print(name[3]) + +def newName(): + a = int(name[3]) + print(a) + if a <=8: + b = a + 1 + b = str(b) + elif a == 9: + b = "A" + elif a == "A": + b = "B" + elif a == "B": + b = "C" + elif a == "C": + b = "D" + elif a == "D": + b = "E" + + newName = name + newName[3] = b + print(newName) + +newName() diff --git a/Partially Saturated Sands BC.py b/Partially Saturated Sands BC.py new file mode 100644 index 0000000..7e7f7a2 --- /dev/null +++ b/Partially Saturated Sands BC.py @@ -0,0 +1,79 @@ +import math + +#USER INPUTS + +moist = 180.4 #moisture in mV +gammad = 15.3 #dry unit weight in kN/m^3 +d60 = .413 # D60 in mm +cu = 1.69 # coefficient of uniformity (unitless) + +emin = 0.468 #minimum void ratio +emax = 0.753 #maximum void ratio +gammadmin = 15.196 #min dry unit weight (kN/m^3) +gammadmax = 18.143 #max dry unit weight (kN/m^3) +gammasea = 10.05 #unit weight of seawater (kN/m^3) +Gs = 2.65 #specific gravity + +#weight volume calculations +def wtvol(): + global S + global gammabulk + global n + global Dr + V = moist/1000 #moisture in volts + theta = 2.832*V**4 - 3.6426*V**3 + 1.3985*V**2 + 0.4112*V - 0.0149 #Calibrated volumetric water content + e = Gs*gammasea/gammad - 1 #void ratio + n = e/(1+e) #porosity + Dr = (emax-e)/(emax-emin) #relative density (%) + S = theta/n*100 #degree of saturation (%) + w = (S*e)/Gs #gravitational water content (%) + gammabulk = gammad*(1+w/100) #bulk unit weight (kN/m^3) + +#matric suction calc +def matric_suction(): + global ms + thetar = 2.1 #residual volumetric water content (%) + c1 = 1.07 + c2 = 12.07 + sal = 30 #salinity (g/kg) + t = 20 #temperature (deg. C) + st1 = 0.073 #surface tension (N/m) + st2 = st1*(1+(3.766*10**-4)*sal+(2.347*10**-6)*sal*t) #surface tension (N/m) + Sr = thetar/n #residual degress of saturation (%) + coefn = (c1/math.log10(cu))+1 + coefa = (c2*st2/(d60/1000))/1000 #kPa + Se = (S-Sr)/(1-Sr/100) #effective degree of saturation (%) + ms = coefa*((Se/100)**(coefn/(1-coefn)) - 1)**(1/coefn) #matric suction (kPa) + +def Duncan():#Duncan Correlation (phi') + global phi + A = 34 + B = 10 + C = 3 + D = 2 + sigN= .05 * gammabulk #normal stress for duncan correlation, assuming D = .05m (kPa) + pa = 1.03 #atmospheric pressure, kPa + phi = A + B* Dr - (C + D*Dr)*math.log10(sigN/pa) + +def Vanapalli(): #Bearing Capacity (Vanapalli & Mohammed) + c = 0 #effective cohesion, kPa + psi = 1 + D = 0.05 + B = 0.0875 + Kp = math.tan(math.radians(45+phi/2))**2 + Nq = math.e**(math.pi*math.tan(math.radians(phi)))*Kp + Nc = (Nq-1)*(1/math.tan(math.radians(phi))) + Ngam = 2*(Nq+1)*math.tan(math.radians(phi)) + Sc = 1+Nq/Nc + Sgam = 0.6 + qu = (c + ms*(S/100**psi)*math.tan(math.radians(phi)))*Nc*Sc+0.5*gammabulk*B*Ngam*Sgam #partially saturated bearing capacity, kPa + print(qu) + + +wtvol() +matric_suction() +Duncan() +Vanapalli() + + +