Mann Morrision Python Software Simulation

import math
import numpy as np
import matplotlib.pyplot as plt

pi = 3.14159265359

f=0.25
Over_samp_factor =100*f

Vm = 5
Phi_v=30

R_perKm =.01
XbyR =10
Line_Length = 100

R = R_perKm * Line_Length
X = R * XbyR
w = 2 * pi * f
L = X / w
Z = math.sqrt(pow(R,2) + pow(X,2))
Im = Vm / Z

Phi_i = (math.atan2(XbyR, 1))*(180/pi)

fsamp = 2* f * Over_samp_factor
dT = 1/ fsamp
T = 1 / f
k=0

v = []
I = []
v_das = []

for t in np.arange(0, T+dT, dT):
    a = Vm * math.sin(2 * pi * f * t + (Phi_v)*(pi/180))
    v.append(a)

no_of_samples = len(v)

Vm_estimate = []
Phi_V_estimate = []
time = []
start = 0

print(v)

for k in range(1, no_of_samples-1):
    v_dash = (v[k+1]-v[k-1])/(2*dT)

    Vm_est= math.sqrt(pow(v[k],2) + pow(v_dash/w,2))

    Vm_estimate.append(round(Vm_est, 8))

    Phi_V_est_rad = (math.atan2(((v[k])/((v_dash)/w)),1))-(w*k*dT);
    Phi_V_est_deg= (Phi_V_est_rad)*(180/pi);

    if(Phi_V_est_deg<0):
        pass
        # Phi_V_est_deg = ((-Phi_V_est_deg)//180+1)*180 + Phi_V_est_deg;
    Phi_V_est_deg = (Phi_V_est_deg)%180;

    Phi_V_estimate.append(Phi_V_est_deg)

    time.append(start)
    start = start + dT

xPoints = np.array(time)
yPoints = np.array(Vm_estimate)

plt.plot(xPoints, yPoints, color='r', label='VPeak')

yPoints = np.array(Phi_V_estimate)
plt.plot(xPoints, yPoints, color='g', label='Phase')

plt.legend(loc='best')

plt.show()