test_cases1.py
221 lines
| 9.2 KiB
| text/x-python
|
PythonLexer
/ SRC / test_cases1.py
Alexis Jeandet
|
r0 | #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) |