BD code - profile overlay.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# SETUP VARIABLES - USER INPUTS\n",
    "BD = 2 #Bluedrop file is from \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 = 'a' #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/08March23/BD2 3.8.23 - Drops Only\")\n",
    "#paste the filepath to an excel file that the analysis results will be printed in"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "set functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "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": 8,
   "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": 13,
   "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)\")\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)\"])\n",
    "    ax1.set_title(\"BD file \"+fileNum)\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": 10,
   "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 for first sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak index: 52569\n",
      "peak index: 52569\n",
      "start of penetration:  52522\n",
      "penend:  52595\n",
      "offset:  -0.11577380812803753\n",
      "[[0.00000000e+00 2.20000000e-01 1.52053084e-05]\n",
      " [1.99912229e-01 3.30017922e-01 4.00037594e-05]\n",
      " [3.99495817e-01 4.39854983e-01 8.21087376e-05]\n",
      " [5.98674146e-01 5.49469017e-01 1.40272516e-04]\n",
      " [7.97335481e-01 6.58798533e-01 2.14064942e-04]\n",
      " [9.95317945e-01 7.67754445e-01 3.03190190e-04]\n",
      " [1.19237563e+00 8.76201423e-01 4.07310090e-04]\n",
      " [1.38814051e+00 9.83936932e-01 5.25952310e-04]\n",
      " [1.58212064e+00 1.09069023e+00 6.58451743e-04]\n",
      " [1.77374864e+00 1.19614909e+00 8.03939924e-04]\n",
      " [1.96245662e+00 1.30000096e+00 9.61384172e-04]\n",
      " [2.14774379e+00 1.40197026e+00 1.12965623e-03]\n",
      " [2.32921072e+00 1.50183716e+00 1.30759779e-03]\n",
      " [2.50655550e+00 1.59943551e+00 1.49405868e-03]\n",
      " [2.67954773e+00 1.69463852e+00 1.68790633e-03]\n",
      " [2.84800107e+00 1.78734364e+00 1.88802091e-03]\n",
      " [3.01175882e+00 1.87746462e+00 2.09329509e-03]\n",
      " [3.17068923e+00 1.96492897e+00 2.30264029e-03]\n",
      " [3.32467829e+00 2.04967395e+00 2.51498711e-03]\n",
      " [3.47362699e+00 2.13164505e+00 2.72928893e-03]\n",
      " [3.61745562e+00 2.21079842e+00 2.94453509e-03]\n",
      " [3.75609998e+00 2.28709873e+00 3.15975203e-03]\n",
      " [3.88950548e+00 2.36051593e+00 3.37399944e-03]\n",
      " [4.01762779e+00 2.43102562e+00 3.58637507e-03]\n",
      " [4.14043070e+00 2.49860789e+00 3.79601466e-03]\n",
      " [4.25788387e+00 2.56324602e+00 4.00209072e-03]\n",
      " [4.36997128e+00 2.62493121e+00 4.20383000e-03]\n",
      " [4.47669855e+00 2.68366655e+00 4.40053078e-03]\n",
      " [4.57808248e+00 2.73946128e+00 4.59154763e-03]\n",
      " [4.67414457e+00 2.79232724e+00 4.77628084e-03]\n",
      " [4.76491813e+00 2.84228276e+00 4.95419091e-03]\n",
      " [4.85044527e+00 2.88935101e+00 5.12479411e-03]\n",
      " [4.93077403e+00 2.93355842e+00 5.28765727e-03]\n",
      " [5.00596825e+00 2.97494014e+00 5.44241827e-03]\n",
      " [5.07611355e+00 3.01354328e+00 5.58879980e-03]\n",
      " [5.14131390e+00 3.04942507e+00 5.72660359e-03]\n",
      " [5.20168508e+00 3.08264921e+00 5.85569679e-03]\n",
      " [5.25735132e+00 3.11328407e+00 5.97600405e-03]\n",
      " [5.30844601e+00 3.14140307e+00 6.08750770e-03]\n",
      " [5.35511105e+00 3.16708429e+00 6.19024484e-03]\n",
      " [5.39749603e+00 3.19041007e+00 6.28430392e-03]\n",
      " [5.43576075e+00 3.21146833e+00 6.36982852e-03]\n",
      " [5.47008109e+00 3.23035589e+00 6.44702876e-03]\n",
      " [5.50065177e+00 3.24717988e+00 6.51618598e-03]\n",
      " [5.52768168e+00 3.26205528e+00 6.57764051e-03]\n",
      " [5.55139414e+00 3.27510499e+00 6.63179003e-03]\n",
      " [5.57202933e+00 3.28646117e+00 6.67909292e-03]\n",
      " [5.58983370e+00 3.29625947e+00 6.72004167e-03]\n",
      " [5.60505753e+00 3.30463762e+00 6.75515459e-03]\n",
      " [5.61796136e+00 3.31173900e+00 6.78498820e-03]\n",
      " [5.62880726e+00 3.31770783e+00 6.81011471e-03]\n",
      " [5.63784755e+00 3.32268299e+00 6.83109370e-03]\n",
      " [5.64531635e+00 3.32679330e+00 6.84845018e-03]\n",
      " [5.65142190e+00 3.33015338e+00 6.86265500e-03]\n",
      " [5.65635092e+00 3.33286597e+00 6.87413333e-03]\n",
      " [5.66027835e+00 3.33502736e+00 6.88328609e-03]\n",
      " [5.66337157e+00 3.33672965e+00 6.89049902e-03]\n",
      " [5.66578670e+00 3.33805877e+00 6.89613339e-03]\n",
      " [5.66766044e+00 3.33908995e+00 6.90050630e-03]\n",
      " [5.66910534e+00 3.33988512e+00 6.90387933e-03]\n",
      " [5.67020945e+00 3.34049275e+00 6.90645736e-03]\n",
      " [5.67104109e+00 3.34095043e+00 6.90839952e-03]\n",
      " [5.67165857e+00 3.34129024e+00 6.90984171e-03]\n",
      " [5.67211111e+00 3.34153929e+00 6.91089878e-03]\n",
      " [5.67243358e+00 3.34171676e+00 6.91165205e-03]\n",
      " [5.67265191e+00 3.34183691e+00 6.91216209e-03]\n",
      " [5.67279043e+00 3.34191314e+00 6.91248570e-03]\n",
      " [5.67287206e+00 3.34195807e+00 6.91267640e-03]\n",
      " [5.67291721e+00 3.34198292e+00 6.91278189e-03]\n",
      " [5.67294281e+00 3.34199700e+00 6.91284169e-03]\n",
      " [5.67295802e+00 3.34200538e+00 6.91287724e-03]\n",
      " [5.67296473e+00 3.34200907e+00 6.91289291e-03]\n",
      " [5.67296600e+00 3.34200977e+00 6.91289587e-03]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\elise\\AppData\\Local\\Temp\\ipykernel_6600\\2952692391.py:296: RuntimeWarning: divide by zero 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.0567\n",
      "Pre-truncation d1          Count  Time (s)  2g (m/s^2)  18g (m/s^2)  50g (m/s^2)  200g (m/s^2)  \\\n",
      "52522  52522.0   26.2615    0.123339    -0.090929    -1.135353      7.583530   \n",
      "52523  52523.0   26.2620    0.589063     0.345147    -0.932839      8.081497   \n",
      "52524  52524.0   26.2625    0.702590     0.445602    -0.669309      8.076366   \n",
      "52525  52525.0   26.2630    0.983377     0.670002    -0.908923      8.360782   \n",
      "52526  52526.0   26.2635    1.305609     1.057197    -0.393516      8.817706   \n",
      "...        ...       ...         ...          ...          ...           ...   \n",
      "52590  52590.0   26.2955   18.296934    18.860279    25.248705     26.429967   \n",
      "52591  52591.0   26.2960   15.236179    15.091169    19.092327     22.685433   \n",
      "52592  52592.0   26.2965   12.507849    12.346727    14.744169     19.967296   \n",
      "52593  52593.0   26.2970   10.511706    10.234559    11.719331     17.858070   \n",
      "52594  52594.0   26.2975    8.852879     8.612671     9.858375     16.261240   \n",
      "\n",
      "       250g (m/s^2)  X55g (m/s^2)  Y55g (m/s^2)  Velocity (m/s)  \\\n",
      "52522     -0.751401      0.761632      1.625056        4.001226   \n",
      "52523     -0.576926      0.699724      1.480714        3.995264   \n",
      "52524     -0.306888      0.932199      1.444090        3.988080   \n",
      "52525     -0.250419      1.303643      1.773275        3.979053   \n",
      "52526      0.370017      1.256890      1.948927        3.967400   \n",
      "...             ...           ...           ...             ...   \n",
      "52590     23.664232      0.664410      1.222623        0.000633   \n",
      "52591     18.143294      0.417782      1.615577        0.000391   \n",
      "52592     14.864469      0.647539      1.608396        0.000218   \n",
      "52593     12.049703      1.313794      1.665845        0.000051   \n",
      "52594      8.607261      1.189121      1.993019        0.000000   \n",
      "\n",
      "       Penetration Depth (m)  Pore Pressure (kPa)  \n",
      "52522               0.000000            -1.429354  \n",
      "52523               0.001999            -1.441340  \n",
      "52524               0.003995            -1.443156  \n",
      "52525               0.005987            -1.436981  \n",
      "52526               0.007973            -1.399932  \n",
      "...                      ...                  ...  \n",
      "52590               0.056729            -1.483112  \n",
      "52591               0.056729            -1.496551  \n",
      "52592               0.056730            -1.483112  \n",
      "52593               0.056730            -1.446062  \n",
      "52594               0.056730            -1.438071  \n",
      "\n",
      "[73 rows x 12 columns]\n"
     ]
    }
   ],
   "source": [
    "fileNum= '014F' # write the bin file number you want to analyze (do not include 'bLog' or '.bin')\n",
    "a = 3 #what drop do you want to analyze? USER INPUT\n",
    "a = a-1\n",
    "\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)\n",
    "\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",
    "peakplot()\n",
    "#overviewplot()\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",
    "\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",
    "\n",
    "peakpick()\n",
    "\n",
    "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",
    "qsbc1 = qsbc\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",
    "\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",
    "\n",
    "duck_dep_qsbc_comboplot(drop1g, drop1, drop1t, acc1Nameg, qdyn1t)\n",
    "\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trial 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak index: 53576\n",
      "start of penetration:  53537\n",
      "[[0.00000000e+00 2.20000000e-01 1.52053084e-05]\n",
      " [2.67703840e-01 3.67325756e-01 5.24517085e-05]\n",
      " [5.34974050e-01 5.14412871e-01 1.19940456e-04]\n",
      " [8.01651552e-01 6.61173801e-01 2.15842204e-04]\n",
      " [1.06747736e+00 8.07466015e-01 3.39534301e-04]\n",
      " [1.33206241e+00 9.53075403e-01 4.90418300e-04]\n",
      " [1.59488199e+00 1.09771320e+00 6.67689633e-04]\n",
      " [1.85533038e+00 1.24104606e+00 8.70280151e-04]\n",
      " [2.11281431e+00 1.38274748e+00 1.09689750e-03]\n",
      " [2.36681864e+00 1.52253397e+00 1.34610152e-03]\n",
      " [2.61693180e+00 1.66017902e+00 1.61637741e-03]\n",
      " [2.86285098e+00 1.79551601e+00 1.90619918e-03]\n",
      " [3.10435842e+00 1.92842507e+00 2.21405621e-03]\n",
      " [3.34128051e+00 2.05881067e+00 2.53844011e-03]\n",
      " [3.57346740e+00 2.18659034e+00 2.87783793e-03]\n",
      " [3.80079196e+00 2.31169412e+00 3.23074425e-03]\n",
      " [4.02314557e+00 2.43406223e+00 3.59566688e-03]\n",
      " [4.24042806e+00 2.55363955e+00 3.97111934e-03]\n",
      " [4.45254383e+00 2.67037346e+00 4.35561946e-03]\n",
      " [4.65940481e+00 2.78421550e+00 4.74769883e-03]\n",
      " [4.86092908e+00 2.89512058e+00 5.14590509e-03]\n",
      " [5.05703541e+00 3.00304399e+00 5.54879464e-03]\n",
      " [5.24764402e+00 3.10794184e+00 5.95493636e-03]\n",
      " [5.43268028e+00 3.20977306e+00 6.36292205e-03]\n",
      " [5.61207525e+00 3.30849969e+00 6.77137139e-03]\n",
      " [5.78576929e+00 3.40408893e+00 7.17894404e-03]\n",
      " [5.95370825e+00 3.49651096e+00 7.58433496e-03]\n",
      " [6.11583464e+00 3.58573416e+00 7.98625533e-03]\n",
      " [6.27209250e+00 3.67172773e+00 8.38344533e-03]\n",
      " [6.42243735e+00 3.75446718e+00 8.77470205e-03]\n",
      " [6.56683974e+00 3.83393631e+00 9.15889370e-03]\n",
      " [6.70527733e+00 3.91012282e+00 9.53494282e-03]\n",
      " [6.83772176e+00 3.98301112e+00 9.90179244e-03]\n",
      " [6.96414222e+00 4.05258423e+00 1.02584150e-02]\n",
      " [7.08451649e+00 4.11882994e+00 1.06038443e-02]\n",
      " [7.19883055e+00 4.18174052e+00 1.09371776e-02]\n",
      " [7.30707557e+00 4.24131113e+00 1.12575681e-02]\n",
      " [7.40924676e+00 4.29753911e+00 1.15642236e-02]\n",
      " [7.50533883e+00 4.35042157e+00 1.18563930e-02]\n",
      " [7.59534603e+00 4.37500000e+00 1.19934265e-02]\n",
      " [7.67927176e+00 4.37500000e+00 1.19934265e-02]\n",
      " [7.75713698e+00 4.37500000e+00 1.19934265e-02]\n",
      " [7.82898389e+00 4.37500000e+00 1.19934265e-02]\n",
      " [7.89487667e+00 4.37500000e+00 1.19934265e-02]\n",
      " [7.95489976e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.00917355e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.05787399e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.10122552e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.13949515e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.17300031e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.20210403e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.22719892e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.24869217e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.26698264e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.28243761e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.29539211e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.30616133e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.31504230e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.32230401e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.32818239e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.33288621e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.33660394e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.33950488e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34173771e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34343045e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34469069e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34560574e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34624695e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34667688e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34695047e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34711504e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34720906e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34725915e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34728267e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34729121e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.34729286e+00 4.37500000e+00 1.19934265e-02]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\elise\\AppData\\Local\\Temp\\ipykernel_6600\\2952692391.py:296: RuntimeWarning: divide by zero 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.0835\n",
      "Pre-truncation d1          Count  Time (s)  2g (m/s^2)  18g (m/s^2)  50g (m/s^2)  200g (m/s^2)  \\\n",
      "53537  53537.0   26.7690   -0.021884    -0.325301    -1.561538     12.106918   \n",
      "53538  53538.0   26.7695   -0.002523    -0.299851    -1.297549     12.115255   \n",
      "53539  53539.0   26.7700    0.408521     0.062016    -0.986648     12.454822   \n",
      "53540  53540.0   26.7705    1.107617     0.777492    -0.579473     13.174357   \n",
      "53541  53541.0   26.7710    2.166898     1.849819     0.150253     14.221275   \n",
      "...        ...       ...         ...          ...          ...           ...   \n",
      "53608  53608.0   26.8045   19.580556    29.610954    46.478073     42.543516   \n",
      "53609  53609.0   26.8050   18.903300    20.800537    31.219584     33.568885   \n",
      "53610  53610.0   26.8055   15.805783    15.601148    21.108273     28.315380   \n",
      "53611  53611.0   26.8060   12.700919    12.385391    15.151038     25.038352   \n",
      "53612  53612.0   26.8065   10.586530    10.319596    11.724236     22.889045   \n",
      "\n",
      "       250g (m/s^2)  X55g (m/s^2)  Y55g (m/s^2)  Velocity (m/s)  \\\n",
      "53537     -2.028470      0.919331      2.899139        5.357827   \n",
      "53538     -1.205325      0.792228      2.695050        5.350327   \n",
      "53539     -0.842620      0.883445      2.951993        5.340481   \n",
      "53540     -1.171299      0.895454      2.760543        5.326619   \n",
      "53541     -0.703619      1.092471      2.430065        5.306413   \n",
      "...             ...           ...           ...             ...   \n",
      "53608     45.897841      5.865627      2.490243        0.001338   \n",
      "53609     29.851214      3.343296      0.332301        0.000665   \n",
      "53610     20.405677      1.789894      0.404113        0.000276   \n",
      "53611     15.071522      1.382278      0.786439        0.000066   \n",
      "53612     12.371866      0.999253      1.018103        0.000000   \n",
      "\n",
      "       Penetration Depth (m)  Pore Pressure (kPa)  \n",
      "53537               0.000000            -2.037040  \n",
      "53538               0.002677            -2.051206  \n",
      "53539               0.005350            -2.067915  \n",
      "53540               0.008017            -2.070458  \n",
      "53541               0.010675            -2.078812  \n",
      "...                      ...                  ...  \n",
      "53608               0.083472            -2.197589  \n",
      "53609               0.083473            -2.206670  \n",
      "53610               0.083473            -2.167077  \n",
      "53611               0.083473            -2.150005  \n",
      "53612               0.083473            -2.178338  \n",
      "\n",
      "[76 rows x 12 columns]\n"
     ]
    }
   ],
   "source": [
    "fileNum= '0182' # write the bin file number you want to analyze (do not include 'bLog' or '.bin')\n",
    "a = 4 #what drop do you want to analyze? USER INPUT\n",
    "a = a-1\n",
    "\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)\n",
    "\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",
    "peakplot()\n",
    "#overviewplot()\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",
    "drop2 = dfCal[dropstart:dropend]  # start and end points of the drop in m/s^2\n",
    "drop2g = dfCalg[dropstart:dropend]  # start and end points of the drop in g\n",
    "drop2 = pd.DataFrame(drop2)  # makes dataframe including all data within the start and end points of the drop\n",
    "drop2g = pd.DataFrame(drop2g)\n",
    "dg = drop2g \n",
    "d = drop2\n",
    "accPick(dg, d)  # chooses what accelerometer to use\n",
    "acc2 = acc\n",
    "acc2g = accg\n",
    "acc2Name = accName\n",
    "acc2Nameg = accNameg\n",
    "#findchangepts() #automatically identifies the moment of penetration\n",
    "\n",
    "peak = int(peaksArray[a])  # count at the ath drop\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",
    "drop2 = dfCal[dropstart:dropend]  # start and end points of the drop in m/s^2\n",
    "drop2g = dfCalg[dropstart:dropend]  # start and end points of the drop in g\n",
    "drop2 = pd.DataFrame(drop2)  # makes dataframe including all data within the start and end points of the drop\n",
    "drop2g = pd.DataFrame(drop2g)\n",
    "dg = drop2g \n",
    "d = drop2\n",
    "accPick(dg, d)  # chooses what accelerometer to use\n",
    "acc2 = acc\n",
    "acc2g = accg\n",
    "acc2Name = accName\n",
    "acc2Nameg = accNameg\n",
    "#findchangepts() #automatically identifies the moment of penetration\n",
    "\n",
    "peakpick()\n",
    "\n",
    "finddropend(a) #automatically find the end of the drop\n",
    "offset = acc2g[penstart] #sets the offset to be the acceleration at the moment of impact\n",
    "#print(drop)\n",
    "dg = dfCalg[penstart:penend] #chops the file to only include penetration \n",
    "d = dfCal[penstart:penend]\n",
    "acc2g = accg[penstart:penend]\n",
    "acc2g = acc2g-offset\n",
    "acc2 = acc[penstart:penend]\n",
    "offsetms = offset * 9.80665\n",
    "acc2 = acc2-offsetms\n",
    "\n",
    "\n",
    "drop2 = d\n",
    "drop2g = dg\n",
    "integration(d, acc2) #double integration to solve for velocity and penetration depth\n",
    "drop2 = d\n",
    "areafind() #area calculations for bearing capacity calculations\n",
    "qdynfun(acc2) #calculation of dynamic bearing capacity\n",
    "qsbcfun(acc2) #calculation of quasi-static bearing capacities\n",
    "qsbc2 = qsbc\n",
    "qdyn2 = qdyn\n",
    "integplot(drop2,acc2Name)\n",
    "#qdynplot(drop1, qdyn1)\n",
    "print(\"Max penetration depth is: \", maxdep)\n",
    "print(\"Pre-truncation d1\", drop2)\n",
    "drop2 = pd.DataFrame(drop2)\n",
    "drop2g = pd.DataFrame(drop2g)\n",
    "acc2 = pd.DataFrame(acc2)\n",
    "qdyn2 = pd.DataFrame(qdyn2)\n",
    "\n",
    "drop2t = drop2 #.truncate(before=trunc_index) #this truncates each df such that only depths below cone height are considered\n",
    "drop2gt = drop2g #.truncate(before=trunc_index)\n",
    "acc2t = acc2 #.truncate(before=trunc_index)\n",
    "acc2t = np.array(acc2t)\n",
    "acc2gt = acc2g #.truncate(before=trunc_index)\n",
    "acc2gt = np.array(acc2gt)\n",
    "drop2t = drop2 #truncate(before=trunc_index)\n",
    "qdyn2t = qdyn2#.truncate(before=trunc_index)\n",
    "\n",
    "duck_dep_qsbc_comboplot(drop2g, drop2, drop2t, acc2Nameg, qdyn2t)\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trial 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak index: 100949\n",
      "start of penetration:  100907\n",
      "[[0.00000000e+00 2.20000000e-01 1.52053084e-05]\n",
      " [2.66738636e-01 3.66794574e-01 5.22609137e-05]\n",
      " [5.33079768e-01 5.13370389e-01 1.19360672e-04]\n",
      " [7.98922656e-01 6.59672005e-01 2.14717646e-04]\n",
      " [1.06412463e+00 8.05620903e-01 3.37800002e-04]\n",
      " [1.32849110e+00 9.51110004e-01 4.88197433e-04]\n",
      " [1.59177893e+00 1.09600549e+00 6.65437420e-04]\n",
      " [1.85371491e+00 1.24015701e+00 8.68940983e-04]\n",
      " [2.11402188e+00 1.38341204e+00 1.09802198e-03]\n",
      " [2.37243880e+00 1.52562692e+00 1.35190351e-03]\n",
      " [2.62872662e+00 1.66667008e+00 1.62973283e-03]\n",
      " [2.88266070e+00 1.80641791e+00 1.93058442e-03]\n",
      " [3.13402669e+00 1.94475244e+00 2.25346318e-03]\n",
      " [3.38262020e+00 2.08156118e+00 2.59731105e-03]\n",
      " [3.62823619e+00 2.21673131e+00 2.96099767e-03]\n",
      " [3.87066541e+00 2.35014765e+00 3.34331656e-03]\n",
      " [4.10970849e+00 2.48170050e+00 3.74300979e-03]\n",
      " [4.34518551e+00 2.61129083e+00 4.15879256e-03]\n",
      " [4.57693287e+00 2.73882862e+00 4.58935892e-03]\n",
      " [4.80480250e+00 2.86423237e+00 5.03339003e-03]\n",
      " [5.02866467e+00 2.98743069e+00 5.48956939e-03]\n",
      " [5.24840262e+00 3.10835932e+00 5.95658140e-03]\n",
      " [5.46390544e+00 3.22695723e+00 6.43310293e-03]\n",
      " [5.67506904e+00 3.34316713e+00 6.91780997e-03]\n",
      " [5.88179617e+00 3.45693551e+00 7.40938230e-03]\n",
      " [6.08399120e+00 3.56820973e+00 7.90649493e-03]\n",
      " [6.28155737e+00 3.67693654e+00 8.40781352e-03]\n",
      " [6.47439809e+00 3.78306279e+00 8.91199894e-03]\n",
      " [6.66242102e+00 3.88653766e+00 9.41772040e-03]\n",
      " [6.84554170e+00 3.98731467e+00 9.92366886e-03]\n",
      " [7.02368141e+00 4.08535050e+00 1.04285555e-02]\n",
      " [7.19676617e+00 4.18060443e+00 1.09311122e-02]\n",
      " [7.36472700e+00 4.27303850e+00 1.14300960e-02]\n",
      " [7.52749694e+00 4.36261587e+00 1.19242822e-02]\n",
      " [7.68501135e+00 4.37500000e+00 1.19934265e-02]\n",
      " [7.83721312e+00 4.37500000e+00 1.19934265e-02]\n",
      " [7.98405598e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.12550275e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.26152304e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.39209275e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.51719014e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.63679240e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.75088012e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.85944394e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.96248245e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.05999483e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.15198734e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.23848664e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.31953630e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.39518590e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.46549063e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.53051457e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.59033671e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.64506419e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.69483550e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.73981356e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.78019337e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.81621422e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.84814974e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.87629067e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.90093717e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.92239281e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.94095590e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.95691594e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.97055634e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.98214808e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.99193615e+00 4.37500000e+00 1.19934265e-02]\n",
      " [1.00001373e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00069478e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00125488e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00171084e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00207796e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00236981e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00259857e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00277532e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00290967e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00300960e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00308146e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00313089e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00316331e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00318357e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00319553e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00320206e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00320512e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00320623e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00320652e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.00320658e+01 4.37500000e+00 1.19934265e-02]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\elise\\AppData\\Local\\Temp\\ipykernel_6600\\2952692391.py:296: RuntimeWarning: divide by zero encountered in log10\n",
      "  srcv = np.log10(vel/0.02) #Velocity portion of the strain rate correction.\n"
     ]
    }
   ],
   "source": [
    "fileNum= '01BE' # write the bin file number you want to analyze (do not include 'bLog' or '.bin')\n",
    "a = 1 #what drop do you want to analyze? USER INPUT\n",
    "a = a-1\n",
    "\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)\n",
    "\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",
    "peakplot()\n",
    "#overviewplot()\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",
    "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 \n",
    "d = drop3\n",
    "accPick(dg, d)  # chooses what accelerometer to use\n",
    "acc3 = acc\n",
    "acc3g = accg\n",
    "acc3Name = accName\n",
    "acc3Nameg = accNameg\n",
    "#findchangepts() #automatically identifies the moment of penetration\n",
    "\n",
    "peak = int(peaksArray[a])  # count at the ath drop\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",
    "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 \n",
    "d = drop3\n",
    "accPick(dg, d)  # chooses what accelerometer to use\n",
    "acc3 = acc\n",
    "acc3g = accg\n",
    "acc3Name = accName\n",
    "acc3Nameg = accNameg\n",
    "#findchangepts() #automatically identifies the moment of penetration\n",
    "\n",
    "peakpick()\n",
    "\n",
    "finddropend(a) #automatically find the end of the drop\n",
    "offset = acc3g[penstart] #sets the offset to be the acceleration at the moment of impact\n",
    "#print(drop)\n",
    "dg = dfCalg[penstart:penend] #chops the file to only include penetration \n",
    "d = dfCal[penstart:penend]\n",
    "acc3g = accg[penstart:penend]\n",
    "acc3g = acc3g-offset\n",
    "acc3 = acc[penstart:penend]\n",
    "offsetms = offset * 9.80665\n",
    "acc3 = acc3-offsetms\n",
    "\n",
    "\n",
    "drop3 = d\n",
    "drop3g = dg\n",
    "integration(d, acc3) #double integration to solve for velocity and penetration depth\n",
    "drop3 = d\n",
    "areafind() #area calculations for bearing capacity calculations\n",
    "qdynfun(acc3) #calculation of dynamic bearing capacity\n",
    "qsbcfun(acc3) #calculation of quasi-static bearing capacities\n",
    "qsbc3 = qsbc\n",
    "qdyn3 = qdyn\n",
    "integplot(drop3,acc3Name)\n",
    "#qdynplot(drop1, qdyn1)\n",
    "drop3 = pd.DataFrame(drop3)\n",
    "drop3g = pd.DataFrame(drop3g)\n",
    "acc3 = pd.DataFrame(acc3)\n",
    "qdyn3 = pd.DataFrame(qdyn3)\n",
    "\n",
    "drop3t = drop3 #.truncate(before=trunc_index) #this truncates each df such that only depths below cone height are considered\n",
    "drop3gt = drop3g #.truncate(before=trunc_index)\n",
    "acc3t = acc3 #.truncate(before=trunc_index)\n",
    "acc3t = np.array(acc3t)\n",
    "acc3gt = acc3g #.truncate(before=trunc_index)\n",
    "acc3gt = np.array(acc3gt)\n",
    "drop3t = drop3 #truncate(before=trunc_index)\n",
    "qdyn3t = qdyn3#.truncate(before=trunc_index)\n",
    "\n",
    "duck_dep_qsbc_comboplot(drop3g, drop3, drop3t, acc3Nameg, qdyn3t)\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "trial 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak index: 101915\n",
      "start of penetration:  101866\n",
      "[[0.00000000e+00 2.20000000e-01 1.52053084e-05]\n",
      " [2.57757345e-01 3.61851890e-01 5.05041996e-05]\n",
      " [5.15175812e-01 5.03517284e-01 1.13951766e-04]\n",
      " [7.72205013e-01 6.44968454e-01 2.03863616e-04]\n",
      " [1.02878502e+00 7.86172420e-01 3.19790098e-04]\n",
      " [1.28483987e+00 9.27087375e-01 4.61459660e-04]\n",
      " [1.54027976e+00 1.06766390e+00 6.28614554e-04]\n",
      " [1.79500552e+00 1.20784741e+00 8.20972270e-04]\n",
      " [2.04891515e+00 1.34758178e+00 1.03821619e-03]\n",
      " [2.30191080e+00 1.48681316e+00 1.27999739e-03]\n",
      " [2.55389868e+00 1.62548993e+00 1.54593523e-03]\n",
      " [2.80478597e+00 1.76356102e+00 1.83561517e-03]\n",
      " [3.05448177e+00 1.90097639e+00 2.14859042e-03]\n",
      " [3.30290078e+00 2.03768911e+00 2.48438824e-03]\n",
      " [3.54996352e+00 2.17365542e+00 2.84251324e-03]\n",
      " [3.79558942e+00 2.30883100e+00 3.22243959e-03]\n",
      " [4.03968391e+00 2.44316380e+00 3.62358903e-03]\n",
      " [4.28213167e+00 2.57659034e+00 4.04531251e-03]\n",
      " [4.52280663e+00 2.70904126e+00 4.48689975e-03]\n",
      " [4.76158383e+00 2.84044779e+00 4.94759838e-03]\n",
      " [4.99833820e+00 2.97074109e+00 5.42661275e-03]\n",
      " [5.23294226e+00 3.09985100e+00 5.92310000e-03]\n",
      " [5.46526946e+00 3.22770789e+00 6.43617743e-03]\n",
      " [5.69519267e+00 3.35424179e+00 6.96492048e-03]\n",
      " [5.92258237e+00 3.47938142e+00 7.50835923e-03]\n",
      " [6.14730952e+00 3.60305577e+00 8.06548587e-03]\n",
      " [6.36924568e+00 3.72519415e+00 8.63525647e-03]\n",
      " [6.58826390e+00 3.84572669e+00 9.21659512e-03]\n",
      " [6.80423851e+00 3.96458424e+00 9.80839533e-03]\n",
      " [7.01704088e+00 4.08169601e+00 1.04095100e-02]\n",
      " [7.22654136e+00 4.19699064e+00 1.10187574e-02]\n",
      " [7.43261706e+00 4.31040051e+00 1.16349459e-02]\n",
      " [7.63515479e+00 4.37500000e+00 1.19934265e-02]\n",
      " [7.83404742e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.02919292e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.22049648e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.40786951e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.59122515e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.77047632e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.94554155e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.11634888e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.28283092e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.44492060e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.60255167e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.75565948e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.90418503e+00 4.37500000e+00 1.19934265e-02]\n",
      " [1.00480813e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.01873087e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.03218299e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.04516102e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.05766176e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.06968266e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.08122269e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.09228205e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.10286119e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11296053e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.12258142e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.13172643e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.14039864e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.14860158e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.15634037e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.16362157e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.17045169e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.17683763e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.18278790e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.18831253e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.19342241e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.19812892e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.20244392e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.20638019e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.20995248e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.21317813e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.21607633e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.21866718e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.22097154e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.22301102e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.22480734e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.22638131e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.22775273e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.22894094e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.22996482e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23084205e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23158903e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23222077e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23275099e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23319234e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23355659e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23385487e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23409718e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23429185e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23444592e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23456579e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23465727e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23472565e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23477568e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23481128e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23483547e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23485090e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23486016e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23486530e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23486759e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23486818e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.23486822e+01 4.37500000e+00 1.19934265e-02]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\elise\\AppData\\Local\\Temp\\ipykernel_6600\\2952692391.py:296: RuntimeWarning: divide by zero encountered in log10\n",
      "  srcv = np.log10(vel/0.02) #Velocity portion of the strain rate correction.\n"
     ]
    }
   ],
   "source": [
    "fileNum= '01F8' # write the bin file number you want to analyze (do not include 'bLog' or '.bin')\n",
    "a = 2 #what drop do you want to analyze? USER INPUT\n",
    "a = a-1\n",
    "\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)\n",
    "\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",
    "peakplot()\n",
    "#overviewplot()\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",
    "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 \n",
    "d = drop4\n",
    "accPick(dg, d)  # chooses what accelerometer to use\n",
    "acc4 = acc\n",
    "acc4g = accg\n",
    "acc4Name = accName\n",
    "acc4Nameg = accNameg\n",
    "#findchangepts() #automatically identifies the moment of penetration\n",
    "\n",
    "peak = int(peaksArray[a])  # count at the ath drop\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",
    "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 \n",
    "d = drop4\n",
    "accPick(dg, d)  # chooses what accelerometer to use\n",
    "acc4 = acc\n",
    "acc4g = accg\n",
    "acc4Name = accName\n",
    "acc4Nameg = accNameg\n",
    "#findchangepts() #automatically identifies the moment of penetration\n",
    "\n",
    "peakpick()\n",
    "\n",
    "finddropend(a) #automatically find the end of the drop\n",
    "offset = acc4g[penstart] #sets the offset to be the acceleration at the moment of impact\n",
    "#print(drop)\n",
    "dg = dfCalg[penstart:penend] #chops the file to only include penetration \n",
    "d = dfCal[penstart:penend]\n",
    "acc4g = accg[penstart:penend]\n",
    "acc4g = acc4g-offset\n",
    "acc4 = acc[penstart:penend]\n",
    "offsetms = offset * 9.80665\n",
    "acc4 = acc4-offsetms\n",
    "\n",
    "\n",
    "drop4 = d\n",
    "drop4g = dg\n",
    "integration(d, acc4) #double integration to solve for velocity and penetration depth\n",
    "drop4 = d\n",
    "areafind() #area calculations for bearing capacity calculations\n",
    "qdynfun(acc4) #calculation of dynamic bearing capacity\n",
    "qsbcfun(acc4) #calculation of quasi-static bearing capacities\n",
    "qsbc4 = qsbc\n",
    "qdyn4 = qdyn\n",
    "integplot(drop4,acc4Name)\n",
    "#qdynplot(drop1, qdyn1)\n",
    "drop4 = pd.DataFrame(drop4)\n",
    "drop4g = pd.DataFrame(drop4g)\n",
    "acc4 = pd.DataFrame(acc4)\n",
    "qdyn4 = pd.DataFrame(qdyn4)\n",
    "\n",
    "drop4t = drop4 #.truncate(before=trunc_index) #this truncates each df such that only depths below cone height are considered\n",
    "drop4gt = drop4g #.truncate(before=trunc_index)\n",
    "acc4t = acc4 #.truncate(before=trunc_index)\n",
    "acc4t = np.array(acc4t)\n",
    "acc4gt = acc4g #.truncate(before=trunc_index)\n",
    "acc4gt = np.array(acc4gt)\n",
    "drop4t = drop4 #truncate(before=trunc_index)\n",
    "qdyn4t = qdyn4#.truncate(before=trunc_index)\n",
    "\n",
    "duck_dep_qsbc_comboplot(drop4g, drop4, drop4t, acc4Nameg, qdyn4t)\n",
    "\n",
    "dep4 = dep\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "trial 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "peak index: 60350\n",
      "start of penetration:  60304\n",
      "[[0.00000000e+00 2.20000000e-01 1.52053084e-05]\n",
      " [2.72580056e-01 3.70009289e-01 5.34215405e-05]\n",
      " [5.44830716e-01 5.19837301e-01 1.22980473e-04]\n",
      " [8.16682537e-01 6.69445820e-01 2.22089324e-04]\n",
      " [1.08804672e+00 8.18785978e-01 3.50271767e-04]\n",
      " [1.35880239e+00 9.67791249e-01 5.07207124e-04]\n",
      " [1.62878662e+00 1.11637197e+00 6.92545557e-04]\n",
      " [1.89779297e+00 1.26441454e+00 9.05849722e-04]\n",
      " [2.16557905e+00 1.41178556e+00 1.14656888e-03]\n",
      " [2.43188175e+00 1.55834022e+00 1.41403305e-03]\n",
      " [2.69643192e+00 1.70393041e+00 1.70745882e-03]\n",
      " [2.95896941e+00 1.84841297e+00 2.02596612e-03]\n",
      " [3.21925780e+00 1.99165777e+00 2.36860452e-03]\n",
      " [3.47708930e+00 2.13355047e+00 2.73437453e-03]\n",
      " [3.73228287e+00 2.27399144e+00 3.12224111e-03]\n",
      " [3.98468320e+00 2.41289519e+00 3.53114727e-03]\n",
      " [4.23415352e+00 2.55018647e+00 3.96001599e-03]\n",
      " [4.48056341e+00 2.68579351e+00 4.40773816e-03]\n",
      " [4.72378579e+00 2.81964635e+00 4.87317132e-03]\n",
      " [4.96370190e+00 2.95167965e+00 5.35515310e-03]\n",
      " [5.20020209e+00 3.08183307e+00 5.85250844e-03]\n",
      " [5.43317879e+00 3.21004740e+00 6.36403945e-03]\n",
      " [5.66252056e+00 3.33626132e+00 6.88851422e-03]\n",
      " [5.88811697e+00 3.46041404e+00 7.42467816e-03]\n",
      " [6.10986577e+00 3.58244931e+00 7.97127423e-03]\n",
      " [6.32767094e+00 3.70231427e+00 8.52704298e-03]\n",
      " [6.54144148e+00 3.81995885e+00 9.09072354e-03]\n",
      " [6.75109402e+00 3.93533717e+00 9.66106435e-03]\n",
      " [6.95654322e+00 4.04840226e+00 1.02368005e-02]\n",
      " [7.15769416e+00 4.15910189e+00 1.08166322e-02]\n",
      " [7.35445265e+00 4.26738420e+00 1.13992528e-02]\n",
      " [7.54673229e+00 4.37320168e+00 1.19833735e-02]\n",
      " [7.73445117e+00 4.37500000e+00 1.19934265e-02]\n",
      " [7.91752929e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.09588962e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.26945916e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.43816957e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.60195563e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.76074847e+00 4.37500000e+00 1.19934265e-02]\n",
      " [8.91447603e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.06307059e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.20647118e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.34462429e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.47748396e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.60500580e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.72714177e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.84384954e+00 4.37500000e+00 1.19934265e-02]\n",
      " [9.95510158e+00 4.37500000e+00 1.19934265e-02]\n",
      " [1.00608812e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.01611827e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.02560161e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.03454028e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.04293663e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.05079327e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.05811478e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.06490877e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.07118474e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.07695340e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.08222733e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.08702155e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.09135356e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.09524322e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.09871261e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.10178522e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.10448653e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.10684484e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.10888992e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11065139e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11215822e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11343834e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11451797e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11542161e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11617203e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11679056e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11729717e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11770958e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11804276e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11830959e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11852133e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11868774e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11881704e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11891574e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11898923e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11904255e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11908037e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11910656e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11912408e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11913552e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11914300e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11914768e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11915010e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11915085e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11915083e+01 4.37500000e+00 1.19934265e-02]\n",
      " [1.11915076e+01 4.37500000e+00 1.19934265e-02]]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\elise\\AppData\\Local\\Temp\\ipykernel_6600\\2952692391.py:296: 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_6600\\2952692391.py:296: RuntimeWarning: invalid value encountered in log10\n",
      "  srcv = np.log10(vel/0.02) #Velocity portion of the strain rate correction.\n"
     ]
    }
   ],
   "source": [
    "fileNum= '023B' # write the bin file number you want to analyze (do not include 'bLog' or '.bin')\n",
    "a = 3 #what drop do you want to analyze? USER INPUT\n",
    "a = a-1\n",
    "\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)\n",
    "\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",
    "peakplot()\n",
    "#overviewplot()\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",
    "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 \n",
    "d = drop5\n",
    "accPick(dg, d)  # chooses what accelerometer to use\n",
    "acc5 = acc\n",
    "acc5g = accg\n",
    "acc5Name = accName\n",
    "acc5Nameg = accNameg\n",
    "#findchangepts() #automatically identifies the moment of penetration\n",
    "\n",
    "peak = int(peaksArray[a])  # count at the ath drop\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",
    "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 \n",
    "d = drop5\n",
    "accPick(dg, d)  # chooses what accelerometer to use\n",
    "acc5 = acc\n",
    "acc5g = accg\n",
    "acc5Name = accName\n",
    "acc5Nameg = accNameg\n",
    "#findchangepts() #automatically identifies the moment of penetration\n",
    "\n",
    "peakpick()\n",
    "\n",
    "finddropend(a) #automatically find the end of the drop\n",
    "offset = acc5g[penstart] #sets the offset to be the acceleration at the moment of impact\n",
    "#print(drop)\n",
    "dg = dfCalg[penstart:penend] #chops the file to only include penetration \n",
    "d = dfCal[penstart:penend]\n",
    "acc5g = accg[penstart:penend]\n",
    "acc5g = acc5g-offset\n",
    "acc5 = acc[penstart:penend]\n",
    "offsetms = offset * 9.80665\n",
    "acc5 = acc5-offsetms\n",
    "\n",
    "\n",
    "drop5 = d\n",
    "drop5g = dg\n",
    "integration(d, acc5) #double integration to solve for velocity and penetration depth\n",
    "drop5 = d\n",
    "areafind() #area calculations for bearing capacity calculations\n",
    "qdynfun(acc5) #calculation of dynamic bearing capacity\n",
    "qsbcfun(acc5) #calculation of quasi-static bearing capacities\n",
    "qsbc5 = qsbc\n",
    "qdyn5 = qdyn\n",
    "integplot(drop5,acc5Name)\n",
    "#qdynplot(drop1, qdyn1)\n",
    "drop5 = pd.DataFrame(drop5)\n",
    "drop5g = pd.DataFrame(drop5g)\n",
    "acc5 = pd.DataFrame(acc5)\n",
    "qdyn5 = pd.DataFrame(qdyn5)\n",
    "\n",
    "drop5t = drop5 #.truncate(before=trunc_index) #this truncates each df such that only depths below cone height are considered\n",
    "drop5gt = drop5g #.truncate(before=trunc_index)\n",
    "acc5t = acc5 #.truncate(before=trunc_index)\n",
    "acc5t = np.array(acc5t)\n",
    "acc5gt = acc5g #.truncate(before=trunc_index)\n",
    "acc5gt = np.array(acc5gt)\n",
    "drop5t = drop5 #truncate(before=trunc_index)\n",
    "qdyn5t = qdyn5#.truncate(before=trunc_index)\n",
    "\n",
    "duck_dep_qsbc_comboplot(drop5g, drop5, drop5t, acc5Nameg, qdyn5t)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'plt' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[1;32mc:\\Users\\elise\\Documents\\GitHub\\BD-Code-Development\\BD code - profile overlay.ipynb Cell 20\u001b[0m in \u001b[0;36m5\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20profile%20overlay.ipynb#X25sZmlsZQ%3D%3D?line=51'>52</a>\u001b[0m         ax2\u001b[39m.\u001b[39mset_title(\u001b[39m\"\u001b[39m\u001b[39mDrop #\u001b[39m\u001b[39m\"\u001b[39m\u001b[39m+\u001b[39m\u001b[39mstr\u001b[39m(a))\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20profile%20overlay.ipynb#X25sZmlsZQ%3D%3D?line=52'>53</a>\u001b[0m         ax2\u001b[39m.\u001b[39mgrid()\n\u001b[1;32m---> <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20profile%20overlay.ipynb#X25sZmlsZQ%3D%3D?line=56'>57</a>\u001b[0m combined_duck_dep_qsbc_comboplot()\n\u001b[0;32m     <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20profile%20overlay.ipynb#X25sZmlsZQ%3D%3D?line=58'>59</a>\u001b[0m plt\u001b[39m.\u001b[39mshow()\n",
      "\u001b[1;32mc:\\Users\\elise\\Documents\\GitHub\\BD-Code-Development\\BD code - profile overlay.ipynb Cell 20\u001b[0m in \u001b[0;36m2\n\u001b[0;32m      <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20profile%20overlay.ipynb#X25sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m \u001b[39mdef\u001b[39;00m \u001b[39mcombined_duck_dep_qsbc_comboplot\u001b[39m():\n\u001b[1;32m----> <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20profile%20overlay.ipynb#X25sZmlsZQ%3D%3D?line=1'>2</a>\u001b[0m     fig, (ax1, ax2) \u001b[39m=\u001b[39m plt\u001b[39m.\u001b[39msubplots(\u001b[39m1\u001b[39m,\u001b[39m2\u001b[39m)\n\u001b[0;32m      <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20profile%20overlay.ipynb#X25sZmlsZQ%3D%3D?line=2'>3</a>\u001b[0m     fig\u001b[39m.\u001b[39mset_size_inches(\u001b[39m7\u001b[39m, \u001b[39m7\u001b[39m)\n\u001b[0;32m      <a href='vscode-notebook-cell:/c%3A/Users/elise/Documents/GitHub/BD-Code-Development/BD%20code%20-%20profile%20overlay.ipynb#X25sZmlsZQ%3D%3D?line=3'>4</a>\u001b[0m     ax1\u001b[39m.\u001b[39mplot(drop1g[acc1Nameg], drop1[\u001b[39m\"\u001b[39m\u001b[39mPenetration Depth (m)\u001b[39m\u001b[39m\"\u001b[39m]\u001b[39m*\u001b[39m\u001b[39m100\u001b[39m, linestyle \u001b[39m=\u001b[39m \u001b[39m\"\u001b[39m\u001b[39msolid\u001b[39m\u001b[39m\"\u001b[39m, color \u001b[39m=\u001b[39m \u001b[39m\"\u001b[39m\u001b[39morange\u001b[39m\u001b[39m\"\u001b[39m) \u001b[39m#marker = 11\u001b[39;00m\n",
      "\u001b[1;31mNameError\u001b[0m: name 'plt' is not defined"
     ]
    }
   ],
   "source": [
    "def combined_duck_dep_qsbc_comboplot():\n",
    "    fig, (ax1, ax2) = plt.subplots(1,2)\n",
    "    fig.set_size_inches(7, 7)\n",
    "    fig.set_dpi(300)\n",
    "    ax1.plot(drop1g[acc1Nameg], drop1[\"Penetration Depth (m)\"]*100, linestyle = \"solid\", color = \"orange\") #marker = 11\n",
    "    ax1.plot(drop1[\"Velocity (m/s)\"], drop1[\"Penetration Depth (m)\"]*100, color = \"orange\", linestyle = \"dashed\")\n",
    "\n",
    "    ax1.plot(drop2g[acc2Nameg], drop2[\"Penetration Depth (m)\"]*100, linestyle = \"solid\", color = \"orangered\") #marker = 11\n",
    "    ax1.plot(drop2[\"Velocity (m/s)\"], drop2[\"Penetration Depth (m)\"]*100, linestyle = \"dashed\", color = \"orangered\")\n",
    "\n",
    "    ax1.plot(drop3g[acc3Nameg], drop3[\"Penetration Depth (m)\"]*100, linestyle = \"solid\", color = \"darkviolet\") #marker = 11\n",
    "    ax1.plot(drop3[\"Velocity (m/s)\"], drop3[\"Penetration Depth (m)\"]*100, linestyle = \"dashed\", color = \"darkviolet\")\n",
    "\n",
    "    ax1.plot(drop4g[acc4Nameg], drop4[\"Penetration Depth (m)\"]*100, linestyle = \"solid\", color = \"dodgerblue\") #marker = 11\n",
    "    ax1.plot(drop4[\"Velocity (m/s)\"], drop4[\"Penetration Depth (m)\"]*100, linestyle = \"dashed\", color = \"dodgerblue\")\n",
    "\n",
    "    ax1.plot(drop5g[acc5Nameg], drop5[\"Penetration Depth (m)\"]*100, linestyle = \"solid\", color = \"limegreen\") #marker = 11\n",
    "    ax1.plot(drop5[\"Velocity (m/s)\"], drop5[\"Penetration Depth (m)\"]*100, linestyle = \"dashed\", color = \"limegreen\")\n",
    "\n",
    "\n",
    "    ax1.set(xlabel=\"Deceleration (g) and Velocity (m/s)\", ylabel=\"Penetration Depth (cm)\")\n",
    "    ax1.set_xlim(left = 0)\n",
    "    ax1.set_ylim(0,dep4.max()*100+.25)\n",
    "    ax1.invert_yaxis()\n",
    "    ax1.legend([\"Deceleration (g)\", \"Velocity (m/s)\"])\n",
    "    ax1.set_title(\"Station Q3\")\n",
    "    ax1.grid()\n",
    "\n",
    "    \n",
    "    if droptype == \"a\": #for partial saturation, only one qsbc curve\n",
    "        y1 = drop1t[\"Penetration Depth (m)\"]*100\n",
    "        ax2.plot(qdyn1t[5:-50], y1[5:-50], label=\"Qdyn\", linestyle = \"dashed\", color = \"orange\") #color = \"k\", marker = 11 #Plots Qdyn, chops off first 40 readings\n",
    "        ax2.plot(qsbc1[5:-50], y1[5:-50], label=\"QSBC\", color = \"orange\") #plots QSBC; chops off last 10 points which go off towards infinity (can adjust as needed)\n",
    "        y2 = drop2t[\"Penetration Depth (m)\"]*100\n",
    "        ax2.plot(qdyn2t[5:-47], y2[5:-47], label=\"Qdyn\", linestyle = \"dashed\", color = \"orangered\") #color = \"k\", marker = 11 #Plots Qdyn, chops off first 40 readings\n",
    "        ax2.plot(qsbc2[5:-47], y2[5:-47], label=\"QSBC\", color = \"orangered\") #plots QSBC; chops off last 10 points which go off towards infinity (can adjust as needed)\n",
    "        y3 = drop3t[\"Penetration Depth (m)\"]*100\n",
    "        ax2.plot(qdyn3t[5:-50], y3[5:-50], label=\"Qdyn\", linestyle = \"dashed\", color = \"darkviolet\") #color = \"k\", marker = 11 #Plots Qdyn, chops off first 40 readings\n",
    "        ax2.plot(qsbc3[5:-50], y3[5:-50], label=\"QSBC\", color = \"darkviolet\") #plots QSBC; chops off last 10 points which go off towards infinity (can adjust as needed)\n",
    "        y4 = drop4t[\"Penetration Depth (m)\"]*100\n",
    "        ax2.plot(qdyn4t[5:-50], y4[5:-50], label=\"Qdyn\", linestyle = \"dashed\", color = \"dodgerblue\") #color = \"k\", marker = 11 #Plots Qdyn, chops off first 40 readings\n",
    "        ax2.plot(qsbc4[5:-50], y4[5:-50], label=\"QSBC\", color = \"dodgerblue\") #plots QSBC; chops off last 10 points which go off towards infinity (can adjust as needed)\n",
    "        y5 = drop5t[\"Penetration Depth (m)\"]*100\n",
    "        ax2.plot(qdyn5t[5:-50], y5[5:-50], label=\"Qdyn\", linestyle = \"dashed\", color = \"limegreen\") #color = \"k\", marker = 11 #Plots Qdyn, chops off first 40 readings\n",
    "        ax2.plot(qsbc5[5:-50], y5[5:-50], label=\"QSBC\", color = \"limegreen\") #plots QSBC; chops off last 10 points which go off towards infinity (can adjust as needed)\n",
    "        \n",
    "        \n",
    "        ax2.set(xlabel=\"Bearing Capacity (kPa)\",)\n",
    "        ax2.set_xlim(left = 0)\n",
    "        ax2.set_ylim(0,dep4.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",
    "    \n",
    "\n",
    "combined_duck_dep_qsbc_comboplot()\n",
    "\n",
    "plt.show()\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
}