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

src/libtoast/src/toast_atm.cpp

@@ -3,9 +3,9 @@
 // All rights reserved.  Use of this source code is governed by
 // a BSD-style license that can be found in the LICENSE file.
 
-// #if !defined(DEBUG)
-// #   define DEBUG
-// #endif
+#if !defined(DEBUG)
+# define DEBUG
+#endif // if !defined(DEBUG)
 
 #include <toast/sys_utils.hpp>
 #include <toast/sys_environment.hpp>
@@ -1021,10 +1021,13 @@ void toast::atm_sim::initialize_kolmogorov() {
     double invkappal = 1 / kappal;     // Optimize
     double kappa0 = 0.75 * kappamin;
     double kappa0sq = kappa0 * kappa0; // Optimize
-    long nkappa = 1000000;             // Number of integration steps needs to
-                                       // be large
-    double upper_limit = 10 * kappamax;
-    double kappastep = upper_limit / (nkappa - 1);
+    long nkappa = 10000;               // Number of integration steps
+    double kappastart = 1e-4;
+    double kappastop = 10 * kappamax;
+
+    // kappa = exp(ikappa * kappascale) * kappastart
+    double kappascale = log(kappastop / kappastart) / (nkappa - 1);
+
     double slope1 = 7. / 6.;
     double slope2 = -11. / 6.;
 
@@ -1048,26 +1051,31 @@ void toast::atm_sim::initialize_kolmogorov() {
     // Precalculate the power spectrum function
 
     std::vector <double> phi(last_kappa - first_kappa);
+    std::vector <double> kappa(last_kappa - first_kappa);
+    # pragma omp parallel for schedule(static, 10)
+    for (long ikappa = first_kappa; ikappa < last_kappa; ++ikappa) {
+        kappa[ikappa - first_kappa] = exp(ikappa * kappascale) * kappastart;
+    }
+
     # pragma omp parallel for schedule(static, 10)
     for (long ikappa = first_kappa; ikappa < last_kappa; ++ikappa) {
-        double kappa = ikappa * kappastep;
-        double kkl = kappa * invkappal;
+        double k = exp(ikappa * kappascale) * kappastart;
+        double kkl = k * invkappal;
         phi[ikappa - first_kappa] =
             (1. + 1.802 * kkl - 0.254 * pow(kkl, slope1))
-            * exp(-kkl * kkl) * pow(kappa * kappa + kappa0sq, slope2);
+            * exp(-kkl * kkl) * pow(k * k + kappa0sq, slope2);
     }
 
-    /*
-       if ( rank == 0 && verbosity > 0) {
-       std::ofstream f;
-       std::ostringstream fname;
-       fname << "kolmogorov_f.txt";
-       f.open( fname.str(), std::ios::out );
-       for ( int ikappa=0; ikappa<nkappa; ++ikappa )
-       f << ikappa*kappastep << " " << phi[ikappa] << std::endl;
-       f.close();
-       }
-     */
+    if ((rank == 0) && (verbosity > 0)) {
+        std::ofstream f;
+        std::ostringstream fname;
+        fname << "kolmogorov_f.txt";
+        f.open(fname.str(), std::ios::out);
+        for (int ikappa = 0; ikappa < nkappa; ++ikappa) {
+            f << kappa[ikappa] << " " << phi[ikappa] << std::endl;
+        }
+        f.close();
+    }
 
     // Newton's method factors, not part of the power spectrum
 
@@ -1086,25 +1094,33 @@ void toast::atm_sim::initialize_kolmogorov() {
     for (long ir = 0; ir < nr; ++ir) {
         double r = rmin_kolmo
                    + (exp(ir * nri * tau) - 1) * enorm * (rmax_kolmo - rmin_kolmo);
+        double rinv = 1 / r;
         double val = 0;
         if (r * kappamax < 1e-2) {
             // special limit r -> 0,
             // sin(kappa.r)/r -> kappa - kappa^3*r^2/3!
-            for (long ikappa = first_kappa; ikappa < last_kappa; ++ikappa) {
-                double kappa = ikappa * kappastep;
-                double kappa2 = kappa * kappa;
+            double r2 = r * r;
+            for (long ikappa = first_kappa; ikappa < last_kappa - 1; ++ikappa) {
+                double k = kappa[ikappa - first_kappa];
+                double kstep = kappa[ikappa + 1 - first_kappa] - k;
+                double kappa2 = k * k;
                 double kappa4 = kappa2 * kappa2;
-                double r2 = r * r;
-                val += phi[ikappa - first_kappa] * (kappa2 - r2 * kappa4 * ifac3);
+                val += phi[ikappa - first_kappa] * (kappa2 - r2 * kappa4 * ifac3) *
+                       kstep;
             }
         } else {
-            for (long ikappa = first_kappa; ikappa < last_kappa; ++ikappa) {
-                double kappa = ikappa * kappastep;
-                val += phi[ikappa - first_kappa] * sin(kappa * r) * kappa;
+            for (long ikappa = first_kappa; ikappa < last_kappa - 1; ++ikappa) {
+                double k1 = kappa[ikappa - first_kappa];
+                double k2 = kappa[ikappa + 1 - first_kappa];
+                double phi1 = phi[ikappa - first_kappa];
+                double phi2 = phi[ikappa + 1 - first_kappa];
+                val += .5 * (phi1 + phi2) * rinv *
+                       (k1 * cos(k1 * r) - k2 * cos(k2 * r)
+                        - rinv * (sin(k1 * r) - sin(k2 * r))
+                       );
             }
             val /= r;
         }
-        val *= kappastep;
         kolmo_x[ir] = r;
         kolmo_y[ir] = val;
     }

src/libtoast/src/toast_math_lapack.cpp

@@ -224,3 +224,24 @@ void toast::lapack_potri(char * UPLO, int * N, double * A, int * LDA,
     #endif // ifdef HAVE_LAPACK
     return;
 }
+
+#define dgelss LAPACK_FUNC(dgelss, DGELSS)
+
+extern "C" void 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);
+
+void toast::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) {
+    #ifdef HAVE_LAPACK
+    dgelss(M, N, NRHS, A, LDA, B, LDB, S, RCOND, RANK, WORK, LWORK, INFO);
+    #else // ifdef HAVE_LAPACK
+    auto here = TOAST_HERE();
+    auto log = toast::Logger::get();
+    std::string msg("TOAST was not compiled with BLAS/LAPACK support.");
+    log.error(msg.c_str(), here);
+    throw std::runtime_error(msg.c_str());
+    #endif // ifdef HAVE_LAPACK
+    return;
+}