notebook_GMOS_a04_make_multiextension_files.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "db057b12",
   "metadata": {},
   "source": [
    "# Create Multi-Extension FITS files from individual extension files\n",
    "\n",
    "Creates multi-extension FITS (MEF) files by combining individual extension files for bias, flatfield, and science images for a particular night of science observations.\n",
    "\n",
    "## Before running notebook:\n",
    "- Install <tt>fpack</tt> and <tt>funpack</tt> if not already installed (https://heasarc.gsfc.nasa.gov/fitsio/fpack/); change paths to local installation location in functions <tt>compress_file_fpack</tt> and <tt>decompress_file_fpack</tt> in this notebook<br><br>\n",
    "\n",
    "- Move or copy all raw science data (MEF files from the Gemini archive) and processed science data (as individual extensions) produced by <tt>notebook_GMOS_a03_12amp_process_science_data</tt> for a particular night of science observations to a single directory.  Also move or copy the processed calibration files (median-combined bias and flatfield files) to the same directory.<br><br>\n",
    "\n",
    "- Run <tt>create_mef_files_all(base_path,data_set)</tt> where <tt>base_path</tt> is the directory containing the data and <tt>data_set</tt> is the designation to be used in the filename of the target data table.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61d36870",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os, sys\n",
    "from astropy.io import fits\n",
    "from astropy.io.fits import getheader\n",
    "from astropy.modeling import models\n",
    "from astropy import nddata\n",
    "from astropy import units as u\n",
    "import astropy.units as u\n",
    "import astropy.coordinates as coord\n",
    "from astropy.time import Time, TimeDelta\n",
    "from astropy.coordinates import SkyCoord,Angle\n",
    "from astroquery.jplhorizons import Horizons\n",
    "import ccdproc as cp\n",
    "from ccdproc import CCDData\n",
    "import numpy as np\n",
    "import math, glob, subprocess\n",
    "import datetime\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2641032",
   "metadata": {},
   "outputs": [],
   "source": [
    "def compress_file_fpack(filename):\n",
    "    fpack_command  = '/Users/hhsieh/Astro/tools/cfitsio/fpack'\n",
    "    fpack_cmd = [fpack_command,filename]\n",
    "    process = subprocess.call(fpack_cmd)\n",
    "    if process == 0:\n",
    "        os.remove(filename)\n",
    "    else:\n",
    "        print('fpack command failed.')\n",
    "    return None\n",
    "\n",
    "def decompress_file_fpack(filename):\n",
    "    funpack_command = '/Users/hhsieh/Astro/tools/cfitsio/funpack'\n",
    "    funpack_cmd = [funpack_command,filename]\n",
    "    process = subprocess.call(funpack_cmd)\n",
    "    if process == 0:\n",
    "        os.remove(filename)\n",
    "    else:\n",
    "        print('funpack command failed.')\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e79cffcb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def output_write_log_entry(log_filepath,log_entry):\n",
    "    print('{:s} - {:s}'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),log_entry))\n",
    "    with open(log_filepath,'a') as of:\n",
    "        of.write('{:s} - {:s}\\n'.format(datetime.datetime.today().strftime('%Y-%m-%d %H:%M:%S'),log_entry))\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7bc444c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def convert_to_horizons_search_term(target_name_from_fits_file):\n",
    "    horizons_search_term = target_name_from_fits_file\n",
    "    if target_name_from_fits_file == '16P':               horizons_search_term = '90000273'\n",
    "    if target_name_from_fits_file == '43P':               horizons_search_term = '90000508'\n",
    "    if target_name_from_fits_file == '46P':               horizons_search_term = '90000544'\n",
    "    if target_name_from_fits_file == '104P':              horizons_search_term = '90000962'\n",
    "    return horizons_search_term\n",
    "\n",
    "def convert_to_long_target_desig(short_desig):\n",
    "    # see https://pdssbn.astro.umd.edu/data_sb/resources/periodic_date.txt\n",
    "    long_desig = short_desig\n",
    "    if short_desig == '16P':               long_desig = '16P/Brooks 2'\n",
    "    if short_desig == '43P':               long_desig = '43P/Wolf-Harrington'\n",
    "    if short_desig == '46P':               long_desig = '46P/Wirtanen'\n",
    "    if short_desig == '104P':              long_desig = '104P/Kowal 2'\n",
    "    if short_desig == '124P':              long_desig = '124P/Mrkos'\n",
    "    if short_desig == '133P':              long_desig = '133P/1996 N2 (Elst-Pizarro 1)'\n",
    "    if short_desig == '162P':              long_desig = '162P/2004 TU12 (Siding Spring 2)'\n",
    "    if short_desig == '169P':              long_desig = '169P/2002 EX12 (NEAT 22)'\n",
    "    if short_desig == '176P':              long_desig = '176P/1999 RE70 (LINEAR 52)'\n",
    "    if short_desig == '182P':              long_desig = '182P/2001 WF2 (LONEOS 6)'\n",
    "    if short_desig == '197P':              long_desig = '197P/2003 KV2 (LINEAR 30)'\n",
    "    if short_desig == '209P':              long_desig = '209P/2004 CB (LINEAR 41)'\n",
    "    if short_desig == '210P':              long_desig = '210P/2003 K2 (Christensen 1)'\n",
    "    if short_desig == '233P':              long_desig = '233P/2005 JR71 (La Sagra 3)'\n",
    "    if short_desig == '238P':              long_desig = '238P/2005 U1 (Read 3)'\n",
    "    if short_desig == '249P':              long_desig = '249P/2006 U1 (LINEAR 53)'\n",
    "    if short_desig == '259P':              long_desig = '259P/2008 R1 (Garradd 4)'\n",
    "    if short_desig == '288P':              long_desig = '288P/2006 VW139 (300163)'\n",
    "    if short_desig == '300163':            long_desig = '288P/2006 VW139 (300163)'\n",
    "    if short_desig == '2006 VW139':        long_desig = '288P/2006 VW139 (300163)'\n",
    "    if short_desig == '304P':              long_desig = '304P/2008 Q2 (Ory 1)'\n",
    "    if short_desig == '311P':              long_desig = '311P/2013 P5 (PanSTARRS 23)'\n",
    "    if short_desig == 'P/2013 P5':         long_desig = '311P/2013 P5 (PanSTARRS 23)'\n",
    "    if short_desig == '313P':              long_desig = '313P/2014 S4 (Gibbs 16)'\n",
    "    if short_desig == 'P/2014 S4':         long_desig = '313P/2014 S4 (Gibbs 16)'\n",
    "    if short_desig == '324P':              long_desig = '324P/2010 R2 (La Sagra 4)'\n",
    "    if short_desig == 'P/2010 R2':         long_desig = '324P/2010 R2 (La Sagra 4)'\n",
    "    if short_desig == '331P':              long_desig = '331P/2012 F5 (Gibbs 14)'\n",
    "    if short_desig == 'P/2012 F5':         long_desig = '331P/2012 F5 (Gibbs 14)'\n",
    "    if short_desig == '348P':              long_desig = '348P/2017 A2 (PanSTARRS 62)'\n",
    "    if short_desig == '358P':              long_desig = '358P/2012 T1 (PanSTARRS 11)'\n",
    "    if short_desig == 'P/2012 T1':         long_desig = '358P/2012 T1 (PanSTARRS 11)'\n",
    "    if short_desig == '364P':              long_desig = '364P/2013 CU129 (PanSTARRS 18)'\n",
    "    if short_desig == '417P':              long_desig = '417P/2015 J3 (NEOWISE 2)'\n",
    "    if short_desig == 'P/2015 J3':         long_desig = '417P/2015 J3 (NEOWISE 2)'\n",
    "    if short_desig == '426P':              long_desig = '426P/2019 A7 (PanSTARRS 90)'\n",
    "    if short_desig == 'P/2019 A7':         long_desig = '426P/2019 A7 (PanSTARRS 90)'\n",
    "    if short_desig == '427P':              long_desig = '427P/2017 S5 (ATLAS 2)'\n",
    "    if short_desig == 'P/2017 S5':         long_desig = '427P/2017 S5 (ATLAS 2)'\n",
    "    if short_desig == '432P':              long_desig = '432P/2021 N4 (PanSTARRS 123)'\n",
    "    if short_desig == '433P':              long_desig = '433P/2005 QN173 (248370)'\n",
    "    if short_desig == '248370':            long_desig = '433P/2005 QN173 (248370)'\n",
    "    if short_desig == '2005 QN173':        long_desig = '433P/2005 QN173 (248370)'\n",
    "    if short_desig == '435P':              long_desig = '435P/2021 T3 (PanSTARRS 125)'\n",
    "    if short_desig == 'P/2021 T3':         long_desig = '435P/2021 T3 (PanSTARRS 125)'\n",
    "    if short_desig == '439P':              long_desig = '439P/2008 WZ96 (LINEAR 56)'\n",
    "    if short_desig == 'P/2008 WZ96':       long_desig = '439P/2008 WZ96 (LINEAR 56)'\n",
    "    if short_desig == '455P':              long_desig = '455P/2017 S9 (PanSTARRS 71)'\n",
    "    if short_desig == 'P/2017 S9':         long_desig = '455P/2017 S9 (PanSTARRS 71)'\n",
    "    if short_desig == '456P':              long_desig = '456P/2021 L4 (PanSTARRS 140)'\n",
    "    if short_desig == 'P/2021 L4':         long_desig = '456P/2021 L4 (PanSTARRS 140)'\n",
    "    if short_desig == '457P':              long_desig = '457P/2020 O1 (Lemmon-PanSTARRS 9)'\n",
    "    if short_desig == 'P/2020 O1':         long_desig = '457P/2020 O1 (Lemmon-PanSTARRS 9)'\n",
    "    if short_desig == '477P':              long_desig = '477P/2018 P3 (PanSTARRS 83)'\n",
    "    if short_desig == 'P/2018 P3':         long_desig = '477P/2018 P3 (PanSTARRS 83)'\n",
    "    if short_desig == 'P/2013 YG46':       long_desig = 'P/2013 YG46 (Spacewatch 7)'\n",
    "    if short_desig == 'P/2015 X6':         long_desig = 'P/2015 X6 (PANSTARRS 52)'\n",
    "    if short_desig == 'P/2016 J1-A':       long_desig = 'P/2016 J1-A (PANSTARRS 59)'\n",
    "    if short_desig == 'P/2016 J1-B':       long_desig = 'P/2016 J1-B (PANSTARRS 59)'\n",
    "    if short_desig == 'P/2016 P1':         long_desig = 'P/2016 P1 (PANSTARRS 57)'\n",
    "    if short_desig == 'P/2019 A3':         long_desig = 'P/2019 A3 (PanSTARRS 88)'\n",
    "    if short_desig == 'P/2019 A4':         long_desig = 'P/2019 A4 (PanSTARRS 89)'\n",
    "    if short_desig == 'P/2021 A5':         long_desig = 'P/2021 A5 (PanSTARRS 117)'\n",
    "    if short_desig == 'P/2022 R5':         long_desig = 'P/2022 R5 (PanSTARRS 134)'\n",
    "    if short_desig == '2000 SY178':        long_desig = '(62412) 2000 SY178'\n",
    "    return long_desig\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "497b59b8",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_geometry_params(object_name,obs_date,obs_time,exp_time):\n",
    "    obs_code  = '500'   # observatory code\n",
    "    time_step = '2d'\n",
    "    date_start = Time('{:s} {:s}'.format(obs_date,obs_time),scale='utc',format='iso')\n",
    "    mid_dt = TimeDelta(exp_time,format='sec') / 2\n",
    "    date_mid = date_start + mid_dt\n",
    "    date_next = date_mid + 1\n",
    "    obj = Horizons(id=convert_to_horizons_search_term(object_name),location=obs_code,epochs=date_mid.jd)\n",
    "    ephems = obj.ephemerides()\n",
    "    #print(ephems.columns)\n",
    "    ephem_params = {}\n",
    "    ephem_params['ra']          = ephems[0]['RA']\n",
    "    ephem_params['dec']         = ephems[0]['DEC']\n",
    "    ra_angle                    = Angle(ephem_params['ra'] * u.deg)\n",
    "    dec_angle                   = Angle(ephem_params['dec'] * u.deg)\n",
    "    ephem_params['ra_hms']      = ra_angle.to_string(unit=u.hour,sep=':',precision=0,pad=True)\n",
    "    ephem_params['dec_dms']     = dec_angle.to_string(unit=u.degree,sep=':',precision=0,pad=True,alwayssign=True)\n",
    "    ephem_params['ra_rate']     = ephems[0]['RA_rate']\n",
    "    ephem_params['dec_rate']    = ephems[0]['DEC_rate']\n",
    "    ephem_params['heliodist']   = ephems[0]['r']\n",
    "    ephem_params['geodist']     = ephems[0]['delta']\n",
    "    ephem_params['phsang']      = ephems[0]['alpha']\n",
    "    ephem_params['lunar_elong'] = ephems[0]['lunar_elong']\n",
    "    ephem_params['lunar_illum'] = ephems[0]['lunar_illum']\n",
    "    ephem_params['antisolarPA'] = ephems[0]['sunTargetPA']\n",
    "    ephem_params['neghelioPA']  = ephems[0]['velocityPA']\n",
    "    ephem_params['orbplang']    = ephems[0]['OrbPlaneAng']\n",
    "    ephem_params['glxlon']      = ephems[0]['GlxLon']\n",
    "    ephem_params['glxlat']      = ephems[0]['GlxLat']\n",
    "    ephem_params['trueanom']    = ephems[0]['true_anom']\n",
    "    \n",
    "    return ephem_params\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c42c9e38",
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_metadata_geminiN(fits_file):\n",
    "    object_name,object_desig,obs_date,obs_time,exp_time,filter_name = '','','','','',''\n",
    "    with fits.open(fits_file) as hdulist:\n",
    "        data = hdulist[0].data\n",
    "        hdr = hdulist[0].header\n",
    "    object_name = hdr['OBJECT']\n",
    "    obs_date    = hdr['DATE-OBS']\n",
    "    obs_time    = hdr['UTSTART'][:10]\n",
    "    exp_time    = hdr['EXPTIME']\n",
    "    filter_name = hdr['FILTER2'][:1]\n",
    "    object_desig = convert_to_long_target_desig(object_name)\n",
    "    return object_name,object_desig,obs_date,obs_time,exp_time,filter_name\n",
    "\n",
    "\n",
    "def get_metadata_geminiS(fits_file):\n",
    "    object_name,object_desig,obs_date,obs_time,exp_time,filter_name = '','','','','',''\n",
    "    with fits.open(fits_file) as hdulist:\n",
    "        data = hdulist[0].data\n",
    "        hdr = hdulist[0].header\n",
    "    object_name = hdr['OBJECT']\n",
    "    obs_date    = hdr['DATE-OBS']\n",
    "    obs_time    = hdr['UTSTART'][:10]\n",
    "    exp_time    = hdr['EXPTIME']\n",
    "    filter_name = hdr['FILTER2'][:1]\n",
    "    object_desig = convert_to_long_target_desig(object_name)\n",
    "    return object_name,object_desig,obs_date,obs_time,exp_time,filter_name\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c5919a71",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_mef_bias(data_path,chip03_filename_fz):\n",
    "    \n",
    "    # Set filenames\n",
    "    chip04_filename_fz = chip03_filename_fz[:-11]+'.04.fits.fz'\n",
    "    chip05_filename_fz = chip03_filename_fz[:-11]+'.05.fits.fz'\n",
    "    chip06_filename_fz = chip03_filename_fz[:-11]+'.06.fits.fz'\n",
    "    chip07_filename_fz = chip03_filename_fz[:-11]+'.07.fits.fz'\n",
    "    chip08_filename_fz = chip03_filename_fz[:-11]+'.08.fits.fz'\n",
    "    chip09_filename_fz = chip03_filename_fz[:-11]+'.09.fits.fz'\n",
    "    chip10_filename_fz = chip03_filename_fz[:-11]+'.10.fits.fz'\n",
    "    \n",
    "    # Unpack fz files\n",
    "    chip03_filename = chip03_filename_fz[:-3]\n",
    "    chip04_filename = chip04_filename_fz[:-3]\n",
    "    chip05_filename = chip05_filename_fz[:-3]\n",
    "    chip06_filename = chip06_filename_fz[:-3]\n",
    "    chip07_filename = chip07_filename_fz[:-3]\n",
    "    chip08_filename = chip08_filename_fz[:-3]\n",
    "    chip09_filename = chip09_filename_fz[:-3]\n",
    "    chip10_filename = chip10_filename_fz[:-3]\n",
    "    decompress_file_fpack(data_path+chip03_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip04_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip05_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip06_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip07_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip08_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip09_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip10_filename_fz)\n",
    "\n",
    "    # Get raw fits header data and Gemini program ID\n",
    "    with fits.open(data_path+chip03_filename) as hdulist:\n",
    "        hdr0 = hdulist[0].header\n",
    "        gemprgid = hdr0['GEMPRGID']\n",
    "        hdu0 = fits.PrimaryHDU(header=hdr0)\n",
    "        output_filename = gemprgid + '_' + chip03_filename[1:-13] + '.gmosn.bias.fits'\n",
    "\n",
    "    # Extract data from processed files\n",
    "    with fits.open(data_path+chip03_filename) as hdu: hdr03,data03 = hdu[0].header,hdu[0].data\n",
    "    with fits.open(data_path+chip04_filename) as hdu: hdr04,data04 = hdu[0].header,hdu[0].data\n",
    "    with fits.open(data_path+chip05_filename) as hdu: hdr05,data05 = hdu[0].header,hdu[0].data\n",
    "    with fits.open(data_path+chip06_filename) as hdu: hdr06,data06 = hdu[0].header,hdu[0].data\n",
    "    with fits.open(data_path+chip07_filename) as hdu: hdr07,data07 = hdu[0].header,hdu[0].data\n",
    "    with fits.open(data_path+chip08_filename) as hdu: hdr08,data08 = hdu[0].header,hdu[0].data\n",
    "    with fits.open(data_path+chip09_filename) as hdu: hdr09,data09 = hdu[0].header,hdu[0].data\n",
    "    with fits.open(data_path+chip10_filename) as hdu: hdr10,data10 = hdu[0].header,hdu[0].data\n",
    "        \n",
    "    # Combine data from processed files\n",
    "    hdu03 = fits.ImageHDU(data03)\n",
    "    hdu04 = fits.ImageHDU(data04)\n",
    "    hdu05 = fits.ImageHDU(data05)\n",
    "    hdu06 = fits.ImageHDU(data06)\n",
    "    hdu07 = fits.ImageHDU(data07)\n",
    "    hdu08 = fits.ImageHDU(data08)\n",
    "    hdu09 = fits.ImageHDU(data09)\n",
    "    hdu10 = fits.ImageHDU(data10)\n",
    "    new_hdul = fits.HDUList([hdu0,hdu03,hdu04,hdu05,hdu06,hdu07,hdu08,hdu09,hdu10])\n",
    "    \n",
    "    # Create headers\n",
    "    new_hdul[1].header             = new_hdul[1].header+hdr03\n",
    "    new_hdul[2].header             = new_hdul[2].header+hdr04\n",
    "    new_hdul[3].header             = new_hdul[3].header+hdr05\n",
    "    new_hdul[4].header             = new_hdul[4].header+hdr06\n",
    "    new_hdul[5].header             = new_hdul[5].header+hdr07\n",
    "    new_hdul[6].header             = new_hdul[6].header+hdr08\n",
    "    new_hdul[7].header             = new_hdul[7].header+hdr09\n",
    "    new_hdul[8].header             = new_hdul[8].header+hdr10\n",
    "    new_hdul[1].header['CHIPFILE'] = (chip03_filename,'Source data filename')\n",
    "    new_hdul[2].header['CHIPFILE'] = (chip04_filename,'Source data filename')\n",
    "    new_hdul[3].header['CHIPFILE'] = (chip05_filename,'Source data filename')\n",
    "    new_hdul[4].header['CHIPFILE'] = (chip06_filename,'Source data filename')\n",
    "    new_hdul[5].header['CHIPFILE'] = (chip07_filename,'Source data filename')\n",
    "    new_hdul[6].header['CHIPFILE'] = (chip08_filename,'Source data filename')\n",
    "    new_hdul[7].header['CHIPFILE'] = (chip09_filename,'Source data filename')\n",
    "    new_hdul[8].header['CHIPFILE'] = (chip10_filename,'Source data filename')\n",
    "\n",
    "    compress_file_fpack(data_path+chip03_filename)\n",
    "    compress_file_fpack(data_path+chip04_filename)\n",
    "    compress_file_fpack(data_path+chip05_filename)\n",
    "    compress_file_fpack(data_path+chip06_filename)\n",
    "    compress_file_fpack(data_path+chip07_filename)\n",
    "    compress_file_fpack(data_path+chip08_filename)\n",
    "    compress_file_fpack(data_path+chip09_filename)\n",
    "    compress_file_fpack(data_path+chip10_filename)\n",
    "\n",
    "    new_hdul.writeto(data_path+output_filename,overwrite=True)\n",
    "    with fits.open(data_path+output_filename) as hdu:\n",
    "        hdu.info()\n",
    "    \n",
    "    if os.path.exists(data_path+output_filename+'.fz'):\n",
    "        os.remove(data_path+output_filename+'.fz')\n",
    "    compress_file_fpack(data_path+output_filename)\n",
    "        \n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "66bb0b41",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_mef_flatfield(data_path,chip1_filename_fz):\n",
    "    \n",
    "    # Set filenames\n",
    "    chip2_filename_fz = chip1_filename_fz[:-14]+'.chip2.fits.fz'\n",
    "    chip3_filename_fz = chip1_filename_fz[:-14]+'.chip3.fits.fz'\n",
    "    \n",
    "    # Unpack fz files\n",
    "    chip1_filename = chip1_filename_fz[:-3]\n",
    "    chip2_filename = chip2_filename_fz[:-3]\n",
    "    chip3_filename = chip3_filename_fz[:-3]\n",
    "    decompress_file_fpack(data_path+chip1_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip2_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip3_filename_fz)\n",
    "\n",
    "    # Get raw fits header data and Gemini program ID\n",
    "    with fits.open(data_path+chip1_filename) as hdulist:\n",
    "        hdr0 = hdulist[0].header\n",
    "        gemprgid = hdr0['GEMPRGID']\n",
    "        hdu0 = fits.PrimaryHDU(header=hdr0)\n",
    "        output_filename = gemprgid + '_' + chip1_filename[1:-11] + '.fits'\n",
    "\n",
    "    # Extract data from processed files\n",
    "    with fits.open(data_path+chip1_filename) as hdu_chip1, fits.open(data_path+chip2_filename) as hdu_chip2, fits.open(data_path+chip3_filename) as hdu_chip3:\n",
    "        hdr1,data1 = hdu_chip1[0].header,hdu_chip1[0].data\n",
    "        hdr2,data2 = hdu_chip2[0].header,hdu_chip2[0].data\n",
    "        hdr3,data3 = hdu_chip3[0].header,hdu_chip3[0].data\n",
    "    \n",
    "    # Combine data from processed files\n",
    "    hdu1 = fits.ImageHDU(data1)\n",
    "    hdu2 = fits.ImageHDU(data2)\n",
    "    hdu3 = fits.ImageHDU(data3)\n",
    "    new_hdul = fits.HDUList([hdu0,hdu1,hdu2,hdu3])\n",
    "    \n",
    "    # Create headers\n",
    "    new_hdul[1].header             = new_hdul[1].header+hdr1\n",
    "    new_hdul[2].header             = new_hdul[2].header+hdr2\n",
    "    new_hdul[3].header             = new_hdul[3].header+hdr3\n",
    "    new_hdul[1].header['CHIPFILE'] = (chip1_filename,'Source data filename')\n",
    "    new_hdul[2].header['CHIPFILE'] = (chip2_filename,'Source data filename')\n",
    "    new_hdul[3].header['CHIPFILE'] = (chip3_filename,'Source data filename')\n",
    "\n",
    "    compress_file_fpack(data_path+chip1_filename)\n",
    "    compress_file_fpack(data_path+chip2_filename)\n",
    "    compress_file_fpack(data_path+chip3_filename)\n",
    "\n",
    "    new_hdul.writeto(data_path+output_filename,overwrite=True)\n",
    "    with fits.open(data_path+output_filename) as hdu:\n",
    "        hdu.info()\n",
    "    \n",
    "    if os.path.exists(data_path+output_filename+'.fz'):\n",
    "        os.remove(data_path+output_filename+'.fz')\n",
    "    compress_file_fpack(data_path+output_filename)\n",
    "        \n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6b8a4f6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_mef_science(data_path,raw_filename_fz,obs_data_table_filepath,log_filepath):\n",
    "\n",
    "    # Set filenames\n",
    "    chip1_filename_fz = 'n'+raw_filename_fz[1:9]+'.'+raw_filename_fz[11:14]+'.chip1.otzfc.fits.fz'\n",
    "    chip2_filename_fz = 'n'+raw_filename_fz[1:9]+'.'+raw_filename_fz[11:14]+'.chip2.otzfc.fits.fz'\n",
    "    chip3_filename_fz = 'n'+raw_filename_fz[1:9]+'.'+raw_filename_fz[11:14]+'.chip3.otzfc.fits.fz'\n",
    "\n",
    "    # Get calib filenames\n",
    "    os.chdir(data_path)\n",
    "    for filename_fz in glob.glob('G?-CAL*.bias.fits.fz'):\n",
    "        bias_filename = filename_fz[:-3]\n",
    "    for filename_fz in glob.glob('G?-CAL*.twiskyflat*.fits.fz'):\n",
    "        flatfield_filename = filename_fz[:-3]\n",
    "    \n",
    "    # Unpack fz files\n",
    "    raw_filename   = raw_filename_fz[:-3]\n",
    "    chip1_filename = chip1_filename_fz[:-3]\n",
    "    chip2_filename = chip2_filename_fz[:-3]\n",
    "    chip3_filename = chip3_filename_fz[:-3]\n",
    "    decompress_file_fpack(data_path+raw_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip1_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip2_filename_fz)\n",
    "    decompress_file_fpack(data_path+chip3_filename_fz)\n",
    "\n",
    "    # Get metadata\n",
    "    if raw_filename_fz[:1] == 'N':\n",
    "        object_name,object_desig,obs_date,obs_time,exp_time,filter_name = get_metadata_geminiN(data_path+chip2_filename)\n",
    "        telescope_name = 'Gemini-N'\n",
    "    elif raw_filename_fz[:1] == 'S':\n",
    "        object_name,object_desig,obs_date,obs_time,exp_time,filter_name = get_metadata_geminiS(data_path+chip2_filename)\n",
    "        telescope_name = 'Gemini-S'\n",
    "        \n",
    "    # Get raw fits header data and Gemini program ID\n",
    "    with fits.open(data_path+raw_filename) as hdulist:\n",
    "        hdr0 = hdulist[0].header\n",
    "        gemprgid = hdr0['GEMPRGID']\n",
    "        ra_obs = hdr0['RA']\n",
    "        dec_obs = hdr0['DEC']\n",
    "        hdu0 = fits.PrimaryHDU(header=hdr0)\n",
    "        output_filename = gemprgid + '_' + raw_filename_fz[1:9]+'.'+raw_filename_fz[11:14]+'.otzfc.fits'\n",
    "\n",
    "    # Retrieve geometry data for objects in images\n",
    "    ephem_params = get_geometry_params(object_name,obs_date,obs_time,exp_time)\n",
    "    \n",
    "    # Write metadata and geometry data to output file\n",
    "    with open(obs_data_table_filepath,'a') as of:\n",
    "        of.write('{:<40s}   {:<10s}   {:<20s}   {:<35s}   {:s}   {:s}   {:>7d}   {:<6s}   {:>8s}   {:>9s}   {:9.5f}   {:9.5f}   {:7.1f}   {:8.1f}   {:6.4f}   {:6.4f}   {:6.1f}   {:8.1f}   {:8.1f}   {:9.1f}   {:9.1f}   {:8.1f}   {:6.1f}   {:6.1f}   {:8.1f}\\n'.format(output_filename,\n",
    "          telescope_name,object_name,object_desig,obs_date,obs_time,\n",
    "          int(round(exp_time)),filter_name,ephem_params['ra_hms'],ephem_params['dec_dms'],\n",
    "          ephem_params['ra'],ephem_params['dec'],ephem_params['ra_rate'],ephem_params['dec_rate'],\n",
    "          ephem_params['heliodist'],ephem_params['geodist'],ephem_params['phsang'],\n",
    "          ephem_params['lunar_elong'],ephem_params['lunar_illum'],\n",
    "          ephem_params['antisolarPA'],ephem_params['neghelioPA'],\n",
    "          ephem_params['orbplang'],ephem_params['glxlon'],ephem_params['glxlat'],\n",
    "          ephem_params['trueanom']))\n",
    "\n",
    "    # Extract data from processed files\n",
    "    with fits.open(data_path+chip1_filename) as hdu_chip1, fits.open(data_path+chip2_filename) as hdu_chip2, fits.open(data_path+chip3_filename) as hdu_chip3:\n",
    "        hdr1,data1 = hdu_chip1[0].header,hdu_chip1[0].data\n",
    "        hdr2,data2 = hdu_chip2[0].header,hdu_chip2[0].data\n",
    "        hdr3,data3 = hdu_chip3[0].header,hdu_chip3[0].data\n",
    "    \n",
    "    # Combine data from processed files\n",
    "    hdu1 = fits.ImageHDU(data1)\n",
    "    hdu2 = fits.ImageHDU(data2)\n",
    "    hdu3 = fits.ImageHDU(data3)\n",
    "    new_hdul = fits.HDUList([hdu0,hdu1,hdu2,hdu3])\n",
    "    \n",
    "    # Create headers\n",
    "    new_hdul[1].header             = new_hdul[1].header+hdr1\n",
    "    new_hdul[2].header             = new_hdul[2].header+hdr2\n",
    "    new_hdul[3].header             = new_hdul[3].header+hdr3\n",
    "    new_hdul[0].header['LONGNAME'] = (object_desig,'Extended target designation')\n",
    "    new_hdul[0].header['TARG_RA']  = (ephem_params['ra'],'Target RA from Horizons (2024-04-06)')\n",
    "    new_hdul[0].header['TARG_DEC'] = (ephem_params['dec'],'Target Dec from Horizons (2024-04-06)')\n",
    "    new_hdul[0].header['RA_HMS']   = (ephem_params['ra_hms'],'Target RA (HH:MM:SS) from Horizons')\n",
    "    new_hdul[0].header['DEC_DMS']  = (ephem_params['dec_dms'],'Target Dec (DD:MM:SS) from Horizons')\n",
    "    new_hdul[0].header['RA_RATE']  = (ephem_params['ra_rate'],'Non-sidereal rate in RA (arcsec/hr)')\n",
    "    new_hdul[0].header['DEC_RATE'] = (ephem_params['dec_rate'],'Non-sidereal rate in Dec (arcsec/hr)')\n",
    "    new_hdul[0].header['R_DIST']   = (ephem_params['heliodist'],'Heliocentric distance (au)')\n",
    "    new_hdul[0].header['D_DIST']   = (ephem_params['geodist'],'Geocentric distance (au)')\n",
    "    new_hdul[0].header['PHSANG']   = (ephem_params['phsang'],'Solar phase angle (deg)')\n",
    "    new_hdul[0].header['LUNELONG'] = (ephem_params['lunar_elong'],'Lunar elongation (deg)')\n",
    "    new_hdul[0].header['LUNILLUM'] = (ephem_params['lunar_illum'],'Lunar illumination (percent)')\n",
    "    new_hdul[0].header['ANTISUNV'] = (ephem_params['antisolarPA'],'Sun-to-target vector PA (deg E of N)')\n",
    "    new_hdul[0].header['NEGHVELV'] = (ephem_params['neghelioPA'],'Neg. heliocen. vel. vector PA (deg E of N)')\n",
    "    new_hdul[0].header['ORBPLANG'] = (ephem_params['orbplang'],'Orbit plane angle (deg)')\n",
    "    new_hdul[0].header['GLXLON']   = (ephem_params['glxlon'],'Galactic longitude (deg)')\n",
    "    new_hdul[0].header['GLXLAT']   = (ephem_params['glxlat'],'Galactic latitude (deg)')\n",
    "    new_hdul[0].header['TRUEANOM'] = (ephem_params['trueanom'],'True anomaly (deg)')\n",
    "    new_hdul[0].header['BIAS_MEF'] = (bias_filename,'Bias MEF file')\n",
    "    new_hdul[0].header['FLAT_MEF'] = (flatfield_filename,'Flatfield MEF file')\n",
    "    new_hdul[1].header['FLATCHIP'] = new_hdul[1].header['FLATUSED']\n",
    "    new_hdul[2].header['FLATCHIP'] = new_hdul[2].header['FLATUSED']\n",
    "    new_hdul[3].header['FLATCHIP'] = new_hdul[3].header['FLATUSED']\n",
    "    del new_hdul[1].header['FLATUSED']\n",
    "    del new_hdul[2].header['FLATUSED']\n",
    "    del new_hdul[3].header['FLATUSED']\n",
    "\n",
    "    compress_file_fpack(data_path+raw_filename)\n",
    "    compress_file_fpack(data_path+chip1_filename)\n",
    "    compress_file_fpack(data_path+chip2_filename)\n",
    "    compress_file_fpack(data_path+chip3_filename)\n",
    "\n",
    "    new_hdul.writeto(data_path+output_filename,overwrite=True)\n",
    "    \n",
    "    if os.path.exists(data_path+output_filename+'.fz'):\n",
    "        os.remove(data_path+output_filename+'.fz')\n",
    "    \n",
    "    compress_file_fpack(data_path+output_filename)\n",
    "    \n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0fe8f35",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_mef_calib_batch(base_path):\n",
    "    os.chdir(base_path)\n",
    "    for data_dirname in sorted(glob.glob('ut*_gemini?')):\n",
    "        print('Processing {:s}...'.format(data_dirname))\n",
    "        calib_data_path = base_path + data_dirname + '/calib/'\n",
    "        os.chdir(calib_data_path)\n",
    "        for chip03_filename_fz in glob.glob('n*.bias.03.fits.fz'):\n",
    "            create_mef_bias(calib_data_path,chip03_filename_fz)\n",
    "        for chip1_filename_fz in glob.glob('n*.gmos?.twiskyflat.*.chip1.fits.fz'):\n",
    "            create_mef_flatfield(calib_data_path,chip1_filename_fz)\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efd7227c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_mef_science_batch(base_path,data_set):\n",
    "    \n",
    "    log_filepath = base_path + 'log_make_multiextension_files_{:s}_{:s}.txt'.format(data_set,datetime.datetime.today().strftime('%Y%m%d_%H%M%S'))\n",
    "    with open(log_filepath,'w') as log_file:\n",
    "        log_file.write('Log: notebook_GMOS_a04_make_multiextension_files.ipynb\\n')\n",
    "\n",
    "    obs_data_table_filename = 'target_data_table_{:s}.txt'.format(data_set)\n",
    "    with open(base_path+obs_data_table_filename,'w') as of:\n",
    "        of.write('File Name                                  Telescope    Target                 LongDesig                             ObsDate      ObsStartUT   ExpTime   Filter   RA (HMS)   Dec (DMS)    RA (deg)   Dec (deg)   RA_rate   Dec_rate        r    delta   PhsAng   LunElong   LunIllum   AntiSolPA   NegHVelPA   OrbPlAng   GlxLon   GlxLat   TrueAnom\\n')\n",
    "    output_write_log_entry(log_filepath,'Target data output file: {:s}'.format(base_path+obs_data_table_filename))\n",
    "    \n",
    "    os.chdir(base_path)\n",
    "    output_write_log_entry(log_filepath,'Processing files in {:s}...'.format(base_path))\n",
    "    for raw_filename_fz in sorted(glob.glob('?20??????S????.fits.fz')):\n",
    "        output_write_log_entry(log_filepath,'Processing {:s}...'.format(raw_filename_fz))\n",
    "        create_mef_science(base_path,raw_filename_fz,base_path+obs_data_table_filename,log_filepath)\n",
    "\n",
    "    output_write_log_entry(log_filepath,'>>> Processing files in {:s} complete.'.format(base_path))\n",
    "    return None\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec870f93",
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_mef_files_all(base_path,data_set):\n",
    "    os.chdir(base_path)\n",
    "    for chip03_filename_fz in glob.glob('n*.bias.03.fits.fz'):\n",
    "        create_mef_bias(base_path,chip03_filename_fz)\n",
    "    for chip1_filename_fz in glob.glob('n*.gmos?.twiskyflat.*.chip1.fits.fz'):\n",
    "        create_mef_flatfield(base_path,chip1_filename_fz)\n",
    "    create_mef_science_batch(base_path,data_set)\n",
    "    return None\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a18eb92c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# data_set: user-specified designation to be used in the filename of the target data table\n",
    "\n",
    "create_mef_files_all(base_path,data_set)\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.7.16"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}