%matplotlib inline

from math import sin,log,sqrt,pi,exp
import numpy as np
from pylab import imshow,show,gray

q1 = 1.0            # Electric charge 1, in C
q2 = 1.0            # Electric charge 2, in C
eps0 = 8.854E-12    # Vacuum permitivity, in F/m**2
dx = 0.001          # sample distance in x, in m
dy = 0.001          # sample distance in y, in m
Dx = 1.0            # grid size in x, in m
Dy = 1.0            # grid size in y, in m
rx1 = -0.05         # position of q1 in the axis x, in m
ry1 = 0.0           # position of q1 in the axis y, in m
rx2 = 0.05          # position of q2 in the axis x, in m
ry2 = 0.0           # position of q2 in the axis y, in m

def phi(q,r):
    return q/(4*pi*eps0*r)
    
Nx = int(Dx/dx)     # number of samples in x
Ny = int(Dy/dy)     # number of samples in y
EP = np.zeros([Nx,Ny])
logEP = np.zeros([Nx,Ny])       # using logarithmic scale makes the plot more eye-appealing... sometimes

ry = Dy/2           # starting position in the axis y: bottom         
for k in range(Ny):
    rx = -Dx/2      # starting position in the axis x: left
    for i in range(Nx):
        r1 = sqrt((rx1-rx)**2+(ry1-ry)**2)   # Distance from charge 1
        r2 = sqrt((rx2-rx)**2+(ry2-ry)**2)   # Distance from charge 2
        if abs(r1) < 1E-10 or abs(r2) < 1E-10:
            EP[k][i] = None
            logEP[k][i] = None
        else:
            EP[k][i] = phi(q1,r1)+phi(q2,r2)
            if abs(EP[k][i]) > 1E-10 and EP[k][i]>0:        # This crap here is to deal with the log of
                logEP[k][i] = log(EP[k][i])                 # negative values of electric potential.
            elif abs(EP[k][i]) > 1E-10 and EP[k][i]<0:      # If you're using different signal 
                logEP[k][i] = -log(abs(EP[k][i]))           # charges, it won't look cool anyway.
            else:                                           # In this case, you might want to use
                logEP[k][i] = None                          # the linear scale instead of logarithmic.
        rx += dx
    ry -= dy
  
imshow(logEP)
gray()
show()

from math import pi
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from numpy import arange, array

# Initial data
pos1 = 5.0*pi/180   # Angular initial position of the pendulum 1, in radians
vel1 = 0.0          # Angular initial velocity of the pendulum 1, in rad/s
pos2 = -5.0*pi/180  # Angular initial position of the pendulum 2, in radians
vel2 = 0.0          # Angular initial velocity of the pendulum 2, in rad/s
pos3 = 5.0*pi/180   # Angular initial position of the pendulum 3, in radians
vel3 = 0.0          # Angular initial velocity of the pendulum 3, in rad/s
m1 = 0.1353         # Mass of the pendulum 1
m2 = 0.1352         # Mass of the pendulum 2
m3 = 0.1349         # Mass of the pendulum 3
l = 41*0.01         # Lenght of the strings (in this simulation, they must be equal)
g = 9.785           # Gravity acceleration
k1 = 2.90           # Spring 1 constant
k2 = 2.39           # Spring 2 constant
deltat = 10.0       # Time interval for the simulation

# Solving the differential equations system
def deriv(y, t):
    return array([y[1], k1*(y[2]-y[0])/m1-(g/l)*y[0], y[3], k2*(y[4]-y[2])/m2-k1*(y[2]-y[0])/m2-(g/l)*y[2], y[5], -k2*(y[4]-y[2])/m3-(g/l)*y[4]])

time = arange(0.0, deltat, 0.01)
yinit = array([pos1, vel1, pos2, vel2, pos3, vel3])
y = odeint(deriv, yinit, time)

# Plotting the graph
ax = plt.axes(ylim = (-0.20, 0.20))
plt.plot(time, y[:,0], label='Mass 1')
plt.plot(time, y[:,2], label='Mass 2')
plt.plot(time, y[:,4], label='Mass 3')
plt.xlabel('t (s)')
plt.ylabel('Amplitude')
plt.title('Graph title')
plt.grid(True)
plt.legend()
plt.show()

import numpy as np
import scipy.integrate as si
import matplotlib.pyplot as plt
from timeit import default_timer as timer

def func(x,m,theta):
    return np.cos(m*theta-x*np.sin(theta))

x0 = 0.0        # start
x1 = 20.0       # end
N = 1000        # number of points
x = np.linspace(x0,x1,N)
theta = np.linspace(0.0,np.pi,N)

start = timer()   # Starting the timer to measure computation time
J0 = np.zeros(N,float)
for i in range(N):
    y0 = np.zeros(N,float)
    for k in range(N):
        y0[k] = func(x[i],0,theta[k])
    J0[i] = si.simps(y=y0,x=theta)    

J1 = np.zeros(N,float)
for i in range(N):
    y1 = np.zeros(N,float)
    for k in range(N):
        y1[k] = func(x[i],1,theta[k])
    J1[i] = si.simps(y=y1,x=theta)
    
J2 = np.zeros(N,float)
for i in range(N):
    y2 = np.zeros(N,float)
    for k in range(N):
        y2[k] = func(x[i],2,theta[k])
    J2[i] = si.simps(y=y2,x=theta)
    
time = timer()-start
print "Calculation time = %f seconds" % time
    
plt.plot(x,J0,label=r'$J_0$')
plt.plot(x,J1,label=r'$J_1$')
plt.plot(x,J2,label=r'$J_2$')
plt.xlabel('x')
plt.ylabel('J')
plt.legend()
plt.show()

import numpy as np
import scipy.integrate as si
import matplotlib.pyplot as plt
from timeit import default_timer as timer

def func(x,m,theta):
    return np.cos(m*theta-x*np.sin(theta))

x0 = 0.0        # start
x1 = 20.0       # end
N = 1000        # number of points
x = np.linspace(x0,x1,N)
theta = np.linspace(0.0,np.pi,N)

start = timer()
J0 = np.array([si.simps(y=np.array([func(xk,0,tk) for tk in theta]),x=theta) for xk in x])
J1 = np.array([si.simps(y=np.array([func(xk,1,tk) for tk in theta]),x=theta) for xk in x]) 
J2 = np.array([si.simps(y=np.array([func(xk,2,tk) for tk in theta]),x=theta) for xk in x]) 
time = timer()-start
print "Calculation time = %f seconds" % time
    
plt.plot(x,J0,label=r'$J_0$')
plt.plot(x,J1,label=r'$J_1$')
plt.plot(x,J2,label=r'$J_2$')
plt.xlabel('x')
plt.ylabel('J')
plt.legend()
plt.show()

import numpy as np
import scipy.integrate as si
import scipy.special as sp
import matplotlib.pyplot as plt
from timeit import default_timer as timer

def func(x,m,theta):
    return np.cos(m*theta-x*np.sin(theta))

x0 = 0.0        # start
x1 = 20.0       # end
N = 1000        # number of points
x = np.linspace(x0,x1,N)
theta = np.linspace(0.0,np.pi,N)

start = timer()
J0 = np.array([sp.jv(0,xk) for xk in x])
J1 = np.array([sp.jv(1,xk) for xk in x])
J2 = np.array([sp.jv(2,xk) for xk in x]) 
time = timer()-start
print "Calculation time = %f seconds" % time
    
plt.plot(x,J0,label=r'$J_0$')
plt.plot(x,J1,label=r'$J_1$')
plt.plot(x,J2,label=r'$J_2$')
plt.xlabel('x')
plt.ylabel('J')
plt.legend()
plt.show()

import numpy as np
import scipy.special as sp
from pylab import imshow,show,gray,hot
import matplotlib.pyplot as plt
from timeit import default_timer as timer

def func(x,m,theta):
    return np.cos(m*theta-x*np.sin(theta))

wl = 500E-9     # wavelength, in meters
size = 2.0      # image size, in micrometers
N = 500         # plot resolution
G = 5E12        # Gigantic constant

start = timer()
dists = np.array([[np.sqrt(float((N/2-i)**2+(N/2-j)**2)) for i in range(N)] for j in range(N)])
dists = 1E-6*np.sqrt(size)*dists/np.max(dists)     # distances matrix

image = np.zeros([N,N],float)
for i in range(N):
    for j in range(N):
        if abs(i-N/2) < 2 and abs(j-N/2)<2:     # This deals with the limit r -> 0
            image[i,j] = 0.5**2
        else:
            image[i,j] = (sp.jv(1,G*2.*np.pi*wl*dists[i,j])/(G*2.*np.pi*wl*dists[i,j]))**2
time = timer()-start
print "Calculation time = %f seconds" % time

imshow(image)
gray()
show()

imshow(image,vmax=0.01)
gray()
show()