In [1]:
from astropy.coordinates import SkyCoord, get_moon, EarthLocation, ICRS, GCRS,AltAz
from astropy.io.fits.hdu.hdulist import HDUList
from datetime import date
from astral import Astral, Location
import astropy.units as u

In [1]:
import pyfits
import pandas as pd
import numpy as np
from itertools import chain
from collections import Counter
from scipy.interpolate import interp1d
from datetime import datetime
from astropy.time import Time

In [2]:
## Finding out which file years to download
all_sky=pd.read_csv('objSKY.csv')
all_sky.head()
print('earliest date=', min(all_sky['MJD']), 'or 12-11-2009')
print('latest date=', max(all_sky['MJD']), 'or 05-12-2016')


earliest date= 55176 or 12-11-2009
latest date= 57520 or 05-12-2016

Data


In [3]:
## To check for reading in the right file
rand_files=('spec-3586-55181-0496.fits','spec-3586-55181-0788.fits','spec-3586-55181-0996.fits',
            'spec-3761-55272-0008.fits','spec-10000-57346-0334.fits','spec-3761-55272-0475.fits',
            'spec-10000-57346-0659.fits','spec-5478-56014-0654.fits','spec-10000-57346-0955.fits',
            'spec-5478-56014-0716.fits','spec-10000-57346-0994.fits')

## To look for patterns
inseq_files=('spec-3663-55176-0010.fits','spec-3663-55176-0012.fits','spec-3663-55176-0020.fits',
             'spec-3663-55176-0024.fits','spec-3663-55176-0036.fits','spec-3663-55176-0038.fits',
             'spec-3663-55176-0048.fits','spec-3663-55176-0052.fits','spec-3663-55176-0056.fits',
             'spec-3663-55176-0068.fits','spec-3663-55176-0075.fits','spec-3663-55176-0078.fits',
             'spec-3663-55176-0090.fits','spec-3663-55176-0094.fits','spec-3663-55176-0096.fits',
             'spec-3663-55176-0108.fits','spec-3663-55176-0112.fits','spec-3663-55176-0114.fits',
             'spec-3663-55176-0128.fits','spec-3663-55176-0134.fits')

solar_files=('g2009.txt','g2010.txt','g2011.txt','g2012.txt','g2013.txt','g2014.txt','g2015.txt','g2016.txt')

Getting MJDs from fits


In [4]:
#a) first load all files and convert all dates into mjds.
def get_mjd(file):
    da=pyfits.open(file)
    return (da[4].header['MJD'])

In [5]:
for i in rand_files:
    print(get_mjd(i))


55181
55181
55181
55272
57345
55272
57345
56014
57345
56014
57345

Master List of Solar File MJDs & Coverage


In [6]:
def merge_subs(lst_of_lsts):
    res = []
    for row in lst_of_lsts:
        for i, resrow in enumerate(res):
            if row[:3]==resrow[:3]:
                res[i] += row[1:]
                break
        else:
            res.append(row)
    return (res)
## Merges lists, but does not remove first 3 cols

In [7]:
def f2(seq): 
    # order preserving
    checked = []
    for e in seq:
        if e not in checked:
            checked.append(e)
    return checked
## removes duplicates

In [8]:
## Enter solar files (the range of years)
def mjds_sol(txt_file, fits_files):
    master_list=[]
    for i in txt_file:
        ## We only need year, month, day of month, and mill. of sol. disk
        ## Cols         1-4    5-6      7-8               31-34
        demo=list(chain.from_iterable((x[:4], x[4:6], x[6:10], x[31:34]) for x in open(i).readlines()))
        ## Organizing cols
        str_list=[demo[i*4:i*4+4] for i in range(int(len(demo)/4))]
        ## Converting to float
        float_dates=[]
        for k in range(len(str_list)):
            float_dates.append([float(i) for i in str_list[k]])
        b=([x[:4] for x in float_dates])
        ve=merge_subs(b)
        new=[]
        summed_mills=[]
        no_mills=[]
        merged_lists=[]
        merged_utc=[]
        for i in range(len(ve)):
            new.append(f2(ve[i]))
            summed_mills.append(sum(new[i][3:]))
            no_mills.append(new[i][:3])
            merged_lists.append(np.append(no_mills[i],summed_mills[i]).tolist())
            ## Conversion to datetime ##
            merged_lists[i][3:3]=[12]
            merged_utc.append([int(n) for n in merged_lists[i]])
        utc_df=pd.DataFrame(merged_utc)
        utc_df.columns=['year','month','day','hour','mill']
        no_mills_mjds=[]
        vt=[]
        merged_mjds=[]
        ## Conversion to MJDs ##
        for i in range(len(utc_df)):
            vt.append(Time(datetime(utc_df['year'][i],utc_df['month'][i],utc_df['day'][i],utc_df['hour'][i]), scale='utc'))
            no_mills_mjds.append(vt[i].mjd)
            merged_mjds.append(np.append(no_mills_mjds[i],summed_mills[i]).tolist())
            master_list.append(merged_mjds[i]) ## lists all info from solar files
        ## preparing to interpolate
    master_df=pd.DataFrame(master_list)
    master_df.columns=['mjd','sun_cov']
    mjd=master_df['mjd']
    sun_cov=master_df['sun_cov']
    ifunc=interp1d(mjd,sun_cov)
    seq_sol=[ifunc(i).tolist() for i in [get_mjd(i) for i in fits_files]]
    return(seq_sol)

