VanVleck.h

Go to the documentation of this file.

// # VanVleck.h: Class of static functions to aid with vanVleck corrections.
// # Copyright (C) 2002
// # Associated Universities, Inc. Washington DC, USA.
// #
// # This library is free software; you can redistribute it and/or modify it
// # under the terms of the GNU Library General Public License as published by
// # the Free Software Foundation; either version 2 of the License, or (at your
// # option) any later version.
// #
// # This library is distributed in the hope that it will be useful, but WITHOUT
// # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
// # FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
// # License for more details.
// #
// # You should have received a copy of the GNU Library General Public License
// # along with this library; if not, write to the Free Software Foundation,
// # Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
// #
// # Correspondence concerning AIPS++ should be addressed as follows:
// #        Internet email: casa-feedback@nrao.edu.
// #        Postal address: AIPS++ Project Office
// #                        National Radio Astronomy Observatory
// #                        520 Edgemont Road
// #                        Charlottesville, VA 22903-2475 USA
 
#ifndef SCIMATH_VANVLECK_H
#define SCIMATH_VANVLECK_H
 
// #! Includes go here
#include <casacore/casa/aips.h>
#include <casacore/casa/Arrays/Matrix.h>
#include <casacore/scimath/Functionals/Interpolate1D.h>
#include <casacore/casa/BasicSL/Constants.h>
 
#include <mutex>
 
namespace casacore {  // # NAMESPACE CASACORE - BEGIN
 
// # Forward Declarations
 
// <summary>
// A class of static functions to aid with vanVleck corrections of lag data.
// </summary>
 
// <use visibility=export>
 
// <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
// </reviewed>
 
// <prerequisite>
//   <li> Familiarity with the issues involved in turning digitally
//        sampled lag data from a correlator into spectral data.
// </prerequisite>
//
// <etymology>
// This provides the functions necessary to determine the van Vleck correction
// for a general n-level by m-level correlator.
// </etymology>
//
// <synopsis>
// This provides the functions necessary to determine the van Vleck correction
// for a general n-level by m-level correlator.
// </synopsis>
//
// <example>
// </example>
//
// <motivation>
// The GBT spectrometer provides the measured auto-correlation and
// cross-correlation lags.  The gbt MeasurementSet filler (gbtmsfiller)
// needs to convert those lags to the spectral domain.  These functions
// allow the filler to calculate the van Vleck correction appropriate
// for each measured set of lags.  They are of general and hence are
// not specific to the GBT spectrometer.
//
// The functions here are static because of the nature of the underlying
// numerical quadrature fortran code used to integrate the
// drbyrho function.
// </motivation>
//
// <thrown>
//    <li>
//    <li>
// </thrown>
//
// <todo asof="2002/07/19">
//   <li> The inverse error functions may be more generally useful.
//        It exists here only as a private member function to be
//        used internally.
// </todo>
 
class VanVleck {
 public:
  // Set the interpolation table size.
  // Should be an odd number.  The default size is 65.
  static void size(uInt npts);
 
  // get the current size.
  static uInt getsize();
 
  // Set the x and y quantization functions.
  // Each matrix should have dimensions (2,n)
  // where n is the number of levels.  The first
  // row (0,...) is the (n-1) threshold levels and
  // the second row is the n quantizations based
  // on those thresholds.  The thresholds may
  // include a DC offset. The (0,(n-1)) element is
  // never used and need not be set.
  static void setQuantization(const Matrix<Double> &qx, const Matrix<Double> &qy);
 
  // Set the x and y quantization levels for the case
  // of equi-spaced levels with a possible non-zero
  // offset.  The total number of levels is given by n,
  // which must be 3 or 9.  If n is not 3 or 9, False
  // will be returned and no quantization will have been
  // set.  For the 3- and 9- level cases a bivarate normal
  // integral calculation will be used.  That is much faster
  // than the more general numerical integration used
  // by setQuantization.
  static Bool setEquiSpaced(Double xlev, Double ylev, Double xmean, Double ymean, Int n);
 
  // Get the data used in setting up the interpolation
  static void getTable(Vector<Double> &rs, Vector<Double> &rhos);
 
  // Given a rho return the corresponding corrected r
  // Returns 0.0 if no quantization has been set yet.
  static Double r(const Double rho);
 
