Line data    Source code

       1             : //# ImageMetaData.cc: Meta information for Images
       2             : //# Copyright (C) 2009
       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             : //# $Id: ImageMetaData.cc 20886 2010-04-29 14:06:56Z gervandiepen $
      27             : 
      28             : #include <components/ComponentModels/GaussianDeconvolver.h>
      29             : #include <casacore/casa/BasicMath/Math.h>
      30             : #include <casacore/scimath/Mathematics/GaussianBeam.h>
      31             : 
      32             : using namespace casacore;
      33             : namespace casa {
      34             : 
      35        2654 : Bool GaussianDeconvolver::deconvolve(
      36             :         Angular2DGaussian& deconvolvedSize,
      37             :         const Angular2DGaussian& convolvedSize,
      38             :         const GaussianBeam& beam
      39             : ) {
      40        2654 :         Unit radians(String("rad"));
      41        2654 :         Unit positionAngleModelUnit = deconvolvedSize.getPA(false).getFullUnit();
      42        2654 :         Unit majorAxisModelUnit = deconvolvedSize.getMajor().getFullUnit();
      43        2654 :         Unit minorAxisModelUnit = deconvolvedSize.getMinor().getFullUnit();
      44             : 
      45             :         // Get values in radians
      46        2654 :         Double majorSource = convolvedSize.getMajor().getValue(radians);
      47        2654 :         Double minorSource = convolvedSize.getMinor().getValue(radians);
      48        2654 :         Double thetaSource = convolvedSize.getPA(true).getValue(radians);
      49        2654 :         Double majorBeam = beam.getMajor().getValue(radians);
      50        2654 :         Double minorBeam = beam.getMinor().getValue(radians);
      51        2654 :         Double thetaBeam = beam.getPA(true).getValue(radians);
      52             :         // Do the sums
      53             : 
      54        2654 :         Double alpha  = square(majorSource*cos(thetaSource)) +
      55        2654 :                 square(minorSource*sin(thetaSource)) -
      56        2654 :                 square(majorBeam*cos(thetaBeam)) -
      57        2654 :                 square(minorBeam*sin(thetaBeam));
      58        2654 :         Double beta   = square(majorSource*sin(thetaSource)) +
      59        2654 :                 square(minorSource*cos(thetaSource)) -
      60        2654 :                 square(majorBeam*sin(thetaBeam)) -
      61        2654 :                 square(minorBeam*cos(thetaBeam));
      62        2654 :         Double gamma  = 2 * ( (square(minorSource)-square(majorSource))*sin(thetaSource)*cos(thetaSource) -
      63        2654 :                 (square(minorBeam)-square(majorBeam))*sin(thetaBeam)*cos(thetaBeam) );
      64             :         // Set result in radians
      65             : 
      66        2654 :         Double s = alpha + beta;
      67        2654 :         Double t = sqrt(square(alpha-beta) + square(gamma));
      68        2654 :         Double limit = min(majorSource,minorSource);
      69        2654 :         limit = min(limit,majorBeam);
      70        2654 :         limit = min(limit,minorBeam);
      71        2654 :         limit = 0.1*limit*limit;
      72             : 
      73        2654 :         if(alpha<0.0 || beta<0.0 || s<t) {
      74          92 :                 if(0.5*(s-t)<limit && alpha>-limit && beta>-limit) {
      75             :                         // Point source. Fill in values of beam
      76         174 :                         deconvolvedSize = GaussianBeam(
      77         174 :                                 Quantity(beam.getMajor().get(majorAxisModelUnit)),
      78         174 :                                 Quantity(beam.getMinor().get(minorAxisModelUnit)),
      79         174 :                                 Quantity(beam.getPA(true).get(positionAngleModelUnit))
      80          87 :                         );
      81             :                         // unwrap
      82          87 :                         deconvolvedSize.setPA(deconvolvedSize.getPA(true));
      83          87 :                         return true;
      84             :                 }
      85             :                 else {
      86           5 :                         throw AipsError("Source may be only (slightly) resolved in one direction");
      87             :                 }
      88             :         }
      89        2562 :         Quantity majax(sqrt(0.5*(s+t)), radians);
      90        2562 :         majax.convert(majorAxisModelUnit);
      91        2562 :         Quantity minax(sqrt(0.5*(s-t)), radians);
      92        2562 :         minax.convert(minorAxisModelUnit);
      93             :         Quantity pa(
      94        2562 :                 abs(gamma)+abs(alpha-beta) == 0.0
      95        2562 :                         ? 0.0
      96        2455 :                         : 0.5*atan2(-gamma,alpha-beta),
      97        2562 :                 radians);
      98        2562 :         pa.convert(positionAngleModelUnit);
      99        2562 :         deconvolvedSize = GaussianBeam(majax, minax, pa);
     100             :         // unwrap
     101        2562 :         deconvolvedSize.setPA(deconvolvedSize.getPA(true));
     102        2562 :         return false;
     103        2669 : }
     104             : 
     105             : 
     106             : 
     107             : } //# NAMESPACE CASA - END
     108             :