Line data    Source code

       1             : //# SpectralEstimate.cc: Get an initial estimate for spectral lines
       2             : //# Copyright (C) 2001,2002,2003,2004
       3             : //# Associated Universities, Inc. Washington DC, USA.
       4             : //#
       5             : //# This library is free software; you can redistribute it and/or modify it
       6             : //# under the terms of the GNU Library General Public License as published by
       7             : //# the Free Software Foundation; either version 2 of the License, or (at your
       8             : //# option) any later version.
       9             : //#
      10             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      11             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      12             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Library General Public
      13             : //# License for more details.
      14             : //#
      15             : //# You should have received a copy of the GNU Library General Public License
      16             : //# along with this library; if not, write to the Free Software Foundation,
      17             : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
      18             : //#
      19             : //# Correspondence concerning AIPS++ should be addressed as follows:
      20             : //#        Internet email: aips2-request@nrao.edu.
      21             : //#        Postal address: AIPS++ Project Office
      22             : //#                        National Radio Astronomy Observatory
      23             : //#                        520 Edgemont Road
      24             : //#                        Charlottesville, VA 22903-2475 USA
      25             : //#
      26             : //# $Id: SpectralEstimate.cc 19933 2007-02-27 05:04:51Z Malte.Marquarding $
      27             : 
      28             : //# Includes
      29             : #include <components/SpectralComponents/SpectralEstimate.h>
      30             : #include <casacore/casa/Arrays/Vector.h>
      31             : #include <casacore/casa/BasicMath/Math.h>
      32             : 
      33             : using namespace casacore;
      34             : namespace casa { //# NAMESPACE CASA - BEGIN
      35             : 
      36             : //# Constructors
      37       16211 : SpectralEstimate::SpectralEstimate(const uInt maxpar) :
      38             :   useWindow_p(false), rms_p(0), cutoff_p(0),
      39             :   windowLow_p(0), windowEnd_p(0),
      40             :   regionLow_p(0), regionEnd_p(0),
      41             :   q_p(2), sigmin_p(0),
      42       16211 :   deriv_p(0), slist_p(maxpar), lprof_p(0) {
      43       16211 :   setQ();
      44       16211 : }
      45             : 
      46           0 : SpectralEstimate::SpectralEstimate(const Double rms,
      47             :                                    const Double cutoff, const Double minsigma,
      48           0 :                                    const uInt maxpar) :
      49             :   useWindow_p(false), rms_p(rms), cutoff_p(cutoff),
      50             :   windowLow_p(0), windowEnd_p(0),
      51             :   regionLow_p(0), regionEnd_p(0),
      52             :   q_p(2), sigmin_p(minsigma),
      53           0 :   deriv_p(0), slist_p(maxpar), lprof_p(0) {
      54           0 :   setQ();
      55           0 : }
      56             : 
      57             : 
      58           0 : SpectralEstimate::SpectralEstimate(const SpectralEstimate &other) :
      59           0 :   useWindow_p(other.useWindow_p), rms_p(other.rms_p), cutoff_p(other.cutoff_p),
      60           0 :   windowLow_p(other.windowLow_p), windowEnd_p(other.windowEnd_p),
      61           0 :   regionLow_p(other.regionLow_p), regionEnd_p(other.regionEnd_p),
      62           0 :   q_p(other.q_p), sigmin_p(other.sigmin_p),
      63           0 :   deriv_p(0), slist_p(other.slist_p), lprof_p(other.lprof_p) {
      64           0 :   setQ(q_p);
      65           0 :   deriv_p = new Double[lprof_p];
      66           0 :   for (uInt i=0; i<lprof_p; i++) deriv_p[i] = other.deriv_p[i];
      67           0 : }
      68             : 
      69       48633 : SpectralEstimate::~SpectralEstimate() {
      70       16211 :   delete [] deriv_p; deriv_p = 0; lprof_p = 0;
      71       16211 : }
      72             : 
      73           0 : SpectralEstimate &SpectralEstimate::operator=(const SpectralEstimate &other) {
      74           0 :   if (this != &other) {
      75           0 :     useWindow_p = other.useWindow_p; 
      76           0 :     rms_p = other.rms_p; 
      77           0 :     cutoff_p = other.cutoff_p;
      78           0 :     windowLow_p = other.windowLow_p;
      79           0 :     windowEnd_p = other.windowEnd_p;
      80           0 :     regionLow_p = other.regionLow_p;
      81           0 :     regionEnd_p = other.regionEnd_p;
      82           0 :     q_p = other.q_p;
      83           0 :     sigmin_p = other.sigmin_p;
      84           0 :     deriv_p = 0;
      85           0 :     slist_p = other.slist_p;
      86           0 :     lprof_p = other.lprof_p;
      87           0 :     setQ(q_p);
      88           0 :     deriv_p = new Double[lprof_p];
      89           0 :     for (uInt i=0; i<lprof_p; i++) deriv_p[i] = other.deriv_p[i];
      90             :   };
      91           0 :   return *this;
      92             : }
      93             : 
      94           0 : void SpectralEstimate::setRMS(const Double rms) {
      95           0 :   rms_p = abs(rms);
      96           0 : }
      97             : 
      98           0 : void SpectralEstimate::setCutoff(const Double cutoff) {
      99           0 :   cutoff_p = max(0.0, cutoff);
     100           0 : }
     101             : 
     102           0 : void SpectralEstimate::setMinSigma(const Double minsigma) {
     103           0 :   sigmin_p = max(0.0, minsigma);
     104           0 : }
     105             : 
     106       32422 : void SpectralEstimate::setQ(const uInt q) {
     107       32422 :   q_p = max(1, Int(q));
     108       32422 :   a_p = 90.0/(q_p*(q_p+1)*(4*q_p*q_p-1)*(2*q_p+3));
     109       32422 :   b_p = (q_p*(q_p+1))/3.0;
     110       32422 : }
     111             : 
     112           0 : void SpectralEstimate::setRegion(const Int lo, const Int hi) {
     113           0 :   regionLow_p = regionEnd_p = 0;
     114           0 :   if (hi > lo) {
     115           0 :     regionLow_p = lo;
     116           0 :     regionEnd_p = hi;
     117             :   };
     118           0 : }
     119             : 
     120           0 : void SpectralEstimate::setWindowing(const Bool win) {
     121           0 :   useWindow_p = win;
     122           0 : }
     123             : 
     124           0 : void SpectralEstimate::setMaxN(const uInt maxpar) {
     125           0 :   slist_p.set(maxpar);
     126           0 : }
     127             : 
     128             : } //# NAMESPACE CASA - END
     129             : 
     130             : 
     131             : //# Cater for Double and Float
     132             : #ifdef AIPS_NO_TEMPLATE_SRC
     133             : #include <components/SpectralComponents/Spectral2Estimate.tcc>
     134             : 
     135             : using namespace casacore;
     136             : namespace casa { //# NAMESPACE CASA - BEGIN
     137             : template SpectralList const & SpectralEstimate::estimate<Float>(Vector<Float> const &, Vector<Float> *);
     138             : template SpectralList const & SpectralEstimate::estimate<Float>(Vector<Float> const &, Vector<Float> const &);
     139             : template SpectralElement SpectralEstimate::convertElement<Float>(Vector<Float> const &,
     140             :                                                                  SpectralElement const &) const;
     141             : template void SpectralEstimate::findga<Float>(Vector<Float> const &); 
     142             : template uInt SpectralEstimate::window<Float>(Vector<Float> const &);
     143             : template void SpectralEstimate::findc2<Float>(Vector<Float> const &); 
     144             : 
     145             : template SpectralList const & SpectralEstimate::estimate<Double>(Vector<Double> const &, Vector<Double> *);
     146             : template SpectralList const & SpectralEstimate::estimate<Double>(Vector<Double> const &, Vector<Double> const &);
     147             : template SpectralElement SpectralEstimate::convertElement<Double>(Vector<Double> const &,
     148             :                                                                  SpectralElement const &) const;
     149             : template void SpectralEstimate::findga<Double>(Vector<Double> const &);
     150             : template uInt SpectralEstimate::window<Double>(Vector<Double> const &);
     151             : template void SpectralEstimate::findc2<Double>(Vector<Double> const &);
     152             : } //# NAMESPACE CASA - END
     153             : #endif