Python Initiation

import os
import numpy as np
import pandas as pd
from astropy.table import Table
from scipy.stats import spearmanr#
from scipy.stats.stats import pearsonr

from astropy.cosmology import FlatLambdaCDM,Planck13,Planck15,z_at_value
from astropy import units as u
from astropy.cosmology import LambdaCDM
cosmo = LambdaCDM(H0=70, Om0=0.27, Ode0=0.73)


from astropy.time import Time
from astropy.time import TimeYearDayTime
from datetime import datetime
import time
from time import strftime,strptime
import calendar
from dateutil.parser import parse

#from adjustText import adjust_text
import matplotlib as mpl
import matplotlib.pyplot as plt
from pylab import cm
from collections import OrderedDict
#from adjustText import adjust_text

%matplotlib inline
%config InlineBackend.figure_format='svg'
# Edit the font, font size, and axes width
mpl.rcParams['font.family'] = 'Times New Roman' #'Avenir'
plt.rcParams['font.size'] = 18
plt.rcParams['axes.linewidth'] = 2

#Install a conda package in the current Jupyter kernel
import sys
!conda install --yes --prefix {sys.prefix} numpy


#Install a pip package in the current Jupyter kernel
import sys
!{sys.executable} -m pip install numpy

#from astropy.time import Time
from astropy.io import fits
from matplotlib.pyplot import MultipleLocator
import matplotlib.dates as mdates

def datetime2mjd(x):
    mjd_ref=59000
    mjd_minus_mdates_num=mdates.date2num(convert_xaxis_time(mjd_ref))-mjd_ref
    
    x=mdates.date2num(x)
    y = x - mjd_minus_mdates_num
    return y

def mjd2datetime(x):
    mjd_ref=59000
    mjd_minus_mdates_num=mdates.date2num(convert_xaxis_time(mjd_ref))-mjd_ref
    y= x + mjd_minus_mdates_num
    y= mdates.num2date(y)
    return y



def datenums2mjd(x):
    #x=mdates.date2num(x)
    mjd_ref=59000
    mjd_minus_mdates_num=mdates.date2num(convert_xaxis_time(mjd_ref))-mjd_ref
    y = x - mjd_minus_mdates_num
    return y

def mjd2numsdate(x):
    mjd_ref=59000
    mjd_minus_mdates_num=mdates.date2num(convert_xaxis_time(mjd_ref))-mjd_ref
    
    y= x + mjd_minus_mdates_num
    #y= mdates.num2date(y)
    return y


def convert_xaxis_mjd(time):
    return Time(time).mjd

def convert_xaxis_time(mjd):
    return Time(mjd,format='mjd').to_datetime()


def date2yday(x):
    """
        x is in matplotlib datenums, so they are floats.
        """
    y = x - mdates.date2num(datetime(2018, 1, 1))
    return y

def yday2date(x):
    """
        return a matplotlib datenum (x is days since start of year of 2018)
        """
    y = x + mdates.date2num(datetime(2018, 1, 1))
    return y


def convert_partial_year(numbers):
    datetimes=[]
    for number in numbers:
        year = int(number)
        d = timedelta(days=(number - year)*(365 + is_leap(year)))
        day_one = datetime(year,1,1)
        date = d + day_one
        datetimes.append(date)
    return datetimes


def is_leap(year):
    if not year%4 and  year%100 or not year%400:
        return True
    return False


def convert_mjd(times):
    timesmjd=[]
    for i in times:
        timesmjd.append(Time(i).mjd)
    return timesmjd


def convert_date(times):
    timesdate=[]
    for i in times:
        timesdate.append(Time(i,format='mjd').datetime)
    return timesdate

def convert_date_single(time):
    timedate=Time(time,format='mjd').datetime
    return timedate
    

import pandas as pd
import numpy as np
import os

def get_obsids(path):
    dirname=os.listdir(path)
    obsids=[]
    for i in dirname:
        if i.isdigit():
            obsids.append(i)
    obsids.sort()        
    return obsids

