Line data    Source code

       1             : //# LimbDarkenedDiskShape.cc: implementation of LimbDarkenedDiskShape.h which defines LimbDarkened Disk shape
       2             : //
       3             : //#  CASA - Common Astronomy Software Applications (http://casa.nrao.edu/)
       4             : //# Copyright (C) 2012
       5             : //# Associated Universities, Inc. Washington DC, USA.
       6             : //#
       7             : //# This library is free software; you can redistribute it and/or modify it
       8             : //# under the terms of the GNU Lesser General Public License as published by
       9             : //# the Free Software Foundation; either version 2.1 of the License, or (at your
      10             : //# option) any later version.
      11             : //#
      12             : //# This library is distributed in the hope that it will be useful, but WITHOUT
      13             : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
      14             : //# FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
      15             : //# License for more details.
      16             : //#
      17             : //# You should have received a copy of the GNU Lesser General Public License
      18             : //# along with this library; if not, write to the Free Software Foundation,
      19             : //# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
      20             : //#
      21             : //#
      22             : //#
      23             : //# $Id$
      24             : 
      25             : #include <components/ComponentModels/LimbDarkenedDiskShape.h>
      26             : //#include <components/ComponentModels/Flux.h>
      27             : #include <casacore/casa/Arrays/Matrix.h>
      28             : #include <casacore/casa/Arrays/Vector.h>
      29             : #include <casacore/casa/Arrays/VectorIter.h>
      30             : #include <casacore/casa/Arrays/ArrayMath.h>
      31             : #include <casacore/casa/Arrays/Array.h>
      32             : #include <casacore/casa/Exceptions/Error.h>
      33             : #include <casacore/casa/Logging/LogIO.h>
      34             : #include <casacore/casa/Logging/LogOrigin.h>
      35             : #include <casacore/casa/BasicSL/Constants.h>
      36             : #include <casacore/casa/BasicMath/Math.h>
      37             : #include <casacore/measures/Measures/MCDirection.h>
      38             : #include <casacore/measures/Measures/MDirection.h>
      39             : #include <casacore/measures/Measures/MeasConvert.h>
      40             : #include <casacore/measures/Measures/MeasRef.h>
      41             : #include <casacore/casa/Quanta/MVAngle.h>
      42             : #include <casacore/casa/Quanta/MVDirection.h>
      43             : #include <casacore/casa/Quanta/Quantum.h>
      44             : #include <casacore/casa/Utilities/Assert.h>
      45             : #include <casacore/casa/BasicSL/String.h>
      46             : 
      47             : #include <gsl/gsl_sf_bessel.h>
      48             : 
      49             : using namespace casacore;
      50             : namespace casa { //# NAMESPACE CASA - BEGIN
      51             : 
      52           0 : LimbDarkenedDiskShape::LimbDarkenedDiskShape()
      53             :   :TwoSidedShape(),
      54           0 :    itsMajValue(Quantity(1,"'").getValue("rad")),
      55           0 :    itsMinValue(Quantity(1,"'").getValue("rad")),
      56           0 :    itsPaValue(Quantity(0,"deg").getValue("rad")),
      57           0 :    itsHeight(1.0/(C::pi*itsMajValue*itsMinValue)),
      58           0 :    itsAttnFactor(0.05)
      59             : {
      60           0 :   DebugAssert(ok(), AipsError);
      61           0 : }
      62             : 
      63           0 : LimbDarkenedDiskShape::LimbDarkenedDiskShape(const MDirection& direction,
      64             :                                          const Quantum<Double>& majorAxis,
      65             :                                          const Quantum<Double>& minorAxis,
      66             :                                          const Quantum<Double>& positionAngle,
      67           0 :                                          const Float& n)
      68           0 :   :TwoSidedShape(direction, majorAxis.getFullUnit(), minorAxis.getFullUnit(), positionAngle.getFullUnit()),
      69           0 :    itsMajValue(majorAxis.getValue("rad")),
      70           0 :    itsMinValue(minorAxis.getValue("rad")),
      71           0 :    itsPaValue(positionAngle.getValue("rad")),
      72           0 :    itsHeight(1.0/(C::pi*itsMajValue*itsMinValue)),
      73           0 :    itsAttnFactor(n)
      74             : {
      75           0 :   DebugAssert(ok(), AipsError);
      76           0 : }
      77             : 
      78             : 
      79           0 : LimbDarkenedDiskShape::LimbDarkenedDiskShape(const MDirection& direction, 
      80             :                     const Quantum<Double>& width,
      81             :                     const Double axialRatio,
      82             :                     const Quantum<Double>& positionAngle, 
      83           0 :                     const Float& n)
      84           0 :   :TwoSidedShape(direction, width.getFullUnit(),
      85           0 :                     width.getFullUnit(), positionAngle.getFullUnit()),
      86           0 :    itsMajValue(width.getValue("rad")),
      87           0 :    itsMinValue(itsMajValue*axialRatio),
      88           0 :    itsPaValue(positionAngle.getValue("rad")),
      89           0 :    itsHeight(1.0/(C::pi*itsMajValue*itsMinValue)),
      90           0 :    itsAttnFactor(n)
      91             : {
      92           0 :   DebugAssert(ok(), AipsError);
      93           0 : }
      94             : 
      95           0 : LimbDarkenedDiskShape::LimbDarkenedDiskShape(const LimbDarkenedDiskShape& other)
      96             :   :TwoSidedShape(other),
      97           0 :    itsMajValue(other.itsMajValue),
      98           0 :    itsMinValue(other.itsMinValue),
      99           0 :    itsPaValue(other.itsPaValue),
     100           0 :    itsHeight(other.itsHeight),
     101           0 :    itsAttnFactor(other.itsAttnFactor)
     102             : {
     103           0 :   DebugAssert(ok(), AipsError);
     104           0 : }
     105             : 
     106             :   // The destructor
     107           0 : LimbDarkenedDiskShape::~LimbDarkenedDiskShape() {
     108           0 :   DebugAssert(ok(), AipsError);
     109           0 : }
     110             :   //#! Operators
     111             :   //The assignment operator
     112           0 : LimbDarkenedDiskShape& LimbDarkenedDiskShape::operator=(const LimbDarkenedDiskShape& other) {
     113           0 :   if (this != &other) {
     114           0 :     TwoSidedShape::operator=(other);
     115           0 :     itsMajValue = other.itsMajValue;
     116           0 :     itsMinValue = other.itsMinValue;
     117           0 :     itsPaValue = other.itsPaValue;
     118           0 :     itsHeight = other.itsHeight;
     119           0 :     itsAttnFactor = other.itsAttnFactor;
     120             :   }
     121           0 :   DebugAssert(ok(), AipsError);
     122           0 :   return *this;
     123             : }
     124             : 
     125             : //#! General Member Functions
     126             : // get the type of the shape (always returns ComponentType::LimbDakenDisk)
     127           0 : ComponentType::Shape LimbDarkenedDiskShape::type() const {
     128           0 :     DebugAssert(ok(), AipsError);
     129           0 :     return ComponentType::LDISK;
     130             : }
     131             : 
     132             : 
     133             : // use diskshape ones?
     134           0 : void LimbDarkenedDiskShape::setWidthInRad(const Double majorAxis,
     135             :                              const Double minorAxis,
     136             :                              const Double positionAngle) {
     137           0 :     itsMajValue = majorAxis;
     138           0 :     itsMinValue = minorAxis;
     139           0 :     itsPaValue = positionAngle;
     140           0 :     AlwaysAssert(itsMajValue > 0 && itsMinValue > 0 && itsMajValue >=itsMinValue,
     141             :                  AipsError);
     142           0 :     itsHeight = 1.0/(C::pi*itsMajValue*itsMinValue);
     143           0 :     DebugAssert(ok(), AipsError);
     144           0 : }
     145             : //
     146           0 : Double LimbDarkenedDiskShape::majorAxisInRad() const {
     147           0 :   DebugAssert(ok(), AipsError);
     148           0 :   return itsMajValue;
     149             : }
     150             : 
     151           0 : Double LimbDarkenedDiskShape::minorAxisInRad() const {
     152           0 :   DebugAssert(ok(), AipsError);
     153           0 :   return itsMinValue;
     154             : }
     155             : 
     156           0 : Double LimbDarkenedDiskShape::positionAngleInRad() const {
     157           0 :   DebugAssert(ok(), AipsError);
     158           0 :   return itsPaValue;
     159             : }
     160             : 
     161           0 : Float LimbDarkenedDiskShape::getAttnFactor() const {
     162           0 :   DebugAssert(ok(), AipsError);
     163           0 :   return itsAttnFactor;
     164             : }
     165             :   //set n factor in darkening equation, I=I0(1-rho^2)^n/2
     166           0 : void LimbDarkenedDiskShape::setAttnFactor(const Float attnFactor) {
     167           0 :   itsAttnFactor=attnFactor;  
     168           0 : }
     169             : 
     170           0 : Vector<Double> LimbDarkenedDiskShape::optParameters() const {
     171           0 :   DebugAssert(ok(), AipsError);
     172           0 :   Vector<Double> optparm(1);
     173           0 :   optparm(0) = (Double)getAttnFactor();
     174           0 :   return optparm;
     175             : }
     176             : 
     177           0 : void LimbDarkenedDiskShape::setOptParameters(const Vector<Double>& newOptParms) {
     178           0 :   DebugAssert(ok(), AipsError);
     179           0 :   setAttnFactor((Float)newOptParms(0));
     180           0 : }
     181             : 
     182             :   // Calculate the proportion of the flux that is in a pixel of specified size
     183             :   // centered in the specified direction. The returned value will always be
     184             :   // between zero and one (inclusive).
     185           0 : Double LimbDarkenedDiskShape::sample(const MDirection& direction,
     186             :                         const MVAngle& pixelLatSize,
     187             :                         const MVAngle& pixelLongSize) const {
     188           0 :   DebugAssert(ok(), AipsError);
     189           0 :   const MDirection& compDir(refDirection());
     190           0 :   const MDirection::Ref& compDirFrame(compDir.getRef());
     191           0 :   const MDirection::MVType* compDirValue = &(compDir.getValue());
     192           0 :   Bool deleteValue = false;
     193             :   // Convert direction to the same frame as the reference direction
     194           0 :   if (ComponentShape::differentRefs(direction.getRef(), compDirFrame)) {
     195           0 :     compDirValue = new MDirection::MVType
     196           0 :       (MDirection::Convert(compDir, direction.getRef())().getValue());
     197           0 :     deleteValue = true;
     198             :   }
     199           0 :   Double retVal = calcSample(*compDirValue, direction.getValue(),
     200           0 :                              itsMajValue/2.0, itsMinValue/2.0,
     201           0 :                              itsHeight*pixelLatSize.radian()*
     202           0 :                              pixelLongSize.radian());
     203           0 :   if (deleteValue) delete compDirValue;
     204           0 :   return retVal;
     205             : }
     206             : 
     207             : 
     208             :   // Same as the previous function except that many directions can be sampled
     209             :   // at once. The reference frame and pixel size must be the same for all the
     210             :   // specified directions.
     211           0 : void LimbDarkenedDiskShape::sample(Vector<Double>& scale,
     212             :                       const Vector<MDirection::MVType>& directions,
     213             :                       const MDirection::Ref& refFrame,
     214             :                       const MVAngle& pixelLatSize,
     215             :                       const MVAngle& pixelLongSize) const {
     216           0 :   DebugAssert(ok(), AipsError);
     217           0 :   const uInt nSamples = directions.nelements();
     218           0 :   DebugAssert(scale.nelements() == nSamples, AipsError);
     219             : 
     220           0 :   const MDirection& compDir(refDirection());
     221           0 :   const MDirection::Ref& compDirFrame(compDir.getRef());
     222           0 :   const MDirection::MVType* compDirValue = &(compDir.getValue());
     223           0 :   Bool deleteValue = false;
     224             :   // Convert direction to the same frame as the reference direction
     225           0 :   if (refFrame != compDirFrame) {
     226           0 :     compDirValue = new MDirection::MVType
     227           0 :       (MDirection::Convert(compDir, refFrame)().getValue());
     228           0 :     deleteValue = true;
     229             :   }
     230           0 :   const Double majRad = itsMajValue/2.0;
     231           0 :   const Double minRad = itsMinValue/2.0;
     232           0 :   const Double pixValue = itsHeight *
     233           0 :     pixelLatSize.radian() * pixelLongSize.radian();
     234           0 :   for (uInt i = 0; i < nSamples; i++) {
     235           0 :     scale(i) = calcSample(*compDirValue, directions(i),
     236             :                           majRad, minRad, pixValue);
     237             :   }
     238           0 :   if (deleteValue) delete compDirValue;
     239           0 : }
     240             : 
     241           0 : DComplex LimbDarkenedDiskShape::visibility(const Vector<Double>& uvw,
     242             :                               const Double& frequency) const {
     243           0 :  DebugAssert(uvw.nelements() == 3, AipsError);
     244           0 :   DebugAssert(frequency > 0, AipsError);
     245           0 :   DebugAssert(ok(), AipsError);
     246           0 :   Double u = uvw(0);
     247           0 :   Double v = uvw(1);
     248           0 :   if (near(u + v, 0.0)) return DComplex(1.0, 0.0);
     249           0 :   if (!nearAbs(itsPaValue, 0.0)) {
     250           0 :     rotateVis(u, v, cos(itsPaValue), sin(itsPaValue));
     251             :   }
     252           0 :   return DComplex(calcVis(u, v, C::pi * frequency/C::c), 0.0);
     253             : }
     254             : 
     255           0 : void LimbDarkenedDiskShape::visibility(Matrix<DComplex>& scale,
     256             :                            const Matrix<Double>& uvw,
     257             :                            const Vector<Double>& frequency) const {
     258             : 
     259           0 :   scale.resize(uvw.ncolumn(), frequency.nelements());
     260             : 
     261           0 :   VectorIterator<DComplex> iter(scale, 0);
     262           0 :   for ( uInt k =0 ; k < frequency.nelements() ; ++k){
     263           0 :     visibility(iter.vector(), uvw, frequency(k));
     264           0 :     iter.next();
     265             :   }
     266           0 : }
     267             : 
     268           0 : void LimbDarkenedDiskShape::visibility(Vector<DComplex>& scale,
     269             :                            const Matrix<Double>& uvw,
     270             :                            const Double& frequency) const {
     271           0 :   DebugAssert(ok(), AipsError);
     272           0 :   const uInt nSamples = scale.nelements();
     273           0 :   DebugAssert(uvw.ncolumn() == nSamples, AipsError);
     274           0 :   DebugAssert(uvw.nrow() == 3, AipsError);
     275           0 :   DebugAssert(frequency > 0, AipsError);
     276             : 
     277           0 :   Bool doRotation = false;
     278           0 :   Double cpa = 1.0, spa = 0.0;
     279           0 :   if (!nearAbs(itsPaValue, 0.0)) {
     280           0 :     doRotation = true;
     281           0 :     cpa = cos(itsPaValue);
     282           0 :     spa = sin(itsPaValue);
     283             :   }
     284             : 
     285           0 :   const Double factor = C::pi * frequency/C::c;
     286             :   Double u, v;
     287           0 :   for (uInt i = 0; i < nSamples; i++) {
     288           0 :     u = uvw(0, i);
     289           0 :     v = uvw(1, i);
     290             :     // DComplex& thisVis = scale(i);
     291             :     ///    thisVis.imag() = 0.0;
     292           0 :     if (near(u + v, 0.0)) {
     293             :       ///      thisVis.real() = 1.0; // avoids dividing by zero in calcVis(...)
     294           0 :       scale[i] = DComplex(1.0, 0.0); // avoids dividing by zero
     295             :       // in calcVis(...)
     296             :     } else {
     297           0 :       if (doRotation) rotateVis(u, v, cpa, spa);
     298             :       ///      thisVis.real() = calcVis(u, v, factor);
     299           0 :       scale[i] = DComplex(calcVis(u, v, factor), 0.0);
     300             :     }
     301             :   }
     302           0 : }
     303             : 
     304             : 
     305           0 : ComponentShape* LimbDarkenedDiskShape::clone() const {
     306           0 :   DebugAssert(ok(), AipsError);
     307           0 :   ComponentShape* tmpPtr = new LimbDarkenedDiskShape(*this);
     308           0 :   AlwaysAssert(tmpPtr != 0, AipsError);
     309           0 :   return tmpPtr;
     310             : }
     311             : 
     312           0 : Bool LimbDarkenedDiskShape::ok() const {
     313             :   // The LogIO class is only constructed if an error is detected for
     314             :   // performance reasons. Both function static and file static variables
     315             :   // where considered and rejected for this purpose.
     316           0 :   if (!TwoSidedShape::ok()) return false;
     317           0 :   if (itsMajValue <= 0) {
     318           0 :     LogIO logErr(LogOrigin("DiskCompRep", "ok()"));
     319             :     logErr << LogIO::SEVERE << "The major axis width is zero or negative"
     320           0 :            << LogIO::POST;
     321           0 :     return false;
     322             :   }
     323           0 :   if (itsMinValue <= 0) {
     324           0 :     LogIO logErr(LogOrigin("DiskCompRep", "ok()"));
     325             :     logErr << LogIO::SEVERE << "The minor axis width is zero or negative"
     326           0 :            << LogIO::POST;
     327           0 :     return false;
     328             :   }
     329           0 :   if (itsMinValue > itsMajValue) {
     330           0 :     LogIO logErr(LogOrigin("DiskCompRep", "ok()"));
     331             :     logErr << LogIO::SEVERE << "The minor axis width is larger than "
     332             :            << "the major axis width"
     333           0 :            << LogIO::POST;
     334           0 :     return false;
     335             :   }
     336           0 :   if (!near(itsHeight, 1.0/(C::pi*itsMajValue*itsMinValue), 2*C::dbl_epsilon)) {
     337           0 :     LogIO logErr(LogOrigin("DiskCompRep", "ok()"));
     338             :     logErr << LogIO::SEVERE << "The disk shape does not have"
     339             :            << " unit area"
     340           0 :            << LogIO::POST;
     341           0 :     return false;
     342             :   }
     343           0 :   return true;
     344             : }
     345             : 
     346             : 
     347             :   // return a pointer to this object.
     348           0 : const ComponentShape* LimbDarkenedDiskShape::getPtr() const {
     349           0 :     return this;
     350             : }
     351             : 
     352           0 : String LimbDarkenedDiskShape::sizeToString() const {
     353             :         return TwoSidedShape::sizeToString(
     354           0 :                 Quantity(itsMajValue, "rad"),
     355           0 :                 Quantity(itsMinValue, "rad"),
     356           0 :                 Quantity(itsPaValue, "rad"),
     357             :                 false
     358           0 :         );
     359             : }
     360             : 
     361             : //need to modify here
     362           0 : Double LimbDarkenedDiskShape::calcSample(const MDirection::MVType& compDirValue,
     363             :                     const MDirection::MVType& dirVal,
     364             :                     const Double majRad, const Double minRad,
     365             :                     const Double pixValue) const {
     366           0 :   const Double separation = compDirValue.separation(dirVal);
     367           0 :   if (separation <= majRad) {
     368           0 :     const Double pa = compDirValue.positionAngle(dirVal) - itsPaValue;
     369           0 :     const Double x = abs(separation*cos(pa));
     370           0 :     const Double y = abs(separation*sin(pa));
     371           0 :     if ((x <= majRad) &&
     372           0 :         (y <= minRad) &&
     373           0 :         (y <= minRad * sqrt(1 - square(x/majRad)))) {
     374           0 :       return pixValue;
     375             :     }
     376             :   }
     377           0 :   return 0.0;
     378             : }
     379             : 
     380           0 : Double LimbDarkenedDiskShape::calcVis(Double u, Double v, const Double factor) const {
     381           0 :   u *= itsMinValue;
     382           0 :   v *= itsMajValue;
     383           0 :   const double r = hypot(u, v) * factor;
     384           0 :   const double eta = 1 + itsAttnFactor/2.0;
     385             :   //return 2.0 * j1(r)/r;
     386             :   // Vi(u,v) for the limb-darkened disk from ALMA memo #594 
     387             :   // assume u, v are != 0.0 (in such case Vi(u,v)=Vo, handled in visibility())
     388           0 :   return pow(C::e,lgamma(eta + 1))*pow(2.0/r,eta)*gsl_sf_bessel_Jnu(eta,r); 
     389             : }
     390             : 
     391             : 
     392           0 : void LimbDarkenedDiskShape::rotateVis(Double& u, Double& v,
     393             :                         const Double cpa, const Double spa) {
     394           0 :   const Double utemp = u;
     395           0 :   u = u * cpa - v * spa;
     396           0 :   v = utemp * spa + v * cpa;
     397           0 : }
     398             : 
     399             : } //# NAMESPACE CASA - END