|
|
#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)
|