def drop_index(data):
    data=data.reset_index(drop=True)
    return data

def get_info(data,label,label_err=None):
    return min(data[label]),max(data[label]),np.mean(data[label])
    

def set_mag_ylim(ax):
    bottom, top = ax.set_ylim()
    if bottom< top:
        ax.set_ylim(top,bottom)
        
def set_mag_xlim(ax):
    bottom, top = ax.set_xlim()
    if bottom< top:
        ax.set_xlim(top,bottom)  

def set_ax_tick(ax):
    ax.xaxis.set_tick_params(which='major', size=10, width=2, direction='in', top='on',)
    ax.xaxis.set_tick_params(which='minor', size=5, width=2, direction='in', top='on')
    ax.yaxis.set_tick_params(which='major', size=10, width=2, direction='in', right='on')
    ax.yaxis.set_tick_params(which='minor', size=5, width=2, direction='in', right='on')

def set_ax_locator(ax,xma=1,xmi=0.2,yma=1,ymi=0.2):
    ax.xaxis.set_major_locator(mpl.ticker.MultipleLocator(xma))
    ax.xaxis.set_minor_locator(mpl.ticker.MultipleLocator(xmi))
    ax.yaxis.set_major_locator(mpl.ticker.MultipleLocator(yma))
    ax.yaxis.set_minor_locator(mpl.ticker.MultipleLocator(ymi))
    
    
def set_ax_legend(ax,bbox_to_anchor=(0.01, 0.99)):
    #ax.xaxis.set_tick_params(which='major', size=10, width=2, rotation=0,)
    handles, labels = ax.get_legend_handles_labels()
    # remove the errorbars
    #hdl = [h[0] for h in handles]
    hdl = handles
    labels_dict=dict(zip(labels, hdl)) #key,values
    by_label=OrderedDict(sorted(labels_dict.items(),key=lambda t:t[0]))
    #by_label = OrderedDict(zip(labels, handles))
    ax.legend(by_label.values(), by_label.keys(), bbox_to_anchor=bbox_to_anchor,
              loc=2, numpoints=1,ncol=1,fontsize=11.)



def plot_secax(ax,mi_interval=365,ma_interval=365*2,rotation=10,):
    secax1 = ax.secondary_xaxis('top', functions=(mjd2numsdate,datenums2mjd))
    secax1.xaxis.set_major_formatter(mdates.DateFormatter("%Y-%m-%d"))
    secax1.xaxis.set_major_locator(mdates.DayLocator(interval=ma_interval))
    secax1.xaxis.set_minor_locator(mdates.DayLocator(interval=mi_interval))
    #secax1.xaxis.set_tick_params(which='major', size=10, width=2, direction='out')
    secax1.xaxis.set_tick_params(which='minor', size=5, width=2, direction='out')
    secax1.xaxis.set_tick_params(which='major', size=10, width=2, direction='out', rotation=rotation,)
 
    
def set_ax_tick_nolabel(ax):
    ax.xaxis.set_tick_params(which='major', size=10, width=2, direction='in', top='off',labelsize=0)
    #ax.yaxis.set_tick_params(which='major', size=10, width=2, direction='in', right='off',labelsize=0)

       
    

from matplotlib.lines import Line2D
for m, func in Line2D.markers.items():
    print(m,func)
colors = cm.get_cmap('tab10', 10)

# Preliminaries: load in python modules.  The uses of these will become clear as we go forward
import numpy as np
from astropy.io import ascii
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import matplotlib.cm as cmx
from scipy.optimize import least_squares, curve_fit
from scipy.stats import f
import emcee
import corner
import os
from timeit import default_timer as timer


def calc_power_law(freq,S0,alpha):
    S = S0 * (freq) ** alpha
    return S

