casa  $Rev:20696$
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
VanVleck.h
Go to the documentation of this file.
00001 //# VanVleck.h: Class of static functions to aid with vanVleck corrections.
00002 //# Copyright (C) 2002
00003 //# Associated Universities, Inc. Washington DC, USA.
00004 //#
00005 //# This library is free software; you can redistribute it and/or modify it
00006 //# under the terms of the GNU Library General Public License as published by
00007 //# the Free Software Foundation; either version 2 of the License, or (at your
00008 //# option) any later version.
00009 //#
00010 //# This library is distributed in the hope that it will be useful, but WITHOUT
00011 //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00012 //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
00013 //# License for more details.
00014 //#
00015 //# You should have received a copy of the GNU Library General Public License
00016 //# along with this library; if not, write to the Free Software Foundation,
00017 //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
00018 //#
00019 //# Correspondence concerning AIPS++ should be addressed as follows:
00020 //#        Internet email: aips2-request@nrao.edu.
00021 //#        Postal address: AIPS++ Project Office
00022 //#                        National Radio Astronomy Observatory
00023 //#                        520 Edgemont Road
00024 //#                        Charlottesville, VA 22903-2475 USA
00025 //#
00026 //# $Id: VanVleck.h 21051 2011-04-20 11:46:29Z gervandiepen $
00027 
00028 #ifndef SCIMATH_VANVLECK_H
00029 #define SCIMATH_VANVLECK_H
00030 
00031 //#! Includes go here
00032 #include <casa/aips.h>
00033 #include <casa/Arrays/Matrix.h>
00034 #include <scimath/Functionals/Interpolate1D.h>
00035 #include <casa/BasicSL/Constants.h>
00036 #include <casa/OS/Mutex.h>
00037 
00038 
00039 namespace casa { //# NAMESPACE CASA - BEGIN
00040 
00041 //# Forward Declarations
00042 
00043 // <summary>
00044 // A class of static functions to aid with vanVleck corrections of lag data.
00045 // </summary>
00046 
00047 // <use visibility=export>
00048 
00049 // <reviewed reviewer="" date="yyyy/mm/dd" tests="" demos="">
00050 // </reviewed>
00051 
00052 // <prerequisite>
00053 //   <li> Familiarity with the issues involved in turning digitally 
00054 //        sampled lag data from a correlator into spectral data.
00055 // </prerequisite>
00056 //
00057 // <etymology>
00058 // This provides the functions necessary to determine the van Vleck correction
00059 // for a general n-level by m-level correlator.
00060 // </etymology>
00061 //
00062 // <synopsis>
00063 // This provides the functions necessary to determine the van Vleck correction
00064 // for a general n-level by m-level correlator.
00065 // </synopsis>
00066 //
00067 // <example>
00068 // </example>
00069 //
00070 // <motivation>
00071 // The GBT spectrometer provides the measured auto-correlation and
00072 // cross-correlation lags.  The gbt MeasurementSet filler (gbtmsfiller)
00073 // needs to convert those lags to the spectral domain.  These functions
00074 // allow the filler to calculate the van Vleck correction appropriate
00075 // for each measured set of lags.  They are of general and hence are
00076 // not specific to the GBT spectrometer.
00077 //
00078 // The functions here are static because of the nature of the underlying
00079 // numerical quadrature fortran code used to integrate the 
00080 // drbyrho function.
00081 // </motivation>
00082 //
00083 // <thrown>
00084 //    <li>
00085 //    <li>
00086 // </thrown>
00087 //
00088 // <todo asof="2002/07/19">
00089 //   <li> The inverse error functions may be more generally useful.
00090 //        It exists here only as a private member function to be
00091 //        used internally.
00092 // </todo>
00093 
00094 class VanVleck
00095 {
00096 public:
00097     // Set the interpolation table size.
00098     // Should be an odd number.  The default size is 65.
00099     static void size(uInt npts);
00100 
00101     // get the current size.
00102     static uInt getsize();
00103     
00104     // Set the x and y quantization functions.
00105     // Each matrix should have dimensions (2,n)
00106     // where n is the number of levels.  The first
00107     // row (0,...) is the (n-1) threshold levels and 
00108     // the second row is the n quantizations based 
00109     // on those thresholds.  The thresholds may 
00110     // include a DC offset. The (0,(n-1)) element is
00111     // never used and need not be set.
00112     static void setQuantization(const Matrix<Double> &qx,
00113                                 const Matrix<Double> &qy);
00114 
00115     // Set the x and y quantization levels for the case
00116     // of equi-spaced levels with a possible non-zero
00117     // offset.  The total number of levels is given by n,
00118     // which must be 3 or 9.  If n is not 3 or 9, False
00119     // will be returned and no quantization will have been
00120     // set.  For the 3- and 9- level cases a bivarate normal
00121     // integral calculation will be used.  That is much faster
00122     // than the more general numerical integration used 
00123     // by setQuantization.
00124     static Bool setEquiSpaced(Double xlev, Double ylev,
00125                               Double xmean, Double ymean,
00126                               Int n);
00127     
00128     // Get the data used in setting up the interpolation
00129     static void getTable(Vector<Double> &rs, Vector<Double> &rhos);
00130     
00131     // Given a rho return the corresponding corrected r
00132     // Returns 0.0 if no quantization has been set yet.
00133     static Double r(const Double rho);
00134     
00135     // Given a measured zero-lag autocorrelation and number of
00136     // levels (n>=3) return the first positive quantizer input
00137     // threshold level.  This can be used to set the up the
00138     // matrix arguments used in setQuantization.
00139     static Double thresh(Int n, Double zerolag)
00140     { return ( (n>3) ? threshNgt3(n,zerolag) : threshN3(zerolag) ); }
00141     
00142     // Predict a given zero-lag given n and a threshold.  This
00143     // is included here to be used as a check against the output
00144     // of thresh.
00145     static Double predict(Int n, Double threshhold)
00146     { return ( (n>3) ? predictNgt3(n,threshhold) : predictN3(threshhold));}
00147 
00148     // Compute an approximation to the mean signal level (DC offset)
00149     // and quantizer threshold setting (both in terms of the r.m.s.
00150     // signal input level) given the observed positive bias (the
00151     // asymptotic limit of the measured autocorrelation at large
00152     // lags) and the zero-lag autocorrelation.
00153     // dcoffset is the mean signal level, threshold is the quantizer
00154     // setting, n is the number of levels, zerolag is the zero-lag
00155     // value and bias is the asymptotic bias.
00156     // Currently, this is only available for the n==3 level case,
00157     // all other cases set the returned dcoffset to 0 and use thresh()
00158     // to set the returned value of threshold.  A return value of F
00159     // indicates that the zerolag and bias values are inconsistent
00160     // and the dcoffset can not be determined.  In that case,
00161     // the returned dcoffset is 0 and thresh() is used to set
00162     // the threshold level.
00163     static Bool dcoff(Double &dcoffset, Double &threshold,
00164                       Int n, Double zerolag, Double bias);
00165                       
00166     
00167 private:
00168     // the number of points to use in setting up the interpolator
00169     static uInt itsSize, itsNx, itsNy;
00170 
00171     static Bool itsEquiSpaced;
00172 
00173     static Double itsXlev, itsYlev, itsXmean, itsYmean;
00174     
00175     // The interpolator
00176     static Interpolate1D<Double, Double> *itsInterp;
00177     
00178     // the quantization functions
00179     static Vector<Double> itsQx0, itsQx1, itsQy0, itsQy1;
00180     
00181     // Useful combinations of the above - to speed up drbydrho
00182     // these are -1/2*(Qx0*Qx0) and -1/2*(Qy0*Qy0)
00183     // These are only used for i up to (itsQx0.nelements() and
00184     // for j up to (itsQy0.nelements()).
00185     static Vector<Double> itsQx0Qx0, itsQy0Qy0;
00186     // This is Qx0[i]*Qy0[j]
00187     static Matrix<Double> itsQx0Qy0;
00188     // This is (Qx1[i+1]-Qx1[i])*(Qy1[j+1]*Qy1[j])
00189     static Matrix<Double> itsQx1Qy1diffs;
00190     // The mutex to make the functions thread-safe.
00191     static Mutex theirMutex;
00192 
00193 
00194     // The fortran numerical integration function will call this.
00195     // For a given rho and quantization functions, this computes,
00196     // via Price's theorem, the value dr/drho of the derivative,
00197     // with respect to rho, of the expected value of the correlator
00198     // output.
00199     static Double drbydrho(Double *rho);
00200     
00201     // For a given rhoi, rhof, this produces a high-accuracy numerical
00202     // approximation to the integral of drbydrho over the range
00203     // rhoi to rhof.  It calls the standard QUADPACK adaptive Gaussian quadrature
00204     // procedure, dqags, to do the numerical integration.
00205     static Double rinc(Double &rhoi, Double &rhof);
00206     
00207     // initialize the interpolator
00208     static void initInterpolator();
00209     
00210     // compute first threshhold for a given zerolag for n>3
00211     static Double threshNgt3(Int n, Double zerolag);
00212     
00213     // compute first threshhold for a given zerolag for n==3
00214     static Double threshN3(Double zerolag)
00215     { return sqrt(2.0)*invErfc(zerolag);}
00216     
00217     // inverse err fn - used by invErfc
00218     static Double invErf(Double x);
00219     
00220     // inverse complementary err fn - used by threshN3
00221     static Double invErfc(Double x);
00222     
00223     // Predict a zero-lag value given the indicated first threshold level
00224     // for n>3.
00225     static Double predictNgt3(Int n, Double threshhold);
00226     
00227     // Predict a zero-lag value given the indicated first threshold level
00228     // for n=3.
00229     static Double predictN3(Double threshhold)
00230     { return ::erfc(threshhold/sqrt(2.0));}
00231 
00232     // implementation of dcoff for the 3-level case
00233     static Bool dcoff3(Double &dcoffset, Double &threshold,
00234                        Double zerolag, Double bias);
00235 };
00236 
00237 
00238 } //# NAMESPACE CASA - END
00239 
00240 #endif