centerofmass.py

import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as integrate

def CubicNatural(x, y):
    m = x.size # m is the number of data points
    n = m-1
    a = np.zeros(m)
    b = np.zeros(n)
    c = np.zeros(m)
    d = np.zeros(n)
    for i in range(m):
        a[i] = y[i]
    h = np.zeros(n)
    for i in range(n):
        h[i] = x[i+1] - x[i]
    u = np.zeros(n)
    u[0] = 0
    for i in range(1, n):
        u[i] = 3*(a[i+1]-a[i])/h[i]-3*(a[i]-a[i-1])/h[i-1]
    s = np.zeros(m)
    z = np.zeros(m)
    t = np.zeros(n)
    s[0] = 1
    z[0] = 0
    t[0] = 0
    for i in range(1, n):
        s[i] = 2*(x[i+1]-x[i-1])-h[i-1]*t[i-1]
        t[i] = h[i]/s[i]
        z[i]=(u[i]-h[i-1]*z[i-1])/s[i]
    s[m-1] = 1
    z[m-1] = 0
    c[m-1] = 0
    for i in np.flip(np.arange(n)):
        c[i] = z[i]-t[i]*c[i+1]
        b[i] = (a[i+1]-a[i])/h[i]-h[i]*(c[i+1]+2*c[i])/3
        d[i] = (c[i+1]-c[i])/(3*h[i])
    return a, b, c, d

def CubicNaturalEval(w, x, coeff):
    m = x.size
    if w<x[0] or w>x[m-1]:
        print('error: spline evaluated outside its domain')
        return
    n = m-1
    p = 0
    for i in range(n):
        if w <= x[i+1]:
            break
        else:
            p += 1
    # p is the number of the subinterval w falls into, i.e., p=i means
    # w falls into the ith subinterval $(x_i,x_{i+1}), and therefore
    # the value of the spline at w is
    # a_i+b_i*(w-x_i)+c_i*(w-x_i)^2+d_i*(w-x_i)^3.
    a = coeff[0]
    b = coeff[1]
    c = coeff[2]
    d = coeff[3]
    return a[p]+b[p]*(w-x[p])+c[p]*(w-x[p])**2+d[p]*(w-x[p])**3

xaxis = np.linspace(0, 14.875, 100)

xi = np.array([0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5, 10, 10.5, 11, 11.5, 12, 12.5, 13, 13.5, 14, 14.5, 14.875])  # Replace with your data points
yi = np.array([0.44, 1.5, 2.875, 4.125, 5.5, 6.875, 8.25, 9.25, 9.375, 9.25, 8.75, 8, 7.375, 6.875, 6.375, 6.375, 6.625, 7.25, 8, 8.25, 7.5, 5.625, 5.25, 5, 5, 5.125, 5.625, 6.375, 7, 7.5, 8.25])  # Replace with your function values at x_data points

for i in CubicNatural(xi,yi):
    print(i)

coeff = CubicNatural(xi, yi)

def f(x):
    return CubicNaturalEval(x, xi, coeff)

def g(x):
    return (f(x)**2)/2

def h(x):
    return (x*f(x))

cuts = 30

#----Trapezoidal Rule Integration----

part = np.linspace(0,14.875,cuts)
fpart = np.zeros(cuts)
gpart = np.zeros(cuts)
hpart = np.zeros(cuts)

for i in range(cuts):
    fpart[i] = f(part[i])
    gpart[i] = g(part[i])
    hpart[i] = h(part[i])


print("Sfdx = ", np.trapz(fpart))
print("Sgdx = ", np.trapz(gpart))
print("Shdx = ", np.trapz(hpart))

print("x = ", np.trapz(hpart)/np.trapz(fpart))
print("y = ", np.trapz(gpart)/np.trapz(fpart))

#------------------------------------


#Alternative version below

# Sf, err = integrate.quad(f, 0, 14.875)
# Sg, err = integrate.quad(g, 0, 14.875)
# Sh, err = integrate.quad(h, 0, 14.875)

# print("Sf(x)dx=", Sf)
# print("Sg(x)dx=", Sg)
# print("Sh(x)dx=", Sh)

# print("x = ", Sh/Sf)
# print("y = ", Sg/Sf)