test_cases2.py
237 lines
| 7.6 KiB
| text/x-python
|
PythonLexer
/ SRC / test_cases2.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 = 'plots2' | |||
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 plots | |||
f = 500 # is the frequency of our signal of test (which gives the w = 2 * pi * f) | |||
# 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 | |||
# 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 | |||
# 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) | |||
# 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) | |||
# 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]) | |||
#This is the number of points used in the FFT | |||
N = 256 | |||
#We may add some noise to the input signals | |||
# b1 = b1 + np.random.random(len(b1))/10 | |||
# b2 = b2 + np.random.random(len(b2))/10 | |||
# b3 = b3 + np.random.random(len(b3))/10 | |||
# e1 = e1 + np.random.random(len(e1))/10 | |||
# e2 = e2 + np.random.random(len(e2))/10 | |||
fm1 = 24576 | |||
fm2 = 4096 | |||
fm3 = 256 | |||
S1, w1 = Spectral_MatriceAvgTime(b1[::4],b2[::4],b3[::4],e1[::4],e2[::4],fm1,N) | |||
S2, w2 = Spectral_MatriceAvgTime(b1[::24],b2[::24],b3[::24],e1[::24],e2[::24],fm2,N) | |||
S3, w3 = Spectral_MatriceAvgTime(b1[::384],b2[::384],b3[::384],e1[::384],e2[::384],fm3,N) | |||
print "S1" | |||
PB1 = PSD_B(S1) | |||
print "Power spectrum density of the Magnetic Field" | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w1/(2*math.pi),PB1) | |||
plt.title("Power spectrum density of the Magnetic Field (fm = %s Hz)" % fm1) | |||
plt.xlabel("frequency [Hz]") | |||
plt.ylabel("PB(w)") | |||
plt.axvline(x = f, color = 'r', linestyle='--') | |||
plt.savefig('plots2/PB1.pdf') | |||
PE1 = PSD_E(S1) | |||
print "Power spectrum density of the Electric Field" | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w1/(2*math.pi),PE1) | |||
plt.title("Power spectrum density of the Electric Field (fm = %s Hz)" % fm1) | |||
plt.xlabel("frequency [Hz]") | |||
plt.ylabel("PE(w)") | |||
plt.axvline(x = f, color = 'r', linestyle='--') | |||
plt.savefig('plots2/PE1.pdf') | |||
e1 = ellipticity(S1) | |||
print "Ellipticity of the electromagnetic wave" | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w1/(2*math.pi),e1) | |||
plt.title("Ellipticity of the electromagnetic wave (fm = %s Hz)" % fm1) | |||
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('plots2/e1.pdf') | |||
d1 = degree_polarization(S1) | |||
print "Degree of polarization of the electromagnetic wave" | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w1/(2*math.pi),d1) | |||
plt.title("Degree of polarization of the electromagnetic wave (fm = %s Hz)" % fm1) | |||
plt.xlabel("frequency [Hz]") | |||
plt.ylabel("d(w)") | |||
plt.axvline(x = f, color = 'r', linestyle='--') | |||
plt.savefig('plots2/d1.pdf') | |||
#[n11,n21,n31] = normal_vector(S1) | |||
print "S2" | |||
PB2 = PSD_B(S2) | |||
print "Power spectrum density of the Magnetic Field" | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w2/(2*math.pi),PB2) | |||
plt.title("Power spectrum density of the Magnetic Field (fm = %s Hz)" % fm2) | |||
plt.xlabel("frequency [Hz]") | |||
plt.ylabel("PB(w)") | |||
plt.axvline(x = f, color = 'r', linestyle='--') | |||
plt.savefig('plots2/PB2.pdf') | |||
plt.close() | |||
PE2 = PSD_E(S2) | |||
print "Power spectrum density of the Electric Field" | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w2/(2*math.pi),PE2) | |||
plt.title("Power spectrum density of the Electric Field (fm = %s Hz)" % fm2) | |||
plt.xlabel("frequency [Hz]") | |||
plt.ylabel("PE(w)") | |||
plt.axvline(x = f, color = 'r', linestyle='--') | |||
plt.savefig('plots2/PE2.pdf') | |||
plt.close() | |||
e2 = ellipticity(S2) | |||
print "Ellipticity of the electromagnetic wave" | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w2/(2*math.pi),e2) | |||
plt.title("Ellipticity of the electromagnetic wave (fm = %s Hz)" % fm2) | |||
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('plots2/e2.pdf') | |||
plt.close() | |||
d2 = degree_polarization(S2) | |||
print "Degree of polarization of the electromagnetic wave" | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w2/(2*math.pi),d2) | |||
plt.title("Degree of polarization of the electromagnetic wave (fm = %s Hz)" % fm2) | |||
plt.xlabel("frequency [Hz]") | |||
plt.ylabel("d(w)") | |||
plt.axvline(x = f, color = 'r', linestyle='--') | |||
plt.savefig('plots2/d2.pdf') | |||
plt.close() | |||
#[n12,n22,n32] = normal_vector(S2) | |||
print "S3" | |||
PB3 = PSD_B(S3) | |||
print "Power spectrum density of the Magnetic Field" | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w3/(2*math.pi),PB3) | |||
plt.title("Power spectrum density of the Magnetic Field (fm = %s Hz)" % fm3) | |||
plt.xlabel("frequency [Hz]") | |||
plt.ylabel("PB(w)") | |||
plt.axvline(x = f, color = 'r', linestyle='--') | |||
plt.savefig('plots2/PB3.pdf') | |||
plt.close() | |||
PE3 = PSD_E(S3) | |||
print "Power spectrum density of the Electric Field " | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w3/(2*math.pi),PE3) | |||
plt.title("Power spectrum density of the Electric Field (fm = %s Hz)" % fm3) | |||
plt.xlabel("frequency [Hz]") | |||
plt.ylabel("PE(w)") | |||
plt.axvline(x = f, color = 'r', linestyle='--') | |||
plt.savefig('plots2/PE3.pdf') | |||
plt.close() | |||
e3 = ellipticity(S3) | |||
print "Ellipticity of the electromagnetic wave" | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w3/(2*math.pi),e3) | |||
plt.title("Ellipticity of the electromagnetic wave (fm = %s Hz)" % fm3) | |||
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('plots2/e3.pdf') | |||
plt.close() | |||
d3 = degree_polarization(S3) | |||
print "Degree of polarization of the electromagnetic wave" | |||
fig = fig + 1 | |||
plt.figure(fig) | |||
plt.plot(w3/(2*math.pi),d3) | |||
plt.title("Degree of polarization of the electromagnetic wave (fm = %s Hz)" % fm3) | |||
plt.xlabel("frequency [Hz]") | |||
plt.ylabel("d(w)") | |||
plt.axvline(x = f, color = 'r', linestyle='--') | |||
plt.savefig('plots2/d3.pdf') | |||
plt.close() | |||
#[n13,n23,n33] = normal_vector(S3) |