#Operations with integers#

import numpy as np
from bin16 import *
from fft import *
import math
import cmath
import matplotlib.pyplot as plt
from basic_parameters_Int import *
from basic_parameters import *

#This is the number of points used for the FFT
N = 256

#These are the limits of our quantization
liminfB = -1
limsupB = +1
liminfE = -10000000
limsupE = +10000000

#############################

#Example of signals taken from a simulation
a = np.random.rand()
b = np.sqrt(1 - a*a)
n = [1,1,1]
E_para = [0,0,0]
f = 1500
fs = 98304
l = 30000

[b1,b2,b3,e1,e2] = SignalsWaveMonochromatic(a,b,n,E_para,f,fs,l)

#############################

#We add noise
b1 = b1 + np.random.uniform(1.0*liminfB/100,1.0*limsupB/100,len(b1))
b2 = b2 + np.random.uniform(1.0*liminfB/100,1.0*limsupB/100,len(b2))
b3 = b3 + np.random.uniform(1.0*liminfB/100,1.0*limsupB/100,len(b3))
e1 = e1 + np.random.uniform(1.0*liminfE/100,1.0*limsupE/100,len(e1))
e2 = e2 + np.random.uniform(1.0*liminfE/100,1.0*limsupE/100,len(e2))

#This will take our float-number signals and quantize them into integers going from -32768 to 32767
#The liminf and limsup is what decides the maximum and minimum values for our quantizations
#quant16 is a function from the bin16 library

b1Q = np.zeros(len(b1))
b2Q = np.zeros(len(b2)) 
b3Q = np.zeros(len(b3)) 
e1Q = np.zeros(len(e1)) 
e2Q = np.zeros(len(e2)) 

for i in range(len(b1)):
	b1Q[i] = quant16(b1[i],liminfB,limsupB)
	b2Q[i] = quant16(b2[i],liminfB,limsupB)
	b3Q[i] = quant16(b3[i],liminfB,limsupB)
	e1Q[i] = quant16(e1[i],liminfE,limsupE)	
	e2Q[i] = quant16(e2[i],liminfE,limsupE)

#We calculate the Spectral Matrices (averaged in time) for the case of float-numbers operations and also for the integer-numbers operations
Si,wi = Spectral_MatriceAvgTime_Int(b1Q[::24],b2Q[::24],b3Q[::24],e1Q[::24],e2Q[::24],4096,N)
Sf,wf = Spectral_MatriceAvgTime(b1[::24],b2[::24],b3[::24],e1[::24],e2[::24],4096,N)

[PBf,PEf,df,ef,nf] = BasicParameters(Sf,wf)
[PBi,PEi,di,ei,ni] = BasicParameters_Int(Si,wi)

mypath = "plotsFloatVsInt"
if not os.path.isdir(mypath):
	os.makedirs(mypath)
	
fm = 4096

alfa = 1.0 * math.pow(2,30)
plt.plot(wf/(2*math.pi),PBf,'g')
plt.scatter(wf/(2*math.pi),1.0 * PBi/alfa)
plt.title("Power spectrum density of the Magnetic Field (fm = %s Hz)" % fm)
plt.xlabel("frequency [Hz]")
plt.ylabel("PB(w)")
plt.savefig('%s/PB.png' % mypath)
plt.close()

alfa = 1.0 * math.pow(2,30)/math.pow(10,14)
plt.plot(wf/(2*math.pi),PEf,'g')
plt.scatter(wf/(2*math.pi),1.0 * PEi/alfa)
plt.title("Power spectrum density of the Electric Field (fm = %s Hz)" % fm)
plt.xlabel("frequency [Hz]")
plt.ylabel("PE(w)")
plt.savefig('%s/PE.png' % mypath)
plt.close()

alfa = 16
plt.plot(wf/(2*math.pi),ef,'g')
plt.scatter(wf/(2*math.pi),1.0 * ei/alfa)
plt.title("Ellipticity of the electromagnetic wave (fm = %s Hz)" % fm)
plt.xlabel("frequency [Hz]")
plt.ylabel("e(w)")
plt.savefig('%s/e.png' % mypath)
plt.close()

alfa = 8
plt.plot(wf/(2*math.pi),df, 'g')
plt.scatter(wf/(2*math.pi),1.0 * di/alfa)
plt.title("Degree of polarization of the electromagnetic wave (fm = %s Hz)" % fm)
plt.xlabel("frequency [Hz]")
plt.ylabel("d(w)")
plt.savefig('%s/d.png' % mypath)
plt.close()