final

scripts/03_gti.py

@@ -0,0 +1,205 @@
+#!/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)
+
+#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)
+
+