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() + + +