# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# +
import os, sys
import scipy
import numpy as np
from astropy.io import fits
import matplotlib.pyplot as plt
from tqdm import tqdm_notebook
from astropy.table import Table
from scipy.ndimage import minimum_filter1d
from scipy.ndimage.filters import percentile_filter

plt.rcParams['font.size'] = 5
# -

d = 23
master_log = master_log = Table.read('/Users/arcticfox/Documents/youngStars/veloce/master_log.tab', format='ascii')
date = '2020-11-{0}'.format(d)
directory = '2011{0}'.format(d)
table = master_log[master_log['ObsDate']==date]
fileformat = '{0}nov3{1:04d}.fits'

table[table['Frame']==110]

# +
files = np.sort([i for i in os.listdir(directory) if i.endswith('.fits')])

science_frames, bias_frames, dark_frames, flat_frames = [], [], [], []

for i in range(len(table)):
    if 'TIC' in table['ObjType'][i]:
        science_frames.append(fileformat.format(d, table['Frame'][i]))
    elif 'BiasFrame' == table['ObjType'][i]:
        bias_frames.append(fileformat.format(d, table['Frame'][i]))
    elif 'DarkFrame' in table['ObjType'][i]:
        dark_frames.append(fileformat.format(d, table['Frame'][i]))
    elif 'FlatField' in table['ObjType'][i]:
        flat_frames.append(fileformat.format(d, table['Frame'][i]))
    else:
        continue
        
#dark_inds = dark_inds[np.argwhere(np.diff(dark_inds)>10)[0][0]+1:]
#bias_inds = bias_inds[4:]
#flat_inds = flat_inds[2:]

#science_frames = np.unique(np.sort([os.path.join(directory, i) for i in science_frames]))
bias_frames = np.unique(np.sort([os.path.join(directory, i) for i in bias_frames]))
dark_frames = np.unique(np.sort([os.path.join(directory, i) for i in dark_frames]))[6:46]#[27:-5]
flat_frames = np.unique(np.sort([os.path.join(directory, i) for i in flat_frames]))[20:]#[23:]
# -

len(bias_frames), len(dark_frames), len(flat_frames)


# ## Creating master frames

def master_file(files, output_fn, fntype='dark'):
    arrs = []
    
    for fn in tqdm_notebook(files):
        hdu = fits.open(fn)
        if hdu[0].data.shape == (4112, 4202):
            arrs.append(hdu[0].data)
        
        hdu.close()
    arrs = np.array(arrs)

    
    if fntype == 'bias' or fntype == 'dark':
        masked = np.copy(arrs) + 0.0
        for i in range(len(arrs)):
            rows, cols = np.where(arrs[i]>1000)
            masked[i][rows,cols] = np.nan
        masked = np.array(masked)

        med = np.nanmedian(masked, axis=0)   
    else:
        med = np.nanmedian(arrs, axis=0)
    np.save(output_fn, med)
    return med

# +
if 'dark_med.npy' not in os.listdir(directory):
    DARK_MED = master_file(dark_frames, 
                           os.path.join(directory, 'dark_med.npy'), 
                           fntype='dark')
else:
    DARK_MED = np.load(os.path.join(directory, 'dark_med.npy'))

if 'bias_med.npy' not in os.listdir(directory):
    BIAS_MED = master_file(bias_frames,
                           os.path.join(directory, 'bias_med.npy'),
                           fntype='bias')
else:
    BIAS_MED = np.load(os.path.join(directory, 'bias_med.npy'))
    
if 'flat_med.npy' not in os.listdir(directory):
    FLAT_MED = master_file(flat_frames,
                           os.path.join(directory, 'flat_med.npy'),
                           fntype='flat')
else:
    FLAT_MED = np.load(os.path.join(directory, 'flat_med.npy'))


# -

# ## Science Frames

def extract_science(files):
    outputfns = []
    arrs = []
    for fn in tqdm_notebook(files):
        hdu = fits.open(os.path.join(directory,fn))
        np.save(fn[:-5]+'.npy', hdu[0].data)
        outputfns.append(fn[:-5]+'.npy')
        arrs.append(hdu[0].data)
    return np.array(arrs), outputfns


science_arrs, science_files = extract_science(science_frames)

directory


# # Creating the dot models

# +
# #%matplotlib inline
def get_outliers(x_value, sigma=0.8, plot=False):
    
    arr = science_arrs[0][x_value:x_value+1][0] + 0.0
    x = np.arange(0,len(arr),1,dtype=int)
    outliers = np.where(arr >= np.nanmedian(arr) + sigma*np.nanstd(arr))[0]
    
    if plot:
        plt.figure(figsize=(1,1))
        plt.plot(x, arr, 'k')
        plt.ylim(800,1800)
        plt.plot(x[outliers], arr[outliers], 'o')
        plt.show()
        
    return outliers

