import numpy as np
import matplotlib.pyplot as plt

c = 1
a = 1
A = 1
B = 1
phi_1 = 0 
phi_2 = 0

t = np.linspace(0,2*np.pi,1000); #Timer Vector
y = A*np.sin(np.sqrt(c/a)*t+phi_1) + B*np.cos(np.sqrt(c/a) + phi_2) #y(t)
d2y_dt2 = np.diff(np.diff(y))/(0.005)**2 #second derivative of y(t)

#Plot results
plt.grid(True)
plt.plot(t,y,label='$y(t)$')
plt.plot(t[2:],d2y_dt2,label='$\ddot{y}(t)$')
plt.legend()
plt.show()


# Coefficients setup
points = 5000
a = 100.
b = 10.
c = 100.

# Time vector for the solution y[n] filled with zeros.
y = np.zeros((points,1))

# Initial conditions and sampling time
dt = 0.01;
y[0] = 50; #y[0]
dydt = 0;  #z[0]

#Simulation
for n in xrange(1,points):
	dydt = ((a/dt)*(dydt) - c*y[n-1]) / ((a/dt) + b + c*dt)
	y[n] = y[n-1] + dt*dydt

#Plot results
plt.grid(True)    
plt.plot(np.linspace(0,dt*points,points),y,label='y(t)')
plt.legend()
plt.show()

# Coefficients setup
points = 5000
a = 100.
b = 10.
c = 100.

# Time vector for the solution y[n] filled with zeros.
y = np.zeros((points,1))

# Initial conditions and sampling time
dt = 0.01;
y[0] = 50; #y[0]
dydt = 0;  #z[0]

# Impulses
f = np.zeros((points,1))
for i in xrange(2):
    f[i*2500] = 50
    
#Simulation
for n in xrange(1,points):
	dydt = f[n] + ((a/dt)*(dydt) - c*y[n-1]) / ((a/dt) + b + c*dt)
	y[n] = y[n-1] + dt*dydt

#Plot results
plt.grid(True)    
plt.plot(np.linspace(0,dt*points,points),y,label='y(t)')
plt.plot(np.linspace(0,dt*points,points),f,label='F(t)')
plt.legend()
plt.show()

# Coefficients setup
points = 5000
a = 100.
b = 10.
c = 100.

# Time vector for the solution y[n] filled with zeros.
y = np.zeros((points,1))

# Initial conditions and sampling time
dt = 0.01;
y[0] = 50; #y[0]
dydt = 0;  #z[0]

# Step
f = np.zeros((points,1))
for i in xrange(2):
    if (i+1)%2  == 0:
        f[i*2500:(i+1)*2500] = 0.5;

#Simulation
for n in xrange(1,points):
	dydt = f[n] + ((a/dt)*(dydt) - c*y[n-1]) / ((a/dt) + b + c*dt)
	y[n] = y[n-1] + dt*dydt

#Plot results
plt.grid(True)    
plt.plot(np.linspace(0,dt*points,points),y,label='y(t)')
plt.plot(np.linspace(0,dt*points,points),f*100,label='F(t)')
plt.legend()
plt.show()

# Coefficients setup
points = 10000 #0.1 seconds with dt = 0.00001
L = 10e-3 #10 mH
R = 30   #30  ohm
C = 10e-6 #10 uF

# Time vector for the solution Q[n] and I[n] filled with zeros.
Q = np.zeros((points,1))
I = np.zeros((points,1))


# Initial conditions and sampling time
dt = 0.00001; # 10 us
Q[0] = 0;  #y[0]
I[0] = 0;  #z[0]

# Source voltage
t = linspace(0,0.1,10000)
f = 500 # natural frequency 500Hz 
Vs = 10*np.cos(2*np.pi*f*t)


#Simulation
for n in xrange(1,points):
	I[n] = (Vs[n] + (L/dt)*(I[n-1]) - (1/C)*Q[n-1]) / ((L/dt) + R + (1/C)*dt)
	Q[n] = Q[n-1] + dt*I[n]

#Voltage on the resistor
VR = R*I

#Plot results
plt.figure()
plt.grid(True)    
plt.ylim(-15,15)
plt.plot(np.linspace(0,dt*points,points),Vs,'b-',label='$V_s(t)$')
plt.legend()
plt.show()

plt.grid(True)    
plt.ylim(-15,15)
plt.plot(np.linspace(0,dt*points,points),VR,'g-',label='$V_R(t)$')
plt.legend()
plt.show()

def gaussian(x, mu, sig):
    return (1/np.sqrt(2*np.pi*sig**2))*np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))

t = linspace(0,4,1000);

plt.grid(True)
plt.ylim(0,150);
width = 0.04;
plt.plot(t,10*gaussian(t,2,width),'g-',linewidth=2);

width = 0.1;
plt.plot(t,10*gaussian(t,2,width),'r-',linewidth=2,label='must be absorbed');

width = 1;
plt.plot(t,10*gaussian(t,2,width),'g-',linewidth=2,label='must not be absorbed');

plt.legend()

# Coefficients setup
points = 1000 #10 seconds with dt = 0.01
m = 30   #10 mH
k = 20   #30  ohm
c = 40   #10 uF

# Time vector for the solution Q[n] and I[n] filled with zeros.
x = np.zeros((points,1))
v = np.zeros((points,1))


# Initial conditions and sampling time
dt = 0.01; # 10 us
x[0] = 0;  #y[0]
v[0] = 0;  #z[0]

# Source voltage
t = linspace(0,10,1000)
width = 0.1
F = 10*gaussian(t,2,width)

#Simulation
for n in xrange(1,points):
	v[n] = (F[n] + (m/dt)*(v[n-1]) - (k)*x[n-1]) / ((m/dt) + c + (k)*dt)
	x[n] = x[n-1] + dt*v[n]


#Plot results
plt.figure()
plt.grid(True)    
plt.ylim(0,150)
plt.plot(np.linspace(0,dt*points,points),F,'b-',label='F(t)$',linewidth=2)
plt.legend()
plt.show()

plt.grid(True)    
#plt.ylim(-15,15)
plt.plot(np.linspace(0,dt*points,points),x,'g-',label='$V_R(t)$')
plt.legend()
plt.show()