{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Importiamo le librerie necessarie."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.signal import find_peaks\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Impariamo a leggere un file in formato CSV. Ci sono molte varianti di questo formato ma il più comune consiste di campi separati da \",\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "spettro = np.genfromtxt('CO.csv', delimiter=',',skip_header=1,names=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternativamente possiamo usare read_csv dalla libreria pandas, che restituisce una tabella"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "sp=pd.read_csv('CO.csv',skiprows=1)\n",
    "sp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Tradizionalmente gli spettri IR hanno l'asse delle ascisse invertito, riportiamolo ad un ordine più adatto alla manipolazione. Inoltre per comodità estraiamo le variabili dalla tabella \"spettro\", la tabella ha le intestazioni 'cm1' e 'A'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "sp_x=spettro['cm1'][::-1] # invertiamo l'ordine del vettore con la notazione [::-1]\n",
    "sp_y=spettro['A'][::-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Soffermiamoci ad osservare lo spettro per definire gli intervalli di interesse"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(sp_x,sp_y)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Andiamo a selezionare una regione interessante dello spettro, in questo caso la banda per la transizione vibrazionale 0→1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "banda_01=(sp_x < 2300) & (sp_x > 2000)\n",
    "banda_01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Abbiamo definito la condizione banda_01 che è un vettore di valori booleani. Adesso possiamo\n",
    "usare la notazione comune a python, matlab, fortran ed altri linguaggi di alto livello che ci\n",
    "consente di selezionare gli elementi di un vettore in corrispondenza dei valori True."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "x_01=sp_x[banda_01]\n",
    "y_01=sp_y[banda_01]\n",
    "\n",
    "plt.plot(x_01,y_01)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Adesso, sfruttando lo spazio per il codice sopra, provate a restringere ulteriormente l'intervallo di definizione dello spettro.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Passiamo all'analisi dello spettro. Troviamo i picchi sfruttando le capacità in scipy.signal.\n",
    "La funzione find_peaks_cwt, non è stata disegnata per analizzare spettri IR, ma funziona comunque abbastanza bene per i nostri scopi."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "peakind, _=find_peaks(y_01, distance=20, height=0.2)\n",
    "plt.plot(x_01,y_01)\n",
    "plt.plot(x_01[peakind],y_01[peakind],'ob')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "find_peaks?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "La funzione find_peaks richiede in ingresso lo spettro ed altri valori opzionali, che servono per trattare il rumore, la presenza di picchi minori e altro. Vengono restituiti gli indici corrispondenti ai massimi trovati nel vettore dello spettro, in questo caso y_01.\n",
    "Abbiamo aggiunto al nostro grafico l'annotazione dei picchi trovati con cerchi blu 'ob'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Adesso andiamo a numerare i picchi, tramite l'indice m=J+1 per le transizioni J→J+1, e m=-J per le transizioni J→J-1. Cominciamo con il creare un indice m fittizio, per aggiustarlo successivamente. Ma prima eliminiamo alcuni picchi erroneamente trovati. Per semplicità andiamo solo a restringere lo spettro."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "m_01=np.arange(len(peakind))\n",
    "p_01=x_01[peakind]\n",
    "plt.plot(m_01,p_01,'o-r')\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Ora bisogna aggiustare l'indice m. Possiamo farlo sistematicamente approfittando della discontinuità intorno ad m=0, visibile anche nel grafico. Intorno ad m=0 la distanza tra i picchi è doppia.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "def trova_00(x):\n",
    "    d=x[1]-x[0]\n",
    "    for i in range(len(x)-1):\n",
    "        di=x[i+1]-x[i]\n",
    "        if di>1.5*d:\n",
    "            return 0.5*(x[i+1]+x[i])\n",
    "\n",
    "def indice_m(x):\n",
    "    \"restituisce l'indice m per una serie di picchi\"\n",
    "    v00=trova_00(x)\n",
    "    nsin=sum(x<v00)\n",
    "    ndex=sum(x>v00)\n",
    "    psin=np.arange(-nsin,0)\n",
    "    pdex=np.arange(ndex)+1\n",
    "    return np.append(psin, pdex)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "m_01=indice_m(p_01)\n",
    "plt.plot(m_01,p_01,'o-r')\n",
    "plt.axvline(0, ls='--')\n",
    "m_01"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Infine andiamo ad effettuare un fit polinomiale delle posizioni dei picchi in funzione di m. Il polinomio, che tiene conto della distorsione centrifuga e dell'accoppiamento rotovibrazionale, è di terzo grado. La libreria numpy consente questo calcolo tramite la funzione polyfit. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "coeff = np.polyfit(m_01, p_01, 3)\n",
    "polinomio=np.poly1d(coeff)\n",
    "valori_fit=polinomio(m_01)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "polyfit ritorna i coefficienti del polinomio. Usate il riquadro sopra per visionarli. Per calcolare i valori del polinomio è conveniente usare la funzione poly1d che permette di trasformare i coefficienti in un oggetto che può essere chiamato come una funzione. Provate per esempio a chiamare polinomio(x) con diversi valori di x.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "Verifichiamo il nostro fit. La grafica è la migliore prima opzione."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "autoscroll": false,
    "ein.tags": "worksheet-0",
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [],
   "source": [
    "plt.plot(m_01,p_01,'o-r')\n",
    "plt.plot(m_01,valori_fit,'+b')\n",
    "plt.show()\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.2"
  },
  "name": "peaks.ipynb"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}