#include <components/ComponentModels/GaussianDeconvolver.h>
#include <casacore/casa/BasicMath/Math.h>
#include <scimath/Mathematics/GaussianBeam.h>
using namespace casacore;
Bool GaussianDeconvolver::deconvolve(
Angular2DGaussian& deconvolvedSize,
const Angular2DGaussian& convolvedSize,
Unit radians(String("rad"));
Unit positionAngleModelUnit = deconvolvedSize.getPA(false).getFullUnit();
Unit majorAxisModelUnit = deconvolvedSize.getMajor().getFullUnit();
Unit minorAxisModelUnit = deconvolvedSize.getMinor().getFullUnit();
Double majorSource = convolvedSize.getMajor().getValue(radians);
Double minorSource = convolvedSize.getMinor().getValue(radians);
Double thetaSource = convolvedSize.getPA(true).getValue(radians);
Double majorBeam = beam.getMajor().getValue(radians);
Double minorBeam = beam.getMinor().getValue(radians);
Double thetaBeam = beam.getPA(true).getValue(radians);
Double alpha = square(majorSource*cos(thetaSource)) +
square(minorSource*sin(thetaSource)) -
square(majorBeam*cos(thetaBeam)) -
square(minorBeam*sin(thetaBeam));
Double beta = square(majorSource*sin(thetaSource)) +
square(minorSource*cos(thetaSource)) -
square(majorBeam*sin(thetaBeam)) -
square(minorBeam*cos(thetaBeam));
Double gamma = 2 * ( (square(minorSource)-square(majorSource))*sin(thetaSource)*cos(thetaSource) -
(square(minorBeam)-square(majorBeam))*sin(thetaBeam)*cos(thetaBeam) );
Double t = sqrt(square(alpha-beta) + square(gamma));
Double limit = min(majorSource,minorSource);
limit = min(limit,majorBeam);
limit = min(limit,minorBeam);
if(alpha<0.0 || beta<0.0 || s<t) {
if(0.5*(s-t)<limit && alpha>-limit && beta>-limit) {
deconvolvedSize = GaussianBeam(
Quantity(beam.getMajor().get(majorAxisModelUnit)),
Quantity(beam.getMinor().get(minorAxisModelUnit)),
Quantity(beam.getPA(true).get(positionAngleModelUnit))
deconvolvedSize.setPA(deconvolvedSize.getPA(true));
throw AipsError("Source may be only (slightly) resolved in one direction");
Quantity majax(sqrt(0.5*(s+t)), radians);
majax.convert(majorAxisModelUnit);
Quantity minax(sqrt(0.5*(s-t)), radians);