<commit_before>"""Perform regular monitoring of the COS FUV and NUV dark rates

"""

from __future__ import print_function, absolute_import, division

import os
import datetime
import numpy as np
import shutil
import glob

from astropy.io import fits
from astropy.time import Time

from .solar import get_solar_data
from .plotting import plot_histogram, plot_time, plot_orbital_rate

from ..utils import corrtag_image
from ..database.db_tables import open_settings, load_connection

base_dir = '/grp/hst/cos/Monitors/Darks/'
web_directory = '/grp/webpages/COS/'

#-------------------------------------------------------------------------------

def get_temp(filename):
    """Get detector temperture during observation from spt filename

    Parameters
    ----------
    filename : str
        FITS file for which the temperature is to be Found

    Returns
    -------
    temperature : float
        Detector temperature at the time of the observation

    """

    with fits.open(filename) as hdu:
        detector = hdu[0].header['DETECTOR']
        segment = hdu[0].header['SEGMENT']
        rootname = hdu[0].header['ROOTNAME']

    if detector == 'FUV' and segment == 'FUVA':
        temp_keyword = 'LDCAMPAT'
    elif detector == 'FUV' and segment == 'FUVB':
        temp_keyword = 'LDCAMPBT'
    elif detector == 'NUV':
        temp_keyword = 'LMMCETMP'
    else:
        raise ValueError('What??? {} {}'.format(detector, segment))

    path, name = os.path.split(filename)
    spt_file = os.path.join(path, rootname + '_spt.fits')

    try:
        temperature = fits.getval(spt_file, temp_keyword, ext=2)
    except IOError:
        temperature = fits.getval(spt_file + '.gz', temp_keyword, ext=2)

    return temperature

#-------------------------------------------------------------------------------

def mjd_to_decyear( time_array ):
    """ pull this out when you get it into astropy.time

    """

    times = Time( time_array, scale='tt', format='mjd' )

    out_times = []
    for value in times:
        year = value.datetime.year
        n_days = (value.datetime - datetime.datetime(value.datetime.year, 1, 1)).total_seconds()
        total_days = (datetime.datetime(value.datetime.year+1, 1, 1) - datetime.datetime(value.datetime.year, 1, 1)).total_seconds()

        fraction = float(n_days) / total_days

        out_times.append(year + fraction)

    return np.array(out_times)

#-------------------------------------------------------------------------------

