High resolution and fix recent

nuclockutils/barycorr.py

@@ -165,7 +165,6 @@ def apply_clock_correction(
                     else:
                         hdu.header[keyname] = corrected_time
 
-
             hdu.header['CREATOR'] = f'NuSTAR Clock Utils - v. {version}'
             hdu.header['DATE'] = Time.now().fits
             hdu.header['PLEPHEM'] = f'JPL-{ephem}'

nuclockutils/clean_clock/app.py

@@ -59,8 +59,8 @@ def recalc(outfile='save_all.pickle'):
         load_and_flag_clock_table(clockfile=CLOCKFILE, shift_non_malindi=True)
 
     table_times = temptable_raw['met']
-    met_start = table_times[0]
-    met_stop = table_times[-1]
+    met_start = clock_offset_table['met'][0]
+    met_stop = clock_offset_table['met'][-1] + 30
     clock_jump_times = np.array([78708320, 79657575, 81043985, 82055671,
                                  293346772])
     clock_jump_times += 30 #  Sum 30 seconds to avoid to exclude these points
@@ -74,14 +74,16 @@ def recalc(outfile='save_all.pickle'):
         time_resolution=10, craig_fit=False, hdf_dump_file='dump.hdf5')
 
     table_new = eliminate_trends_in_residuals(
-        table_new, clock_offset_table_corr, gtis)
+        table_new, clock_offset_table_corr, gtis,
+        fixed_control_points=np.arange(291e6, 295e6, 86400))
 
     mets = np.array(table_new['met'])
     start = mets[0]
     stop = mets[-1]
 
     good_mets = clock_offset_table['met'] < stop
     clock_offset_table = clock_offset_table[good_mets]
+    # print(clock_offset_table['met'][-10:], table_new['met'][-1])
 
     clock_mets = clock_offset_table['met']
     clock_mjds = clock_offset_table['mjd']

nuclockutils/diagnostics/bary_and_fold_all.py

@@ -9,6 +9,7 @@
 import astropy
 from astropy import log
 from astropy.io import fits
+import pint
 
 from .get_crab_ephemeris import get_crab_ephemeris
 from .fold_to_ephemeris import get_events_from_fits, \
@@ -96,6 +97,13 @@ def main(args=None):
     parser.add_argument('--plot-phaseogram', default=False,
                         action='store_true',
                         help='Plot the phaseogram (requires PINT)')
+    parser.add_argument("--expocorr", default=False,
+                        action='store_true',
+                        help="Calculate the exposure from (NuSTAR only, the"
+                             "event file must be unfiltered)")
+    parser.add_argument("--use-standard-barycorr", default=False,
+                        action='store_true',
+                        help="Use standard barycorr instead of nuclockutils")
 
     args = parser.parse_args(args)
 
@@ -135,10 +143,20 @@ def main(args=None):
             mkdir(outdir)
 
             if not os.path.exists(bary_file):
-                cmd = f'{fname} {attorb_file} -p {parfile} ' \
-                      f'-c {args.clockfile} -o {bary_file}'
-
-                nubarycorr(cmd.split())
+                if args.use_standard_barycorr:
+                    model = pint.models.get_model(parfile)
+                    ra = model.RAJ.to("degree").value
+                    dec = model.DECJ.to("degree").value
+                    cmd = f'barycorr {fname} {bary_file} ' \
+                          f'orbitfiles={attorb_file} ' \
+                          f'clockfile={args.clockfile} ' \
+                          f'refframe=ICRS ra={ra} dec={dec}'
+                    sp.check_call(cmd.split())
+                else:
+                    cmd = f'{fname} {attorb_file} -p {parfile} ' \
+                          f'-c {args.clockfile} -o {bary_file}'
+
+                    nubarycorr(cmd.split())
 
             if args.plot_phaseogram:
                 cmd = f'photonphase {bary_file} {parfile} ' \

nuclockutils/diagnostics/compare_pulses.py

@@ -78,13 +78,14 @@ def main(args=None):
     mjds = []
     plt.figure(figsize=(6, 9))
     ref_profile = None
+    phase_time_plot_template, prof_plot_template = None, None
     if args.template is not None and os.path.exists(args.template):
         mjd, phase, prof, _, _, period_ms = \
             format_profile_and_get_phase(args.template, template=None)
-        phase_time_plot, prof_plot = format_for_plotting(phase, prof, period_ms)
+        phase_time_plot_template, prof_plot_template = format_for_plotting(phase, prof, period_ms)
         local_max = phase[np.argmax(prof)] * period_ms
 
-        plt.plot(phase_time_plot, prof_plot,
+        plt.plot(phase_time_plot_template, prof_plot_template,
             drawstyle='steps-mid', label=args.template, color='k')
         ref_profile = prof
         ref_max = local_max
@@ -100,13 +101,18 @@ def main(args=None):
 
         local_max = phase[np.argmax(prof)] * period_ms
 
-        phases.append(phase_res * period_ms)
+        local_phase_time = phase_res * period_ms
+        phases.append(local_phase_time)
         phase_errs.append(phase_res_err * period_ms)
         mjds.append(mjd)
         maxs.append(local_max)
 
         plt.plot(phase_time_plot, (i + 1) * 0.2 + prof_plot,
                  drawstyle='steps-mid', label=f, alpha=0.5, color='grey')
+        if phase_time_plot_template is not None:
+            plt.plot(phase_time_plot_template + local_phase_time,
+                     (i + 1) * 0.2 + prof_plot_template,
+                     drawstyle='steps-mid', label=f, alpha=0.5, color='grey')
         if len(files) < 2:
             continue
         for plot_shift in [0, period_ms, 2 * period_ms]: