<jupyter_start><jupyter_text>RPS Nanopore Particle Sizing Algorithm  <NAME>  * Importing Modules<jupyter_code>import numpy as np
import pandas as pd
import scipy
import scipy.integrate as integrate
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d<jupyter_output><empty_output><jupyter_text>* Defining a test sigma value<jupyter_code>test_sig = 1.264<jupyter_output><empty_output><jupyter_text>*Part 1 - Calculating Resistance Changes*  * Function for calculating B Values<jupyter_code>def calc_b_values(rtip, xpart, rpart, a):
    """Calculates B values for Resistance Change calculations from 
    Gyurcsanyi Analytical Chem 2014, 86(10), 4688-4697 SI
    
    inputs: pipet tip radius, position of particle center, 
    particle radius, neccessary angle
    
    outputs: B values in an array (b3 is squared)
    -------------------------------------------------------------------------"""
    # Degrees to Rad
    angle = a * np.pi/180
    # Calculating B values 
    b1 = (1-np.cos(angle)**2)*xpart**2 - np.sin(2*angle)*rtip*xpart+np.cos(angle)**2 * rtip**2 - rpart**2
    b2 = (np.cos(angle)**2 - 1)*xpart + np.cos(angle)*np.sin(angle)*rtip
    b3_sq = (np.tan(angle)**2+1)*rpart**2 - (np.tan(angle)*xpart - rtip)**2
    b4 = np.tan(angle)**2*xpart - np.tan(angle)*rtip - (np.tan(angle)**2 + 1)*rpart
    b5 = xpart = np.tan(angle)*rtip
    return [b1, b2, b3_sq, b4, b5]
calc_b_values(110e-9, -100e-9, 50e-9, 8)<jupyter_output><empty_output><jupyter_text>* Function for calculating delta R-partial<jupyter_code>def calc_delta_r_part(rtip, xpart, rpart, angle, sigma):
    """Calculating delta r for the particle being partially inside of the pipette
    
    inputs: pipet tip radius, position of particle center, 
    particle radius, cone angle, bulk conductivity value
    
    output: delta R-part value"""
    ang = angle * np.pi/180
    beta= 54.8 * np.pi/180 
    
    # Negating xpart for pipette calculations 
    xpp = -1 * xpart
    
    # Calculating B Values (access)
    bvals = calc_b_values(rtip, xpart, rpart, beta)
    b1 = bvals[0]
    b2 = bvals[1]
    b3_sqr = bvals[2]
    b3 = np.sqrt(b3_sqr)
    b4 = bvals[-2]
    b5 = bvals[-1]
    
    # Calculating B Values (pipette)
    bvals_p = calc_b_values(rtip, xpp, rpart, angle)
    b1p = bvals_p[0]
    b2p = bvals_p[1]
    b3_sqr_p = bvals_p[2]
    b3p = np.sqrt(b3_sqr)
    b4p = bvals_p[-2]
    b5p = bvals_p[-1]
    
    # Calculating Ro access
    num1 = rpart - xpart
    denom1 = sigma * np.pi * rtip * (np.tan(beta)*(num1) + rtip)
    Roa = num1 / denom1
    
    # Calculating Ro pipette
    num2 = rpart - xpp
    denom2 = sigma * np.pi * rtip * (np.tan(ang)*(num2) + rtip)
    Rop = num2 / denom2
    
    # Calculating R access
    if b3_sqr < 0:
        fx1 = np.cos(beta)*xpart
        fx2 = np.sin(beta)*rtip
        fx3 = np.arctan((fx1+fx2)/np.sqrt(b1))
        fx4 = np.cos(beta)*np.sqrt(b1) 
        fx5 = np.arctan((b2+rpart)/fx4)
        num = -1 *np.cos(beta)*(fx3 - fx5)
        denom = np.pi * sigma * np.sqrt(b1)
        Ra = num / denom 
    if b3_sqr > 0: 
        fx1 = np.abs(b3 + b4)
        fx2 = np.abs(b3 - b4)
        fx3 = np.abs(b3 + b5)
        fx4 = np.abs(b3 - b5)
        num = np.log((fx1*fx3)/(fx2*fx4))
        denom = 2 * sigma * np.pi * b3
        Ra = num / denom
    else:
        pass
    
    # Calculating R pipette
    if b3_sqr_p < 0:
        fx1 = np.cos(ang)* xpp
        fx2 = np.sin(ang)*rtip
        fx3 = np.arctan((fx1+fx2)/np.sqrt(b1p))
        fx4 = np.cos(ang)*np.sqrt(b1p) 
        fx5 = np.arctan((b2p+rpart)/fx4)
        num = -1 *np.cos(ang)*(fx3 - fx5)
        denom = np.pi * sigma * np.sqrt(b1p)
        Rp = num / denom 
    if b3_sqr_p > 0: 
        fx1 = np.abs(b3p + b4p)
        fx2 = np.abs(b3p - b4p)
        fx3 = np.abs(b3p + b5p)
        fx4 = np.abs(b3p - b5p)
        num = np.log((fx1*fx3)/(fx2*fx4))
        denom = 2 * sigma * np.pi * b3p
        Rp = num / denom
    else:
        pass
    
    # Calculating delta R
    Ro = Roa + Rop
    R = Ra + Rp
    delta_R = R - Ro
    return delta_R
