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