def alpha_calc(data):    
    #Get lightcurve values
    freqs = data['frequency']
    flux = data['flux']
    flux_errs = data['rms']
    
    #Use the scipy curve_fit algorithm to calculate the best fit value
    popt, pcov = curve_fit(calc_power_law, freqs, flux ,sigma=flux_errs, p0=(50,-0.61),absolute_sigma=True)
    
    alpha = popt[1] #Best-fit spectral index
    alpha_err = np.sqrt(np.diag(pcov))[1] #Uncertainty in alpha
    
    return alpha, alpha_err
 
def scale_data(data, alpha, alpha_err, ref_freq=5.0):
    #calculate a scaling factor for the flux density and uncertainty
    f_scale = (ref_freq/data['frequency'])**alpha
    rms_scale = np.abs(f_scale*np.log(ref_freq/data['frequency'])*alpha_err)
    
    #scale the flux and uncertainty - don't forget to add errors in quadrature
    scaled_flux = data['flux'] * f_scale
    scaled_rms = np.abs(scaled_flux) * np.sqrt((data['rms']/data['flux'])**2 + (rms_scale/f_scale)**2)
    
    #Add two new columns to the data
    data['scaled_flux'] = scaled_flux
    data['scaled_rms'] = scaled_rms
    
    return data 
    
#Select the data at the ~162 day epoch and print the spectral index + uncertainty
sel_data = data[data['delta_t'] == 162.89]
alpha, alpha_err = alpha_calc(sel_data)
print("alpha = %.1f+/-%.1f"%(alpha, alpha_err))    
data = scale_data(data, alpha, alpha_err)   

figure_n= 4
fig = plt.figure(figsize=(6,figure_n*3))
fig.subplots_adjust(hspace=0.0, wspace = 0.0)
ax = fig.add_subplot(figure_n,1,1)



dataplot=

for i in range(len(dataplot)):
    x=dataplot.iloc[i]['mjd']
    y=dataplot.iloc[i]['scaled_flux']
    xerr=dataplot.iloc[i]['mjderr']
    yerr=dataplot.iloc[i]['scaled_rms']
               
    ax.errorbar(x=x,y=y,
           xerr=xerr,
           yerr=yerr,
           marker='o',ms=11., mew=1, capsize=0,mec=color,ecolor=color,
           elinewidth=2,fmt='o',ls='',fillstyle='none',label='$F_\mathrm{5\,GHz}$')
            
            

ax_x=ax.twinx()
ax2 = fig.add_subplot(figure_n,1,2,sharex=ax)
ax3 = fig.add_subplot(figure_n,1,3,sharex=ax)
ax4 = fig.add_subplot(figure_n,1,4,sharex=ax)


#set_mag_ylim(ax)


ax.set_title('title')
ax.set_ylabel(r'$\mathrm{Mag}$')
ax_x.set_ylabel(r'$\mathrm{Mag}$')
ax.set_xlabel(r'$\mathrm{MJD(d)}$')


set_ax_tick(ax) 
set_ax_locator(ax,xma=1,xmi=0.2,yma=1,ymi=0.2)
set_ax_legend(ax,bbox_to_anchor=(0.01, 0.99))
ax.text(0.0, 0.1, "LogLocator(base=10, numticks=15)",fontsize=15, transform=ax.transAxes)


plot_secax(ax,mi_interval=365,ma_interval=365*2,rotation=10,)
ax.xaxis.set_tick_params(which='major', size=10, width=2, direction='in', top='off',labelsize=0)


ax.axhspan(1,2, facecolor='#2ca02c', alpha=0.5)    
ax.axvspan(56293,57033, facecolor='#2ca02c', alpha=0.05)

ax.grid(alpha=0.1,which='major', linestyle='--', linewidth=1)

range_l,range_r=(54930+1,54970-1)
ax.set_xlim(range_l,range_r)
ax.set_ylim(0.1,9.9)


#ax.semilogx() 
ax.semilogy() 
#ax.set_xscale('log')
ax.xaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=15))
ax.yaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=15))




save_lc_img_path='/Users/lyubing/path/lc_%d_%d.png'%(range_l,range_r)
plt.savefig(save_lc_img_path,dpi=400, transparent=False, bbox_inches='tight')