def pull_orbital_info(dataset, step=25):
    """ Pull second by second orbital information from the dataset

    """

    SECOND_PER_MJD = 1.15741e-5

    info = {}
    info['obsname'] = os.path.split(dataset)[-1]

    hdu = fits.open(dataset)
    try:
        timeline = hdu['timeline'].data
        segment = hdu[0].header['segment']
    except KeyError:
        return

    if segment == 'N/A':
        segment = 'NUV'
        xlim = (0, 1024)
        ylim = (0, 1204)
        pha = (-1, 1)
    elif segment == 'FUVA':
        xlim = (1200, 15099)
        ylim = (380, 680)
        pha = (2, 23)
    elif segment == 'FUVB':
        xlim = (950, 15049)
        ylim = (440, 720)
        pha = (2, 23)
    else:
        raise ValueError('What segment is this? {}'.format(segment))

    info['rootname'] = hdu[0].header['rootname']
    info['detector'] = segment
    info['temp'] = get_temp(dataset)

    times = timeline['time'][::step].copy()
    lat = timeline['latitude'][:-1][::step].copy().astype(np.float64)
    lon = timeline['longitude'][:-1][::step].copy().astype(np.float64)
    try:
        sun_lat = timeline['sun_lat'][:-1][::step].copy().astype(np.float64)
        sun_lon = timeline['sun_lon'][:-1][::step].copy().astype(np.float64)
    except KeyError:
        sun_lat = lat.copy() * 0
        sun_lon = lat.copy() * 0


    mjd = hdu[1].header['EXPSTART'] + \
        times.copy()[:-1].astype(np.float64) * \
        SECOND_PER_MJD

    decyear = mjd_to_decyear(mjd)

    if not len(times):
        blank = np.array([0])
        return


    events = hdu['events'].data
    filtered_index = np.where((events['PHA'] > pha[0]) &
                              (events['PHA'] < pha[1]) &
                              (events['XCORR'] > xlim[0]) &
                              (events['XCORR'] < xlim[1]) &
                              (events['YCORR'] > ylim[0]) &
                              (events['YCORR'] < ylim[1]))

    ta_index = np.where((events['XCORR'] > xlim[0]) &
                        (events['XCORR'] < xlim[1]) &
                        (events['YCORR'] > ylim[0]) &
                        (events['YCORR'] < ylim[1]))


    counts = np.histogram(events[filtered_index]['time'], bins=times)[0]
    ta_counts = np.histogram(events[ta_index]['time'], bins=times)[0]

    npix = float((xlim[1] - xlim[0]) * (ylim[1] - ylim[0]))
    counts = counts / npix / step
    ta_counts = ta_counts / npix / step

    if not len( lat ) == len(counts):
        lat = lat[:-1]
        lon = lon[:-1]
        sun_lat = sun_lat[:-1]
        sun_lon = sun_lon[:-1]

    assert len(lat) == len(counts), \
        'Arrays are not equal in length {}:{}'.format(len(lat), len(counts))


    for i in range(len(counts)):
        ### - better solution than round?
        info['date'] = round(decyear[i], 3)
        info['dark'] = round(counts[i], 7)
        info['ta_dark'] = round(ta_counts[i], 7)
        info['latitude'] = round(lat[i], 7)
        info['longitude'] = round(lon[i], 7)
        info['sun_lat'] = round(sun_lat[i], 7)
        info['sun_lon'] = round(sun_lon[i], 7)

        yield info

#-------------------------------------------------------------------------------

def compile_phd():
<commit_msg>pull_orbital_info needed conditions to yield for corner-case datasets<commit_after>"""Perform regular monitoring of the COS FUV and NUV dark rates

"""

from __future__ import print_function, absolute_import, division

import os
import datetime
import numpy as np
import shutil
import glob

from astropy.io import fits
from astropy.time import Time

from .solar import get_solar_data
from .plotting import plot_histogram, plot_time, plot_orbital_rate

from ..utils import corrtag_image
from ..database.db_tables import open_settings, load_connection

base_dir = '/grp/hst/cos/Monitors/Darks/'
web_directory = '/grp/webpages/COS/'

#-------------------------------------------------------------------------------

def get_temp(filename):
    """Get detector temperture during observation from spt filename

    Parameters
    ----------
    filename : str
        FITS file for which the temperature is to be Found

    Returns
    -------
    temperature : float
        Detector temperature at the time of the observation

    """

    with fits.open(filename) as hdu:
        detector = hdu[0].header['DETECTOR']
        segment = hdu[0].header['SEGMENT']
        rootname = hdu[0].header['ROOTNAME']

    if detector == 'FUV' and segment == 'FUVA':
        temp_keyword = 'LDCAMPAT'
    elif detector == 'FUV' and segment == 'FUVB':
        temp_keyword = 'LDCAMPBT'
    elif detector == 'NUV':
        temp_keyword = 'LMMCETMP'
    else:
        raise ValueError('What??? {} {}'.format(detector, segment))

    path, name = os.path.split(filename)
    spt_file = os.path.join(path, rootname + '_spt.fits')

    try:
        temperature = fits.getval(spt_file, temp_keyword, ext=2)
    except IOError:
        temperature = fits.getval(spt_file + '.gz', temp_keyword, ext=2)

    return temperature

#-------------------------------------------------------------------------------

def mjd_to_decyear( time_array ):
    """ pull this out when you get it into astropy.time

    """

    times = Time( time_array, scale='tt', format='mjd' )

    out_times = []
    for value in times:
        year = value.datetime.year
        n_days = (value.datetime - datetime.datetime(value.datetime.year, 1, 1)).total_seconds()
        total_days = (datetime.datetime(value.datetime.year+1, 1, 1) - datetime.datetime(value.datetime.year, 1, 1)).total_seconds()

        fraction = float(n_days) / total_days

        out_times.append(year + fraction)

    return np.array(out_times)