def group_inds(values, sep):
    results = []

    for i, v in enumerate(values):
        if i == 0:
            mini = maxi = v
            temp = [v]
        else:
            # SETS 4 CADENCE LIMIT                                                                                                                     
            if (np.abs(v-maxi) <= sep):
                temp.append(v)
                if v > maxi:
                    maxi = v
                if v < mini:
                    mini = v
            else:
                results.append(int(np.nanmin(temp)))
                mini = maxi = v
                temp = [v]

            # GETS THE LAST GROUP                                                                                                                      
            if i == len(values)-1:
                results.append(int(np.nanmin(temp)))

    return np.array(results)


# +
rows, cols = np.where(science_arrs[0] > 1100)
mask = np.zeros(science_arrs[0].shape)

mask[rows,cols] = 1
# -

plt.imshow(DARK_MED[3000:3500,3000:3500], vmin=300, vmax=2000)

# %matplotlib notebook
plt.figure(figsize=(1,1))
plt.imshow(mask, cmap='Greys_r', vmin=0, vmax=1)#, alpha=0.9)
plt.plot(cols, rows, '.')
plt.colorbar()
plt.show()

# +
plt.figure(figsize=(1,1))
sargs = np.argsort(cols)
cols, rows = cols[sargs]+0, rows[sargs]+0
starts = np.where((rows<=795) & (rows>=770))[0]
ends = np.where((rows<=3330) & (rows>=3290))[0]

starts = group_inds(starts, sep=100)
ends = group_inds(ends, sep=100)

plt.plot(cols, rows, 'k.', ms=1)
plt.plot(cols[ends], rows[ends], '.', ms=1)
plt.plot(cols[starts], rows[starts], 'r.', ms=1)
# -

starts = np.delete(starts, [0, 35])
ends = np.delete(ends, [17, 39, 40])

len(starts), len(ends)

# +
mid_ends = np.where((rows>=2892) & (rows<=2908))[0]
mid_starts = np.where((rows>=1160) & (rows<=1180))[0]
mid = np.where((rows>=1995) & (rows<=2010))[0]

mid_starts = group_inds(mid_starts, sep=100)
mid_ends = group_inds(mid_ends, sep=100)
mid = group_inds(mid, sep=100)

plt.figure(figsize=(1,1))
ms = 3
plt.plot(cols, rows, 'k.', ms=1)
plt.plot(cols[mid_ends], rows[mid_ends], 'b.', ms=ms)
plt.plot(cols[starts], rows[starts], 'r.', ms=ms)
for i in range(len(mid_starts)):
    plt.plot(cols[mid_starts[i]], rows[mid_starts[i]], '.', ms=ms)
plt.plot(cols[ends], rows[ends], 'g.', ms=ms)
plt.plot(cols[mid], rows[mid], 'y.', ms=ms)
# -

len(starts), len(mid_starts), len(mid), len(mid_ends), len(ends)

#starts = np.delete(starts, [23])
mid_starts = np.delete(mid_starts, [27, 29, 33, 35, 38, 42, 43])
mid = np.delete(mid, [17, 38, 39])
mid_ends = np.delete(mid_ends, [24, 30, 37, 40, 41, 42])
#ends = np.delete(ends, [-1])

len(starts), len(mid_starts), len(mid), len(mid_ends), len(ends)

plt.figure(figsize=(1,1))
plt.plot(cols, rows, 'k.', ms=1)
plt.plot(cols[ends], rows[ends], 'b.', ms=1)
plt.plot(cols[starts], rows[starts], 'r.', ms=1)
plt.plot(cols[mid_starts], rows[mid_starts], '.', c='darkorange', ms=1)
plt.plot(cols[mid_ends], rows[mid_ends], 'g.', ms=1)
plt.plot(cols[mid], rows[mid], 'y.', ms=1)

dot_array = np.array([starts, mid_starts, mid, mid_ends, ends])
fit_x = np.arange(300, 4000, 1)

plt.figure(figsize=(1,1))
plt.plot(rows, cols, 'k.', ms=1)
models = np.zeros((len(mid), len(fit_x)))
for i in range(len(mid)):
    plt.plot(rows[dot_array[:,i]], cols[dot_array[:,i]], '.', ms=1)
    fit = np.polyfit(rows[dot_array[:,i]], cols[dot_array[:,i]], deg=2)
    model = np.poly1d(fit)
    plt.plot(fit_x, model(fit_x), lw=1)
    models[i] = model(fit_x)

np.save('./{0}/models.npy'.format(directory), models)

# ## Discretize model gap fits

discrete = np.zeros(models.shape, dtype=int)
for i in range(len(models)):
    discrete[i] = np.round(models[i])

np.save('./{0}/discrete_models.npy'.format(directory), discrete)

StarcoderdataJupyterScriptsDedupFiltered