calc_delta_r_part(110e-9, -50e-9, 50e-9, 8, test_sig)<jupyter_output>/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:19: RuntimeWarning: invalid value encountered in sqrt
/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in sqrt<jupyter_text>* Function for calculating Delta R-out<jupyter_code>def calc_delta_r_out(rtip, xpart, rpart, angle, sigma):
    """Calculating delta r for the particle being completely outside/inside of the pipette
    
    inputs: pipet tip radius, position of particle center, 
    particle radius, cone angle, bulk conductivity value
    
    output: delta R-out value"""
    ang = angle * np.pi/180
    
    # Calculating B Values
    bvals = calc_b_values(rtip, xpart, rpart, angle)
    b1 = bvals[0]
    b2 = bvals[1]
    
    # Setting limits of integration for Ro
    low_lim = xpart - rpart
    up_lim = xpart + rpart
    
    # Calculating Ro
    a = rtip
    b = np.tan(ang)
    def f(x):
        return (np.pi * sigma * (a - b*x)**2)**-1
    Ro, err = integrate.quad(f, low_lim, up_lim)
    
    # Calculating R
    fx1 = np.arctan((b2-rpart)/(np.cos(ang)*np.sqrt(b1)))
    fx2 = np.arctan((b2+rpart)/(np.cos(ang)*np.sqrt(b1)))
    num = np.cos(ang) * (fx1 - fx2)
    den = sigma * np.pi * np.sqrt(b1)
    R = - num / den
    
    # Calculating delta R
    delta_R = R - Ro
    return delta_R
calc_delta_r_out(110e-9, -100e-9, 50e-9, 54.8, test_sig)<jupyter_output><empty_output><jupyter_text>* Function for calculating Delta R based on a given particle size<jupyter_code>def calc_delta_r(rtip, xpart, rpart, alpha, sigma): 
    """Uses B values from previous function to calculate delta R for a given set of parameters 
    
    inputs: pipet tip radius, position of particle center, 
    particle radius, cone angle, bulk conductivity value
    
    output: delta R value
    -------------------------------------------------------------------------------"""
    beta = 54.8
    if xpart <= -rpart:
        delta_R = calc_delta_r_out(rtip, xpart, rpart, beta, sigma) 
    if -rpart <= xpart <= rpart:
        delta_R = calc_delta_r_part(rtip, xpart, rpart, alpha, sigma)
    if xpart >= rpart:
        delta_R = calc_delta_r_out(rtip, -1 * xpart, rpart, alpha, sigma)
    else:
        pass
    return delta_R
calc_delta_r(110e-9, 49e-9, 50e-9, 8, test_sig)<jupyter_output>/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:19: RuntimeWarning: invalid value encountered in sqrt
/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in sqrt<jupyter_text>* Function to convert Delta R values calculated above into R-peak*<jupyter_code>def calc_R_peak(rtip, xpart, rpart, alpha, sigma):
    """Converts delta R value into a "delta I over I" value than can be
    compared to an experimental RPS response
    
    inputs: pipet tip radius, position of particle center, 
    particle radius, cone angle, bulk conductivity value
    
    output: delta I over I value"""
    angle = alpha * np.pi/180
    
    # Calculating Delta R
    dR = calc_delta_r(rtip, xpart, rpart, alpha, sigma)
    
    # Calculating R 
    dtip = 2 * rtip 
    term1 = 2 / (np.pi * sigma * dtip * np.tan(angle))
    term2 = 1 / (2 * dtip * sigma)
    R = term1 + term2
    
    # Calculating Current Blockage 
    CB = - dR / (R + dR)
    
    # Calculating R-peak
    R_peak = - R * CB / (CB + 1)
    
    # Calculating R-peak*
    R_peak_f = sigma * dtip * R_peak
    
    return R_peak_f
calc_R_peak(110e-9, 50e-9, 100e-9, 12, test_sig)<jupyter_output>/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:19: RuntimeWarning: invalid value encountered in sqrt
/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in sqrt<jupyter_text>* Function to vary xpart and find the maximum Rpeak value for the translocation event<jupyter_code>def find_max_R_peak(rtip, rpart, alpha, sigma):
    """Varies xpart argument in calc_R_peak 
    and exports the maximum R-peak* value
    
    Inputs: tip radius, particle radius, 
    cone angle, bulk conductivity value
    
    Output: maximum R-peak* value 
    -----------------------------------------------------------"""
    # Setting range for xpart
    lower_x = -2 * rpart
    upper_x = 2 * rpart 
    xrange = np.arange(start = lower_x, stop = upper_x, step = 1e-9)
    
    # Creating empty array for Rpeak* values
    rpeak_arr = []
    
    # Iterating through xrange to calculate Rpeak* values 
    # at different particle positions
    for i in xrange:
        rpeak_val = calc_R_peak(rtip, i, rpart, alpha, sigma)
        rpeak_arr.append(rpeak_val)
    
    rp_f = np.nanmax(rpeak_arr)
    return rp_f
