generated from erosita/uds
213 lines
6.8 KiB
Python
Executable File
213 lines
6.8 KiB
Python
Executable File
#!/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)
|
||
|
||
|