#-------------------------------------------------------------------------------

def pull_orbital_info(dataset, step=25):
    """ Pull second by second orbital information from the dataset

    """

    SECOND_PER_MJD = 1.15741e-5

    info = {}
    info['obsname'] = os.path.split(dataset)[-1]

    hdu = fits.open(dataset)
    try:
        timeline = hdu['timeline'].data
        segment = hdu[0].header['segment']
    except KeyError:
        yield info
        raise StopIteration

    if segment == 'N/A':
        segment = 'NUV'
        xlim = (0, 1024)
        ylim = (0, 1204)
        pha = (-1, 1)
    elif segment == 'FUVA':
        xlim = (1200, 15099)
        ylim = (380, 680)
        pha = (2, 23)
    elif segment == 'FUVB':
        xlim = (950, 15049)
        ylim = (440, 720)
        pha = (2, 23)
    else:
        raise ValueError('What segment is this? {}'.format(segment))

    info['rootname'] = hdu[0].header['rootname']
    info['detector'] = segment
    info['temp'] = get_temp(dataset)

    times = timeline['time'][::step].copy()
    lat = timeline['latitude'][:-1][::step].copy().astype(np.float64)
    lon = timeline['longitude'][:-1][::step].copy().astype(np.float64)
    try:
        sun_lat = timeline['sun_lat'][:-1][::step].copy().astype(np.float64)
        sun_lon = timeline['sun_lon'][:-1][::step].copy().astype(np.float64)
    except KeyError:
        sun_lat = lat.copy() * 0
        sun_lon = lat.copy() * 0


    mjd = hdu[1].header['EXPSTART'] + \
        times.copy()[:-1].astype(np.float64) * \
        SECOND_PER_MJD

    decyear = mjd_to_decyear(mjd)

    if not len(times):
        blank = np.array([0])
        yield info
        raise StopIteration

    events = hdu['events'].data
    filtered_index = np.where((events['PHA'] > pha[0]) &
                              (events['PHA'] < pha[1]) &
                              (events['XCORR'] > xlim[0]) &
                              (events['XCORR'] < xlim[1]) &
                              (events['YCORR'] > ylim[0]) &
                              (events['YCORR'] < ylim[1]))

    ta_index = np.where((events['XCORR'] > xlim[0]) &
                        (events['XCORR'] < xlim[1]) &
                        (events['YCORR'] > ylim[0]) &
                        (events['YCORR'] < ylim[1]))


    counts = np.histogram(events[filtered_index]['time'], bins=times)[0]
    ta_counts = np.histogram(events[ta_index]['time'], bins=times)[0]

    npix = float((xlim[1] - xlim[0]) * (ylim[1] - ylim[0]))
    counts = counts / npix / step
    ta_counts = ta_counts / npix / step

    if not len(lat) == len(counts):
        lat = lat[:-1]
        lon = lon[:-1]
        sun_lat = sun_lat[:-1]
        sun_lon = sun_lon[:-1]

    assert len(lat) == len(counts), \
        'Arrays are not equal in length {}:{}'.format(len(lat), len(counts))

    if not len(counts):
        yield info
    else:
        for i in range(len(counts)):
            ### - better solution than round?
            info['date'] = round(decyear[i], 3)
            info['dark'] = round(counts[i], 7)
            info['ta_dark'] = round(ta_counts[i], 7)
            info['latitude'] = round(lat[i], 7)
            info['longitude'] = round(lon[i], 7)
            info['sun_lat'] = round(sun_lat[i], 7)
            info['sun_lon'] = round(sun_lon[i], 7)

            yield info

#-------------------------------------------------------------------------------

def compile_phd():

{'fraction_non_alphanumeric': 0.19811077715757835, 'fraction_numerical': 0.020094461142121083, 'mean_word_length': 1.9915232468533264, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '5604485', 'n_tokens_mistral': 3888, 'n_tokens_neox': 3659, 'n_words': 824}