BD code - single drop analysis.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "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",
    "import math\n",
    "from math import pi\n",
    "from math import sqrt\n",
    "import statistics\n",
    "import scipy.integrate\n",
    "from scipy.signal import find_peaks\n",
    "from scipy.signal import argrelmin\n",
    "from numpy import trapz\n",
    "from scipy.integrate import cumtrapz"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "user inputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SETUP VARIABLES - USER INPUTS\n",
    "BD = 2 #Bluedrop file is from \n",
    "fileNum = '0021' # write the bin file number you want to analyze (do not include 'bLog' or '.bin')\n",
    "soiltype = \"s\" #s = sand, c=clay, m=mixed, u=unknown\n",
    "atype = 'm'  # m = mantle area (best for sands), p = projected area (best for clays)\n",
    "tiptype = 'c'  # c = cone, p = parabolic, b = blunt\n",
    "#offset = 1 # this value is subtracted from the accelerometer readings\n",
    "droptype = 'w' #w = water, #a = air\n",
    "sign = \"uk\" #enter an effective unit weight value in kg/m^3 or \"uk\" if unknown \n",
    "# paste the filepath to the folder where the BD data is stored\n",
    "binFilepath = Path(\"H:/.shortcut-targets-by-id/1aF9t2aiRGWTftJMZFAOBixqvQniFBjnb/Duck  2023/Data/Intertidal/BlueDrop, Samples & Moisture Gage/07March23/BD2 3.7.23 - Drops Only\")\n",
    "#paste the filepath to an excel file that the analysis results will be printed in\n",
    "outputPath = Path(\"H:/.shortcut-targets-by-id/1aF9t2aiRGWTftJMZFAOBixqvQniFBjnb/Duck  2023/Data/Intertidal/BlueDrop, Samples & Moisture Gage/07March23/BlueDrop Data Processing March 7 2023.xlsx\") #  Path to pre-existing Excel File\n",
    "plotPath = Path(\"H:/.shortcut-targets-by-id/1aF9t2aiRGWTftJMZFAOBixqvQniFBjnb/Duck  2023/Data/Intertidal/BlueDrop, Samples & Moisture Gage/07March23/BD2 3.8.23 - Drops Only/Analysis Figures\")\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "set functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# FUNCTIONS FOR USE IN ANALYSIS\n",
    "def masslength(tiptype): #sets the mass and length of the pentrometer based on the tip\n",
    "    global mass\n",
    "    global length\n",
    "    if tiptype == 'c':\n",
    "        mass = 7.71 #kg\n",
    "        length = 8.48-0.93 #originally 7.87, 7.57 for perfect 60 deg, 8.48 measured - .93 measured 90 deg\n",
    "    elif tiptype == 'p':\n",
    "        mass = 9.15\n",
    "        length = 8.26\n",
    "    elif tiptype == 'b':\n",
    "        mass = 10.30\n",
    "        length = 8.57\n",
    "\n",
    "def dropstartend(peak): #after locating the peaks, this function chops the minute long file into a smaller segment immediately before and after the peak\n",
    "    global dropstart\n",
    "    global dropend\n",
    "    if peak <= 1500:\n",
    "        dropstart = 1\n",
    "        dropend = peak + 500\n",
    "    elif peak > 119500:\n",
    "        dropstart = peak - 1500\n",
    "        dropend = 120000\n",
    "    else:\n",
    "        dropstart = peak - 1500\n",
    "        dropend = peak + 500\n",
    "\n",
    "def accPick(dg, d): #this function picks the smallest accelerometer that's not maxed out to perform the integration on\n",
    "    maxAcc = dg[\"250g (g)\"].max()\n",
    "    global acc\n",
    "    global accName\n",
    "    global accNameg\n",
    "    global accg\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 = dfCal[\"2g (m/s^2)\"]\n",
    "            accg = dfCalg[\"2g (g)\"]\n",
    "            accName = \"2g (m/s^2)\"\n",
    "            accNameg = \"2g (g)\"\n",
    "        else:\n",
    "            acc = dfCal[\"18g (m/s^2)\"]\n",
    "            accg = dfCalg[\"18g (g)\"]\n",
    "            accName = \"18g (m/s^2)\"\n",
    "            accNameg = \"18g (g)\"\n",
    "    elif maxAcc < 18: #- offset\n",
    "        acc = dfCal[\"18g (m/s^2)\"]\n",
    "        accg = dfCalg[\"18g (g)\"]\n",
    "        accName = \"18g (m/s^2)\"\n",
    "        accNameg = \"18g (g)\"\n",
    "    elif maxAcc < 50: #- offset:\n",
    "        acc = dfCal[\"50g (m/s^2)\"]\n",
    "        accg = dfCalg[\"50g (g)\"]\n",
    "        accName = \"50g (m/s^2)\"\n",
    "        accNameg = \"50g (g)\"\n",
    "    else:\n",
    "        acc = dfCal[\"250g (m/s^2)\"]\n",
    "        accg = dfCalg[\"50g (g)\"]\n",
    "        accName = \"250g (m/s^2)\"\n",
    "        accNameg = \"250g (g)\"\n",
    "\n",
    "def findchangepts(): #This function picks the moment that the BD impacts the ground (not in use currently)\n",
    "    global drop\n",
    "    jlist = list()\n",
    "    global jindex\n",
    "    print(\"finding start of drop...\")\n",
    "    for i in range(4,len(accg)-4):\n",
    "        p1 = 1\n",
    "        #print(p1)\n",
    "        p2 = i\n",
    "        #print(p2)\n",
    "        p3 = len(accg)\n",
    "        #print(p3)\n",
    "        sample1 = list(accg[p1:p2-1])\n",
    "        #print(sample1)\n",
    "        sample2 = list(accg[p2:p3])\n",
    "        #print(sample2)\n",
    "        stat1 = math.log(statistics.variance(sample1))\n",
    "        stat2 = math.log(statistics.variance(sample2))\n",
    "        #print(stat1)\n",
    "\n",
    "        j1 = (i-1)*stat1\n",
    "        j2 = ((len(accg)-1)-i+1)*stat2\n",
    "        j = j1+j2\n",
    "        #print(j)\n",
    "        jlist.append(j)\n",
    "\n",
    "    drop = min(jlist)\n",
    "    #print(\"drop is\")\n",
    "    #print(drop)\n",
    "    jlist = np.array(jlist)\n",
    "    #print(jlist)\n",
    "    #print(jlist.size)\n",
    "    jlist = np.insert(jlist, 0, (0,0,0,0)) #reshape to match up with dataframe d\n",
    "    jlist = np.append(jlist, (0,0,0,0)) #reshape to match up with dataframe d\n",
    "    #print(jlist.size) #should be 2000\n",
    "    jindex = np.where(jlist==drop) #finds the index of the drop start\n",
    "    jindex = int(jindex[0]) #converts the index into a number from a tuple\n",
    "    print(\"Jlist: \", jlist)\n",
    "\n",
    "def finddropend(n): #finds the location after the peak where the deceleration is 1 \n",
    "    global penend\n",
    "    below0list = list()\n",
    "    #for i in range(dropstart+jindex, dropend, 1):\n",
    "    for i in range(peak, dropend, 1):\n",
    "        if accg[i] < 1: #the offset is applied later in the code\n",
    "            penend = i \n",
    "            below0list = np.append(below0list, penend)\n",
    "    penend=int(min(below0list))\n",
    "\n",
    "def integration(d, accint): #integrates the deceleration data to solve for velocity and penetration depth\n",
    "    global vel\n",
    "    global maxvel\n",
    "    global dep\n",
    "    global maxdep\n",
    "    #accint = acc[penstart:penend]\n",
    "    #print(len(d))\n",
    "    #print(len(accint))\n",
    "    vel = scipy.integrate.cumtrapz(accint, 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",
    "    maxvel = round(maxvel,1)\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",
    "    maxdep = round(maxdep,4)\n",
    "    d.insert(9, \"Velocity (m/s)\", vel)\n",
    "    d.insert(10, \"Penetration Depth (m)\", dep)\n",
    "\n",
    "def peakpick():\n",
    "    global penstart\n",
    "    fig, (ax1) = plt.subplots(1)\n",
    "    plt.plot(dfCalg['Count'], dfCalg[accNameg], label=accNameg)\n",
    "    plt.plot(dfCalg['Count'], dfCalg['2g (g)'], label='2g')\n",
    "\n",
    "    ax1.set_xlim(left=dropstart, right=dropend)\n",
    "    ax1.legend()\n",
    "    ax1.set(ylabel=\"Deceleration (g)\")\n",
    "    ax1.set(xlabel=\"Time\")\n",
    "    ax1.set_title('Double-click start of drop #' + str(a))\n",
    "    startendpt = plt.ginput(2, 0)\n",
    "    pentimes = []\n",
    "    for t in startendpt:\n",
    "        for x in t:\n",
    "            pentimes.append(x)\n",
    "\n",
    "    penstart = int(pentimes[0])\n",
    "    print(\"start of penetration: \", penstart)\n",
    "\n",
    "def CheckingFunction(): # (to be used with auto peak finding) This function is added to check the values of the start and end of the drop if the autofinding function doesn't work\n",
    "\n",
    "    global jindex\n",
    "    global num1\n",
    "    \n",
    "    xvalues=np.linspace(0,len(dg)-1,len(dg))\n",
    "    plt.figure(num=1)\n",
    "    temp_dg=dg['250g (g)']\n",
    "    plt.plot(xvalues,temp_dg)\n",
    "    plt.grid(visible=True,which='both')\n",
    "    plt.show()\n",
    "    print ('The start of the drop is', jindex)\n",
    "    print ('The end of the drop is', num1)\n",
    "    print('Are the start and end of the drop correct?(y/n)')\n",
    "    drop_input=input()\n",
    "    if drop_input=='n':   \n",
    "        peakpick()\n",
    "        jindex = penstart\n",
    "        num1 = penend\n",
    "\n",
    "def areafind(): #finds the embedded area based on the penetration depth, area type, and the tip\n",
    "    global area\n",
    "    global trunc_index\n",
    "    a1 = list() #placeholder for the penetrated area at that specific depth\n",
    "    d1 = dep*100 #penetration depth array, in cm\n",
    "    #print(len(d1))\n",
    "    r = list() #placeholder for the radius at that specific depth\n",
    "    if tiptype == 'c':\n",
    "        if atype == 'm':\n",
    "            for k in range(0,len(d1)):\n",
    "                if d1[k]<length:\n",
    "                    r.append(.22+d1[k]*((4.375-.22)/7.55))\n",
    "                    #r.append(d1[k]*np.tan(30*pi/180)) this is from the original matlabcode\n",
    "                    a1.append(pi*r[k]*(sqrt((r[k]**2)+(d1[k]**2))))\n",
    "                elif d1[k]>=length:\n",
    "                    r.append(4.375)\n",
    "                    a1.append(pi*r[k]*(sqrt((r[k]**2)+(length**2))))\n",
    "                a1[k] = a1[k]/10000\n",
    "            area = a1\n",
    "            #print(r)\n",
    "            print(np.column_stack((d1,r, area)))\n",
    "        elif atype == 'p':\n",
    "            for k in range(0,len(d1)):\n",
    "                if d1[k]<length:\n",
    "                    r.append(.22+d1[k]*((4.375-.22)/7.55))\n",
    "                    # r.append(d1[k]*np.tan(30*pi/180))\n",
    "                    a1.append(pi*r[k]**2)\n",
    "                elif d1[k]>=length:\n",
    "                    r.append(4.375)\n",
    "                    if r[k-1]<r[k]:\n",
    "                        trunc_index = r.index(r[k-1])+jindex+dropstart #performs analysis only on readings beyond depth of cone height\n",
    "                    a1.append(pi*(r[k])**2)\n",
    "                a1[k] = a1[k]/10000\n",
    "            area = a1\n",
    "\n",
    "    elif tiptype == 'b':\n",
    "        if atype =='m':\n",
    "            for k in range(0,len(d1)):\n",
    "                if d1[k]<length:\n",
    "                    r.append(4.375)\n",
    "                    a1.append(pi*r[k]**2 + 2*pi*r[k]*d1[k])\n",
    "                if d1(k)>=length:\n",
    "                    r.append(4.375)\n",
    "                    a1.append(pi*r[k]**2 + 2*pi*r[k]*length)\n",
    "                a1[k]=a1[k]/10000\n",
    "            area = a1\n",
    "        elif atype == 'p':\n",
    "            for k in range(0,len(d1)):\n",
    "                a1.append(pi*4.375**2)\n",
    "                a1[k]=a1[k]/10000\n",
    "            area = a1\n",
    "\n",
    "    '''elif tiptype == \"p\":\n",
    "        if atype == 'm':\n",
    "            for k in range(1,len(d)):\n",
    "                if d1[k]<length:\n",
    "                     r[k]=sqrt(2.4184*d1[k])\n",
    "                           polarfun = @(theta,r) r.*sqrt(0.745*r.^2 + 1);\n",
    "                           A1(k)= integral2(polarfun,0,2*pi,0,r(k));\n",
    "            \n",
    "        elif atype == 'p':\n",
    "            for k in range(1,len(d)):'''\n",
    "\n",
    "def qdynfun(acc): # calculates dynamic bearing capacity\n",
    "    global qdyn\n",
    "    global maxbcdep\n",
    "    global maxqdyn\n",
    "    global qdyntable\n",
    "\n",
    "    buoy = 1020*0.002473\n",
    "    if droptype == \"w\": #water drops\n",
    "        Fbe = (mass-buoy)*acc #drop force\n",
    "\n",
    "    elif droptype ==\"a\": #air drops\n",
    "        Fbe = mass*acc\n",
    "\n",
    "    qdyn = (Fbe/area)/1000 #Dynamic bearing capacity (kPa)\n",
    "    qdyntable = pd.DataFrame()\n",
    "    qdyntable[\"Qdyn\"] = qdyn\n",
    "    qdyntable[\"Depth\"] = dep\n",
    "    qdyntablet = qdyntable#.truncate(before=trunc_index)\n",
    "    qdynt = qdyntablet[\"Qdyn\"]\n",
    "    #print(qdyntable)\n",
    "\n",
    "    maxqdyn = qdynt.max() #gives maximum qdyn\n",
    "    maxbcdep = qdyntablet.loc[qdyntable[\"Qdyn\"]==maxqdyn, 'Depth'].item()\n",
    "\n",
    "    maxqdyn = round(maxqdyn,1)\n",
    "    maxbcdep = round(maxbcdep,4)\n",
    "    #print(\"maxbcdep = \", maxbcdep) #finds the penetration depth corresponding to the maximum qdyn\n",
    "\n",
    "\n",
    "    \"\"\"bcmin = argrelmin(qdyn)\n",
    "    peaksArray = np.array(bcmin)\n",
    "    print(\"Peaks Array = \", peaksArray)\"\"\"\n",
    "\n",
    "def qsbcfun(acc): #calculates quasi-static bearing capacity\n",
    "    global bctable\n",
    "    global qsbc\n",
    "    global maxHqsbc\n",
    "    global maxLqsbc\n",
    "    global maxAqsbc\n",
    "\n",
    "    if droptype == \"w\":\n",
    "        srcv = np.log10(vel/0.02) #Velocity portion of the strain rate correction. \n",
    "        srfk = [0.2, 0.4, 1, 1.5] #list of strain rate factors to run; must be in increasing order\n",
    "\n",
    "        bctable = pd.DataFrame()\n",
    "        for i in range(0,len(srfk)):\n",
    "            fsr = 1+srfk[i]*srcv\n",
    "            qsbc = qdyn/fsr\n",
    "            bctable.insert(i, \"qsbc for srf = \"+str(srfk[i]), qsbc)\n",
    "\n",
    "        bctable_avgs = pd.DataFrame()\n",
    "        bctable_avgs[\"High Average qsbc\"]= (bctable.iloc[:,0]+bctable.iloc[:,1])/2 #Average of the lowest 2 strain rate factors\n",
    "        bctable_avgs[\"Low Average qsbc\"] = (bctable.iloc[:,int(len(srfk))-1]+bctable.iloc[:,int(len(srfk)-2)])/2 #average of the highest 2 strain rate factors\n",
    "        maxHqsbc = bctable_avgs[\"High Average qsbc\"].max()\n",
    "        maxHqsbc = round(maxHqsbc, 1)\n",
    "        maxLqsbc = bctable_avgs[\"Low Average qsbc\"].max()\n",
    "        maxLqsbc = round(maxLqsbc,1)\n",
    "        maxAqsbc = \"-\"\n",
    "\n",
    "        #print(bctable)\n",
    "\n",
    "    elif droptype == \"a\":\n",
    "        srcv = np.log10(vel/0.02) #Velocity portion of the strain rate correction. \n",
    "        srfk = 0.31*maxvel\n",
    "        fsr = 1+srfk*srcv\n",
    "        qsbc = qdyn/fsr\n",
    "        maxHqsbc = \"-\"\n",
    "        maxLqsbc = \"-\"\n",
    "        maxAqsbc = qsbc.max()\n",
    "        maxAqsbc = round(maxAqsbc,1)\n",
    "\n",
    "def dr(accg): #Albatal 2019 method for calculating relative density (Duck, NC specific!)\n",
    "    global Dr\n",
    "    maxacc = accg.max() #maximum deceleration, in g\n",
    "    #print(\"Max deceleration = \", maxacc)\n",
    "    Dr = -2.18*10**-4*maxacc**3+1.29*10**-2*maxacc**2+1.61*maxacc-13.09 #Albatal 2019\n",
    "    Dr = round(Dr, 1)\n",
    "    #print(\"Relative density is \", dr)\n",
    "\n",
    "def ignore_drop():\n",
    "    global keep_drop\n",
    "    print('Keep Drop? (y/n)')\n",
    "    keep_drop=input()\n",
    "\n",
    "def duncan_correlation(): #duncan correlation for caluclating friction angle\n",
    "    global phi\n",
    "    global sign\n",
    "    #coefficients A, B, C, and D used in the Duncan Correlation. Duncan[34,10,3,2], Albatal[34,10,2,5]\n",
    "    A = 34\n",
    "    B = 10\n",
    "    C = 3\n",
    "    D = 2\n",
    "    pa = 101.325 #kPa\n",
    "    print(\"max penetration depth = \", maxdep)\n",
    "    if sign == \"uk\":\n",
    "        sign = 1120*maxdep #1120 = 70pcf in kg/m^3\n",
    "\n",
    "    phi = A + B*Dr/100-(C+D*Dr/100)*math.log(sign/pa)\n",
    "    phi = round(phi, 1)\n",
    "    print(\"phi =\", phi)\n",
    "\n",
    "def firmnessfactor(acc, drop):\n",
    "    global maxacc\n",
    "    global ff\n",
    "    maxacc = acc.max()\n",
    "    g = 9.80665\n",
    "    tp = len(drop)/2000\n",
    "\n",
    "    ff = maxacc/(g*tp*maxvel)\n",
    "    ff = round(ff,1)\n",
    "    print(\"Firmness factor = \", ff)\n",
    "\n",
    "def bctruncate(dropg, drop, acc, qdyn):\n",
    "    dropg = pd.DataFrame(dropg)\n",
    "    drop = pd.DataFrame(drop)\n",
    "    qdyn = pd.DataFrame(qdyn)\n",
    "    acc = pd.DataFrame(acc)\n",
    "    dropg = dropg.truncate(before=trunc_index)\n",
    "    drop = drop.truncate(before=trunc_index)\n",
    "    acc = acc.truncate(before=trunc_index)\n",
    "    qdyn = qdyn.truncate(before=trunc_index)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plots for exporting\n",
    "def overviewplot(): #Plot showing all accellerometers and pore pressure readings\n",
    "    fig, (ax1, ax2, ax3) = plt.subplots(3)\n",
    "    fig.set_size_inches(14, 7)\n",
    "    #plt.tight_layout()\n",
    "    ax1.plot(time, g250g, label=\"250g\")\n",
    "    ax1.plot(time, g50g, label=\"50g\")\n",
    "    ax1.plot(time, g18g, label=\"18g\")\n",
    "    ax1.plot(time, g2g, label=\"2g\")\n",
    "    ax1.legend(loc = \"upper right\")\n",
    "    ax1.set(ylabel=\"Deceleration (g)\")\n",
    "    ax1.set_title(\"BD file \"+fileNum)\n",
    "    ax1.grid()\n",
    "\n",
    "    ax2.plot(time, ppm, label=\"Pore Pressure\")\n",
    "    ax2.set(ylabel=\"Pore Pressure (kPa)\")\n",
    "    ax2.grid()\n",
    "\n",
    "    ax3.plot(time, gX55g, label=\"X 55g\")\n",
    "    ax3.plot(time, gY55g, label=\"Y 55g\")\n",
    "    ax3.legend(loc = \"upper right\")\n",
    "    ax3.set(ylabel=\"Deceleration (g)\")\n",
    "    ax3.set(xlabel=\"Time (s)\")\n",
    "    ax3.grid()\n",
    "\n",
    "    fig.subplots_adjust(bottom=.1, left = .1)\n",
    "    plotName = fileNum+\" Overview.png\"\n",
    "    #plt.savefig(plotPath / plotName)\n",
    "    plt.show()\n",
    "\n",
    "def integplot(drop, accName): #Deceleration,Velocity,and penetration depth vs time plots\n",
    "    fig, (ax1, ax2, ax3) = plt.subplots(3)\n",
    "    fig.set_size_inches(7, 7)\n",
    "    #resize to 7,7, 7,14\n",
    "    #plt.tight_layout()\n",
    "    timereset = drop[\"Time (s)\"]- drop[\"Time (s)\"].iloc[0]\n",
    "    ax1.plot(timereset, drop[accName], color = \"k\" ) #,marker = 11\n",
    "    ax1.set(ylabel=\"Deceleration (m/s^2)\")\n",
    "    ax1.set_title(\"BD file \"+fileNum+ \" Drop #\"+str(a))\n",
    "    ax1.grid()\n",
    "\n",
    "    ax2.plot(timereset, drop['Velocity (m/s)'] , color = \"k\" )\n",
    "    ax2.set(ylabel=\"Velocity (m/s)\")\n",
    "    ax2.grid()\n",
    "\n",
    "    ax3.plot(timereset, drop[\"Penetration Depth (m)\"] , color = \"k\")\n",
    "    ax3.set(ylabel=\"Penetration Depth (m)\", xlabel=\"Time(s)\")\n",
    "    ax3.grid()\n",
    "\n",
    "    fig.subplots_adjust(bottom=.1, left = .1)\n",
    "    plotName = fileNum+\" \"+str(a)+\" Integ\"\n",
    "    #plt.savefig(plotPath / plotName)\n",
    "    plt.show()\n",
    "\n",
    "def dep_qsbc_comboplot(dropg, drop, accNameg, qdyn):\n",
    "    fig, (ax1, ax2) = plt.subplots(1,2)\n",
    "    fig.set_size_inches(7, 7)\n",
    "    ax1.plot(dropg[accNameg], drop[\"Penetration Depth (m)\"]*100, color = \"k\", linestyle = \"solid\") #marker = 11\n",
    "    ax1.plot(drop[\"Velocity (m/s)\"], drop[\"Penetration Depth (m)\"]*100, color = \"k\", linestyle = \"dashed\")\n",
    "    ax1.plot(drop[\"X55g (m/s^2)\"], drop[\"Penetration Depth (m)\"]*100, color = \"b\", linestyle = \"dashed\")\n",
    "    ax1.plot(drop[\"Y55g (m/s^2)\"], drop[\"Penetration Depth (m)\"]*100, color = \"r\", linestyle = \"dashed\")\n",
    "    ax1.set(xlabel=\"Deceleration (g) and Velocity (m/s)\", ylabel=\"Penetration Depth (cm)\")\n",
    "    ax1.invert_yaxis()\n",
    "    ax1.legend([\"Deceleration (g)\", \"Velocity (m/s)\", \"X tilt (g)\", \"Y tilt (g)\"], loc = \"upper right\")\n",
    "    ax1.set_title(\"BD file \"+fileNum)\n",
    "    ax1.grid()\n",
    "\n",
    "    if droptype == \"a\": #for partial saturation, only one qsbc curve\n",
    "        ax2.plot(qdyn, drop[\"Penetration Depth (m)\"]*100, label=\"Qdyn\") #color = \"k\", marker = 11 #Plots Qdyn\n",
    "        ax2.plot(qsbc, drop[\"Penetration Depth (m)\"]*100, label=\"QSBC\") #plots QSBC\n",
    "        ax2.set(xlabel=\"Bearing Capacity (kPa)\")\n",
    "        ax2.set_xlim(0,800)\n",
    "        ax2.invert_yaxis()\n",
    "        ax2.legend([\"Qdyn\", \"QSBC\"], loc = \"upper right\")\n",
    "        ax2.set_title(\"Drop #\"+str(a))\n",
    "        ax2.grid()\n",
    "\n",
    "    else: #for saturated drops, use a range of strain-rate factors\n",
    "        ax2.plot(qdyn, drop[\"Penetration Depth (m)\"]*100, label=\"Qdyn\") #color = \"k\", marker = 11 #Plots Qdyn\n",
    "        ax2.plot(bctable.iloc[:,0], drop[\"Penetration Depth (m)\"]*100, label=str(bctable.columns[0])) #plots QSBC with smallest srf\n",
    "        ax2.plot(bctable.iloc[:,len(bctable.columns)-1], drop[\"Penetration Depth (m)\"]*100, label=str(bctable.columns[len(bctable.columns)-1])) #plots QSBC with largest srf\n",
    "        for i in range(1,len(bctable.columns)-1): #plots qsbc for all other srfs\n",
    "            ax2.plot(bctable.iloc[:,i], drop[\"Penetration Depth (m)\"]*100, label=str(bctable.columns[i]), color = '.6')\n",
    "        ax2.set(xlabel=\"Bearing Capacity (kPa)\")\n",
    "        ax2.set_xlim(0,)\n",
    "        ax2.invert_yaxis()\n",
    "        ax2.legend([\"Qdyn\", str(bctable.columns[0]), str(bctable.columns[len(bctable.columns)-1])], loc = \"upper right\")\n",
    "        ax2.set_title(\"Drop #\"+str(a))\n",
    "        ax2.grid()\n",
    "    \n",
    "    fig.subplots_adjust(bottom=.1, left = .1)\n",
    "    \n",
    "\n",
    "\n",
    "\n",
    "#plots for troubleshooting \n",
    "def depthplot(dropg, drop, accNameg): #Velocity and develeration vs. penetration depth \n",
    "    fig, (ax1) = plt.subplots(1)\n",
    "    ax1.plot(dropg[accNameg], drop[\"Penetration Depth (m)\"]*100, color = \"k\", linestyle = \"solid\") #marker = 11\n",
    "    ax1.plot(drop[\"Velocity (m/s)\"], drop[\"Penetration Depth (m)\"]*100, color = \"k\", linestyle = \"dashed\")\n",
    "    ax1.set(xlabel=\"Deceleration (g) and Velocity (m/s)\", ylabel=\"Penetration Depth (cm)\")\n",
    "    ax1.invert_yaxis()\n",
    "    ax1.legend([\"Deceleration (g)\", \"Velocity (m/s)\"])\n",
    "    ax1.set_title(\"BD file \"+fileNum)\n",
    "    plt.show()\n",
    "\n",
    "def peakplot(): # Plot showing peak deceleration\n",
    "    peakplot = plt.plot(x)\n",
    "    peakplot = plt.plot(peaks, x[peaks], \"x\")\n",
    "    plt.show()\n",
    "\n",
    "def qdynplot(drop, qdyn): #Plot showing dynamic bearing capacity vs depth\n",
    "    #print(len(drop))\n",
    "    #print(len(qdyn))\n",
    "    #print((bctable.iloc[:,0]))\n",
    "    #print(len(bctable.iloc[:,len(bctable.columns)-1]))\n",
    "    fig, (ax1) = plt.subplots(1)\n",
    "\n",
    "    if droptype == \"a\": #for partial saturation, only one qsbc curve\n",
    "        ax1.plot(qdyn, drop[\"Penetration Depth (m)\"]*100, label=\"Qdyn\") #color = \"k\", marker = 11 #Plots Qdyn\n",
    "        ax1.plot(qsbc, drop[\"Penetration Depth (m)\"]*100, label=\"QSBC\") #plots QSBC\n",
    "        ax1.set(xlabel=\"Bearing Capacity (kPa)\")\n",
    "        ax1.set_xlim(0,)\n",
    "        ax1.invert_yaxis()\n",
    "        ax1.legend([\"Qdyn\", \"QSBC\"])\n",
    "        ax1.set_title(\"Bearing Capacity- \"+fileNum+ \" \"+str(a))\n",
    "        ax1.grid()\n",
    "\n",
    "    else:\n",
    "        ax1.plot(qdyn, drop[\"Penetration Depth (m)\"]*100, label=\"Qdyn\") #color = \"k\", marker = 11 #Plots Qdyn\n",
    "        ax1.plot(bctable.iloc[:,0], drop[\"Penetration Depth (m)\"]*100, label=str(bctable.columns[0])) #plots QSBC with smallest srf\n",
    "        ax1.plot(bctable.iloc[:,len(bctable.columns)-1], drop[\"Penetration Depth (m)\"]*100, label=str(bctable.columns[len(bctable.columns)-1])) #plots QSBC with largest srf\n",
    "        for i in range(1,len(bctable.columns)-1): #plots qsbc for all other srfs\n",
    "            ax1.plot(bctable.iloc[:,i], drop[\"Penetration Depth (m)\"]*100, label=str(bctable.columns[i]), color = \"k\")\n",
    "        ax1.set(xlabel=\"Bearing Capacity (kPa)\", ylabel=\"Penetration Depth (cm)\")\n",
    "        ax1.set_xlim(0,)\n",
    "        ax1.invert_yaxis()\n",
    "        ax1.legend([\"Qdyn\", str(bctable.columns[0]), str(bctable.columns[len(bctable.columns)-1])])\n",
    "        ax1.set_title(\"Bearing Capacity- \"+fileNum+ \" \"+str(a))\n",
    "    \n",
    "    plt.show()\n",
    "\n",
    "def tiltplot():\n",
    "    fig, (ax1) = plt.subplots(1)\n",
    "    plt.plot(time, gX55g, label=\"2g\")\n",
    "    plt.plot(time, gY55g, label=\"2g\")\n",
    "    ax1.legend()\n",
    "    ax1.set(ylabel=\"Deceleration (g)\")\n",
    "    ax1.set(xlabel=\"Time (s)\")\n",
    "    ax1.set_title(\"BD file \"+fileNum)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib.font_manager import FontProperties\n",
    "font = FontProperties()\n",
    "font.set_name('Calibri')\n",
    "\n",
    "def duck_dep_qsbc_comboplot(dropg, drop, dropt, accNameg, qdynt):\n",
    "    font = {'fontname': 'Calibri'}\n",
    "    fig, (ax1, ax2) = plt.subplots(1,2)\n",
    "    fig.set_size_inches(7, 7)\n",
    "    ax1.plot(dropg[accNameg], drop[\"Penetration Depth (m)\"]*100, color = \"k\", linestyle = \"solid\") #marker = 11\n",
    "    ax1.plot(drop[\"Velocity (m/s)\"], drop[\"Penetration Depth (m)\"]*100, color = \"k\", linestyle = \"dashed\")\n",
    "    #ax1.plot(drop[\"X55g (m/s^2)\"], drop[\"Penetration Depth (m)\"]*100, color = \"b\", linestyle = \"dashed\")\n",
    "    #ax1.plot(drop[\"Y55g (m/s^2)\"], drop[\"Penetration Depth (m)\"]*100, color = \"r\", linestyle = \"dashed\")\n",
    "    ax1.set(xlabel=\"Deceleration (g) and Velocity (m/s)\", ylabel=\"Penetration Depth (cm)\", FontProperties = font)\n",
    "    ax1.set_xlim(left = 0)\n",
    "    ax1.set_ylim(0,dep.max()*100+.25)\n",
    "    ax1.invert_yaxis()\n",
    "    ax1.legend([\"Deceleration (g)\", \"Velocity (m/s)\"], FontProperties = font)\n",
    "    ax1.set_title(\"BD file \"+fileNum, FontProperties = font)\n",
    "    ax1.grid()\n",
    "\n",
    "    \n",
    "    if droptype == \"a\": #for partial saturation, only one qsbc curve\n",
    "        y = dropt[\"Penetration Depth (m)\"]*100\n",
    "        ax2.plot(qdynt[10:-50], y[10:-50], label=\"Qdyn\") #color = \"k\", marker = 11 #Plots Qdyn, chops off first 40 readings\n",
    "        ax2.plot(qsbc[10:-50], y[10:-50], label=\"QSBC\") #plots QSBC; chops off last 10 points which go off towards infinity (can adjust as needed)\n",
    "        ax2.set(xlabel=\"Bearing Capacity (kPa)\",)\n",
    "        ax2.set_xlim(0,qdyn.max()+10)\n",
    "        ax2.set_ylim(0,dep.max()*100+.25)\n",
    "        ax2.invert_yaxis()\n",
    "        ax2.legend([\"Qdyn\", \"QSBC\"], loc = \"upper right\")\n",
    "        ax2.set_title(\"Drop #\"+str(a))\n",
    "        ax2.grid()\n",
    "\n",
    "    else: #for saturated drops, use a range of strain-rate factors\n",
    "        y = list(dropt[\"Penetration Depth (m)\"]*100)\n",
    "        ax2.plot(qdynt[10:-20], y[10:-20], label=\"Qdyn\") #color = \"k\", marker = 11 #Plots Qdyn\n",
    "        ax2.plot(bctablet.iloc[10:-20,0], y[10:-20], label=str(bctablet.columns[0])) #plots QSBC with smallest srf\n",
    "        ax2.plot(bctablet.iloc[10:-20,len(bctable.columns)-1], y[10:-20], label=str(bctable.columns[len(bctable.columns)-1])) #plots QSBC with largest srf\n",
    "        for i in range(1,len(bctablet.columns)-1): #plots qsbc for all other srfs\n",
    "            ax2.plot(bctablet.iloc[10:-20,i], y[10:-20], label=str(bctablet.columns[i]), color = '.6')\n",
    "            #bctablet.iloc[:,i].pop()\n",
    "        ax2.set(xlabel=\"Bearing Capacity (kPa)\")\n",
    "        ax2.set_xlim(0,)\n",
    "        ax2.invert_yaxis()\n",
    "        ax2.legend([\"Qdyn\", str(bctablet.columns[0]), str(bctablet.columns[len(bctablet.columns)-1])], loc = \"upper right\")\n",
    "        ax2.set_title(\"Drop #\"+str(a))\n",
    "        ax2.grid()\n",
    "\n",
    "    #fig.subplots_adjust(bottom=.1, left = .1)\n",
    "    plotName = fileNum+\" \"+str(a)+\" BC\"\n",
    "    #plt.savefig(plotPath / plotName)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Exporting Functions\n",
    "def exporttoexcel():\n",
    "    with pd.ExcelWriter(outputPath, mode=\"a\", if_sheet_exists='replace') as writer:\n",
    "        output_table.to_excel(writer, sheet_name = fileNum, index=False)\n",
    "\n",
    "def troubleshooting_export():\n",
    "        with pd.ExcelWriter(troubleshootingPath, mode=\"a\", if_sheet_exists = 'new') as writer:\n",
    "            drop1.to_excel(writer, sheet_name = str(a), index = False)\n",
    "            qdyntable.to_excel(writer, sheet_name = \"qdyn\", index = False)\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Bring in data and calibrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#Set the penetrometer mass and length\n",
    "masslength(tiptype)\n",
    "\n",
    "# READ BD DATA IN\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",
    "\n",
    "# APPLY CALIBRATION FACTORS\n",
    "if BD == 3:  # calibration factors from March 2023\n",
    "    g2g = (df['g2g']-33570.1)/1614577.9 #- offset# accelerometers are in g\n",
    "    g18g = (df['g18g']+13495)/163387.2 #- offset\n",
    "    g50g = (df['g50g']-238817.4)/63779.5 #- offset\n",
    "    ppm = ((df['ppm']-139040.1)/20705) * 6.89475729 # converts to kPa\n",
    "    g200g = (df['g200g'] -262332.4)/38888.7 #- offset\n",
    "    gX55g = (df['gX55g']-70406)/59754.3\n",
    "    gY55g = (df['gY55g']-69421.1)/141871.5\n",
    "    g250g = (df['g250g']-39077.4)/13746.9 #- offset\n",
    "\n",
    "if BD == 2: # calibration factors from March 2023\n",
    "    g2g = (df['g2g']+36597.6)/1637627.3 #- offset# accelerometers are in g\n",
    "    g18g = (df['g18g']-26185.3)/160297.2 #- offset\n",
    "    g50g = (df['g50g']-212256.9)/63968.7 #- offset\n",
    "    ppm = ((df['ppm']+55518.9)/18981.7) * 6.89475729 # converts to kPa\n",
    "    g200g = (df['g200g']-175499.4)/30583.8 #- offset\n",
    "    gX55g = (df['gX55g']-53629.9)/68590.9\n",
    "    gY55g = (df['gY55g']-43694.3)/68280.3\n",
    "    g250g = (df['g250g']-39619.9)/13545.8 #- offset\n",
    "\n",
    "if BD == 1: # calibration factors from March 2023\n",
    "    g2g = (df['g2g']-43727.6)/1625064 #- offset # accelerometers are in g\n",
    "    g18g = (df['g18g']-45085.6)/160925.7 #- offset\n",
    "    g50g = (df['g50g']-173493.4)/63944.6 #- offset\n",
    "    ppm = ((df['ppm']+31776.1)/20679.7) * 6.89475729 # this is kPa\n",
    "    g200g = (df['g200g'] -731734.3)/32192.4  #- offset\n",
    "    gX55g = (df['gX55g'] -68837.9)/52137.3\n",
    "    gY55g = (df['gY55g']-68015.9)/28074.9\n",
    "    g250g = (df['g250g']+10580)/13688.7 #- offset\n",
    "\n",
    "time = (df['Count']-df['Count'].iloc[0]+1)/2000 # gives time in s\n",
    "count = df[\"Count\"]-df['Count'].iloc[0]\n",
    "\n",
    "# make a new dataframe of the calibrated values in units of g\n",
    "dfCalg = pd.DataFrame([count, time, g2g, g18g, g50g, g200g, g250g, gX55g, gY55g, ppm])\n",
    "dfCalg = dfCalg.T\n",
    "dfCalg.columns = ['Count', '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([count, time, g2g, g18g, g50g, g200g, g250g, gX55g, gY55g, ppm])\n",
    "dfCal = dfCal.T\n",
    "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\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)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "locate drops"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "#Locate the drops\n",
    "x = np.array(g250g)  # what accelerometer to get the peaks from - use 250 because it will never be maxed out\n",
    "max250 = max(g250g)\n",
    "\n",
    "peaks, _ = find_peaks(x, height = 2, distance=1000, prominence=3)  # finds the largest peaks more than 2g spaced at least 10000 counts apart\n",
    "\n",
    "peaksArray = np.array(peaks)  # prints a list of the count where the peaks occur\n",
    "#peaksArray = filter(filterpeak, peaksArray)\n",
    "\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",
    "# For each drop, find the start and end points and integrate to solve for velocity and acceleration\n",
    "#peakplot()\n",
    "\n",
    "\n",
    "\n",
    "peakplot()\n",
    "#overviewplot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak index: 28061\n"
     ]
    }
   ],
   "source": [
    "a = 1 #what drop do you want to analyze? USER INPUT\n",
    "a = a-1\n",
    "\n",
    "peak = int(peaksArray[a])  # count at the ath drop\n",
    "print(\"peak index:\" , peak)\n",
    "dropstartend(peak) #zooms in the drop file to only consider 500 counts before and 1500 counts after the peak deceleration\n",
    "#print(dropstart, dropend)\n",
    "\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 \n",
    "d = drop1\n",
    "accPick(dg, d)  # chooses what accelerometer to use\n",
    "acc1 = acc\n",
    "acc1g = accg\n",
    "acc1Name = accName\n",
    "acc1Nameg = accNameg\n",
    "#findchangepts() #automatically identifies the moment of penetration\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "start of penetration:  28024\n"
     ]
    }
   ],
   "source": [
    "\n",
    "peakpick()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "penend:  64989\n",
      "offset:  0.0\n",
      "[[0.00000000e+00 2.20000000e-01 1.52053084e-05]\n",
      " [1.47248589e-01 3.01035482e-01 3.16932046e-05]\n",
      " [2.93532681e-01 3.81540171e-01 5.77012654e-05]\n",
      " [4.38084895e-01 4.61091753e-01 9.21317591e-05]\n",
      " [5.80000723e-01 5.39192451e-01 1.34144303e-04]\n",
      " [7.18323948e-01 6.15316027e-01 1.82836669e-04]\n",
      " [8.52172394e-01 6.88976993e-01 2.37195082e-04]\n",
      " [9.80844224e-01 7.59789106e-01 2.96148468e-04]\n",
      " [1.10382817e+00 8.27471000e-01 3.58623365e-04]\n",
      " [1.22072637e+00 8.91803719e-01 4.23553725e-04]\n",
      " [1.33120009e+00 9.52600846e-01 4.89881597e-04]\n",
      " [1.43497379e+00 1.00971074e+00 5.56580655e-04]\n",
      " [1.53185183e+00 1.06302574e+00 6.22688185e-04]\n",
      " [1.62171550e+00 1.11248052e+00 6.87324157e-04]\n",
      " [1.70450408e+00 1.15804165e+00 7.49694332e-04]\n",
      " [1.78020420e+00 1.19970178e+00 8.09093652e-04]\n",
      " [1.84884215e+00 1.23747538e+00 8.64907878e-04]\n",
      " [1.91047506e+00 1.27139389e+00 9.16611196e-04]\n",
      " [1.96519393e+00 1.30150739e+00 9.63771454e-04]\n",
      " [2.01313002e+00 1.32788811e+00 1.00605745e-03]\n",
      " [2.05445293e+00 1.35062939e+00 1.04323796e-03]\n",
      " [2.08935689e+00 1.36983813e+00 1.07516827e-03]\n",
      " [2.11804690e+00 1.38562713e+00 1.10177422e-03]\n",
      " [2.14074175e+00 1.39811682e+00 1.12305074e-03]\n",
      " [2.15768381e+00 1.40744056e+00 1.13906656e-03]\n",
      " [2.16914172e+00 1.41374621e+00 1.14996230e-03]\n",
      " [2.17541393e+00 1.41719800e+00 1.15594872e-03]\n",
      " [2.17682859e+00 1.41797653e+00 1.15730107e-03]\n",
      " [2.17373135e+00 1.41627202e+00 1.15434128e-03]\n",
      " [2.16647849e+00 1.41228055e+00 1.14742511e-03]\n",
      " [2.15544204e+00 1.40620684e+00 1.13694085e-03]\n",
      " [2.14101961e+00 1.39826973e+00 1.12331249e-03]\n",
      " [2.12363316e+00 1.38870143e+00 1.10699250e-03]\n",
      " [2.10371222e+00 1.37773832e+00 1.08844020e-03]\n",
      " [2.08168302e+00 1.36561496e+00 1.06810692e-03]\n",
      " [2.05796348e+00 1.35256136e+00 1.04642767e-03]\n",
      " [2.03296492e+00 1.33880388e+00 1.02381987e-03]\n",
      " [2.00709296e+00 1.32456573e+00 1.00068203e-03]\n",
      " [1.98072513e+00 1.31005469e+00 9.77372669e-04]\n",
      " [1.95418527e+00 1.29544898e+00 9.54188470e-04]\n",
      " [1.92773594e+00 1.28089309e+00 9.31360081e-04]\n",
      " [1.90157749e+00 1.26649728e+00 9.09054454e-04]\n",
      " [1.87585341e+00 1.25234052e+00 8.87382737e-04]\n",
      " [1.85066482e+00 1.23847846e+00 8.66415376e-04]\n",
      " [1.82608112e+00 1.22494928e+00 8.46193134e-04]\n",
      " [1.80214698e+00 1.21177758e+00 8.26734510e-04]\n",
      " [1.77889355e+00 1.19898049e+00 8.08045976e-04]\n",
      " [1.75635003e+00 1.18657409e+00 7.90131870e-04]\n",
      " [1.73454604e+00 1.17457468e+00 7.72996373e-04]\n",
      " [1.71350615e+00 1.16299577e+00 7.56639383e-04]\n",
      " [1.69324813e+00 1.15184715e+00 7.41055457e-04]\n",
      " [1.67378617e+00 1.14113663e+00 7.26236581e-04]\n",
      " [1.65513123e+00 1.13087023e+00 7.12172624e-04]\n",
      " [1.63728776e+00 1.12105042e+00 6.98849061e-04]\n",
      " [1.62025774e+00 1.11167827e+00 6.86250184e-04]\n",
      " [1.60404827e+00 1.10275769e+00 6.74364764e-04]\n",
      " [1.58866893e+00 1.09429396e+00 6.63183989e-04]\n",
      " [1.57412640e+00 1.08629075e+00 6.52697517e-04]\n",
      " [1.56042531e+00 1.07875062e+00 6.42894237e-04]\n",
      " [1.54756884e+00 1.07167530e+00 6.33762728e-04]\n",
      " [1.53556033e+00 1.06506664e+00 6.25292477e-04]\n",
      " [1.52440953e+00 1.05893001e+00 6.17478233e-04]\n",
      " [1.51412948e+00 1.05327258e+00 6.10317705e-04]\n",
      " [1.50472698e+00 1.04809809e+00 6.03804987e-04]\n",
      " [1.49620435e+00 1.04340783e+00 5.97931905e-04]\n",
      " [1.48856131e+00 1.03920162e+00 5.92689364e-04]\n",
      " [1.48179550e+00 1.03547818e+00 5.88067794e-04]\n",
      " [1.47590859e+00 1.03223844e+00 5.84061300e-04]\n",
      " [1.47090237e+00 1.02948336e+00 5.80664956e-04]\n",
      " [1.46677724e+00 1.02721317e+00 5.77873800e-04]\n",
      " [1.46354243e+00 1.02543295e+00 5.75689763e-04]\n",
      " [1.46121534e+00 1.02415228e+00 5.74121148e-04]\n",
      " [1.45981206e+00 1.02338002e+00 5.73176276e-04]\n",
      " [1.45934307e+00 1.02312192e+00 5.72860662e-04]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\elise\\AppData\\Local\\Temp\\ipykernel_26304\\2952692391.py:275: RuntimeWarning: divide by zero encountered in log10\n",
      "  srcv = np.log10(vel/0.02) #Velocity portion of the strain rate correction.\n",
      "C:\\Users\\elise\\AppData\\Local\\Temp\\ipykernel_26304\\2952692391.py:275: RuntimeWarning: invalid value encountered in log10\n",
      "  srcv = np.log10(vel/0.02) #Velocity portion of the strain rate correction.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max penetration depth is:  0.0218\n",
      "Pre-truncation d1          Count  Time (s)  2g (m/s^2)  18g (m/s^2)  50g (m/s^2)  200g (m/s^2)  \\\n",
      "64915  64915.0   32.4580   10.364800    10.014808   -37.606040    -44.381730   \n",
      "64916  64916.0   32.4585   10.651060    10.286805   -37.432653    -44.095712   \n",
      "64917  64917.0   32.4590   11.566406    11.226008   -37.014594    -43.465958   \n",
      "64918  64918.0   32.4595   12.300198    11.865991   -36.347262    -42.915405   \n",
      "64919  64919.0   32.4600   12.220535    11.808422   -35.526473    -42.829150   \n",
      "...        ...       ...         ...          ...          ...           ...   \n",
      "64984  64984.0   32.4925   19.781746   147.322909   144.427755     70.443785   \n",
      "64985  64985.0   32.4930   19.781297   109.785541   104.877040     39.510823   \n",
      "64986  64986.0   32.4935   19.781746    79.014942    67.639834     13.834893   \n",
      "64987  64987.0   32.4940   19.781332    58.460689    36.943951     -5.796812   \n",
      "64988  64988.0   32.4945   19.781129    38.658312    12.791691    -19.787348   \n",
      "\n",
      "       250g (m/s^2)  X55g (m/s^2)  Y55g (m/s^2)  Velocity (m/s)  \\\n",
      "64915    -31.076014    -22.352249     -8.096526        2.951189   \n",
      "64916    -31.221530    -22.736132     -8.867066        2.938755   \n",
      "64917    -30.614850    -22.562848     -8.741826        2.912609   \n",
      "64918    -28.786121    -22.444038     -8.562728        2.869480   \n",
      "64919    -27.682078    -22.441321     -8.687824        2.807153   \n",
      "...             ...           ...           ...             ...   \n",
      "64984    155.175623      1.238733      6.838586       -0.073680   \n",
      "64985    113.945240      0.277811      3.244123       -0.055712   \n",
      "64986     74.769462     -8.853878     -2.856426       -0.037371   \n",
      "64987     44.224040    -13.018254     -5.955675       -0.018760   \n",
      "64988     20.078446    -17.082119     -8.859597        0.000000   \n",
      "\n",
      "       Penetration Depth (m)  Pore Pressure (kPa)  \n",
      "64915               0.000000           102.421527  \n",
      "64916               0.001472           102.461846  \n",
      "64917               0.002935           102.497080  \n",
      "64918               0.004381           102.495627  \n",
      "64919               0.005800           102.501438  \n",
      "...                      ...                  ...  \n",
      "64984               0.014668           112.754020  \n",
      "64985               0.014635           110.921878  \n",
      "64986               0.014612           108.595742  \n",
      "64987               0.014598           105.846077  \n",
      "64988               0.014593           102.788028  \n",
      "\n",
      "[74 rows x 12 columns]\n",
      "Firmness factor =  260.9\n",
      "max penetration depth =  0.0218\n",
      "phi = 41.0\n",
      "Keep Drop? (y/n)\n"
     ]
    },
    {
     "ename": "FileNotFoundError",
     "evalue": "[Errno 2] No such file or directory: 'H:\\\\.shortcut-targets-by-id\\\\1aF9t2aiRGWTftJMZFAOBixqvQniFBjnb\\\\Duck  2023\\\\Data\\\\Intertidal\\\\BlueDrop, Samples & Moisture Gage\\\\07March23\\\\BD2 3.8.23 - Drops Only\\\\Analysis Figures\\\\004D 0 BC - REJECTED.png'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mFileNotFoundError\u001b[0m                         Traceback (most recent call last)",
      "\u001b[1;32mc:\\Users\\elise\\Documents\\GitHub\\BD-Code-Development\\BD code - single drop analysis.ipynb Cell 16\u001b[0m in \u001b[0;36m7\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20single%20drop%20analysis.ipynb#X25sZmlsZQ%3D%3D?line=67'>68</a>\u001b[0m dep_qsbc_comboplot(drop1g,drop1,acc1Nameg, qdyn1)\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20single%20drop%20analysis.ipynb#X25sZmlsZQ%3D%3D?line=68'>69</a>\u001b[0m plotName \u001b[39m=\u001b[39m fileNum\u001b[39m+\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m \u001b[39m\u001b[39m\"\u001b[39m\u001b[39m+\u001b[39m\u001b[39mstr\u001b[39m(a)\u001b[39m+\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m BC - REJECTED\u001b[39m\u001b[39m\"\u001b[39m\n\u001b[1;32m---> <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20single%20drop%20analysis.ipynb#X25sZmlsZQ%3D%3D?line=69'>70</a>\u001b[0m plt\u001b[39m.\u001b[39;49msavefig(plotPath \u001b[39m/\u001b[39;49m plotName)\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20single%20drop%20analysis.ipynb#X25sZmlsZQ%3D%3D?line=70'>71</a>\u001b[0m output_table1 \u001b[39m=\u001b[39m pd\u001b[39m.\u001b[39mDataFrame([BD, fileNum, a, soiltype, atype, tiptype, droptype,\u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m, \u001b[39m'\u001b[39m\u001b[39mNaN\u001b[39m\u001b[39m'\u001b[39m])\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20single%20drop%20analysis.ipynb#X25sZmlsZQ%3D%3D?line=71'>72</a>\u001b[0m output_table1 \u001b[39m=\u001b[39m output_table1\u001b[39m.\u001b[39mT\n",
      "File \u001b[1;32mc:\\Users\\elise\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\matplotlib\\pyplot.py:977\u001b[0m, in \u001b[0;36msavefig\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m    974\u001b[0m \u001b[39m@_copy_docstring_and_deprecators\u001b[39m(Figure\u001b[39m.\u001b[39msavefig)\n\u001b[0;32m    975\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39msavefig\u001b[39m(\u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs):\n\u001b[0;32m    976\u001b[0m     fig \u001b[39m=\u001b[39m gcf()\n\u001b[1;32m--> 977\u001b[0m     res \u001b[39m=\u001b[39m fig\u001b[39m.\u001b[39msavefig(\u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs)\n\u001b[0;32m    978\u001b[0m     fig\u001b[39m.\u001b[39mcanvas\u001b[39m.\u001b[39mdraw_idle()   \u001b[39m# need this if 'transparent=True' to reset colors\u001b[39;00m\n\u001b[0;32m    979\u001b[0m     \u001b[39mreturn\u001b[39;00m res\n",
      "File \u001b[1;32mc:\\Users\\elise\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\matplotlib\\figure.py:3058\u001b[0m, in \u001b[0;36mFigure.savefig\u001b[1;34m(self, fname, transparent, **kwargs)\u001b[0m\n\u001b[0;32m   3054\u001b[0m     \u001b[39mfor\u001b[39;00m ax \u001b[39min\u001b[39;00m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39maxes:\n\u001b[0;32m   3055\u001b[0m         stack\u001b[39m.\u001b[39menter_context(\n\u001b[0;32m   3056\u001b[0m             ax\u001b[39m.\u001b[39mpatch\u001b[39m.\u001b[39m_cm_set(facecolor\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mnone\u001b[39m\u001b[39m'\u001b[39m, edgecolor\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39mnone\u001b[39m\u001b[39m'\u001b[39m))\n\u001b[1;32m-> 3058\u001b[0m \u001b[39mself\u001b[39m\u001b[39m.\u001b[39mcanvas\u001b[39m.\u001b[39mprint_figure(fname, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs)\n",
      "File \u001b[1;32mc:\\Users\\elise\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\matplotlib\\backend_bases.py:2319\u001b[0m, in \u001b[0;36mFigureCanvasBase.print_figure\u001b[1;34m(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)\u001b[0m\n\u001b[0;32m   2315\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m   2316\u001b[0m     \u001b[39m# _get_renderer may change the figure dpi (as vector formats\u001b[39;00m\n\u001b[0;32m   2317\u001b[0m     \u001b[39m# force the figure dpi to 72), so we need to set it again here.\u001b[39;00m\n\u001b[0;32m   2318\u001b[0m     \u001b[39mwith\u001b[39;00m cbook\u001b[39m.\u001b[39m_setattr_cm(\u001b[39mself\u001b[39m\u001b[39m.\u001b[39mfigure, dpi\u001b[39m=\u001b[39mdpi):\n\u001b[1;32m-> 2319\u001b[0m         result \u001b[39m=\u001b[39m print_method(\n\u001b[0;32m   2320\u001b[0m             filename,\n\u001b[0;32m   2321\u001b[0m             facecolor\u001b[39m=\u001b[39mfacecolor,\n\u001b[0;32m   2322\u001b[0m             edgecolor\u001b[39m=\u001b[39medgecolor,\n\u001b[0;32m   2323\u001b[0m             orientation\u001b[39m=\u001b[39morientation,\n\u001b[0;32m   2324\u001b[0m             bbox_inches_restore\u001b[39m=\u001b[39m_bbox_inches_restore,\n\u001b[0;32m   2325\u001b[0m             \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs)\n\u001b[0;32m   2326\u001b[0m \u001b[39mfinally\u001b[39;00m:\n\u001b[0;32m   2327\u001b[0m     \u001b[39mif\u001b[39;00m bbox_inches \u001b[39mand\u001b[39;00m restore_bbox:\n",
      "File \u001b[1;32mc:\\Users\\elise\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\matplotlib\\backend_bases.py:1648\u001b[0m, in \u001b[0;36m_check_savefig_extra_args.<locals>.wrapper\u001b[1;34m(*args, **kwargs)\u001b[0m\n\u001b[0;32m   1640\u001b[0m     _api\u001b[39m.\u001b[39mwarn_deprecated(\n\u001b[0;32m   1641\u001b[0m         \u001b[39m'\u001b[39m\u001b[39m3.3\u001b[39m\u001b[39m'\u001b[39m, name\u001b[39m=\u001b[39mname, removal\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39m3.6\u001b[39m\u001b[39m'\u001b[39m,\n\u001b[0;32m   1642\u001b[0m         message\u001b[39m=\u001b[39m\u001b[39m'\u001b[39m\u001b[39m%(name)s\u001b[39;00m\u001b[39m() got unexpected keyword argument \u001b[39m\u001b[39m\"\u001b[39m\u001b[39m'\u001b[39m\n\u001b[0;32m   1643\u001b[0m                 \u001b[39m+\u001b[39m arg \u001b[39m+\u001b[39m \u001b[39m'\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m which is no longer supported as of \u001b[39m\u001b[39m'\u001b[39m\n\u001b[0;32m   1644\u001b[0m                 \u001b[39m'\u001b[39m\u001b[39m%(since)s\u001b[39;00m\u001b[39m and will become an error \u001b[39m\u001b[39m'\u001b[39m\n\u001b[0;32m   1645\u001b[0m                 \u001b[39m'\u001b[39m\u001b[39m%(removal)s\u001b[39;00m\u001b[39m'\u001b[39m)\n\u001b[0;32m   1646\u001b[0m     kwargs\u001b[39m.\u001b[39mpop(arg)\n\u001b[1;32m-> 1648\u001b[0m \u001b[39mreturn\u001b[39;00m func(\u001b[39m*\u001b[39margs, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs)\n",
      "File \u001b[1;32mc:\\Users\\elise\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\matplotlib\\_api\\deprecation.py:415\u001b[0m, in \u001b[0;36mdelete_parameter.<locals>.wrapper\u001b[1;34m(*inner_args, **inner_kwargs)\u001b[0m\n\u001b[0;32m    405\u001b[0m     deprecation_addendum \u001b[39m=\u001b[39m (\n\u001b[0;32m    406\u001b[0m         \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mIf any parameter follows \u001b[39m\u001b[39m{\u001b[39;00mname\u001b[39m!r}\u001b[39;00m\u001b[39m, they should be passed as \u001b[39m\u001b[39m\"\u001b[39m\n\u001b[0;32m    407\u001b[0m         \u001b[39mf\u001b[39m\u001b[39m\"\u001b[39m\u001b[39mkeyword, not positionally.\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[0;32m    408\u001b[0m     warn_deprecated(\n\u001b[0;32m    409\u001b[0m         since,\n\u001b[0;32m    410\u001b[0m         name\u001b[39m=\u001b[39m\u001b[39mrepr\u001b[39m(name),\n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m    413\u001b[0m                  \u001b[39melse\u001b[39;00m deprecation_addendum,\n\u001b[0;32m    414\u001b[0m         \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mkwargs)\n\u001b[1;32m--> 415\u001b[0m \u001b[39mreturn\u001b[39;00m func(\u001b[39m*\u001b[39minner_args, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39minner_kwargs)\n",
      "File \u001b[1;32mc:\\Users\\elise\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\matplotlib\\backends\\backend_agg.py:541\u001b[0m, in \u001b[0;36mFigureCanvasAgg.print_png\u001b[1;34m(self, filename_or_obj, metadata, pil_kwargs, *args)\u001b[0m\n\u001b[0;32m    494\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[0;32m    495\u001b[0m \u001b[39mWrite the figure to a PNG file.\u001b[39;00m\n\u001b[0;32m    496\u001b[0m \n\u001b[1;32m   (...)\u001b[0m\n\u001b[0;32m    538\u001b[0m \u001b[39m    *metadata*, including the default 'Software' key.\u001b[39;00m\n\u001b[0;32m    539\u001b[0m \u001b[39m\"\"\"\u001b[39;00m\n\u001b[0;32m    540\u001b[0m FigureCanvasAgg\u001b[39m.\u001b[39mdraw(\u001b[39mself\u001b[39m)\n\u001b[1;32m--> 541\u001b[0m mpl\u001b[39m.\u001b[39;49mimage\u001b[39m.\u001b[39;49mimsave(\n\u001b[0;32m    542\u001b[0m     filename_or_obj, \u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mbuffer_rgba(), \u001b[39mformat\u001b[39;49m\u001b[39m=\u001b[39;49m\u001b[39m\"\u001b[39;49m\u001b[39mpng\u001b[39;49m\u001b[39m\"\u001b[39;49m, origin\u001b[39m=\u001b[39;49m\u001b[39m\"\u001b[39;49m\u001b[39mupper\u001b[39;49m\u001b[39m\"\u001b[39;49m,\n\u001b[0;32m    543\u001b[0m     dpi\u001b[39m=\u001b[39;49m\u001b[39mself\u001b[39;49m\u001b[39m.\u001b[39;49mfigure\u001b[39m.\u001b[39;49mdpi, metadata\u001b[39m=\u001b[39;49mmetadata, pil_kwargs\u001b[39m=\u001b[39;49mpil_kwargs)\n",
      "File \u001b[1;32mc:\\Users\\elise\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\matplotlib\\image.py:1675\u001b[0m, in \u001b[0;36mimsave\u001b[1;34m(fname, arr, vmin, vmax, cmap, format, origin, dpi, metadata, pil_kwargs)\u001b[0m\n\u001b[0;32m   1673\u001b[0m pil_kwargs\u001b[39m.\u001b[39msetdefault(\u001b[39m\"\u001b[39m\u001b[39mformat\u001b[39m\u001b[39m\"\u001b[39m, \u001b[39mformat\u001b[39m)\n\u001b[0;32m   1674\u001b[0m pil_kwargs\u001b[39m.\u001b[39msetdefault(\u001b[39m\"\u001b[39m\u001b[39mdpi\u001b[39m\u001b[39m\"\u001b[39m, (dpi, dpi))\n\u001b[1;32m-> 1675\u001b[0m image\u001b[39m.\u001b[39msave(fname, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mpil_kwargs)\n",
      "File \u001b[1;32mc:\\Users\\elise\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\PIL\\Image.py:2317\u001b[0m, in \u001b[0;36mImage.save\u001b[1;34m(self, fp, format, **params)\u001b[0m\n\u001b[0;32m   2315\u001b[0m         fp \u001b[39m=\u001b[39m builtins\u001b[39m.\u001b[39mopen(filename, \u001b[39m\"\u001b[39m\u001b[39mr+b\u001b[39m\u001b[39m\"\u001b[39m)\n\u001b[0;32m   2316\u001b[0m     \u001b[39melse\u001b[39;00m:\n\u001b[1;32m-> 2317\u001b[0m         fp \u001b[39m=\u001b[39m builtins\u001b[39m.\u001b[39;49mopen(filename, \u001b[39m\"\u001b[39;49m\u001b[39mw+b\u001b[39;49m\u001b[39m\"\u001b[39;49m)\n\u001b[0;32m   2319\u001b[0m \u001b[39mtry\u001b[39;00m:\n\u001b[0;32m   2320\u001b[0m     save_handler(\u001b[39mself\u001b[39m, fp, filename)\n",
      "\u001b[1;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'H:\\\\.shortcut-targets-by-id\\\\1aF9t2aiRGWTftJMZFAOBixqvQniFBjnb\\\\Duck  2023\\\\Data\\\\Intertidal\\\\BlueDrop, Samples & Moisture Gage\\\\07March23\\\\BD2 3.8.23 - Drops Only\\\\Analysis Figures\\\\004D 0 BC - REJECTED.png'"
     ]
    }
   ],
   "source": [
    "finddropend(a) #automatically find the end of the drop\n",
    "print(\"penend: \", penend)\n",
    "offset = acc1g[penstart] #sets the offset to be the acceleration at the moment of impact\n",
    "print(\"offset: \", offset)\n",
    "#print(drop)\n",
    "dg = dfCalg[penstart:penend] #chops the file to only include penetration \n",
    "d = dfCal[penstart:penend]\n",
    "acc1g = accg[penstart:penend]\n",
    "acc1g = acc1g-offset\n",
    "acc1 = acc[penstart:penend]\n",
    "offsetms = offset * 9.80665\n",
    "acc1 = acc1-offsetms\n",
    "\n",
    "\n",
    "drop1 = d\n",
    "drop1g = dg\n",
    "integration(d, acc1) #double integration to solve for velocity and penetration depth\n",
    "drop1 = d\n",
    "areafind() #area calculations for bearing capacity calculations\n",
    "qdynfun(acc1) #calculation of dynamic bearing capacity\n",
    "qsbcfun(acc1) #calculation of quasi-static bearing capacities\n",
    "qdyn1 = qdyn\n",
    "integplot(drop1,acc1Name)\n",
    "#qdynplot(drop1, qdyn1)\n",
    "print(\"Max penetration depth is: \", maxdep)\n",
    "print(\"Pre-truncation d1\", drop1)\n",
    "drop1 = pd.DataFrame(drop1)\n",
    "drop1g = pd.DataFrame(drop1g)\n",
    "acc1 = pd.DataFrame(acc1)\n",
    "qdyn1 = pd.DataFrame(qdyn1)\n",
    "if droptype == \"w\":\n",
    "    bctable = pd.DataFrame(bctable)\n",
    "\n",
    "drop1t = drop1 #.truncate(before=trunc_index) #this truncates each df such that only depths below cone height are considered\n",
    "drop1gt = drop1g #.truncate(before=trunc_index)\n",
    "acc1t = acc1 #.truncate(before=trunc_index)\n",
    "acc1t = np.array(acc1t)\n",
    "acc1gt = acc1g #.truncate(before=trunc_index)\n",
    "acc1gt = np.array(acc1gt)\n",
    "drop1t = drop1 #truncate(before=trunc_index)\n",
    "qdyn1t = qdyn1#.truncate(before=trunc_index)\n",
    "if droptype == \"w\":\n",
    "    bctablet = bctable#.truncate(before=trunc_index)\n",
    "\n",
    "#bctruncate(drop1g, drop1, acc1, qdyn1)\n",
    "firmnessfactor(acc1t, drop1) #calculates the firmness factor\n",
    "if soiltype == \"s\":\n",
    "    dr(acc1gt) #calculates relative density\n",
    "    duncan_correlation() #calculates friction angle\n",
    "elif soiltype == \"u\":\n",
    "    if maxdep <= .2:\n",
    "        dr(acc1gt)\n",
    "        duncan_correlation()\n",
    "pentime=(penend - penstart)/2 #ms #calculates the time of penetration\n",
    "\n",
    "dep_qsbc_comboplot(drop1g,drop1,acc1Nameg, qdyn1)\n",
    "plt.show()\n",
    "ignore_drop()\n",
    "if keep_drop == \"y\":\n",
    "    duck_dep_qsbc_comboplot(drop1g, drop1, drop1t, acc1Nameg, qdyn1t)\n",
    "    print(\"pentime= \", pentime)\n",
    "\n",
    "    output_table1 = pd.DataFrame([BD, fileNum, a, soiltype, atype, tiptype, droptype, maxdep*100, maxacc, maxvel, maxbcdep*100, maxqdyn, maxHqsbc, maxLqsbc, maxAqsbc, ff, pentime, Dr, phi])\n",
    "    output_table1 = output_table1.T\n",
    "    output_table = output_table1\n",
    "    output_table.columns = [\"BlueDrop\", \"BD File\", \"Drop No.\", \"Soil Type\", \"Area Type\", \"Tip Type\", \"Drop Type\", \"Penetration Depth (cm)\", \"Maximum Deceleration (m/s^2)\", \"Impact Velocity (m/s)\", \"Depth of Max QSBC (cm)\", \"Max Qdyn (kPa)\", \"Max QSBC-Upper Bounds (kPa)\", \"Max QSBC-Lower Bounds (kPa)\", \"Max QSBC-Average (kPa)\", \"Firmness Factor (m-1)\", \"Penetration Time (ms)\", \"Relative Density (%)\", \"Friction Angle (deg.)\"]  # Creates an empty Dataframe in pandas from the output data\n",
    "else:\n",
    "    dep_qsbc_comboplot(drop1g,drop1,acc1Nameg, qdyn1)\n",
    "    plotName = fileNum+\" \"+str(a)+\" BC - REJECTED\"\n",
    "    plt.savefig(plotPath / plotName)\n",
    "    output_table1 = pd.DataFrame([BD, fileNum, a, soiltype, atype, tiptype, droptype,'NaN', 'NaN', 'NaN', 'NaN', 'NaN', 'NaN', 'NaN', 'NaN', 'NaN', 'NaN', 'NaN', 'NaN'])\n",
    "    output_table1 = output_table1.T\n",
    "    output_table = output_table1\n",
    "    output_table.columns = [\"BlueDrop\", \"BD File\", \"Drop No.\", \"Soil Type\", \"Area Type\", \"Tip Type\", \"Drop Type\", \"Penetration Depth (cm)\", \"Maximum Deceleration (m/s^2)\", \"Impact Velocity (m/s)\", \"Depth of Max QSBC (cm)\", \"Max Qdyn (kPa)\", \"Max QSBC-Upper Bounds (kPa)\", \"Max QSBC-Lower Bounds (kPa)\", \"Max QSBC-Average (kPa)\", \"Firmness Factor (m-1)\", \"Penetration Time (ms)\", \"Relative Density (%)\", \"Friction Angle (deg.)\"]  # Creates an empty Dataframe in pandas from the output data\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}