short_version_calibration_silver_bullet.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "79ac4c4c-c91a-45f6-a960-c248f143182d",
   "metadata": {},
   "source": [
    "# Calibration silver bullet spectrometer"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "c492cf4f-dcd5-489e-b883-bb849b1312d6",
   "metadata": {},
   "source": [
    "The original script has been written by Raphael Deirmendjian in 2022 (see Github: https://github.com/Rafale05/Silver-Bullet-calibration). The plotting routines have been updated regarding to the original file and it has been made into a Jupyter notebook.\n",
    "\n",
    "In the cell below, the calibration data (both wavelength, background and absolute calibration) from the silver bullet spectrometer is processed (for more information see this Youtube video: https://youtu.be/WP0HZCWyedU). In addition, 5 plotting functions are defined, which you can choose after running the cell below.\n",
    "\n",
    "The relevant functions are:\n",
    "* Plot_raw\n",
    "* Outlier_report_peirce\n",
    "* Plot_norm\n",
    "* Plot_std\n",
    "* Plot_std_variation\n",
    "\n",
    "If you want to plot the data for specific years (can be 1 or several), input the year as a list variable (e.g. Plot_raw([2021,2022]) ).\n",
    "\n",
    "Otherwise, if you want to plot all available years, you can use: Plot_raw(all_years)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "40ff2ebf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import statistics\n",
    "import scipy.special\n",
    "from sklearn.metrics import r2_score\n",
    "import time\n",
    "\n",
    "##FUNCTIONS\n",
    "## Read_Sorted returns a 381 wl containing X measurements of the same wavelength\n",
    "def Read_Sorted(path,year):        \n",
    "    path_final=path\n",
    "    f=open(path_final,\"r\")\n",
    "    data1=[]\n",
    "    data2=[]\n",
    "    for ligne in f:\n",
    "        data1.append(ligne.strip().split())\n",
    "    data1.pop(0)\n",
    "    data1.pop(0)\n",
    "    data1.pop(0)\n",
    "    data1.pop(0)\n",
    "    data1.pop(0)\n",
    "    data1.pop(0)\n",
    "    data1.pop(0)\n",
    "    data1.pop(0)\n",
    "    realdata=[]\n",
    "    data2=[]\n",
    "    for i in data1:\n",
    "        if i[0]!=str(year):\n",
    "            for k in i:\n",
    "                data2.append(k)\n",
    "        if i[0]==str(year):\n",
    "            realdata.append(data2)\n",
    "            data2=[]\n",
    "    realdata.append(data2)\n",
    "    realdata.pop(0)\n",
    "    wl=np.linspace(8285,8706,381,True)\n",
    "    values_sorted=[]\n",
    "    for i in range (0,len(realdata[0])):\n",
    "        values_sorted.append([int(realdata[0][i])])\n",
    "    for i in range(1,len(realdata)):\n",
    "        for k in range (0,len(realdata[i])):\n",
    "            values_sorted[k].append(int(realdata[i][k]))\n",
    "    return(values_sorted)\n",
    "\n",
    "#Returns the adjusted wavelength tab depending on the .Cal file\n",
    "def Adjustment(path):\n",
    "    file=open(path)\n",
    "    data_loc=[]\n",
    "    for k in file :\n",
    "        data_loc.append(k.strip())\n",
    "    x0=float(data_loc[1])\n",
    "    x1=float(data_loc[2])\n",
    "    x2=float(data_loc[3])\n",
    "    x3=float(data_loc[4])\n",
    "    wl=np.linspace(0,381,381,True)\n",
    "    adjusted_wl=[]\n",
    "    for (i) in wl:\n",
    "        adjusted_wl.append((x0+float(i)*x1+float(i)**2*x2+float(i)**3*x3))\n",
    "        \n",
    "    return(adjusted_wl)\n",
    "\n",
    "##Calculate Peirce's coef (related to Peirce's criterion outlier treatment)\n",
    "def peirce_dev(N: int, n: int, m: int) -> float:\n",
    "    \"\"\"Peirce's criterion\n",
    "    \n",
    "    Returns the squared threshold error deviation for outlier identification\n",
    "    using Peirce's criterion based on Gould's methodology.\n",
    "    \n",
    "    Arguments:\n",
    "        - int, total number of observations (N)\n",
    "        - int, number of outliers to be removed (n)\n",
    "        - int, number of model unknowns (m)\n",
    "    Returns:\n",
    "        float, squared error threshold (x2)\n",
    "    \"\"\"\n",
    "    # Assign floats to input variables:\n",
    "    N = float(N)\n",
    "    n = float(n)\n",
    "    m = float(m)\n",
    "\n",
    "    # Check number of observations:\n",
    "    if N > 1:\n",
    "        # Calculate Q (Nth root of Gould's equation B):\n",
    "        Q = (n ** (n / N) * (N - n) ** ((N - n) / N)) / N\n",
    "        #\n",
    "        # Initialize R values (as floats)\n",
    "        r_new = 1.0\n",
    "        r_old = 0.0  # <- Necessary to prompt while loop\n",
    "        #\n",
    "        # Start iteration to converge on R:\n",
    "        while abs(r_new - r_old) > (N * 2.0e-16):\n",
    "            # Calculate Lamda\n",
    "            # (1/(N-n)th root of Gould's equation A'):\n",
    "            ldiv = r_new ** n\n",
    "            if ldiv == 0:\n",
    "                ldiv = 1.0e-6\n",
    "            Lamda = ((Q ** N) / (ldiv)) ** (1.0 / (N - n))\n",
    "            # Calculate x-squared (Gould's equation C):\n",
    "            x2 = 1.0 + (N - m - n) / n * (1.0 - Lamda ** 2.0)\n",
    "            # If x2 goes negative, return 0:\n",
    "            if x2 < 0:\n",
    "                x2 = 0.0\n",
    "                r_old = r_new\n",
    "            else:\n",
    "                # Use x-squared to update R (Gould's equation D):\n",
    "                r_old = r_new\n",
    "                r_new = np.exp((x2 - 1) / 2.0) * scipy.special.erfc(\n",
    "                    np.sqrt(x2) / np.sqrt(2.0)\n",
    "                )\n",
    "    else:\n",
    "        x2 = 0.0\n",
    "    return np.sqrt(x2)\n",
    "\n",
    "#Eliminates outliers using Peirce's criterion and returns additionnaly the amount\n",
    "#of outliers and their standard deviation\n",
    "def Peirce_crit(dataset):\n",
    "    mean=statistics.mean(dataset)\n",
    "    std=statistics.stdev(dataset)\n",
    "    rang=1\n",
    "    count=0\n",
    "    std_count=[]\n",
    "    newdata=dataset[:]\n",
    "    while rang !=0:\n",
    "        to_delete=[]\n",
    "        R=peirce_dev(len(newdata),rang,1)\n",
    "        dev_max=R*std\n",
    "        for i in range (0,len(newdata)):\n",
    "            dev=abs(newdata[i]-mean)\n",
    "            if dev>dev_max:\n",
    "                to_delete.append(i)\n",
    "                std_count.append(abs((newdata[i]-mean)/std))\n",
    "        if len(to_delete)==0:\n",
    "            rang=0\n",
    "        else:\n",
    "            for k in reversed(to_delete):\n",
    "                del newdata[k]\n",
    "            rang=len(to_delete)+1\n",
    "            count+=len(to_delete)\n",
    "    return(newdata,count,std_count)\n",
    "\n",
    "#Return an average tab \n",
    "def Average(dataset):\n",
    "    average=[]\n",
    "    for i in range (0,len(dataset)):\n",
    "        average.append(statistics.mean(dataset[i]))\n",
    "    return(average)\n",
    "\n",
    "#Subtracts the bg from the abs cal\n",
    "def Bg_subtraction(bg,abscal):\n",
    "    result=[]\n",
    "    for i in range(0,len(bg)):\n",
    "        result.append(abscal[i]-bg[i])\n",
    "    return(result)\n",
    "#Calculates the known brightness \n",
    "def Cal_fact(r,R,rho,wl,alpha,cert):\n",
    "        alpha_rad=(alpha*np.pi)/180\n",
    "        return(rho*cert*((r/R)**2)*np.cos(alpha_rad))\n",
    "    \n",
    "#Plots the calibration factor\n",
    "def Plot_raw(year):\n",
    "    count=0\n",
    "    snr=[]\n",
    "    for i in year:\n",
    "        path=dico[\"a\"+str(i)]\n",
    "        wl=Adjustment(path[2])\n",
    "        bg=Read_Sorted(dico[\"a\"+str(i)][0],int(i)) #raw file\n",
    "        abscal=Read_Sorted(dico[\"a\"+str(i)][1],int(i)) #raw file\n",
    "        #outliers sorting with Peirce's crit\n",
    "        bg_out=[]\n",
    "        abscal_out=[]\n",
    "        for ki in range(0,len(wl)):\n",
    "            bg_out.append(Peirce_crit(bg[ki])[0])\n",
    "            abscal_out.append(Peirce_crit(abscal[ki])[0])\n",
    "        #averaging\n",
    "        bg_m=[]\n",
    "        abscal_m=[]\n",
    "        for k in range(0,len(wl)):\n",
    "            bg_m.append(statistics.mean(bg_out[k]))\n",
    "            abscal_m.append(statistics.mean(abscal_out[k]))\n",
    "        #truncate\n",
    "        bot=8374\n",
    "        top=8525\n",
    "        bg_c=[]\n",
    "        abscal_c=[]\n",
    "        wl_c=[]\n",
    "        for k in range(0,len(wl)):\n",
    "            if wl[k]>=bot and wl[k]<=top:\n",
    "                bg_c.append(bg_m[k])\n",
    "                abscal_c.append(abscal_m[k])\n",
    "                wl_c.append(wl[k])\n",
    "        #Calculating the known brightness\n",
    "        B=[]\n",
    "        for k in range (0,len(wl_c)):\n",
    "            B.append(Cal_fact(r,dicop[\"p\"+str(i)][0],rho,wl[k],\n",
    "                              dicop[\"p\"+str(i)][1],synthetic_fit(wl_c[k])))\n",
    "        k=[] #k=calibration factor\n",
    "        for j in range (0,len(B)):\n",
    "            k.append(B[j]/(abscal_c[j]-bg_c[j]))\n",
    "        \n",
    "        plt.scatter(wl_c,k,marker='+') #real plot\n",
    "        #Interpolates\n",
    "        fit=np.polyfit(wl_c,k,1,rcond=None,full=True)\n",
    "        k_int=[]\n",
    "        for l in wl_c:\n",
    "            k_int.append(interpol(l,fit))\n",
    "        r_2=r2_score(k,k_int)\n",
    "        plt.plot(wl_c,k_int,\n",
    "                 label=str(i)+' linear regression : '+str('%.3E' %fit[0][0])+\n",
    "                 'x+' +str('%.3E' %fit[0][1] )+\n",
    "                 ', r²='+str('%.2E' %r_2))\n",
    "        count+=1\n",
    "        snr.append(statistics.mean(abscal_c)/statistics.mean(bg_c))\n",
    "    plt.grid()\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')\n",
    "    #plt.title(\"Values filtered with Peirce's criterion\")\n",
    "    plt.xlabel('Wavelength [Å]')\n",
    "    plt.ylabel('C [R/Å/CTS]')\n",
    "    #plt.savefig('NAME.png')\n",
    "    #plt.figure(2)\n",
    "    #plt.scatter(year,snr,color='black')\n",
    "    #plt.xlabel(\"years\")\n",
    "    #plt.ylabel('SNR estimate')\n",
    "    #plt.show()\n",
    "    #return(\"SNR values:\", snr)\n",
    "\n",
    "#Gives a report on the amount of outliers deleted\n",
    "def Outlier_report_peirce(year):\n",
    "    for i in year:\n",
    "        path=dico[\"a\"+str(i)]\n",
    "        wl=Adjustment(path[2])\n",
    "        bg=Read_Sorted(dico[\"a\"+str(i)][0],int(i)) #raw file\n",
    "        abscal=Read_Sorted(dico[\"a\"+str(i)][1],int(i)) #raw file\n",
    "        #outliers sorting with Peirce's crit\n",
    "        bg_out=[]\n",
    "        abscal_out=[]\n",
    "        nb_bg=[]\n",
    "        nb_abs=[]\n",
    "        for ki in range(0,len(wl)):\n",
    "            bg_out.append(Peirce_crit(bg[ki])[0])\n",
    "            abscal_out.append(Peirce_crit(abscal[ki])[0])\n",
    "            nb_bg.append(Peirce_crit(bg[ki])[1])\n",
    "            nb_abs.append(Peirce_crit(abscal[ki])[1])\n",
    "        plt.scatter(wl,nb_bg,label=str(i)+', background. '+\n",
    "                    str(sum(nb_bg))+'/'+\n",
    "                    str(len(bg[0] )*381)+' datapoints removed ('+ str(round(sum(nb_bg)/(len(bg[0])*381)*100,1))+'%)')\n",
    "        plt.scatter(wl,nb_abs,label=str(i)+', absolute calibration. '+\n",
    "                    str(sum(nb_abs))+'/'+\n",
    "                    str(len(abscal[0] )*381)+' datapoints removed ('+ str(round(sum(nb_abs)/(len(abscal[0])*381)*100,1))+'%)')\n",
    "    plt.grid()\n",
    "    plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')\n",
    "    plt.title(\"Outlier deletion report\")\n",
    "    plt.xlabel('Wavelength [Å]')\n",
    "    plt.ylabel('Deleted outliers [CTS]')\n",
    "    plt.show()\n",
    "\n",
    "#Gives a report on the standard deviation of the deleted outliers\n",
    "def Outlier_report_peirce_std(year):\n",
    "    for i in year:\n",
    "        path=dico[\"a\"+str(i)]\n",
    "        wl=Adjustment(path[2])\n",
    "        bg=Read_Sorted(dico[\"a\"+str(i)][0],int(i)) #raw file\n",
    "        abscal=Read_Sorted(dico[\"a\"+str(i)][1],int(i)) #raw file\n",
    "        #outliers sorting with Peirce's crit\n",
    "        bg_out=[]\n",
    "        abscal_out=[]\n",
    "        dev_bg=[]\n",
    "        dev_abs=[]\n",
    "        for ki in range(0,len(wl)):\n",
    "            for w in (Peirce_crit(bg[ki])[2]):\n",
    "                dev_bg.append(w)\n",
    "            for w in (Peirce_crit(abscal[ki])[2]):\n",
    "                dev_abs.append(w)\n",
    "        plt.hist(dev_bg,label=str(i)+', background ')\n",
    "        plt.hist(dev_abs,label=str(i)+', absolute calibration ')\n",
    "    plt.grid()\n",
    "    plt.legend()\n",
    "    plt.title(\"Outlier deletion report\")\n",
    "    plt.xlabel('Number of STD from mean [CTS]')\n",
    "    plt.ylabel('Deleted outliers [CTS]')\n",
    "    plt.show()  \n",
    "      \n",
    "#Plots the calibration factor, normalized at 8430 \n",
    "def Plot_norm(year):\n",
    "    norm=8430\n",
    "    for i in year:\n",
    "        path=dico[\"a\"+str(i)]\n",
    "        wl=Adjustment(path[2])\n",
    "        bg=Read_Sorted(dico[\"a\"+str(i)][0],int(i)) #raw file\n",
    "        abscal=Read_Sorted(dico[\"a\"+str(i)][1],int(i)) #raw file\n",
    "        #outliers sorting with Peirce's crit\n",
    "        bg_out=[]\n",
    "        abscal_out=[]\n",
    "        for ki in range(0,len(wl)):\n",
    "            bg_out.append(Peirce_crit(bg[ki])[0]) \n",
    "            abscal_out.append(Peirce_crit(abscal[ki])[0])\n",
    "        #averaging\n",
    "        bg_m=[]\n",
    "        abscal_m=[]\n",
    "        for k in range(0,len(wl)):\n",
    "            bg_m.append(statistics.mean(bg_out[k]))\n",
    "            abscal_m.append(statistics.mean(abscal_out[k]))\n",
    "        #truncate\n",
    "        bot=8374\n",
    "        top=8525\n",
    "        bg_c=[]\n",
    "        abscal_c=[]\n",
    "        wl_c=[]\n",
    "        for k in range(0,len(wl)):\n",
    "            if wl[k]>=bot and wl[k]<=top:\n",
    "                bg_c.append(bg_m[k])\n",
    "                abscal_c.append(abscal_m[k])\n",
    "                wl_c.append(wl[k])\n",
    "        #Calculating the known brightness\n",
    "        B=[]\n",
    "        for k in range (0,len(wl_c)):\n",
    "            B.append(Cal_fact(r,dicop[\"p\"+str(i)][0],rho,wl[k],\n",
    "                              dicop[\"p\"+str(i)][1],synthetic_fit(wl_c[k])))\n",
    "        k=[] #k=calibration factor\n",
    "        for j in range (0,len(B)):\n",
    "            k.append(B[j]/(abscal_c[j]-bg_c[j]))\n",
    "        #Interpolates\n",
    "        fit=np.polyfit(wl_c,k,1,rcond=None,full=True)\n",
    "        # print(fit) #to print the fit values for the different years (the first value gives the slope)\n",
    "        k_int=[]\n",
    "        wl_b=np.linspace(bot,top,2)\n",
    "        for l in wl_b:\n",
    "            k_int.append(l*fit[0][0]+(1-(norm*fit[0][0])))\n",
    "        plt.plot(wl_b,k_int,label=i)\n",
    "        \n",
    "    plt.grid()\n",
    "    plt.legend()\n",
    "    plt.xlabel('Wavelength [Å]')\n",
    "    plt.ylabel('C [R/Å/CTS]')\n",
    "    #plt.savefig('norm_calibration_allyears.png')\n",
    "    plt.show()\n",
    "    \n",
    "#Plots the std \n",
    "def Plot_std(year):\n",
    "    bg_std=[]\n",
    "    abs_std=[]\n",
    "    count=0\n",
    "    fig,axs=plt.subplots(2,2)\n",
    "    avg_std_bg=[]\n",
    "    avg_std_abs=[]\n",
    "    for i in year:\n",
    "        bg_std=[]\n",
    "        abs_std=[]\n",
    "        path=dico[\"a\"+str(i)]\n",
    "        wl=Adjustment(path[2])\n",
    "        bg=Read_Sorted(dico[\"a\"+str(i)][0],int(i)) #raw file\n",
    "        abscal=Read_Sorted(dico[\"a\"+str(i)][1],int(i)) #raw file\n",
    "        #outliers sorting with Peirce's crit\n",
    "        bg_out=[]\n",
    "        abscal_out=[]\n",
    "        bg_std_out=[]\n",
    "        abs_std_out=[]\n",
    "        for ki in range(0,len(wl)):\n",
    "            bg_out.append(Peirce_crit(bg[ki])[0])\n",
    "            abscal_out.append(Peirce_crit(abscal[ki])[0])\n",
    "        #averaging\n",
    "        bg_m=[]\n",
    "        abscal_m=[]\n",
    "        bg_out_m=[]\n",
    "        abscal_out_m=[]\n",
    "        for k in range(0,len(wl)):\n",
    "            bg_m.append(statistics.mean(bg[k]))\n",
    "            abscal_m.append(statistics.mean(abscal[k]))\n",
    "            bg_std.append(statistics.stdev(bg[k]))\n",
    "            abs_std.append(statistics.stdev(abscal[k]))\n",
    "            bg_std_out.append(statistics.stdev(bg_out[k]))\n",
    "            abs_std_out.append(statistics.stdev(abscal_out[k]))\n",
    "        #truncate\n",
    "        bot=8374\n",
    "        top=8525\n",
    "        bg_std_c=[]\n",
    "        abs_std_c=[]\n",
    "        wl_c=[]\n",
    "        bg_std_out_c=[]\n",
    "        abs_std_out_c=[]\n",
    "        for k in range(0,len(wl)):\n",
    "            if wl[k]>=bot and wl[k]<=top:\n",
    "                bg_std_c.append(bg_std[k])\n",
    "                abs_std_c.append(abs_std[k])\n",
    "                abs_std_out_c.append(abs_std_out[k])\n",
    "                bg_std_out_c.append(bg_std_out[k])\n",
    "                wl_c.append(wl[k])\n",
    "\n",
    "        avg_std_bg.append(statistics.mean(bg_std_c))\n",
    "        avg_std_abs.append(statistics.mean(abs_std_c))\n",
    "        \n",
    "        axs[0,0].plot(wl_c,bg_std_c,label=str(i))\n",
    "        axs[0,0].set_title('Background, outliers not treated')\n",
    "\n",
    "        axs[1,0].plot(wl_c,bg_std_out_c,label=str(i))\n",
    "        axs[1,0].set_title('Background, outliers treated')\n",
    "\n",
    "        axs[0,1].plot(wl_c,abs_std_c,label=str(i))\n",
    "        axs[0,1].set_title('Absolute calibration, outliers not treated')\n",
    "\n",
    "        axs[1,1].plot(wl_c,abs_std_out_c,label=str(i))\n",
    "        axs[1,1].set_title('Absolute calibration, outliers treated')\n",
    "       \n",
    "        count+=1\n",
    "    for ax in axs.flat:\n",
    "        ax.set(xlabel='Wavelength [Å]', ylabel='Standard deviation')\n",
    "        ax.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')\n",
    "    \n",
    "    fig.tight_layout(pad=0.5)\n",
    "    fig.set_size_inches(18.5, 10.5)\n",
    "\n",
    "    plt.figure(2)\n",
    "    plt.plot(year,avg_std_bg,label='Background')\n",
    "    plt.plot(year,avg_std_abs,label='Absolute calibration')\n",
    "    plt.ylabel('Standard deviation')\n",
    "    plt.title('Outliers not treated')\n",
    "    plt.xlabel(\"Years\")\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "#Plots the relative standard deviation\n",
    "def Plot_std_variation(year):\n",
    "    measurements_bg=[35,20,32,34,42,38,34,26,21,29,33]\n",
    "    measurements_abs=[30,23,32,41,47,38,39,24,18,30,38]\n",
    "    bg_std=[]\n",
    "    abs_std=[]\n",
    "    count=0\n",
    "    fig,axs=plt.subplots(2,2)\n",
    "    avg_std_bg=[]\n",
    "    avg_std_abs=[]\n",
    "    avg_std_out_bg=[]\n",
    "    avg_std_out_abs=[]\n",
    "    for i in year:\n",
    "        bg_std=[]\n",
    "        abs_std=[]\n",
    "        path=dico[\"a\"+str(i)]\n",
    "        wl=Adjustment(path[2])\n",
    "        bg=Read_Sorted(dico[\"a\"+str(i)][0],int(i)) #raw file\n",
    "        abscal=Read_Sorted(dico[\"a\"+str(i)][1],int(i)) #raw file\n",
    "        #outliers sorting with Peirce's crit\n",
    "        bg_out=[]\n",
    "        abscal_out=[]\n",
    "        bg_std_out=[]\n",
    "        abs_std_out=[]\n",
    "        for ki in range(0,len(wl)):\n",
    "            bg_out.append(Peirce_crit(bg[ki])[0])\n",
    "            abscal_out.append(Peirce_crit(abscal[ki])[0])\n",
    "        #averaging\n",
    "        bg_m=[]\n",
    "        abscal_m=[]\n",
    "        bg_out_m=[]\n",
    "        abscal_out_m=[]\n",
    "        for k in range(0,len(wl)):\n",
    "            bg_m.append(statistics.mean(bg[k]))\n",
    "            abscal_m.append(statistics.mean(abscal[k]))\n",
    "            bg_out_m.append(statistics.mean(bg_out[k]))\n",
    "            abscal_out_m.append(statistics.mean(abscal_out[k]))\n",
    "            bg_std.append(statistics.stdev(bg[k]))\n",
    "            abs_std.append(statistics.stdev(abscal[k]))\n",
    "            bg_std_out.append(statistics.stdev(bg_out[k]))\n",
    "            abs_std_out.append(statistics.stdev(abscal_out[k]))\n",
    "            \n",
    "        #truncate\n",
    "        bot=8374\n",
    "        top=8525\n",
    "        bg_std_c=[]\n",
    "        abs_std_c=[]\n",
    "        wl_c=[]\n",
    "        bg_std_out_c=[]\n",
    "        abs_std_out_c=[]\n",
    "        for k in range(0,len(wl)):\n",
    "            if wl[k]>=bot and wl[k]<=top:\n",
    "                bg_std_c.append(bg_std[k]/bg_m[k])\n",
    "                abs_std_c.append(abs_std[k]/abscal_m[k])\n",
    "                abs_std_out_c.append(abs_std_out[k]/abscal_out_m[k])\n",
    "                bg_std_out_c.append(bg_std_out[k]/bg_out_m[k])\n",
    "                wl_c.append(wl[k])\n",
    "\n",
    "        avg_std_bg.append(statistics.mean(bg_std_c))\n",
    "        avg_std_abs.append(statistics.mean(abs_std_c))\n",
    "\n",
    "        avg_std_out_bg.append(statistics.mean(bg_std_out_c))\n",
    "        avg_std_out_abs.append(statistics.mean(abs_std_out_c))\n",
    "        \n",
    "        axs[0,0].plot(wl_c,bg_std_c,label=str(i))\n",
    "        axs[0,0].set_title('Background, outliers not treated')\n",
    "\n",
    "        axs[1,0].plot(wl_c,bg_std_out_c,label=str(i))\n",
    "        axs[1,0].set_title('Background, outliers treated')\n",
    "\n",
    "        axs[0,1].plot(wl_c,abs_std_c,label=str(i))\n",
    "        axs[0,1].set_title('Absolute calibration, outliers not treated')\n",
    "\n",
    "        axs[1,1].plot(wl_c,abs_std_out_c,label=str(i))\n",
    "        axs[1,1].set_title('Absolute calibration, outliers treated')\n",
    "\n",
    "        count+=1\n",
    "    for ax in axs.flat:\n",
    "        ax.set(xlabel='Wavelength [Å]', ylabel='Relative standard deviation')\n",
    "        ax.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')\n",
    "    \n",
    "    fig.tight_layout(pad=0.5)\n",
    "    fig.set_size_inches(18.5, 10.5)\n",
    "        \n",
    "    plt.figure(2)\n",
    "    plt.plot(year,avg_std_bg,label='Background')\n",
    "    plt.plot(year,avg_std_abs,label='Absolute calibration')\n",
    "    plt.ylabel('Relative standard deviation')\n",
    "    plt.title('Ouliers not removed')\n",
    "    plt.xlabel(\"Years\")\n",
    "    plt.legend()\n",
    "    plt.figure(3)\n",
    "    plt.plot(year,avg_std_out_bg,label='Background')\n",
    "    plt.plot(year,avg_std_out_abs,label='Absolute calibration')\n",
    "    plt.title(\"Outliers removed\")\n",
    "    plt.ylabel('Relative standard deviation')\n",
    "    plt.xlabel(\"Years\")\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    #return(bg_std_c,abs_std_c)\n",
    "\n",
    "#Interpolation fuction to plot the linear interpolation of the calibration factor\n",
    "def interpol(x,FIT):\n",
    "        return(FIT[0][0]*x+FIT[0][1])\n",
    "\n",
    "##Lamp\n",
    "Lamp_200W_B_range_of_interest=[456129,465687,473653,481425,492220,498682]\n",
    "Lamp_200W_wavelength_range_of_interest=np.arange(8200,8800,100)\n",
    "fiti=np.polyfit(Lamp_200W_wavelength_range_of_interest,Lamp_200W_B_range_of_interest,1,rcond=None, full=True)\n",
    "def synthetic_fit(x):\n",
    "\treturn(fiti[0][0]*x+fiti[0][1])\n",
    "score=r2_score(Lamp_200W_B_range_of_interest,\n",
    "               synthetic_fit(Lamp_200W_wavelength_range_of_interest))   \n",
    "    \n",
    "##INPUTS AND VARIABLES\n",
    "##In this section, input the files and parameters\n",
    "scan_start=8285\n",
    "scan_stop=8705\n",
    "nb_scan=381\n",
    "r=2\n",
    "rho=0.98\n",
    "scanning_wl=np.linspace(scan_start,scan_stop,nb_scan,True)\n",
    "\n",
    "path_data_2011_bg=(\"./calibration_files/CalData/s0202112_background_for_abs.lyr\")\n",
    "path_data_2011_abs=(\"./calibration_files/CalData/s0202112_abscal.lyr\")\n",
    "path_cal_2011=(\"./calibration_files/WavelengthCalibrations/Wavelength-Feb2011.Cal\")\n",
    "R_2011=51.85 #m\n",
    "a_2011=44 #deg\n",
    "\n",
    "path_data_2013_bg=(\"./calibration_files/CalData/s0602132Silver_BG.lyr\")\n",
    "path_data_2013_abs=(\"./calibration_files/CalData/s0602132Silver_ABS.lyr\")\n",
    "path_cal_2013=(\"./calibration_files/WavelengthCalibrations/Wavelength-Feb2013.Cal\")\n",
    "R_2013=52.05 #m\n",
    "a_2013=43.5 #deg\n",
    "\n",
    "path_data_2014_bg=(\"./calibration_files/CalData/silver_Background_0502142.lyr\")\n",
    "path_data_2014_abs=(\"./calibration_files/CalData/silver_ABS_CAL_0502142.lyr\")\n",
    "path_cal_2014=(\"./calibration_files/WavelengthCalibrations/Wavelength-Feb2014.Cal\")\n",
    "R_2014=51.87 #m\n",
    "a_2014=45 #deg\n",
    "\n",
    "path_data_2015_bg=(\"./calibration_files/CalData/s1202152BackgroundScans.lyr\")\n",
    "path_data_2015_abs=(\"./calibration_files/CalData/s1202152ABS_CAL.lyr\")\n",
    "path_cal_2015=(\"./calibration_files/WavelengthCalibrations/Wavelength-Feb2015.Cal\")\n",
    "R_2015=51.88 #m\n",
    "a_2015=46 #deg\n",
    "\n",
    "path_data_2016_bg=(\"./calibration_files/CalData/Back_10022016_silver_s1002161.lyr\")\n",
    "path_data_2016_abs=(\"./calibration_files/CalData/ABS_CAL_10022016_silver_s1002161.lyr\")\n",
    "path_cal_2016=(\"./calibration_files/WavelengthCalibrations/Wavelength-Feb2016.Cal\")\n",
    "R_2016=51.93 #m\n",
    "a_2016=45 #deg\n",
    "\n",
    "path_data_2017_bg=(\"./calibration_files/CalData/s0102172_BACKGROUND.lyr\")\n",
    "path_data_2017_abs=(\"./calibration_files/CalData/s0102172_ABS_CAL.lyr\")\n",
    "path_cal_2017=(\"./calibration_files/WavelengthCalibrations/Wavelength-Feb2017.Cal\")\n",
    "R_2017=51.974 #m\n",
    "a_2017=44 #deg\n",
    "\n",
    "path_data_2018_bg=(\"./calibration_files/CalData/s0802182_BACKGROUND.lyr\")\n",
    "path_data_2018_abs=(\"./calibration_files/CalData/s0802182_ABS_CAL.lyr\")\n",
    "path_cal_2018=(\"./calibration_files/WavelengthCalibrations/Wavelength-Feb2018.Cal\")\n",
    "R_2018=51.92 #m\n",
    "a_2018=44 #deg\n",
    "\n",
    "path_data_2019_bg=(\"./calibration_files/CalData/DARK_SCAN_ABS_CAL_s0702192.lyr\")\n",
    "path_data_2019_abs=(\"./calibration_files/CalData/ABS_CALs0702192.lyr\")\n",
    "path_cal_2019=(\"./calibration_files/WavelengthCalibrations/Wavelength-Feb2019.Cal\")\n",
    "R_2019=52.02 #m\n",
    "a_2019=45.2 #deg\n",
    "\n",
    "path_data_2020_bg=(\"./calibration_files/CalData/s2002202_CAL_BACKGROUND.lyr\")\n",
    "path_data_2020_abs=(\"./calibration_files/CalData/s2002202_ABS_CAL.lyr\")\n",
    "path_cal_2020=(\"./calibration_files/WavelengthCalibrations/Wavelength-Jan2020.Cal\")\n",
    "R_2020=51.98 #m\n",
    "a_2020=45 #deg\n",
    "\n",
    "path_data_2021_bg=(\"./calibration_files/CalData/s1802212_BG_AbsCal.lyr\")\n",
    "path_data_2021_abs=(\"./calibration_files/CalData/s1802212_ABS_cal.lyr\")\n",
    "path_cal_2021=(\"./calibration_files/WavelengthCalibrations/Wavelength-Feb2021.Cal\")\n",
    "R_2021=51.95 #m\n",
    "a_2021=45.5 #deg\n",
    "\n",
    "path_data_2022_bg=(\"./calibration_files/CalData/s1002222_BG_AbsCal.lyr\")\n",
    "path_data_2022_abs=(\"./calibration_files/CalData/s1002222_ABS_cal.lyr\")\n",
    "path_cal_2022=(\"./calibration_files/WavelengthCalibrations/Wavelength-Feb2022.Cal\")\n",
    "R_2022=51.96 #m\n",
    "a_2022=45  #deg\n",
    "\n",
    "path_data_2023_bg=(\"./calibration_files/CalData/s1502232_BG_ABScal.lyr\")\n",
    "path_data_2023_abs=(\"./calibration_files/CalData/s1502232_ABScal.lyr\")\n",
    "path_cal_2023=(\"./calibration_files/WavelengthCalibrations/Wavelength-Dec2022.Cal\")\n",
    "R_2023=52.08 #m\n",
    "a_2023=45  #deg\n",
    "\n",
    "path_data_2024_bg=(\"./calibration_files/CalData/s0802242_BG_ABScal.lyr\")\n",
    "path_data_2024_abs=(\"./calibration_files/CalData/s0802242_ABScal.lyr\")\n",
    "path_cal_2024=(\"./calibration_files/WavelengthCalibrations/Wavelength-Dec2023.Cal\")\n",
    "R_2024=52.09 #m\n",
    "a_2024=45  #deg\n",
    "\n",
    "dico=dict(a2011=[path_data_2011_bg,path_data_2011_abs,path_cal_2011],\n",
    "          a2013=[path_data_2013_bg,path_data_2013_abs,path_cal_2013],\n",
    "          a2014=[path_data_2014_bg,path_data_2014_abs,path_cal_2014],\n",
    "          a2015=[path_data_2015_bg,path_data_2015_abs,path_cal_2015],\n",
    "          a2016=[path_data_2016_bg,path_data_2016_abs,path_cal_2016],\n",
    "          a2017=[path_data_2017_bg,path_data_2017_abs,path_cal_2017],\n",
    "          a2018=[path_data_2018_bg,path_data_2018_abs,path_cal_2018],\n",
    "          a2019=[path_data_2019_bg,path_data_2019_abs,path_cal_2019],\n",
    "          a2020=[path_data_2020_bg,path_data_2020_abs,path_cal_2020],\n",
    "          a2021=[path_data_2021_bg,path_data_2021_abs,path_cal_2021],\n",
    "          a2022=[path_data_2022_bg,path_data_2022_abs,path_cal_2022],\n",
    "          a2023=[path_data_2023_bg,path_data_2023_abs,path_cal_2023],\n",
    "          a2024=[path_data_2024_bg,path_data_2024_abs,path_cal_2024],\n",
    "          )\n",
    "dicop=dict(p2011=[R_2011,a_2011],p2013=[R_2013,a_2013],p2014=[R_2014,a_2014],\n",
    "           p2015=[R_2015,a_2015],p2016=[R_2016,a_2016],p2017=[R_2017,a_2017],\n",
    "           p2018=[R_2018,a_2018],p2019=[R_2019,a_2019],p2020=[R_2020,a_2020],\n",
    "           p2021=[R_2021,a_2021],p2022=[R_2022,a_2022],p2023=[R_2023,a_2023],\n",
    "           p2024=[R_2024,a_2024])\n",
    "all_years=[2013,2015,2017,2018,2022,2023,2024]"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "eaf81054",
   "metadata": {},
   "source": [
    "A few examples and explanations about each plot: \n",
    "\n",
    "Plot_raw shows the non-normalized linear regression of the calibration factor for which outliers have been removed. The second plot gives an estimate of the Signal to Noise Ratio (SNR) (SNR = mean value of absolute calibration/mean value background measurement).\n",
    "\n",
    "Outlier_report_peirce shows how many outliers have been removed per wavelength for both the background and absolute calibration.\n",
    "\n",
    "Plot_norm shows the normalized linear regressions of the calibration factors, which makes it easier to compare the slopes for stability analysis of the instrument.\n",
    "\n",
    "Plot_std gives the standard deviation as a function of wavelength for the background and the absolute calibration both before and after the outliers have been treated. The last plot shows the average standard deviation per year for both the background and the absolute calibration before the outliers have been treated.\n",
    "\n",
    "Plot_std_variation shows the relative standard deviation (e.g. dividing the standard deviation by the measurement for each wavelength) as a function of wavelength for the background and the absolute calibration both before and after the outliers have been treated. The last two plots show the average relative standard deviation per year for both the background and the absolute calibration before and after the outliers have been treated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d08415e9-d658-4eb8-a04f-4b7928b8c449",
   "metadata": {},
   "outputs": [],
   "source": [
    "Plot_raw(all_years)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "faf4c83e",
   "metadata": {},
   "outputs": [],
   "source": [
    "Outlier_report_peirce(all_years)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b6f7cd89",
   "metadata": {},
   "outputs": [],
   "source": [
    "Plot_norm(all_years)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "fieldwork_silverbullet",
   "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.8.15"
  },
  "vscode": {
   "interpreter": {
    "hash": "95b7d4c420825ca60cfd7b9696c2affa304b32dffa3f4f643c507fa7d53b5f1d"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}