#!/usr/bin/env python
"""
НАЗВАНИЕ:

 gti.py
 
 
НАЗНАЧЕНИЕ:

   Cкрипт для создания интервалов хорошего времени по данным из файла, сбинированного по 60 минут программой create_binned_filed.py (на основе оригинального файла lc_orig_60-120.csv, выгруженного с СРГ L2 монитора https://monitor.srg.cosmos.ru/ )
     Скрипт аппроксимирует кривую блеска полиномом, чтобы убрать синусоидальный тренд,  ищет времена вспышек и возвращает интервалы времени без вспышек в офрмате секунд с начала 1 января 2000 года (Москва) (файл gti_sec.csv) и в формате дестичного года (файл gti_dyear.csv).

ВЫЗОВ:

   conda activate
   ./gti.py


УПРАВЛЕНИЕ:

   Название файла надо вставить внутрь скрипта (ищите lc_bin_60min_60-120.csv)
 
ПАРАМЕТРЫ:

   N/A


ВЫВОД:
   Файл с графиком monitor.png   записывается в текущую директорию.
   Файл с интервалами времени gti.csv записывается в текущую директорию.
  

ИСТОРИЯ:
	   Филиппова Екатерина, ИКИ РАН, kate@cosmos.ru
   Сентябрь, 2025
   
"""

import matplotlib.pyplot as plt
import numpy as np
from matplotlib import ticker
import pandas as pd
import datetime
import dateutil
from os.path import dirname
import inspect

from scipy.optimize import curve_fit
from scipy.ndimage import gaussian_filter1d
from scipy import signal
from scipy.stats import median_abs_deviation
import astropy.io.fits as pf
from astropy.time import Time


import monitor
from monitor.config import *
from monitor.utils import *

""" find root dir """
root_path=dirname(dirname(dirname(inspect.getfile(monitor))))
print("Monitor root path: {}".format(root_path))


products_dir=root_path+'/products/'
data_dir=root_path+'/data/'

create_folder(products_dir)

fig, ax = plt.subplots(figsize=(9, 5), dpi=100)
ax.grid(visible=True,linestyle='dotted', )
ax.xaxis.set_minor_locator(ticker.MultipleLocator(1/12))
ax.yaxis.set_minor_locator(ticker.MultipleLocator(0.2))
ax.tick_params(axis="both", width=1, labelsize=14)

#geminga_dt = datetime.date.fromisoformat('2023-04-16')
#day_of_year = geminga_dt.timetuple().tm_yday  # returns 1 for January 1st
#geminga_tm=geminga_dt.year+(day_of_year/365)
#plt.axvline(x = geminga_tm, color = 'b', linewidth=3 , linestyle='dashed', label = 'Geminga scan')



#plt.xlim([2020,2025])
#plt.ylim([2,12])
plt.ylabel('Count rate (counts s$^{-1}$)',fontsize=14, fontweight='normal')
plt.xlabel('Year',fontsize=14, fontweight='normal')

for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(1)

cl='black'

input_file='default.csv'

# Сюда надо вставить название файла, скачанного с https://monitor.srg.cosmos.ru/
df = pd.read_csv(data_dir+input_file)

tm=[]
rate=[]
error=[]
for index, row in df.iterrows():
    dt = dateutil.parser.parse(row['timestamp'])  
    dyear = get_decimal_year(dt)
    tm.append(dyear)
    rate.append(float(row['value_60.0-120.0']))
    error.append(float(row['error_60.0-120.0']))

rate=np.array(rate)
tm=np.array(tm)
tmp=np.where(rate>0,1,0)
rate=np.compress(tmp,rate)
tm=np.compress(tmp,tm)
error=np.compress(tmp,error)

tmp=np.where(rate < 4.3,1,0)
tm_fit=np.compress(tmp,tm)
rate_fit=np.compress(tmp,rate)