In [9]:
## enter solar files used and the fits files used, return mill of sol disk
mjds_sol(solar_files, rand_files)
#solar_files[0]


Out[9]:
[224.5, 224.5, 224.5, 138.0, 123.0, 138.0, 123.0, 550.5, 123.0, 550.5, 123.0]

IGNORE EVERYTHING BELOW HERE


In [12]:
master_df=pd.DataFrame(master_list)
master_df.columns=['mjd','sun_cov']

In [13]:
master_df[:8]


Out[13]:
mjd sun_cov
0 54840.5 30.0
1 54841.5 51.0
2 54842.5 92.0
3 54843.5 18.0
4 54844.5 17.0
5 54850.5 18.0
6 54873.5 9.0
7 54874.5 13.0

In [37]:
mjd=master_df['mjd']
sun_cov=master_df['sun_cov']
#y=merged_lists[5]
ifunc=interp1d(mjd,sun_cov)
ifunc(54849)


Out[37]:
array(17.75)

In [44]:
seq_files_mjd=[]
rand_files_mjd=[]
for i in inseq_files:
    seq_files_mjd.append(get_mjd(i))
for i in rand_files:
    rand_files_mjd.append(get_mjd(i))

In [51]:
seq_files_mjd=[get_mjd(i) for i in inseq_files]
rand_files_mjd=[get_mjd(i) for i in rand_files]

In [52]:
seq_sol=[ifunc(i).tolist() for i in seq_files_mjd]
rand_sol=[ifunc(i).tolist() for i in rand_files_mjd]

In [53]:
rand_files_mjd


Out[53]:
[55181, 55181, 55181, 55272, 57345, 55272, 57345, 56014, 57345, 56014, 57345]

In [54]:
rand_sol


Out[54]:
[224.5, 224.5, 224.5, 138.0, 123.0, 138.0, 123.0, 550.5, 123.0, 550.5, 123.0]

Trial codes


In [107]:
#b) Then set up one massive interpolated function from mjds to mill of sol. disk
## Pulling out cols of interest
master_list=[]
for i in solar_files:
    demo=list(chain.from_iterable((x[:4], x[4:6], x[6:10], x[31:34]) for x in open(i).readlines()))
    ## Organizing cols
    str_list=[demo[i*4:i*4+4] for i in range(int(len(demo)/4))]
    ## Converting to float
    float_dates=[]
    for k in range(len(str_list)):
        float_dates.append([float(i) for i in str_list[k]])
    b=([x[:4] for x in float_dates])
    ve=merge_subs(b)
    new=[]
    summed_mills=[]
    no_mills=[]
    merged_lists=[]
    merged_utc=[]
    for i in range(len(ve)):
        new.append(f2(ve[i]))
        summed_mills.append(sum(new[i][3:]))
        no_mills.append(new[i][:3])
        merged_lists.append(np.append(no_mills[i],summed_mills[i]).tolist())
        ## Conversion to datetime ##
        merged_lists[i][3:3]=[12]
        merged_utc.append([int(n) for n in merged_lists[i]])
    utc_df=pd.DataFrame(merged_utc)
    utc_df.columns=['year','month','day','hour','mill']
    no_mills_mjds=[]
    vt=[]
    merged_mjds=[]
    ## Conversion to MJDs ##
    for i in range(len(utc_df)):
        vt.append(Time(datetime(utc_df['year'][i],utc_df['month'][i],utc_df['day'][i],utc_df['hour'][i]), scale='utc'))
        no_mills_mjds.append(vt[i].mjd)
        merged_mjds.append(np.append(no_mills_mjds[i],summed_mills[i]).tolist())
        master_list.append(merged_mjds[i]) ## lists all info from solar files
    ## preparing to interpolate
master_df=pd.DataFrame(master_list)
master_df.columns=['mjd','sun_cov']
mjd=master_df['mjd']
sun_cov=master_df['sun_cov']
ifunc=interp1d(mjd,sun_cov)
seq_sol=[ifunc(i).tolist() for i in [get_mjd(i) for i in rand_files]] #fits_files]]
seq_sol


Out[107]:
[224.5, 224.5, 224.5, 138.0, 123.0, 138.0, 123.0, 550.5, 123.0, 550.5, 123.0]

In [329]:
#list_example=[[1,1,2016,212], [1,2,2016,170], [1,3,2016,150], [1,5,2016,96], [1,6,2016,150], [1,8,2016,321]]

In [330]:
#zip(list_example)

In [ ]:


In [ ]: