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