#!/usr/bin/python3
import numpy as np
import scipy.signal as sig
import matplotlib.pyplot as plt

spettro = np.genfromtxt('CO.csv', delimiter=',',skip_header=1,names=True)
spettro = spettro[:][::-1]  # invertiamo gli assi
sp_x=spettro['cm1']
sp_y=spettro['A']

banda_01=(sp_x < 2300) & (sp_x > 2000)

x_01=sp_x[banda_01]
y_01=sp_y[banda_01]

plt.plot(x_01,y_01)
plt.show()

peakind= sig.find_peaks_cwt(y_01, np.arange(5,20))

px=x_01[peakind]; py=y_01[peakind] #posizioni e altezze dei picchi

plt.plot(x_01,y_01)
plt.plot(px, py, 'ob', markersize=3)
plt.show()

px=px[py>0.1]; py=py[py>0.1] # prendo solo i picchi più alti di 0.1

plt.plot(x_01,y_01)
plt.plot(px, py, 'ob')
for i in range(len(px)):
  pxi=px[i]; pyi=py[i]
  plt.annotate('%.1f'%(pxi), xy=(pxi,pyi), xytext=(pxi-2, pyi+0.2), arrowprops=dict(
    facecolor='black', width=0.2, headwidth=3, headlength=5, shrink=0.1))
plt.xlim((px[0]-10,px[-1]+10))
plt.show()

m=np.arange(len(px))

int_medio=(px[-1]-px[0])/len(px)
intervalli=px[1:]-px[0:-1]
pos= np.argmax(intervalli>int_medio*1.6)

m[pos+1:]=m[pos+1:]-pos
m[:pos+1]=m[:pos+1]-pos-1
plt.plot(m,px,'o-r')
plt.show()

coeff = np.polyfit(m, px, 3)
polinomio=np.poly1d(coeff)
valori_fit=polinomio(m)

plt.plot(m,px,'o-r')
plt.plot(m,valori_fit,'+b')
plt.show()

print(polinomio)