{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "cell-0",
   "metadata": {
    "id": "08e72b31"
   },
   "source": [
    "### Notebook basics\n",
    "- Μπορείτε να εκτελέσετε τον κώδικα κελί-κελί\n",
    "- για να εκτελέσετε το κελί, πατάτε shift+enter\n",
    "- με escape+a εισάγετε νέο κελί από πάνω, escape+b νέο κελί από κάτω\n",
    "- για να τυπώσετε μια μεταβλητη x γράφετε print(x) ή απλά x (όπου x το όνομα της μεταβλητής)\n",
    "- στην python δεν δηλώνουμε τύπο μεταβλητής, αποδίδεται κατά την εξίωση\n",
    "- πακέτο numpy (short for Numerical Python) για πράξεις πινάκων"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tα δεδομένα που θα χρησιμοποιήσουμε είναι από το πακέτο PyHammer, που προσδιορίζει αυτόματα το φασματικό τύπο αστέρων.\n",
    "# Για το σκοπό αυτό χρησιμοποιεί πρότυπα φάσματα (templates) που θα χρησιμοποιήσουμε και σε αυτήν την άσκηση.\n",
    "# Ολόκληρο το πακέτο μπορεί να βρεθεί στο https://github.com/BU-hammerTeam/PyHammer\n",
    "# Κατεβάζουμε μόνο τα 13 αρχεία FITS που χρειαζόμαστε, απευθείας από το GitHub.\n",
    "#https://github.com/BU-hammerTeam/PyHammer/tree/master/resources\n",
    "# Προσοχή: το link \"github.com/.../blob/...\" δίνει ιστοσελίδα (HTML) και όχι το αρχείο. Χρειάζεται το \"raw.githubusercontent.com/...\".\n",
    "\n",
    "import os\n",
    "import urllib.request\n",
    "\n",
    "URL = \"https://raw.githubusercontent.com/BU-hammerTeam/PyHammer/master/resources/templates/\"\n",
    "DATA_DIR = \"pyHammer_templates\"          # φάκελος όπου αποθηκεύονται τα αρχεία\n",
    "os.makedirs(DATA_DIR, exist_ok=True)\n",
    "\n",
    "# Λεξικό: όνομα αστέρα -> όνομα αρχείου στο GitHub\n",
    "# (τα O, B και A0 δεν έχουν μεταλλικότητα στο όνομα· τα υπόλοιπα είναι νάνοι με [Fe/H] = +0.0)\n",
    "files = {'O5': 'O5.fits',\n",
    "         'B0': 'B0.fits',\n",
    "         'B5': 'B5.fits',\n",
    "         'A0': 'A0.fits',\n",
    "         'A6': 'A6_+0.0_Dwarf.fits',\n",
    "         'F0': 'F0_+0.0_Dwarf.fits',\n",
    "         'F5': 'F5_+0.0_Dwarf.fits',\n",
    "         'G0': 'G0_+0.0_Dwarf.fits',\n",
    "         'G5': 'G5_+0.0_Dwarf.fits',\n",
    "         'K0': 'K0_+0.0_Dwarf.fits',\n",
    "         'K5': 'K5_+0.0_Dwarf.fits',\n",
    "         'M0': 'M0_+0.0_Dwarf.fits',\n",
    "         'M5': 'M5_+0.0_Dwarf.fits'}\n",
    "\n",
    "for star in files:\n",
    "    path = os.path.join(DATA_DIR, files[star])\n",
    "    if not os.path.isfile(path):                      # κατεβάζουμε μόνο αν δεν υπάρχει ήδη\n",
    "        urllib.request.urlretrieve(URL + files[star], path)\n",
    "    print(star, \"->\", path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-2",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Τι αρχεία κατεβάσαμε στο φάκελο pyHammer_templates;\n",
    "print(sorted(os.listdir(DATA_DIR)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-3",
   "metadata": {
    "id": "e36f1c67"
   },
   "source": [
    "Σημείωση για μονάδες: μήκος κύματος [ Angstrom ] και ροή [ 10$^{-17}$ erg cm$^{-2}$ s$^{-1}$ Angstrom$^{-1}$ ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Τα αρχεία FITS διαβάζονται με τη βιβλιοθήκη astropy. Ένα αρχείο FITS έχει ένα ή περισσότερα μέρη (HDU).\n",
    "from astropy.io import fits\n",
    "import numpy as np   # ονοματίζουμε τη βιβλιοθήκη για ευκολία χρήσης\n",
    "\n",
    "# hdul = fits.open(\"pyHammer_templates/B0.fits\") # Header Data Unit list\n",
    "hdul = fits.open(\"pyHammer_templates/B5.fits\")\n",
    "hdul.info()        # 3 μέρη: 0 = σχόλια, 1 = φάσμα, 2 = ιδιότητες του προτύπου\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94721ac1-9b1c-49fb-ba1e-6ae20060fcc0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Το HDU 0 έχει μόνο σχόλια: από πόσα φάσματα SDSS/BOSS φτιάχτηκε το πρότυπο\n",
    "# print(repr(hdul[0].header))\n",
    "# print(repr(hdul[1].header))\n",
    "print(repr(hdul[2].header))\n",
    "# print(repr(hdul[2].data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# HDU 1: το φάσμα. Στήλες: LogLam = log10(λ / Å), Flux, PropErr και Std (σφάλματα της ροής)\n",
    "spec = hdul[1].data\n",
    "print(spec.columns.names)\n",
    "\n",
    "wavelength1 = 10**spec[\"LogLam\"]     # μήκος κύματος σε Å (το αρχείο δίνει log10 του λ)\n",
    "flux1 = spec[\"Flux\"]                 # ροή σε 1e-17 erg/s/cm^2/Å\n",
    "PropErr = spec[\"PropErr\"]\n",
    "Std = spec[\"Std\"]\n",
    "\n",
    "print(wavelength1[:5])\n",
    "print(flux1[:5])\n",
    "print(PropErr[:5])\n",
    "print(Std[:5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# HDU 2: ιδιότητες του προτύπου (χρώματα SDSS, και για κάποια πρότυπα Teff και log g)\n",
    "# Κάθε στήλη έχει 3 τιμές: [μέση τιμή, σφάλμα φωτομετρίας, τυπική απόκλιση των τιμών που προστέθηκαν]\n",
    "prop = hdul[2].data\n",
    "print(prop.columns.names)\n",
    "print(prop)\n",
    "\n",
    "# Το B0 δεν έχει Teff. Τα πρότυπα F, G, K έχουν:\n",
    "prop_F0 = fits.open(\"pyHammer_templates/F0_+0.0_Dwarf.fits\")[2].data\n",
    "print(prop_F0.columns.names)\n",
    "print(\"F0: Teff =\", prop_F0[\"Teff\"][0], \"K, log g =\", prop_F0[\"Log_g\"][0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-7",
   "metadata": {},
   "outputs": [],
   "source": [
    "stars_b = {'O5': 50000}\n",
    "stars_b['G5']=5654.9395 #8000\n",
    "stars_b['G5']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Θέλουμε να κάνουμε το ίδιο πράγμα αυτόματα - όχι n φορές για n αστέρια.\n",
    "# Tα ονόματα των αστέρων, όπως και τη θερμοκρασία τους, που θα μας χρειαστεί παρακάτω, τα οργανώνουμε σε λεξικό.\n",
    "# Λεξικό = μία λίστα που τα στοιχεία της τα καλούμε με το όνομα.\n",
    "\n",
    "# Λιστα:\n",
    "# stars = [\"O5\",\"B0\",\"B5\",\"A0\",\"A6\",\"F0\",\"F5\",\"G0\",\"G5\",\"K0\",\"K5\",\"M0\",\"M5\"]\n",
    "\n",
    "# Λεξικό: θερμοκρασίες (K) αστέρων της κύριας ακολουθίας (νάνοι, τάξη V)\n",
    "# από τον πίνακα των Pecaut & Mamajek 2013, ApJS 208, 9 (έκδοση 2022).\n",
    "#https://ui.adsabs.harvard.edu/abs/2013ApJS..208....9P/abstract\n",
    "# Για τα πρότυπα F, G, K θα τις αντικαταστήσουμε παρακάτω με την Teff που δίνει το ίδιο το αρχείο FITS.\n",
    "stars = {'O5': 41400, 'B0': 31400, 'B5': 15700, 'A0': 9700, 'A6': 7910, 'F0': 7220, 'F5': 6550,\n",
    "         'G0': 5930, 'G5': 5660, 'K0': 5270, 'K5': 4440, 'M0': 3850, 'M5': 3060}\n",
    "# Από τη θεωρία, γενικά εύρη θερμοκρασίων (Κ):\n",
    "# O > 30000, B 30000 – 10000, A 10000 – 7500, F 7500 – 6000, G 6000 – 5200, K 5200 – 3700, M 3700 – 2400\n",
    "\n",
    "# Παραδειγμα\n",
    "for i in stars: print(i, stars[i])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Τώρα μπορούμε να φτιαξουμε ένα λεξικό από διδιάστατους πίνακες,\n",
    "# ώστε να χρησιμοποιούμε το όνομα του κάθε αστεριού για να βρίσκουμε τα δεδομένα του.\n",
    "# Σε ένα λεξικό, μπορούμε να προσθέσουμε στοιχεία ως spectra['NewQuantity']=value\n",
    "\n",
    "spectra = {}  # κενό λεξικό\n",
    "\n",
    "for i in stars:\n",
    "    hdul = fits.open(os.path.join(DATA_DIR, files[i]))   # το πλήρες path: pyHammer_templates/....fits\n",
    "    spec = hdul[1].data\n",
    "    spectra[i] = (10**spec[\"LogLam\"], spec[\"Flux\"])     # [wavelength, flux] για κάθε αστέρι\n",
    "\n",
    "    # Αν το αρχείο έχει Teff (πρότυπα F, G, K), τη χρησιμοποιούμε αντί για την τιμή του πίνακα\n",
    "    prop = hdul[2].data\n",
    "    if \"Teff\" in prop.columns.names:\n",
    "        print(f\"{i}: Teff από το FITS = {prop['Teff'][0]:.0f} K  (πίνακας: {stars[i]} K)\")\n",
    "        stars[i] = round(float(prop[\"Teff\"][0]))\n",
    "    else:\n",
    "        print(f\"{i}: χωρίς Teff στο FITS, κρατάμε την τιμή του πίνακα {stars[i]} K\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "49713d04-18bb-4156-a1d8-7a58d44f3246",
   "metadata": {},
   "outputs": [],
   "source": [
    "stars"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-10",
   "metadata": {
    "id": "4c374896"
   },
   "outputs": [],
   "source": [
    "# Τωρα το εκάστοτε spectra[i] είναι ο διδιαστατος πίνακας των δεδομένων ενός αστεριού, με μία κολώνα για το λ, και μία κολώνα για τη ροή\n",
    "print(spectra[\"G0\"][0])\n",
    "print(spectra[\"G0\"][1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-11",
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 565
    },
    "id": "74693b79",
    "outputId": "35cd7b7c-80aa-4821-ab81-def69ec935df"
   },
   "outputs": [],
   "source": [
    "# Ένα πρώτο plot να δουμε τα φάσματα της βιβλιοθήκης συγκρίτικά το ένα με το άλλο\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "plt.figure(figsize=(10,6))\n",
    "\n",
    "for i in stars:\n",
    "    plt.plot(spectra[i][0], spectra[i][1]/np.mean(spectra[i][1]), label=i)  # Mπορούμε να τα πλοτάρουμε κανονικοποιημένα / μέση τιμη\n",
    "\n",
    "plt.xlabel(\"Wavelength [Å]\")\n",
    "plt.ylabel(\"Flux / mean flux\")\n",
    "plt.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-12",
   "metadata": {
    "id": "7e9cbd1f"
   },
   "outputs": [],
   "source": [
    "# Θέλουμε να εστιάσουμε σε κάποια γραμμή απορρόφησης ανά αστερι. Εδώ δίνονται τα κεντρικά μήκη κύματος βασικων φασματικών γραμμών.\n",
    "# Μπορείτε να προσθέσετε γραμμές κατά το δοκούν.\n",
    "\n",
    "# Γραμμή απορρόφησης <-> μήκος κύματος (Α) + ευρος συνεχούς πάλι εν είδει λεξικού\n",
    "lines = {\n",
    "    \"Ha\":   {\"lam\": 6562.8,  \"cont_left\": [6505,6530],  \"cont_right\": [6620,6645]},      # Balmer a (2→3)\n",
    "    \"Hb\":   {\"lam\": 4861.3,  \"cont_left\": [4780,4820],  \"cont_right\": [4920,4950]},      # Balmer a (2→4)\n",
    "    \"Hg\":   {\"lam\": 4340.5,  \"cont_left\": [4250,4290],  \"cont_right\": [4390,4440]},      # Balmer a (2→5)\n",
    "    \"Hd\":   {\"lam\": 4101.7,  \"cont_left\": [4020,4060],  \"cont_right\": [4150,4190]},      # Balmer a (2→6)\n",
    "    \"Pa7\":  {\"lam\": 10049,   \"cont_left\": [10025,10030],\"cont_right\": [10075,10120]},    # Paschen δ (3→7)\n",
    "   }\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-13",
   "metadata": {
    "id": "891655c7",
    "outputId": "b9453ebe-6247-42bc-a500-7de0bfd69146"
   },
   "outputs": [],
   "source": [
    "# Ας εστιάσουμε σε κάποια γραμμή απορρόφησης για κάποιο αστερι\n",
    "# θα χρειαστεί να τρέξετε τον κώδικα για πολλές γραμμές κι αστέρια.\n",
    "\n",
    "selected_star=\"G0\"\n",
    "selected_line = \"Ha\"\n",
    "\n",
    "wavelength_init=spectra[selected_star][0]\n",
    "flux_init=spectra[selected_star][1]\n",
    "\n",
    "central_l= lines[selected_line][\"lam\"] # επιλεχθέν μήκος κύματος (διαλέγετε από παραπάνω λεξικό ή απλά βάλετε την τιμή με το χέρι σε Α)\n",
    "dl = 100                               # βρείτε κατάλληλο εύρος του κομμενου φάσματος (σε A) γύρω από τη γραμμή που να ξεπερνάει (αρκετα)\n",
    "                                       # τη διακριτική ικανότητα του οργάνου και να δίνει καλή μετρηση του συνεχούς.\n",
    "\n",
    "# πινακας με True / False για τη συνθηκη\n",
    "selected_range = (wavelength_init >= central_l-dl) & (wavelength_init <= central_l+dl)\n",
    "\n",
    "# αντικατάσταση των πινάκων με νεους πίνακες, που εχουν μονο τις τιμές εκεί όπου η συνθήκη είναι True.\n",
    "wavelength = wavelength_init[selected_range]\n",
    "flux = flux_init[selected_range]\n",
    "\n",
    "plt.plot(wavelength, flux)\n",
    "plt.axvline(central_l,color='lightgrey')\n",
    "plt.xlabel(\"Wavelength [Å]\")\n",
    "plt.ylabel(\"Flux [ 10**(-17) erg / (s * cm^2 * Å) ]\")\n",
    "#plt.ylim(min(flux),max(flux))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-14",
   "metadata": {
    "id": "fc79b29e",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Τώρα ας ξεχωρίσουμε κι ας φιτάρουμε το συνεχές, χρησιμοποιώντας μια ευθεία, κι ας ορίσουμε το ευρος της γραμμής\n",
    "\n",
    "# Η μάσκα είναι σε Α. Θα πρέπει να αλλαζει για κάθε φασματική γραμμή.\n",
    "\n",
    "mask_left =lines[selected_line][\"cont_left\"]  # η έτοιμη τιμή ή απλά βάζετε ένα δικο σας πίνακα, π.χ., [6505, 6530]\n",
    "mask_right=lines[selected_line][\"cont_right\"]\n",
    "\n",
    "wavelength_left = np.mean(mask_left)    # μέση τιμή του αντίστοιχου εύρους\n",
    "wavelength_right = np.mean(mask_right)\n",
    "\n",
    "cont_left  = (wavelength > mask_left[0])  & (wavelength < mask_left[1] ) # πίνακας με True/False εαν τηρείται η συνθηκη.\n",
    "cont_right = (wavelength > mask_right[0]) & (wavelength < mask_right[1] )\n",
    "cont = cont_left | cont_right # αριστερό ή δεξί παραθυρο = συνολο συνεχους\n",
    "\n",
    "flux_left = np.mean( flux[ cont_left ] ) # μέση τιμή της ροής, για pixel στα οποία ισχύει η συνθήκη του μήκους κύματος\n",
    "flux_right = np.mean( flux[ cont_right] )\n",
    "\n",
    "slope= (flux_right-flux_left)/(wavelength_right-wavelength_left)\n",
    "continuum = slope*(wavelength-wavelength_left) + flux_left\n",
    "\n",
    "plt.plot(wavelength, flux, label='data')\n",
    "plt.plot(wavelength, continuum, color='orange',label='continuum fit')\n",
    "\n",
    "in_line = (wavelength > mask_left[1]) & (wavelength < mask_right[0])\n",
    "plt.plot( wavelength[cont_left], flux[cont_left],color='orange', label='continuum range')\n",
    "plt.plot( wavelength[cont_right], flux[cont_right],color='orange')\n",
    "plt.plot( wavelength[in_line], flux[in_line], color='blue',label='line range')\n",
    "plt.axvline(mask_left[1],color='lightgrey')\n",
    "plt.axvline(mask_right[0],color='lightgrey')\n",
    "\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-15",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Τώρα η ροή της γραμμής: μετά την αφαίρεση του συνεχούς, αθροίζουμε το εμβαδον y*dx σε όλα τα pixel εντός των ορίων.\n",
    "# Το βήμα dλ δεν είναι σταθερό (τα pixel είναι ισαπέχοντα σε log λ), οπότε κρατάμε ένα dλ για κάθε pixel.\n",
    "\n",
    "dlambda = np.gradient(wavelength)\n",
    "\n",
    "line_flux = np.abs( np.sum((flux[in_line] - continuum[in_line]) * dlambda[in_line]) )  # 1e-17 erg / s /cm^2\n",
    "\n",
    "print(\"Flux\",line_flux, \"*1e-17 erg/s/cm^2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-16",
   "metadata": {},
   "outputs": [],
   "source": [
    "# έλεγχος εαν η ανίχνευση ειναι σημαντική (πάνω από το θόρυβο)\n",
    "\n",
    "# θόρυβος ανά pixel <-> standard deviation\n",
    "noise_per_pixel = np.std((flux[cont] - continuum[cont]))*np.mean(dlambda[cont])\n",
    "n_pixels = len( wavelength[in_line] )\n",
    "noise = np.sqrt(n_pixels)*noise_per_pixel\n",
    "\n",
    "print(\"Signal-to-Noise ratio: \",line_flux/noise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-17",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Yπολογισμός equivalent width (σε Å).\n",
    "\n",
    "EW = np.sum((1.0 - flux[in_line] / continuum[in_line]) * dlambda[in_line])\n",
    "\n",
    "print(\"Equivalent width:\",EW, \"Å\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-18",
   "metadata": {},
   "source": [
    "## Από τα βήματα σε συναρτήσεις\n",
    "Τα παραπάνω κελιά δείχνουν βήμα-βήμα πώς μετράμε μια γραμμή: κόβουμε το φάσμα, βρίσκουμε το συνεχές, μετράμε ροή, S/N και equivalent width.\n",
    "Για να μην αντιγράφουμε όλα τα κελιά για κάθε αστέρα και κάθε γραμμή, βάζουμε **τον ίδιο κώδικα** μέσα σε μία συνάρτηση, `measure_line(star, line)`.\n",
    "Από εδώ και πέρα αλλάζετε μόνο τα ορίσματα. Για νέα γραμμή, προσθέστε μία καταχώρηση στο λεξικό `lines`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-19",
   "metadata": {},
   "outputs": [],
   "source": [
    "def measure_line(star, line, dl=100, cont_left=None, cont_right=None, plot=True):\n",
    "    # star: όνομα αστέρα, π.χ. \"G0\"      line: όνομα γραμμής από το λεξικό lines, π.χ. \"Ha\"\n",
    "    # dl: μισό εύρος (Å) του κομμένου φάσματος γύρω από τη γραμμή\n",
    "    # cont_left, cont_right: παράθυρα συνεχούς [Å, Å]. Αν δεν δοθούν, παίρνουμε αυτά του λεξικού lines.\n",
    "    if cont_left is None:\n",
    "        cont_left = lines[line][\"cont_left\"]\n",
    "    if cont_right is None:\n",
    "        cont_right = lines[line][\"cont_right\"]\n",
    "    central_l = lines[line][\"lam\"]\n",
    "\n",
    "    # 1. κόβουμε το φάσμα γύρω από τη γραμμή\n",
    "    wavelength_init, flux_init = spectra[star]\n",
    "    selected_range = (wavelength_init >= central_l-dl) & (wavelength_init <= central_l+dl)\n",
    "    wavelength = wavelength_init[selected_range]\n",
    "    flux = flux_init[selected_range]\n",
    "\n",
    "    # 2. συνεχές: ευθεία από τη μέση ροή στα δύο παράθυρα\n",
    "    in_left  = (wavelength > cont_left[0])  & (wavelength < cont_left[1])\n",
    "    in_right = (wavelength > cont_right[0]) & (wavelength < cont_right[1])\n",
    "    cont = in_left | in_right\n",
    "    slope = (np.mean(flux[in_right]) - np.mean(flux[in_left])) / (np.mean(cont_right) - np.mean(cont_left))\n",
    "    continuum = slope*(wavelength - np.mean(cont_left)) + np.mean(flux[in_left])\n",
    "\n",
    "    # 3. ροή γραμμής, θόρυβος, S/N και equivalent width\n",
    "    in_line = (wavelength > cont_left[1]) & (wavelength < cont_right[0])\n",
    "    dlambda = np.gradient(wavelength)\n",
    "    line_flux = np.abs(np.sum((flux[in_line] - continuum[in_line]) * dlambda[in_line]))\n",
    "    noise = np.sqrt(np.sum(in_line)) * np.std(flux[cont] - continuum[cont]) * np.mean(dlambda[cont])\n",
    "    EW = np.sum((1.0 - flux[in_line] / continuum[in_line]) * dlambda[in_line])\n",
    "\n",
    "    # 4. γράφημα, ίδιο με το βήμα-βήμα παράδειγμα\n",
    "    if plot:\n",
    "        plt.figure(figsize=(8,4))\n",
    "        plt.plot(wavelength, flux, label='data')\n",
    "        plt.plot(wavelength, continuum, color='orange', label='continuum fit')\n",
    "        plt.plot(wavelength[in_left],  flux[in_left],  color='orange', lw=3, label='continuum range')\n",
    "        plt.plot(wavelength[in_right], flux[in_right], color='orange', lw=3)\n",
    "        plt.axvline(cont_left[1], color='lightgrey')\n",
    "        plt.axvline(cont_right[0], color='lightgrey')\n",
    "        plt.xlabel(\"Wavelength [Å]\")\n",
    "        plt.ylabel(\"Flux [ 10**(-17) erg / (s * cm^2 * Å) ]\")\n",
    "        plt.title(f\"{star}, {line}:  EW = {EW:.2f} Å,  S/N = {line_flux/noise:.0f}\")\n",
    "        plt.legend()\n",
    "        plt.show()\n",
    "\n",
    "    print(f\"{star} {line}: flux = {line_flux:.3g} *1e-17 erg/s/cm^2,  S/N = {line_flux/noise:.1f},  EW = {EW:.2f} Å\")\n",
    "    return EW"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-20",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Έλεγχος: πρέπει να δώσει το ίδιο EW με τα βήμα-βήμα κελιά παραπάνω\n",
    "measure_line(\"G0\", \"Ha\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-21",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Αλλάζετε μόνο τα ορίσματα. Παραδείγματα:\n",
    "measure_line(\"A0\", \"Hb\")\n",
    "\n",
    "# measure_line(\"A0\", \"Hb\", dl=150)                                          # μεγαλύτερο κομμάτι φάσματος\n",
    "# measure_line(\"A0\", \"Hb\", cont_left=[4740,4780], cont_right=[4940,4980])   # δικά σας παράθυρα συνεχούς"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-22",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Η ίδια μέτρηση για όλα τα αστέρια, με ένα βρόχο. plot=False για να μη βγουν 13 γραφήματα\n",
    "# (βάλτε plot=True για να ελέγξετε τα παράθυρα του συνεχούς σε κάθε αστέρα)\n",
    "EW_Ha = {}\n",
    "for i in stars:\n",
    "    EW_Ha[i] = measure_line(i, \"Ha\", plot=False)\n",
    "\n",
    "plt.figure(figsize=(10,4))\n",
    "plt.bar(list(EW_Ha.keys()), list(EW_Ha.values()))\n",
    "plt.xlabel(\"Spectral type  (θερμότερα αριστερά)\")\n",
    "plt.ylabel(\"EW(Hα) [Å]\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-23",
   "metadata": {
    "id": "43a33256-0bdd-4cec-b2c7-0a6374e8aa48"
   },
   "outputs": [],
   "source": [
    "# Ερωτηματα άσκησης\n",
    "\n",
    "# 1) κάντε την μέτρηση του EW της Ηα για όλα τα αστέρια και βρείτε σε ποιό είναι μέγιστη.\n",
    "# 2) για τον αστέρα F0, χρησιμοποιείστε τα EW της Hδ (n=2->6) και της Paschen 3->7 για να υπολογίσετε τη θερμοκρασία του αστέρα.\n",
    "#     Συγκρίνετε τη θερμοκρασία που βγάλατε με αυτή του μελανου σώματος που περιγράφει τη συνεχή εκπομπή στο παρακάτω σχημα."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-24",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ερώτημα 2: μία κλήση της συνάρτησης για κάθε γραμμή\n",
    "EW_Hd = measure_line(\"F0\", \"Hd\")\n",
    "# EW_Pa7 = measure_line(\"F0\", \"Pa7\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-25",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Σύγκριση με θεωρητική τιμή BB\n",
    "\n",
    "def black_body_luminosity(wavelength,T=1000): # wavelength in Angstrom\n",
    "    h, c, kb = 6.626e-27, 2.9979e18, 1.381e-16    # erg s, Å/s, erg/K\n",
    "    B_lambda = 2*h*c**2/ wavelength**5/(np.exp(h*c/(wavelength*kb*T))-1)\n",
    "    return B_lambda/np.mean(B_lambda) # κανονικοποιημένο σχήμα (οι μονάδες δεν μας ενδιαφέρουν εδώ)\n",
    "\n",
    "def plot_bb(star, T=None):\n",
    "    # star: όνομα αστέρα      T: θερμοκρασία BB σε K. Αν δεν δοθεί, παίρνουμε τη θερμοκρασία αναφοράς stars[star]\n",
    "    if T is None:\n",
    "        T = stars[star]\n",
    "    x = spectra[star][0]\n",
    "    y = spectra[star][1]/np.mean(spectra[star][1])\n",
    "    plt.figure(figsize=(8,4))\n",
    "    plt.plot(x, y, label=star+\", data\")\n",
    "    plt.plot(x, black_body_luminosity(x, T), label=f\"BB, T = {T:.0f} K\")\n",
    "    plt.xlabel(\"Wavelength [Å]\")\n",
    "    plt.ylabel(\"Flux / mean flux\")\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "plot_bb(\"F0\")             # με τη θερμοκρασία αναφοράς\n",
    "# plot_bb(\"F0\", T=7000)   # καντε overplot ένα μελαν σωμα με την Τ που βρήκατε από τα EW των γραμμών του Η"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cell-26",
   "metadata": {},
   "source": [
    "### EXTRA (προαιρετικό): θερμοκρασία από προσαρμογή μέλανος σώματος\n",
    "Βρίσκουμε τη θερμοκρασία *T* του μέλανος σώματος που ταιριάζει καλύτερα στο συνεχές (4000–9000 Å) κάθε αστέρα και τη συγκρίνουμε με τη θερμοκρασία αναφοράς.\n",
    "Ερώτηση: γιατί για τους αστέρες O και B η προσαρμογή δίνει πολύ μεγάλες (και αβέβαιες) θερμοκρασίες;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cell-27",
   "metadata": {},
   "outputs": [],
   "source": [
    "T_grid = np.logspace(np.log10(2000), np.log10(50000), 500)   # δοκιμαστικές θερμοκρασίες (K)\n",
    "\n",
    "for i in stars:\n",
    "    x, y = spectra[i]\n",
    "    m = (x > 4000) & (x < 9000)            # αφήνουμε έξω το άλμα Balmer (< 3650-4000 Å)\n",
    "    chi2 = []\n",
    "    for T in T_grid:\n",
    "        B = black_body_luminosity(x[m], T)\n",
    "        a = np.sum(y[m]*B)/np.sum(B*B)     # καλύτερη κανονικοποίηση για αυτή την T\n",
    "        chi2.append(np.sum((y[m]-a*B)**2))\n",
    "    T_fit = T_grid[np.argmin(chi2)]\n",
    "    print(f\"{i}:  T_BB = {T_fit:7.0f} K   T_ref = {stars[i]:6d} K\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39c298cc-1fc8-4e20-b3f5-7da2e0743ddf",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "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.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}