bsproc/get_multicam_lc.py at main · NGTS/bsproc · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Fri May  7 14:55:44 2021

@author: ed
"""

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import sys
import argparse as ap

def ParseArgs():
    parser = ap.ArgumentParser()
    parser.add_argument('--name', type=str)
    parser.add_argument('--night', type=str)
    parser.add_argument('--actions', type=int, nargs='*')
    parser.add_argument('--apers', type=float, nargs='*', default=None)
    parser.add_argument('--output', type=str, default='./bsproc_outputs/')
    return parser.parse_args()

def lb(time, flux, err, bin_width):
    '''
    Function to bin the data into bins of a given width. time and bin_width
    must have the same units

    Input - time, flux, err, bin_width
    Return - time_bin, flux_bin, err_bin
    '''

    edges = np.arange(np.min(time), np.max(time), bin_width)
    dig = np.digitize(time, edges)
    time_binned = (edges[1:] + edges[:-1]) / 2
    flux_binned = np.array([np.nan if len(flux[dig == i]) == 0 else flux[dig == i].mean() for i in range(1, len(edges))])
    err_binned = np.array([np.nan if len(flux[dig == i]) == 0 else np.sqrt(np.sum(err[dig==i]**2))/len(err[dig==i]) for i in range(1, len(edges))])
    time_bin = time_binned[~np.isnan(err_binned)]
    err_bin = err_binned[~np.isnan(err_binned)]
    flux_bin = flux_binned[~np.isnan(err_binned)]

    return time_bin, flux_bin, err_bin

if __name__ == "__main__":
    args = ParseArgs()
    name = args.name
 #   root_dir = '/ngts/scratch/brightstars/PAOPhot2/'+name+'/'
    root_dir = args.output+'/'+name+'/'
    filedir = root_dir+'analyse_outputs/'+args.night+'/data_files/'
    opdir = root_dir+'analyse_outputs/'+args.night+'/'
    bjd, flux, err = np.array([]), np.array([]), np.array([])
    flux0, err0 = np.array([]), np.array([])
    actions = np.array([])
    airmass = np.array([])
    fwhm_sep = np.array([])
    fwhm_tl = np.array([])
    fwhm_rgw = np.array([])
    skybg = np.array([])
    if len(args.apers) == 1:
        apers = [args.apers[0] for ac in args.actions]
        file_name = opdir+name+f'_NGTS_'+args.night+f'_A{args.apers[0]}_bsproc_lc.dat'
    else:
        apers = args.apers
        file_name = opdir+name+f'_NGTS_'+args.night+f'_multiap_bsproc_lc.dat'
    for ac, rap in zip(args.actions, apers):
        df = pd.read_csv(filedir+f'action{ac}_bsproc_dat.csv',
                         index_col='NExposure')
        t = np.array(df.BJD)
        f = np.array(df.loc[:, f'FluxNormA{rap}'])
        e = np.array(df.loc[:, f'FluxNormErrA{rap}'])
        f0 = np.array(df.loc[:, f'FluxA{rap}'])
        e0 = np.array(df.loc[:, f'FluxErrA{rap}'])

        bjd = np.append(bjd, t)
        flux = np.append(flux, f)
        err = np.append(err, e)
        flux0 = np.append(flux0, f0)
        err0 = np.append(err0, e0)
        actions = np.append(actions, np.zeros_like(t)+ac)

        am = np.array(df.Airmass)
        fw_sep = np.array(df.FWHM_SEP)
        fw_tl = np.array(df.FWHM_TL)
        fw_rgw = np.array(df.FWHM_RGW)
        bg = np.array(df.loc[:, f'SkyBgA{rap}'])

        airmass = np.append(airmass, am)
        fwhm_sep = np.append(fwhm_sep, fw_sep)
        fwhm_tl = np.append(fwhm_tl, fw_tl)
        fwhm_rgw = np.append(fwhm_rgw, fw_rgw)
        skybg = np.append(skybg, bg)

    sig = np.std(flux)
    idbin = ((flux > 1-6.*sig) & (flux < 1+6.*sig))
    tb, fb, eb = lb(bjd[idbin], flux[idbin], err[idbin], 5/1440.)

    plt.plot(bjd, flux, '.k', alpha=0.2)
    plt.errorbar(tb, fb, yerr=eb, fmt='bo')
    plt.xlabel('Time (BJD)')
    plt.ylabel('Norm Flux')

    plt.savefig(opdir+name+f'_NGTS_'+args.night+f'_A{rap}_bsproc_lc.png')
    plt.show(block=False)
    save = input('Save over autosaved plot? [y/n] :  ')
    if save == 'y':
        plt.savefig(opdir+name+f'_NGTS_'+args.night+f'_A{rap}_bsproc_lc.png')
    plt.close()

    op = np.column_stack((actions, bjd, airmass, flux, err, flux0, err0,
                          skybg, fwhm_sep, fwhm_tl, fwhm_rgw))

    header =  'Object: '+args.name + \
              '\n Night: '+args.night + \
             f'\n Actions: {args.actions}' + \
             f'\n Aperture Radii: {apers} pixels' + \
              '\n ActionID   BJD   Airmass   FluxNorm   FluxNormErr   Flux   FluxErr  SkyBg   FWHM_SEP  FWHM_TL   FWHM_RGW'

    np.savetxt(file_name,
               op, header=header,
               fmt='%i %.8f %.6f %.8f %.8f %.8f %.8f %.3f %.3f %.3f %.3f', delimiter=' ')