import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

def divide(x,y,window=250):
    """
    Takes a data vector and a window to produce a vector of the running average ratio.
    The window specifies the number of points on either side of the central point, so
    the total number of points in the average is 2*window+1. Dividing is 
    done by the sum-division is defined by the equation below.
        sum(x)/n
    f = --------
        sum(y)/n
    """
    x = np.array(x)
    y = np.array(y)
    
    # Check that x and y are the same length
    if len(x) != len(y): 
        print ("Error: x and y must be the same length")
        return 0     
    
    N = len(x) # Number of points in the dataset
    fraction = np.ones(N) # Make array for slopes
    
    # Pad data with window number of points NaN on either side
    x_padded = np.empty(2*window+N)
    x_padded[0:window] = 0
    x_padded[window:N+window] = x
    x_padded[N+window:2*window+N] = 0
    
    y_padded = np.empty(2*window+N)
    y_padded[0:window] = 0
    y_padded[window:N+window] = y
    y_padded[N+window:2*window+N] = 0
    
    sum_x    = np.sum(x_padded[0:2*window+1])
    sum_y    = np.sum(y_padded[0:2*window+1])

    n = np.empty(N)
    n[0:window] = np.arange(window+1,2*window+1)
    n[window:N-window] = window*2+1
    n[N-window:N] = np.arange(2*window,window,-1)
    
    fraction[0] = (sum_x/n[0])/(sum_y/n[0])
    
    for i in range(1,N):
        sum_x    = sum_x - x_padded[i-1] + x_padded[2*window+i]
        sum_y    = sum_y - y_padded[i-1] + y_padded[2*window+i]
        fraction[i] = (sum_x/n[i])/(sum_y/n[i])
    return fraction
import os
filepath = os.getcwd()
data = pd.read_csv(os.path.join(filepath, "plate.txt"), sep=r"\s+", engine="python")

data["lx"] = 0.3 * (data["time"] - data["time"].iloc[0])

window_physical = 0.000001
dx = np.mean(np.diff(data["lx"]))
window_points = max(1, int(window_physical / dx))

data["mu"] = divide(-data["sf"], data["nf"])

plt.figure()
ax = plt.gca()

ax.plot(data["lx"], data["mu"], label="mu")

ax.set_xlabel("lx")
ax.set_ylabel("mu")

lx_max = data["lx"].max()
line_pos = 0.5 * lx_max

ax.axvline(x=line_pos, color="red", linestyle="--", alpha=0.6, linewidth=2, label="9/10")

ax.legend()

plt.savefig("mu_plot.png", dpi=300)

plt.figure()
plt.plot(data["lx"], data["thk"])
plt.xlabel("lx")
plt.ylabel("thk")

lx_max = data["lx"].max()
line_pos = 0.5 * lx_max

plt.axvline(x=line_pos, color="red", linestyle="--", alpha=0.6, linewidth=2)

plt.savefig("thk_plot.png", dpi=300)