Line data    Source code

       1             : //# FluxCalcLogFreqPolynomial.cc: Implementation of FluxCalcLogFreqPolynomial.h
       2             : //# Copyright (C) 2010
       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: casa-feedback@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             : #include <components/ComponentModels/FluxCalcLogFreqPolynomial.h>
      27             : #include <casacore/casa/BasicSL/String.h>
      28             : #include <casacore/measures/Measures/MFrequency.h>
      29             : 
      30             : // Handy for passing anonymous arrays to functions.
      31             : #include <casacore/scimath/Mathematics/RigidVector.h>
      32             : 
      33             : #include <map>
      34             : 
      35             : #include <casacore/casa/Logging/LogIO.h>
      36             : 
      37             : using namespace casacore;
      38             : namespace casa { //# NAMESPACE CASA - BEGIN
      39             : 
      40       12104 : Bool FluxCalcLogFreqPolynomial::operator()(Flux<Double>& value,
      41             :                                            Flux<Double>& error,
      42             :                                            const MFrequency& mfreq,
      43             :                                            const Bool updatecoeffs)
      44             : {
      45       24208 :   LogIO os(LogOrigin("FluxCalcLogFreqPolynomial", "()", WHERE));
      46       12104 :   Double dt = log10(mfreq.get(freqUnit_p).getValue());
      47       12104 :   if (updatecoeffs || coeffs_p(0).nelements()==0) {
      48       12094 :     coeffs_p(0).resize();
      49       12094 :     coeffs_p(1).resize();
      50       12094 :     coeffs_p=getCurrentCoeffs();
      51             :   }
      52             :   {
      53       12104 :      validfrange_p=getValidFreqRange();
      54       12104 :      Double ghzinfreq = mfreq.get("GHz").getValue();
      55       12104 :      Double minfreq = validfrange_p(0).get("GHz").getValue();
      56       12104 :      Double maxfreq = validfrange_p(1).get("GHz").getValue();
      57       12104 :      if (!(minfreq==0.0 && maxfreq ==0.0)) {
      58             :        //cerr<<" range:"<<minfreq<<","<<maxfreq<<"  infreq="<<ghzinfreq<<endl;
      59        7776 :         if (minfreq > ghzinfreq || maxfreq < ghzinfreq) {
      60             :           //cerr<<"Input "<<String::toString(mfreq)<<" is out of the valid frequency range of the flux model: "<<String::toString(validfrange_p)<<endl;
      61           0 :           os<<LogIO::WARN<<"Input "<<String::toString(mfreq) 
      62             :                          <<" (Hz) is out of the valid frequency range of the flux model: "
      63           0 :                          <<String::toString(validfrange_p)<<" (Hz)"<<LogIO::POST;
      64             :         }
      65             :     }
      66             :   }
      67       12104 :   Double fluxCoeff = coeffs_p(0)[coeffs_p(0).nelements() - 1];
      68       58640 :   for(Int order = coeffs_p(0).nelements() - 2; order >= 0; --order)
      69       46536 :     fluxCoeff = fluxCoeff * dt + coeffs_p(0)[order];
      70             :   
      71       12104 :   Double coeffErr = 0.0;
      72       12104 :   Double dtpow = 1.0;
      73       70704 :   for(uInt order = 0; order < coeffs_p(1).nelements(); ++order){
      74       58600 :     coeffErr += square(coeffs_p(1)[order] * dtpow);
      75       58600 :     dtpow *= dt;
      76             :   }
      77             : 
      78       12104 :   Double fluxDensity = pow(10.0, fluxCoeff);
      79       12104 :   Double fluxError = coeffErr > 0.0 ? C::ln10 * fluxDensity * sqrt(coeffErr) : 0.0;
      80             : 
      81             :   //cerr<<"FluxDensity=="<<fluxDensity<<endl;
      82       12104 :   value.setValue(fluxDensity);
      83       12104 :   error.setValue(fluxError);
      84             : 
      85             :   // In this case the hard part of matching the std and src has already been done.
      86       12104 :   return true;
      87       12104 : }
      88             : 
      89         103 : Bool FluxCalcLogFreqPolynomial::setSource(const String& sourceName, const MDirection& sourceDir)
      90             : {
      91         103 :   Bool success = FluxCalcVQS::setSource(sourceName, sourceDir);
      92             : 
      93         103 :   if(success)
      94         103 :     success = setSourceCoeffs();
      95         103 :   return success;
      96             : }
      97             : 
      98         103 : void FluxCalcLogFreqPolynomial::setFreqUnit(const String& freqUnit)
      99             : {
     100         103 :   freqUnit_p = freqUnit;
     101         103 : }
     102             : 
     103           0 : void FluxCalcLogFreqPolynomial::fill_coeffs(const Vector<Float>& lfv)
     104             : {
     105           0 :   coeffs_p(0).resize();
     106           0 :   coeffs_p(1).resize();
     107           0 :   coeffs_p(0) = lfv;
     108           0 : }
     109             : 
     110             : 
     111          10 : FluxCalcLogFreqBrokenPolynomial::FluxCalcLogFreqBrokenPolynomial() :
     112             :   FluxCalcLogFreqPolynomial::FluxCalcLogFreqPolynomial(),
     113          10 :   break_freq_p(Quantity(0.0, "Hz"))
     114             : {
     115          10 : }
     116             : 
     117          10 : Bool FluxCalcLogFreqBrokenPolynomial::operator()(Flux<Double>& value,
     118             :                                                  Flux<Double>& error,
     119             :                                                  const MFrequency& mfreq,
     120             :                                                  const Bool updatecoeffs)
     121             : {
     122             :   if (updatecoeffs) {
     123             :     //do nothing for now
     124             :     ;
     125             :   }
     126          10 :   Double break_freq_in_Hz = break_freq_p.get("Hz").getValue();
     127             :   
     128          10 :   if(break_freq_in_Hz > 0.0){
     129           0 :     if(mfreq.get("Hz").getValue() > break_freq_in_Hz){
     130           0 :       if(in_low_state_p)
     131           0 :         fill_coeffs(high_coeffs_p);
     132             :     }
     133           0 :     else if(!in_low_state_p)
     134           0 :       fill_coeffs(low_coeffs_p);
     135             :   }
     136             : 
     137             :   // Proceed...
     138          10 :   return FluxCalcLogFreqPolynomial::operator()(value, error, mfreq);
     139             : }
     140             : 
     141           0 : Bool FluxCalcLogFreqPolynomialSH :: operator()( Flux<Double>& value,
     142             :                                                 Flux<Double>& error,
     143             :                                                 const MFrequency& mfreq,
     144             :                                                 const Bool /* updatecoeffs */)
     145             : {
     146           0 :   Double S = 0.;
     147           0 :   Double dS2 = 0.;
     148           0 :   coeffs_p(0).resize();
     149           0 :   coeffs_p(1).resize();
     150           0 :   coeffs_p=getCurrentCoeffs();
     151           0 :   if ( coeffs_p( 0 ).nelements() > 0 ) {
     152           0 :     Double logF = log10( mfreq.get( freqUnit_p ).getValue() / 0.150 );
     153           0 :     Double logS = 0.;
     154           0 :     for ( uInt order = coeffs_p( 0 ).nelements() - 1; order >= 1; --order)
     155           0 :       logS = ( logS + coeffs_p( 0 )[ order ] ) * logF;
     156           0 :     S = coeffs_p( 0 )[ 0 ] * pow( 10.0, logS );
     157           0 :     if ( coeffs_p( 1 ).nelements() > 0 ) {
     158           0 :       for ( uInt order = coeffs_p( 1 ).nelements() - 1; order >= 1; --order)
     159           0 :         dS2 = ( dS2 + square( coeffs_p( 1 )[ order ] ) ) * square( logF );
     160           0 :       dS2 = square( S * coeffs_p( 1 )[ 0 ] / coeffs_p( 0 )[ 0 ] ) +
     161           0 :             square( C::ln10 * S ) * dS2;
     162             :     }
     163             :   }
     164             : 
     165           0 :   Double dS = ( dS2 > 0.0 ) ? sqrt( dS2 ) : 0.0;
     166           0 :   value.setValue( S );
     167           0 :   error.setValue( dS );
     168             : 
     169           0 :   return true;
     170             : }
     171           0 : Bool FluxCalcLogFreqPolynomialSH::setSource(const String& sourceName, const MDirection& sourceDir)
     172             : {
     173           0 :   Bool success = FluxCalcVQS::setSource(sourceName, sourceDir);
     174             : 
     175           0 :   if(success)
     176           0 :     success = setSourceCoeffs();
     177           0 :   return success;
     178             : }
     179             : 
     180           0 : void FluxCalcLogFreqPolynomialSH::setFreqUnit(const String& freqUnit)
     181             : {
     182           0 :   freqUnit_p = freqUnit;
     183           0 : }
     184             : } //# NAMESPACE CASA - END