Simulation tweaks

.travis.yml

@@ -196,7 +196,7 @@ before_install:
     - pip install --upgrade pip
     - pip install --upgrade setuptools
     - pip install -q numpy scipy matplotlib cmake cython astropy ephem healpy numba toml
-    - pip install -q https://github.com/healpy/pysm/archive/master.tar.gz
+    - pip install -q pysm3
     # Install all compiled dependencies from a pre-cached location
     - echo "Fetching dependencies from ${DEPSURL}"
     - wget ${DEPSURL}

pipelines/toast_ground_schedule.py

@@ -653,6 +653,7 @@ def attempt_scan(
     fout,
     fout_fmt,
     ods,
+    boresight_angle,
 ):
     """ Attempt scanning the visible patches in order until success.
     """
@@ -661,7 +662,16 @@ def attempt_scan(
         if isinstance(patch, CoolerCyclePatch):
             # Cycle the cooler
             t = add_cooler_cycle(
-                args, t, stop_timestamp, observer, sun, moon, fout, fout_fmt, patch
+                args,
+                t,
+                stop_timestamp,
+                observer,
+                sun,
+                moon,
+                fout,
+                fout_fmt,
+                patch,
+                boresight_angle,
             )
             success = True
             break
@@ -705,6 +715,7 @@ def attempt_scan(
                         patch,
                         el,
                         ods,
+                        boresight_angle,
                     )
                     patch.step_azel()
                     break
@@ -789,6 +800,7 @@ def attempt_scan_pole(
     fout,
     fout_fmt,
     ods,
+    boresight_angle,
 ):
     """ Attempt scanning the visible patches in order until success.
     """
@@ -800,7 +812,16 @@ def attempt_scan_pole(
         if isinstance(patch, CoolerCyclePatch):
             # Cycle the cooler
             t = add_cooler_cycle(
-                args, tstart, stop_timestamp, observer, sun, moon, fout, fout_fmt, patch
+                args,
+                tstart,
+                stop_timestamp,
+                observer,
+                sun,
+                moon,
+                fout,
+                fout_fmt,
+                patch,
+                boresight_angle,
             )
             success = True
             break
@@ -848,6 +869,7 @@ def attempt_scan_pole(
                         patch,
                         el,
                         ods,
+                        boresight_angle,
                         subscan=subscan,
                     )
                     el += np.radians(args.pole_el_step)
@@ -1258,6 +1280,7 @@ def add_scan(
     patch,
     el,
     ods,
+    boresight_angle,
     subscan=-1,
 ):
     """ Make an entry for a CES in the schedule file.
@@ -1374,6 +1397,7 @@ def add_scan(
             to_UTC(t2),
             to_MJD(t1),
             to_MJD(t2),
+            boresight_angle,
             patch.name,
             (azmin + args.boresight_offset_az_deg) % 360,
             (azmax + args.boresight_offset_az_deg) % 360,
@@ -1418,7 +1442,9 @@ def add_scan(
 
 
 @function_timer
-def add_cooler_cycle(args, tstart, tstop, observer, sun, moon, fout, fout_fmt, patch):
+def add_cooler_cycle(
+    args, tstart, tstop, observer, sun, moon, fout, fout_fmt, patch, boresight_angle
+):
     """ Make an entry for a cooler cycle in the schedule file.
     """
     log = Logger.get()
@@ -1447,6 +1473,7 @@ def add_cooler_cycle(args, tstart, tstop, observer, sun, moon, fout, fout_fmt, p
         to_UTC(t2),
         to_MJD(t1),
         to_MJD(t2),
+        boresight_angle,
         patch.name,
         az,
         az,
@@ -1539,6 +1566,17 @@ def get_visible(args, observer, sun, moon, patches, el_min):
     return visible, not_visible
 
 
+@function_timer
+def get_boresight_angle(args, t, t0=0):
+    """ Return the scheduled boresight angle at time t.
+    """
+    if args.boresight_angle_step == 0 or args.boresight_angle_time == 0:
+        return 0
+
+    istep = int((t - t0) / 60 / args.boresight_angle_time)
+    return (args.boresight_angle_step * istep) % 360
+
+
 @function_timer
 def apply_blockouts(args, t_in):
     """ Check if `t` is inside a blockout period.
@@ -1640,15 +1678,15 @@ def build_schedule(args, start_timestamp, stop_timestamp, patches, observer, sun
     )
 
     fout_fmt0 = (
-        "#{:20} {:20} {:14} {:14} "
-        "{:35} {:8} {:8} {:8} {:5} "
-        "{:8} {:8} {:8} {:8} "
-        "{:8} {:8} {:8} {:8} {:5} "
-        "{:5} {:3}\n"
+        "#{:>20} {:>20} {:>14} {:>14} {:>8} "
+        "{:35} {:>8} {:>8} {:>8} {:>5} "
+        "{:>8} {:>8} {:>8} {:>8} "
+        "{:>8} {:>8} {:>8} {:>8} {:>5} "
+        "{:>5} {:>3}\n"
     )
 
     fout_fmt = (
-        " {:20} {:20} {:14.6f} {:14.6f} "
+        " {:20} {:20} {:14.6f} {:14.6f} {:8.2f} "
         "{:35} {:8.2f} {:8.2f} {:8.2f} {:5} "
         "{:8.2f} {:8.2f} {:8.2f} {:8.2f} "
         "{:8.2f} {:8.2f} {:8.2f} {:8.2f} {:5.2f} "
@@ -1661,6 +1699,7 @@ def build_schedule(args, start_timestamp, stop_timestamp, patches, observer, sun
             "Stop time UTC",
             "Start MJD",
             "Stop MJD",
+            "Rotation",
             "Patch name",
             "Az min",
             "Az max",
@@ -1689,6 +1728,7 @@ def build_schedule(args, start_timestamp, stop_timestamp, patches, observer, sun
     last_successful = t
     while True:
         t, blocked = apply_blockouts(args, t)
+        boresight_angle = get_boresight_angle(args, t)
         if t > stop_timestamp:
             break
         if t - last_successful > 86400 or blocked:
@@ -1767,6 +1807,7 @@ def build_schedule(args, start_timestamp, stop_timestamp, patches, observer, sun
                 fout,
                 fout_fmt,
                 ods,
+                boresight_angle,
             )
         else:
             success, t = attempt_scan(
@@ -1785,6 +1826,7 @@ def build_schedule(args, start_timestamp, stop_timestamp, patches, observer, sun
                 fout,
                 fout_fmt,
                 ods,
+                boresight_angle,
             )
 
         if args.operational_days and len(ods) > args.operational_days:
@@ -1977,6 +2019,20 @@ def parse_args():
         type=np.float,
         help="Maximum allowed sun elevation [deg]",
     )
+    parser.add_argument(
+        "--boresight-angle-step",
+        required=False,
+        default=0,
+        type=np.float,
+        help="Boresight rotation step size [deg]",
+    )
+    parser.add_argument(
+        "--boresight-angle-time",
+        required=False,
+        default=0,
+        type=np.float,
+        help="Boresight rotation step interval [minutes]",
+    )
     parser.add_argument(
         "--start",
         required=False,

pipelines/toast_ground_sim.py

@@ -225,6 +225,7 @@ def create_observation(args, comm, telescope, ces, verbose=True):
             focalplane.detquats,
             totsamples,
             detranks=ndetrank,
+            boresight_angle=ces.boresight_angle,
             firsttime=ces.start_time,
             rate=args.sample_rate,
             site_lon=site.lon,
@@ -235,7 +236,7 @@ def create_observation(args, comm, telescope, ces, verbose=True):
             el=ces.el,
             scanrate=args.scan_rate,
             scan_accel=args.scan_accel,
-            sinc_modulation=args.scan_sinc_modulate,
+            cosecant_modulation=args.scan_cosecant_modulate,
             CES_start=None,
             CES_stop=None,
             sun_angle_min=args.sun_angle_min,

pipelines/toast_ground_sim_simple.py

@@ -280,7 +280,7 @@ def create_observations(args, comm, schedule):
                 el=ces.el,
                 scanrate=args.scan_rate,
                 scan_accel=args.scan_accel,
-                sinc_modulation=args.scan_sinc_modulate,
+                cosecant_modulation=args.scan_cosecant_modulate,
                 CES_start=None,
                 CES_stop=None,
                 sun_angle_min=args.sun_angle_min,

src/libtoast/include/toast/math_lapack.hpp

@@ -39,6 +39,10 @@ void lapack_pocon(char * UPLO, int * N, double * A, int * LDA, double * ANORM,
                   double * RCOND, double * WORK, int * IWORK, int * INFO);
 
 void lapack_potri(char * UPLO, int * N, double * A, int * LDA, int * INFO);
+
+void lapack_dgelss(int * M, int * N, int * NRHS, double * A, int * LDA,
+                   double * B, int * LDB, double * S, double * RCOND,
+                   int * RANK, double * WORK, int * LWORK, int * INFO);
 }
 
 #endif // ifndef TOAST_LAPACK_HPP

src/libtoast/include/toast/tod_filter.hpp

@@ -10,6 +10,12 @@ namespace toast {
 void filter_polynomial(int64_t order, size_t nsignal, uint8_t * flags,
                        std::vector <double *> const & signals, size_t nscan,
                        int64_t const * starts, int64_t const * stops);
+void bin_templates(double * signal, double * templates, uint8_t * good,
+                   double * invcov, double * proj, size_t nsample, size_t ntemplate);
+void chebyshev(double * x, double * templates, size_t start_order, size_t stop_order,
+               size_t nsample);
+void add_templates(double * signal, double * templates, double * coeff, size_t nsample,
+                   size_t ntemplate);
 }
 
 #endif // ifndef TOAST_TOD_FILTER_HPP