|
|
#Operations with float numbers#
|
|
|
|
|
|
import numpy as np
|
|
|
from filters import *
|
|
|
from bin16 import *
|
|
|
from fft import *
|
|
|
import math
|
|
|
import cmath
|
|
|
import matplotlib.pyplot as plt
|
|
|
from basic_parameters import *
|
|
|
import os
|
|
|
|
|
|
################################################# Simulations ##########################################################
|
|
|
|
|
|
#This creates a new folder to save the plots
|
|
|
mypath = 'plots1'
|
|
|
if not os.path.isdir(mypath):
|
|
|
os.makedirs(mypath)
|
|
|
|
|
|
fs = 98304 # is the sampling frequency of our signals
|
|
|
time_step = 1.0/fs
|
|
|
t_ini = 0
|
|
|
t_fin = 10
|
|
|
time_vec = np.arange(t_ini,t_fin,time_step)
|
|
|
fig = 0 #index for the plots1
|
|
|
|
|
|
f = 2500 # is the frequency of our signal of test (which gives the w = 2 * pi * f)
|
|
|
|
|
|
print "These are the results obtained for the basic parameters on a monochromatic electromagnetic wave"
|
|
|
print "The rounding errors of the LFR are still not taken into account - we're operating with float numbers \n"
|
|
|
|
|
|
# we must give the vector n, which has the information about the direction of k
|
|
|
n = np.array([1,1,1]) # for now we take this example
|
|
|
n = n/np.linalg.norm(n)
|
|
|
l = 3000 # this is the wavelength of our wave (which is not necessarily determined by c = l * f)
|
|
|
k = 2.0 * math.pi / l
|
|
|
|
|
|
print "The n vector is:", n
|
|
|
print "The wavelength l is:", l, "[m]"
|
|
|
print "The wave number k is then:", k, "[m]^-1 \n"
|
|
|
|
|
|
# the values of 'a' and 'b' tell us the intensity of the magnetic field
|
|
|
a = np.random.random() # we take a value at random between 0 and 1
|
|
|
b = np.sqrt(1 - a*a) # we want to force a^2 + b^2 = 1
|
|
|
phi = 0
|
|
|
|
|
|
print "The parameters who describe the magnetic field B are 'a' and 'b', who shall form a complex phasor given by: [a -b*1j 0] * exp(1j * (w*t + phi))"
|
|
|
print "We're considering phi equal to ", phi
|
|
|
print "The vector for the present case is: [", a, ", -j", b, ", 0]"
|
|
|
print "('a' and 'b' were chosen randomly between 0 and 1) \n"
|
|
|
|
|
|
# the matrix R gives us an example of transformation to the SCM referencial
|
|
|
# it was determined by the use of the normalized vector n and then two random vectors who form an orthogonal basis with it
|
|
|
R = Matrix_SCM(n)
|
|
|
|
|
|
# now we transform the magnetic field to the SCM referencial
|
|
|
B = MagneticComponents(a,b,2.0 * math.pi * f,R)
|
|
|
|
|
|
print "From the n vector we can construct an Orthonormal Matrix R capable of taking the components of B and taking them into the SCM reference"
|
|
|
print "R = ", R
|
|
|
print "The values of the components of B in the SCM reference are then: "
|
|
|
print "B = ", B
|
|
|
print " "
|
|
|
|
|
|
# now we have the caracteristics for our functions b1(t), b2(t) and b3(t) who shall enter in the LFR
|
|
|
# b1(t) = cmath.polar(B[0])[0] * cos(w*t + cmath.polar(B[0])[1] + phi)
|
|
|
# b2(t) = cmath.polar(B[1])[0] * cos(w*t + cmath.polar(B[1])[1] + phi)
|
|
|
# b3(t) = cmath.polar(B[2])[0] * cos(w*t + cmath.polar(B[2])[1] + phi)
|
|
|
|
|
|
b1 = cmath.polar(B[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[0,0])[1] + phi)
|
|
|
b2 = cmath.polar(B[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[1,0])[1] + phi)
|
|
|
b3 = cmath.polar(B[2,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(B[2,0])[1] + phi)
|
|
|
|
|
|
print "Now we finally have the functions b1(t), b2(t) and b3(t) who are going to serve as examples for inputs in the LFR analyzer"
|
|
|
print "b1(t) = ", cmath.polar(B[0,0])[0], "cos(w*t + ", cmath.polar(B[0,0])[1], " + ", phi, ")"
|
|
|
print "b2(t) = ", cmath.polar(B[1,0])[0], "cos(w*t + ", cmath.polar(B[1,0])[1], " + ", phi, ")"
|
|
|
print "b3(t) = ", cmath.polar(B[2,0])[0], "cos(w*t + ", cmath.polar(B[2,0])[1], " + ", phi, ")"
|
|
|
print "(the value of w is,", 2.0 * math.pi * f, ")"
|
|
|
print " "
|
|
|
|
|
|
graph = raw_input("Do you wish to plot the graphs for the magnetic field signals? (Y/N): ")
|
|
|
if graph == 'Y':
|
|
|
print "Next we plot the graph for the three magnetic field signals: b1(t) [red], b2(t) [blue], b3(t) [green]"
|
|
|
fig = fig + 1
|
|
|
plt.figure(fig)
|
|
|
p1, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b1[0:len(time_vec)/2000],'r')
|
|
|
p2, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b2[0:len(time_vec)/2000],'b')
|
|
|
p3, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b3[0:len(time_vec)/2000],'g')
|
|
|
plt.title("The three input signals for the Magnetic Field")
|
|
|
plt.xlabel("time [ms]")
|
|
|
plt.ylabel("b(t)")
|
|
|
plt.legend([p1, p2, p3], ["b1(t)", "b2(t)", "b3(t)"])
|
|
|
plt.show()
|
|
|
print " "
|
|
|
|
|
|
# the component of E who's parallel to n is of our choice. It represents the influence of the presence of charges where the electromagnetic
|
|
|
# wave is travelling. However, the component of E who's perpendicular to n is determined by B and n as follows
|
|
|
|
|
|
E_orth = -1.0 * 2.0 * math.pi * f/k * np.cross(n,np.array(B.transpose()))
|
|
|
E_para = np.zeros(3) # for now, we'll consider that E_para is zero
|
|
|
E = E_orth + E_para
|
|
|
E = E.transpose()
|
|
|
|
|
|
e1 = cmath.polar(E[0,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[0,0])[1])
|
|
|
e2 = cmath.polar(E[1,0])[0] * np.cos(2.0 * math.pi * f * time_vec + cmath.polar(E[1,0])[1])
|
|
|
|
|
|
print "The electric field has two components: one parallel to the n vector (E_para) and one orthogonal (E_orth)"
|
|
|
print "E_para represents the degree of polarization of the medium where the wave is propagating"
|
|
|
print "E_para = ", E_para
|
|
|
print "E_orth is fully determined by the B vector representing the magnetic field"
|
|
|
print "E_orth = -w/k * (n x B)"
|
|
|
print "E_orth = ", E_orth
|
|
|
print " "
|
|
|
|
|
|
print "We have then the functions e1(t), e2(t)who are going to serve as examples for inputs in the LFR analyzer"
|
|
|
print "e1(t) = ", cmath.polar(E[0,0])[0], "cos(w*t + ", cmath.polar(E[0,0])[1], " + ", phi, ")"
|
|
|
print "e2(t) = ", cmath.polar(E[1,0])[0], "cos(w*t + ", cmath.polar(E[1,0])[1], " + ", phi, ")"
|
|
|
print "(the value of w is,", 2.0 * math.pi * f, ")"
|
|
|
print " "
|
|
|
|
|
|
graph = raw_input("Do you wish to plot the graphs for the electric field signals? (Y/N): ")
|
|
|
if graph == 'Y':
|
|
|
print "Next we plot the graph for two of the electric field signals: e1(t) [red], e2(t) [blue]"
|
|
|
print "(The electric field is being considered in the SCM referencial, which is not totally right, but should work for now)"
|
|
|
fig = fig + 1
|
|
|
plt.figure(fig)
|
|
|
p1, = plt.plot(1000*time_vec[0:len(time_vec)/2000],e1[0:len(time_vec)/2000],'r')
|
|
|
p2, = plt.plot(1000*time_vec[0:len(time_vec)/2000],e2[0:len(time_vec)/2000],'b')
|
|
|
plt.title("The two input signals for the Electric Field")
|
|
|
plt.xlabel("time [ms]")
|
|
|
plt.ylabel("e(t)")
|
|
|
plt.legend([p1, p2], ["e1(t)", "e2(t)"])
|
|
|
plt.show()
|
|
|
print " "
|
|
|
|
|
|
#This is the number of points used in the FFT
|
|
|
N = 256
|
|
|
|
|
|
#We may add some noise to the input signals
|
|
|
b1 = b1 + 10.0 * np.random.random(len(b1))/100
|
|
|
b2 = b2 + 10.0 * np.random.random(len(b2))/100
|
|
|
b3 = b3 + 10.0 * np.random.random(len(b3))/100
|
|
|
e1 = e1 + 10.0 * 10000000.0 * np.random.random(len(e1))/100
|
|
|
e2 = e2 + 10.0 * 10000000.0 * np.random.random(len(e2))/100
|
|
|
|
|
|
graph = raw_input("Do you wish to plot the graphs for the noised magnetic and electric field signals? (Y/N): ")
|
|
|
if graph == 'Y':
|
|
|
print "Next we plot the graph for the three magnetic field signals: b1(t) [red], b2(t) [blue], b3(t) [green]"
|
|
|
print "As well as for the two electric field signals: e1(t) [red], e2(t) [blue]"
|
|
|
fig = fig + 1
|
|
|
plt.figure(fig)
|
|
|
p1, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b1[0:len(time_vec)/2000],'r')
|
|
|
p2, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b2[0:len(time_vec)/2000],'b')
|
|
|
p3, = plt.plot(1000*time_vec[0:len(time_vec)/2000],b3[0:len(time_vec)/2000],'g')
|
|
|
plt.title("The three input signals for the Magnetic Field")
|
|
|
plt.xlabel("time [ms]")
|
|
|
plt.ylabel("b(t)")
|
|
|
plt.legend([p1, p2, p3], ["b1(t)", "b2(t)", "b3(t)"])
|
|
|
plt.show()
|
|
|
|
|
|
fig = fig + 1
|
|
|
plt.figure(fig)
|
|
|
p1, = plt.plot(1000*time_vec[0:len(time_vec)/2000],e1[0:len(time_vec)/2000],'r')
|
|
|
p2, = plt.plot(1000*time_vec[0:len(time_vec)/2000],e2[0:len(time_vec)/2000],'b')
|
|
|
plt.title("The two input signals for the Electric Field")
|
|
|
plt.xlabel("time [ms]")
|
|
|
plt.ylabel("e(t)")
|
|
|
plt.legend([p1, p2], ["e1(t)", "e2(t)"])
|
|
|
plt.show()
|
|
|
print " "
|
|
|
|
|
|
S, w = Spectral_Matrice(b1,b2,b3,e1,e2,fs,N,offset = 0)
|
|
|
|
|
|
print "Now we calculate the Spectral Matrice for the ensemble of signals entering the LFR."
|
|
|
print "(For now, we're just taking one spectral matrice, for all the values of frequencies. The next step will be to average them in time and frequence)"
|
|
|
print "Then, from this Spectral Matrice we calculate the basic parameters:"
|
|
|
|
|
|
PB = PSD_B(S)
|
|
|
print "Power spectrum density of the Magnetic Field"
|
|
|
fig = fig + 1
|
|
|
plt.figure(fig)
|
|
|
plt.plot(w/(2*math.pi),PB)
|
|
|
plt.title("Power spectrum density of the Magnetic Field")
|
|
|
plt.xlabel("frequency [Hz]")
|
|
|
plt.ylabel("PB(w)")
|
|
|
plt.axvline(x = f, color = 'r', linestyle='--')
|
|
|
plt.savefig('plots1/PB.png')
|
|
|
|
|
|
PE = PSD_E(S)
|
|
|
print "Power spectrum density of the Electric Field"
|
|
|
fig = fig + 1
|
|
|
plt.figure(fig)
|
|
|
plt.plot(w/(2*math.pi),PE)
|
|
|
plt.title("Power spectrum density of the Electric Field")
|
|
|
plt.xlabel("frequency [Hz]")
|
|
|
plt.ylabel("PE(w)")
|
|
|
plt.axvline(x = f, color = 'r', linestyle='--')
|
|
|
plt.savefig('plots1/PE.png')
|
|
|
|
|
|
e = ellipticity(S)
|
|
|
print "Ellipticity of the electromagnetic wave"
|
|
|
fig = fig + 1
|
|
|
plt.figure(fig)
|
|
|
plt.plot(w/(2*math.pi),e)
|
|
|
plt.title("Ellipticity of the electromagnetic wave")
|
|
|
plt.xlabel("frequency [Hz]")
|
|
|
plt.ylabel("e(w)")
|
|
|
plt.axvline(x = f, color = 'r', linestyle = '--')
|
|
|
plt.axhline(y = 2.0/(a/b + b/a), color = 'r', linestyle='--')
|
|
|
plt.savefig('plots1/e.png')
|
|
|
|
|
|
d = degree_polarization(S)
|
|
|
print "Degree of polarization of the electromagnetic wave"
|
|
|
fig = fig + 1
|
|
|
plt.figure(fig)
|
|
|
plt.plot(w/(2*math.pi),d)
|
|
|
plt.title("Degree of polarization of the electromagnetic wave")
|
|
|
plt.xlabel("frequency [Hz]")
|
|
|
plt.ylabel("d(w)")
|
|
|
plt.savefig('plots1/d.png')
|
|
|
|
|
|
#[n1,n2,n3] = normal_vector(S)
|