  // Given a measured zero-lag autocorrelation and number of
  // levels (n>=3) return the first positive quantizer input
  // threshold level.  This can be used to set the up the
  // matrix arguments used in setQuantization.
  static Double thresh(Int n, Double zerolag) {
    return ((n > 3) ? threshNgt3(n, zerolag) : threshN3(zerolag));
  }
 
  // Predict a given zero-lag given n and a threshold.  This
  // is included here to be used as a check against the output
  // of thresh.
  static Double predict(Int n, Double threshhold) {
    return ((n > 3) ? predictNgt3(n, threshhold) : predictN3(threshhold));
  }
 
  // Compute an approximation to the mean signal level (DC offset)
  // and quantizer threshold setting (both in terms of the r.m.s.
  // signal input level) given the observed positive bias (the
  // asymptotic limit of the measured autocorrelation at large
  // lags) and the zero-lag autocorrelation.
  // dcoffset is the mean signal level, threshold is the quantizer
  // setting, n is the number of levels, zerolag is the zero-lag
  // value and bias is the asymptotic bias.
  // Currently, this is only available for the n==3 level case,
  // all other cases set the returned dcoffset to 0 and use thresh()
  // to set the returned value of threshold.  A return value of F
  // indicates that the zerolag and bias values are inconsistent
  // and the dcoffset can not be determined.  In that case,
  // the returned dcoffset is 0 and thresh() is used to set
  // the threshold level.
  static Bool dcoff(Double &dcoffset, Double &threshold, Int n, Double zerolag, Double bias);
 
 private:
  // the number of points to use in setting up the interpolator
  static uInt itsSize, itsNx, itsNy;
 
  static Bool itsEquiSpaced;
 
  static Double itsXlev, itsYlev, itsXmean, itsYmean;
 
  // The interpolator
  static Interpolate1D<Double, Double> *itsInterp;
 
  // the quantization functions
  static Vector<Double> itsQx0, itsQx1, itsQy0, itsQy1;
 
  // Useful combinations of the above - to speed up drbydrho
  // these are -1/2*(Qx0*Qx0) and -1/2*(Qy0*Qy0)
  // These are only used for i up to (itsQx0.nelements() and
  // for j up to (itsQy0.nelements()).
  static Vector<Double> itsQx0Qx0, itsQy0Qy0;
  // This is Qx0[i]*Qy0[j]
  static Matrix<Double> itsQx0Qy0;
  // This is (Qx1[i+1]-Qx1[i])*(Qy1[j+1]*Qy1[j])
  static Matrix<Double> itsQx1Qy1diffs;
  // The mutex to make the functions thread-safe.
  static std::mutex theirMutex;
 
  // The fortran numerical integration function will call this.
  // For a given rho and quantization functions, this computes,
  // via Price's theorem, the value dr/drho of the derivative,
  // with respect to rho, of the expected value of the correlator
  // output.
  static Double drbydrho(Double *rho);
 
  // For a given rhoi, rhof, this produces a high-accuracy numerical
  // approximation to the integral of drbydrho over the range
  // rhoi to rhof.  It calls the standard QUADPACK adaptive Gaussian quadrature
  // procedure, dqags, to do the numerical integration.
  static Double rinc(Double &rhoi, Double &rhof);
 
  // initialize the interpolator
  static void initInterpolator();
 
  // compute first threshhold for a given zerolag for n>3
  static Double threshNgt3(Int n, Double zerolag);
 
  // compute first threshhold for a given zerolag for n==3
  static Double threshN3(Double zerolag) { return sqrt(2.0) * invErfc(zerolag); }
 
  // inverse err fn - used by invErfc
  static Double invErf(Double x);
 
  // inverse complementary err fn - used by threshN3
  static Double invErfc(Double x);
 
  // Predict a zero-lag value given the indicated first threshold level
  // for n>3.
  static Double predictNgt3(Int n, Double threshhold);
 
  // Predict a zero-lag value given the indicated first threshold level
  // for n=3.
  static Double predictN3(Double threshhold) { return ::erfc(threshhold / sqrt(2.0)); }
 
  // implementation of dcoff for the 3-level case
  static Bool dcoff3(Double &dcoffset, Double &threshold, Double zerolag, Double bias);
};
 
}  // namespace casacore
 
#endif