# fit polinom of degree 
degree = 6
#initial_guess = [1] * (degree + 1)  # Initial guess for coefficients
#initial_guess = [-3.67265411e+08, 5.44675290e+05, -2.69261063e+02, 4.43698213e-02]
initial_guess = [ -6.10911582e+07,  3.02005656e+04, 1.49297834e+01, 6.53394516e-09, -3.64859359e-06, -1.80368300e-09, 8.91651441e-13] 

params, covariance = curve_fit(polynomial_func, tm_fit, rate_fit, p0=initial_guess)
print("curve_fit:",params)

# Generate fitted curve (original bin)
y_fit = polynomial_func(tm, *params)

# Delaem interpolyaciyu v 100 raz
year_min = min(tm)
year_max = max(tm)
year_step = (year_max - year_min)/len(tm)
year_sim = np.arange(year_min, year_max, year_step/100)
rate_interp = np.interp(year_sim, tm, rate)

# Generate fitted curve (interpolated)
y_fit_interp = polynomial_func(year_sim, *params)

plt.plot(tm, rate, color=cl, linewidth=1, linestyle='solid')
plt.plot(tm, y_fit, color="red", linewidth=2, linestyle='solid', label=f'Degree {degree} fit')
plt.ylabel('Count rate (counts s$^{-1}$)',fontsize=14, fontweight='normal')
plt.xlabel('Year',fontsize=14, fontweight='normal')
fig.savefig(products_dir+"bestfit.png", bbox_inches='tight')
plt.show()
fig.clear() # Clears the specific figure object

fig, ax = plt.subplots(figsize=(9, 5), dpi=100)
ax.grid(visible=True,linestyle='dotted', )
ax.xaxis.set_minor_locator(ticker.MultipleLocator(1/12))
ax.yaxis.set_minor_locator(ticker.MultipleLocator(0.2))
ax.tick_params(axis="both", width=1, labelsize=14)

#plt.xlim([2020,2025])
#plt.ylim([0,200])

for axis in ['top','bottom','left','right']:
    ax.spines[axis].set_linewidth(1)

    
resid=rate/y_fit
resid_interp=rate_interp/y_fit_interp
#smoothed_resid = gaussian_filter1d(resid, sigma=2)
outbursts = detect_outbursts(year_sim, resid_interp, threshold=threshold, min_duration=1)


plt.plot(tm, resid, color="black", linewidth=1, linestyle='solid',label="Residuals")
#plt.plot(tm, smoothed_resid, color="red", linewidth=1, linestyle='solid',label="Residuals")
plt.axhline(y = threshold, color = 'b', linewidth=1 , linestyle='dashed', label = 'Threshold')
plt.ylabel('Ratio',fontsize=14, fontweight='normal')
plt.xlabel('Year',fontsize=14, fontweight='normal')
fig.savefig(products_dir+"resid.png", bbox_inches='tight')


# Testing
"""
for start,end,peak,peak_flux in outbursts:
    print(start,end,peak,peak_flux)
    plt.axvline(x = start, color = 'r', linewidth=1 , linestyle='solid')
    plt.axvline(x = end, color = 'g', linewidth=1 , linestyle='solid')
    plt.axvline(x = peak, color = 'b', linewidth=1 , linestyle='solid')
"""
plt.show()


save_tex_table(outbursts,output=products_dir+'flares.tex')

tstart=[tm[0]]
tstop=[]
for start,end,peak,peak_flux in outbursts:
    if end != tm[-1]:
        tstop.append(start)
        tstart.append(end)
    else:
        tstop.append(start)

if outbursts[-1][1] != tm[-1]:
    tstop.append(tm[-1])

df = pd.DataFrame({'tstart': tstart,'tstop': tstop})
df.to_csv(products_dir+'gti_dyear.csv', index=False)

tstart_obj = Time(tstart, format='decimalyear')
tstop_obj = Time(tstop,format='decimalyear')
tstart_sec=(tstart_obj.mjd-51543.75)*86400
tstop_sec=(tstop_obj.mjd-51543.75)*86400

df = pd.DataFrame({'tstart': tstart_sec,'tstop': tstop_sec})
df.to_csv(products_dir+'gti_sec.csv', index=False)