find_max_R_peak(110e-9, 100e-9, 12, test_sig)<jupyter_output>/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:19: RuntimeWarning: invalid value encountered in sqrt
/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in sqrt<jupyter_text>*Part 2 - Simulation*<jupyter_code>def simulate_variable_particle_size(rtip, alpha, sigma): 
    """Finds the max Rpeak* values for a range of particle sizes"""
    # Creating axis 
    ratios = np.arange(start = 0.01, stop = 0.9, step = 0.01)
    rparts = ratios * rtip 
    
    # Creating empty list to store values 
    rpeak_list = []
    
    # Generating Axis values
    for i in rparts:
        rpeak_val = find_max_R_peak(rtip, i, alpha, sigma)
        rpeak_list.append(rpeak_val)
    
    # Turning list of Rpeak values into an array
    rpeaks = np.array(rpeak_list)
    return rpeaks
simulate_variable_particle_size(110e-9, 3, 1.0)
def simulate_variable_cone_angle(rtip, sigma):
    """Finds Rpeak vs dpart/dtip relationship for a range of cone angles"""
    # Creating angle array
    ang_arr = np.arange(start = 3, stop = 12, step = 0.5)
    
    # Creating ratio axis
    ratios = np.arange(start = 0.01, stop = 0.9, step = 0.01)
    
    # Creating empty list to store results
    rpeaks_list = []
    
    for ang in ang_arr: 
        rp_vals = simulate_variable_particle_size(110e-9, ang, test_sig)
        rpeaks_list.append(rp_vals)
        
    # Turning List into an array
    rpeaks_arr = np.array(rpeaks_list)
    return ang_arr, ratios, rpeaks_arr
ang, size, rpeak = simulate_variable_cone_angle(110e-9, test_sig)<jupyter_output>/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:19: RuntimeWarning: invalid value encountered in sqrt
/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in sqrt<jupyter_text>*Part 3 - Plotting*Starts off by checking axis dimensions before doing actual plotting function<jupyter_code>def plot_Rpeak_surface(rtip, sigma):
    """Uses simulate_variable_angle to perform full simulation, 
    stores output arrays as local variabls used for plotting
    
    inputs: tip radius, bulk conductivity value
    
    outputs: plotted delta R surface"""
    # Simulating
    ang, size, rpeak = simulate_variable_cone_angle(rtip, sigma)
    
    # Creating 2-D coordinate grid using angle and size as axes
    X, Y = np.meshgrid(size, ang)
    
    # Plotting
    # Y and X are flipped to match axes in Gyurcsanyi SI
    fig = plt.figure(figsize=(12,8), dpi= 200)
    ax = plt.axes(projection='3d')
    ax.plot_surface(Y, X, rpeak, cmap = 'viridis')
    ax.set_xlabel('Cone Angle (Degrees)')
    ax.set_ylabel('rpart/rtip')
    ax.set_zlabel('Rpeak*')
    ax.set_title("Simulated Rpeak* values (rtip = "+ 
                 str(rtip) + ", sigma = " + str(sigma) + ")")
    
    return ang, size, rpeak
plot_Rpeak_surface(110e-9, test_sig)<jupyter_output>/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:19: RuntimeWarning: invalid value encountered in sqrt
/home/mcalli75/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in sqrt<jupyter_text>*Part 4 - Curve Fitting*<jupyter_code>X, Y = np.meshgrid(size, ang)
Z = rpeak

X = X.flatten()
Y = Y.flatten()

A = np.array([X**0*Y**0, X, Y, X**2, X*Y, Y**2, X**3, X**2*Y, X*Y**2, 
             Y**3, X**4, X**3*Y, X**2*Y**2, X*Y**3, Y**4, X**5, X**4*Y,
             X**3*Y**2, X**2*Y**3, X*Y**4, Y**5, X**6, X**5*Y, X**4*Y**2,
             X**3*Y**3, X**2*Y**4, X*Y**5, Y**6, X**7, X**6*Y, X**5*Y**2,
             X**4*Y**3, X**3*Y**4, X**2*Y**5, X*Y**6, Y**7, X**8,
             X**7*Y, X**6*Y**2, X**5*Y**3, X**4*Y**4, X**3*Y**5,
             X**2*Y**6, X*Y**7, Y**8]).T
B = Z.flatten()

coeff, r, rank, s = np.linalg.lstsq(A, B)
coeff<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1376367445108232, 'fraction_numerical': 0.03669795950034913, 'mean_word_length': 2.7778429073856974, 'pattern_counts': {'":': 0, '<': 45, '<?xml version=': 0, '>': 43, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '11', 'dir': 'jupyter-structured-clean-dedup', 'id': '618394', 'license': 'MIT', 'n_tokens_mistral': 5191, 'n_tokens_neox': 4822, 'n_words': 1172, 'path': 'RPSParticleSizing/Real Time RNP Particle Sizing Algorithm Development Notebook - 9.3.20.ipynb', 'repo_name': 'cmcalli716/RPSParticleSizing', 'size': '12899'}