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: 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: 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        2666 : Bool GaussianDeconvolver::deconvolve(
      36             :         Angular2DGaussian& deconvolvedSize,
      37             :         const Angular2DGaussian& convolvedSize,
      38             :         const GaussianBeam& beam
      39             : ) {
      40        5332 :         Unit radians(String("rad"));
      41        5332 :         Unit positionAngleModelUnit = deconvolvedSize.getPA(false).getFullUnit();
      42        5332 :         Unit majorAxisModelUnit = deconvolvedSize.getMajor().getFullUnit();
      43        5332 :         Unit minorAxisModelUnit = deconvolvedSize.getMinor().getFullUnit();
      44             : 
      45             :         // Get values in radians
      46        2666 :         Double majorSource = convolvedSize.getMajor().getValue(radians);
      47        2666 :         Double minorSource = convolvedSize.getMinor().getValue(radians);
      48        2666 :         Double thetaSource = convolvedSize.getPA(true).getValue(radians);
      49        2666 :         Double majorBeam = beam.getMajor().getValue(radians);
      50        2666 :         Double minorBeam = beam.getMinor().getValue(radians);
      51        2666 :         Double thetaBeam = beam.getPA(true).getValue(radians);
      52             :         // Do the sums
      53             : 
      54        2666 :         Double alpha  = square(majorSource*cos(thetaSource)) +
      55        2666 :                 square(minorSource*sin(thetaSource)) -
      56        2666 :                 square(majorBeam*cos(thetaBeam)) -
      57        2666 :                 square(minorBeam*sin(thetaBeam));
      58        2666 :         Double beta   = square(majorSource*sin(thetaSource)) +
      59        2666 :                 square(minorSource*cos(thetaSource)) -
      60        2666 :                 square(majorBeam*sin(thetaBeam)) -
      61        2666 :                 square(minorBeam*cos(thetaBeam));
      62        2666 :         Double gamma  = 2 * ( (square(minorSource)-square(majorSource))*sin(thetaSource)*cos(thetaSource) -
      63        2666 :                 (square(minorBeam)-square(majorBeam))*sin(thetaBeam)*cos(thetaBeam) );
      64             :         // Set result in radians
      65             : 
      66        2666 :         Double s = alpha + beta;
      67        2666 :         Double t = sqrt(square(alpha-beta) + square(gamma));
      68        2666 :         Double limit = min(majorSource,minorSource);
      69        2666 :         limit = min(limit,majorBeam);
      70        2666 :         limit = min(limit,minorBeam);
      71        2666 :         limit = 0.1*limit*limit;
      72             : 
      73        2666 :         if(alpha<0.0 || beta<0.0 || s<t) {
      74          98 :                 if(0.5*(s-t)<limit && alpha>-limit && beta>-limit) {
      75             :                         // Point source. Fill in values of beam
      76         186 :                         deconvolvedSize = GaussianBeam(
      77         186 :                                 Quantity(beam.getMajor().get(majorAxisModelUnit)),
      78         186 :                                 Quantity(beam.getMinor().get(minorAxisModelUnit)),
      79         186 :                                 Quantity(beam.getPA(true).get(positionAngleModelUnit))
      80          93 :                         );
      81             :                         // unwrap
      82          93 :                         deconvolvedSize.setPA(deconvolvedSize.getPA(true));
      83          93 :                         return true;
      84             :                 }
      85             :                 else {
      86           5 :                         throw AipsError("Source may be only (slightly) resolved in one direction");
      87             :                 }
      88             :         }
      89        5136 :         Quantity majax(sqrt(0.5*(s+t)), radians);
      90        2568 :         majax.convert(majorAxisModelUnit);
      91        5136 :         Quantity minax(sqrt(0.5*(s-t)), radians);
      92        2568 :         minax.convert(minorAxisModelUnit);
      93             :         Quantity pa(
      94        2568 :                 abs(gamma)+abs(alpha-beta) == 0.0
      95        2568 :                         ? 0.0
      96        2461 :                         : 0.5*atan2(-gamma,alpha-beta),
      97        2568 :                 radians);
      98        2568 :         pa.convert(positionAngleModelUnit);
      99        2568 :         deconvolvedSize = GaussianBeam(majax, minax, pa);
     100             :         // unwrap
     101        2568 :         deconvolvedSize.setPA(deconvolvedSize.getPA(true));
     102        2568 :         return false;
     103             : }
     104             : 
     105             : 
     106             : 
     107             : } //# NAMESPACE CASA - END
     108             :