Line data Source code
1 : //# SkyComponent.cc: this defines SkyComponent
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2003
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: SkyComponent.cc 21292 2012-11-28 14:58:19Z gervandiepen $
27 :
28 : #include <casacore/casa/Quanta/QMath.h>
29 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
30 : #include <components/ComponentModels/SkyComponent.h>
31 : #include <components/ComponentModels/ComponentShape.h>
32 : #include <components/ComponentModels/Flux.h>
33 : #include <components/ComponentModels/SkyCompRep.h>
34 : #include <components/ComponentModels/SpectralModel.h>
35 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
36 : #include <casacore/casa/Arrays/Vector.h>
37 : #include <casacore/casa/Containers/RecordInterface.h>
38 : #include <casacore/casa/Exceptions/Error.h>
39 : #include <iomanip>
40 : #include <casacore/casa/Logging/LogIO.h>
41 : #include <casacore/casa/Logging/LogOrigin.h>
42 : #include <casacore/casa/Quanta/MVTime.h>
43 : #include <casacore/casa/Utilities/Precision.h>
44 :
45 : #include <casacore/measures/Measures/MDirection.h>
46 : #include <casacore/measures/Measures/MFrequency.h>
47 :
48 : #include <casacore/casa/Quanta/MVAngle.h>
49 : #include <casacore/casa/Utilities/Assert.h>
50 : #include <casacore/casa/BasicSL/String.h>
51 :
52 : #include <iomanip>
53 :
54 : using namespace casacore;
55 : namespace casa { //# NAMESPACE CASA - BEGIN
56 :
57 401 : SkyComponent::SkyComponent()
58 401 : :itsCompPtr(new SkyCompRep)
59 401 : {}
60 :
61 1 : SkyComponent::SkyComponent(const ComponentType::Shape& shape)
62 1 : :itsCompPtr(new SkyCompRep(shape))
63 1 : {}
64 :
65 1 : SkyComponent::SkyComponent(const ComponentType::Shape& shape,
66 1 : const ComponentType::SpectralShape& spectralModel)
67 1 : :itsCompPtr(new SkyCompRep(shape, spectralModel))
68 1 : {}
69 :
70 1 : SkyComponent::SkyComponent(const Flux<Double>& flux,
71 : const ComponentShape& shape,
72 1 : const SpectralModel& spectrum)
73 1 : :itsCompPtr(new SkyCompRep(flux, shape, spectrum))
74 1 : {}
75 :
76 401 : SkyComponent::SkyComponent(const SkyComponent& other)
77 : :SkyCompBase(other),
78 401 : itsCompPtr (other.itsCompPtr)
79 401 : {}
80 :
81 803 : SkyComponent::~SkyComponent() {}
82 :
83 201 : SkyComponent& SkyComponent::operator=(const SkyComponent& other) {
84 201 : if (this != &other)
85 201 : itsCompPtr = other.itsCompPtr;
86 201 : return *this;
87 : }
88 :
89 402 : Flux<Double>& SkyComponent::flux() {
90 402 : return itsCompPtr->flux();
91 : }
92 :
93 0 : const Flux<Double>& SkyComponent::flux() const {
94 0 : return itsCompPtr->flux();
95 : }
96 :
97 200 : const ComponentShape& SkyComponent::shape() const {
98 200 : return itsCompPtr->shape();
99 : }
100 :
101 405 : ComponentShape& SkyComponent::shape() {
102 405 : return itsCompPtr->shape();
103 : }
104 :
105 0 : void SkyComponent::setShape(const ComponentShape& newShape) {
106 0 : itsCompPtr->setShape(newShape);
107 0 : }
108 :
109 0 : const SpectralModel& SkyComponent::spectrum() const {
110 0 : return itsCompPtr->spectrum();
111 : }
112 :
113 0 : SpectralModel& SkyComponent::spectrum() {
114 0 : return itsCompPtr->spectrum();
115 : }
116 :
117 0 : void SkyComponent::setSpectrum(const SpectralModel& newSpectrum) {
118 0 : itsCompPtr->setSpectrum(newSpectrum);
119 0 : }
120 :
121 0 : String& SkyComponent::label() {
122 0 : return itsCompPtr->label();
123 : }
124 :
125 0 : const String& SkyComponent::label() const {
126 0 : return itsCompPtr->label();
127 : }
128 :
129 0 : Vector<Double>& SkyComponent::optionalParameters() {
130 0 : return itsCompPtr->optionalParameters();
131 : }
132 :
133 0 : const Vector<Double>& SkyComponent::optionalParameters() const {
134 0 : return itsCompPtr->optionalParameters();
135 : }
136 :
137 0 : Bool SkyComponent::isPhysical() const {
138 0 : return itsCompPtr->isPhysical();
139 : }
140 :
141 0 : Flux<Double> SkyComponent::sample(const MDirection& direction,
142 : const MVAngle& pixelLatSize,
143 : const MVAngle& pixelLongSize,
144 : const MFrequency& centerFrequency) const {
145 0 : return itsCompPtr->sample(direction, pixelLatSize, pixelLongSize,
146 0 : centerFrequency);
147 : }
148 :
149 0 : void SkyComponent::sample(Cube<Double>& samples, const Unit& reqUnit,
150 : const Vector<MVDirection>& directions,
151 : const MeasRef<MDirection>& dirRef,
152 : const MVAngle& pixelLatSize,
153 : const MVAngle& pixelLongSize,
154 : const Vector<MVFrequency>& frequencies,
155 : const MeasRef<MFrequency>& freqRef) const {
156 0 : itsCompPtr->sample(samples, reqUnit,
157 : directions, dirRef, pixelLatSize, pixelLongSize,
158 0 : frequencies, freqRef);
159 0 : }
160 :
161 0 : Flux<Double> SkyComponent::visibility(const Vector<Double>& uvw,
162 : const Double& frequency) const {
163 0 : return itsCompPtr->visibility(uvw, frequency);
164 : }
165 :
166 200 : void SkyComponent::visibility(Cube<DComplex>& visibilities,
167 : const Matrix<Double>& uvws,
168 : const Vector<Double>& frequencies) const {
169 200 : itsCompPtr->visibility(visibilities, uvws, frequencies);
170 200 : }
171 :
172 200 : Bool SkyComponent::fromRecord(String& errorMessage,
173 : const RecordInterface& record) {
174 200 : return itsCompPtr->fromRecord(errorMessage, record);
175 : }
176 :
177 200 : Bool SkyComponent::toRecord(String& errorMessage,
178 : RecordInterface& record) const {
179 200 : return itsCompPtr->toRecord(errorMessage, record);
180 : }
181 :
182 0 : SkyComponent SkyComponent::copy() const {
183 0 : SkyComponent newComp(flux().copy(), shape(), spectrum());
184 0 : newComp.label() = label();
185 0 : newComp.optionalParameters() = optionalParameters();
186 0 : return newComp;
187 : }
188 :
189 0 : void SkyComponent::fromPixel (Double& fluxRatio, const Vector<Double>& parameters,
190 : const Unit& brightnessUnitIn,
191 : const GaussianBeam& restoringBeam,
192 : const CoordinateSystem& cSys,
193 : ComponentType::Shape componentShape,
194 : Stokes::StokesTypes stokes)
195 : {
196 0 : itsCompPtr->fromPixel(fluxRatio, parameters, brightnessUnitIn, restoringBeam,
197 : cSys, componentShape, stokes);
198 0 : }
199 :
200 0 : Vector<Double> SkyComponent::toPixel (const Unit& brightnessUnitIn,
201 : const GaussianBeam& restoringBeam,
202 : const CoordinateSystem& cSys,
203 : Stokes::StokesTypes stokes) const
204 : {
205 : return itsCompPtr->toPixel(brightnessUnitIn, restoringBeam,
206 0 : cSys, stokes);
207 : }
208 :
209 :
210 :
211 :
212 0 : Bool SkyComponent::ok() const {
213 0 : if (itsCompPtr.null()) {
214 0 : LogIO logErr(LogOrigin("SkyComponent", "ok()"));
215 : logErr << LogIO::SEVERE << "Internal pointer is not pointing to anything"
216 0 : << LogIO::POST;
217 0 : return false;
218 : }
219 0 : if (!itsCompPtr->ok()) {
220 0 : LogIO logErr(LogOrigin("SkyComponent", "ok()"));
221 : logErr << LogIO::SEVERE << "Component representation is not ok"
222 0 : << LogIO::POST;
223 0 : return false;
224 : }
225 0 : return true;
226 : }
227 :
228 0 : String SkyComponent::summarize(
229 : const DirectionCoordinate *const &dc, Bool longErrOnGreatCircle
230 : ) const {
231 0 : ostringstream ldpar;
232 0 : if (shape().type()==ComponentType::LDISK) {
233 0 : ldpar << " (limb-darkening exponent: "<<optionalParameters()(0) <<" )"<<endl;
234 : }
235 0 : ostringstream summary;
236 0 : summary << "SUMMARY OF COMPONENT " << label() << endl;
237 0 : summary << "Shape: " << shape().ident() << ldpar.str() <<endl;
238 0 : const Flux<Double> myFlux = flux();
239 0 : Quantum<Vector<std::complex<double> > > fluxValue;
240 0 : myFlux.value(fluxValue);
241 0 : summary << "Flux density: " << fluxValue << " +/- " << myFlux.errors() << endl;
242 0 : summary << "Spectral model: " << spectrum().ident() << endl;
243 0 : std::shared_ptr<Vector<Double>> pixCoords;
244 0 : summary << "Position: " << positionToString(pixCoords, dc, longErrOnGreatCircle) << endl;
245 0 : summary << "Size: " << endl << shape().sizeToString() << endl;
246 0 : return summary.str();
247 : }
248 :
249 0 : String SkyComponent::positionToString(
250 : std::shared_ptr<casacore::Vector<casacore::Double>>& pixelCoords,
251 : const DirectionCoordinate *const &dc ,Bool longErrOnGreatCircle
252 : ) const {
253 0 : pixelCoords.reset();
254 : // FIXME essentially cut and paste of Gareth's python code. Needs work.
255 0 : ostringstream position;
256 0 : MDirection mdir = shape().refDirection();
257 :
258 0 : Quantity lat = mdir.getValue().getLat("rad");
259 0 : String latString = MVAngle(lat).string(MVAngle::ANGLE_CLEAN, 8);
260 :
261 0 : Quantity longitude = mdir.getValue().getLong("rad");
262 0 : String longString = MVTime(longitude).string(MVTime::TIME, 9);
263 :
264 0 : Quantity dLat = shape().refDirectionErrorLat();
265 0 : dLat.convert("rad");
266 :
267 0 : Quantity dLong = shape().refDirectionErrorLong();
268 0 : dLong.convert("rad");
269 :
270 : // Add error estimates to ra/dec strings if an error is given (either >0)
271 :
272 0 : uInt precision = 1;
273 0 : if ( dLong.getValue() != 0 || dLat.getValue() != 0 ) {
274 0 : dLong.convert("s");
275 0 : dLat.convert("arcsec");
276 0 : Double drasec = roundDouble(dLong.getValue());
277 0 : Double ddecarcsec = roundDouble(dLat.getValue());
278 0 : Vector<Double> dravec(2), ddecvec(2);
279 0 : dravec.set(drasec);
280 0 : ddecvec.set(ddecarcsec);
281 0 : precision = precisionForValueErrorPairs(dravec,ddecvec);
282 0 : longString = MVTime(longitude).string(MVTime::TIME, 6+precision);
283 0 : latString = MVAngle(lat).string(MVAngle::ANGLE, 6+precision);
284 : }
285 0 : std::pair<String, String> labels = _axisLabels(dc);
286 0 : position << "Position ---" << endl;
287 0 : position << " --- " << labels.first << " " << longString;
288 0 : if (dLong.getValue() == 0) {
289 0 : position << " (fixed)" << endl;
290 : }
291 : else {
292 0 : Quantity timeError = longErrOnGreatCircle ? dLong/cos(lat) : dLong;
293 0 : position << " +/- " << std::fixed
294 0 : << std::setprecision(precision) << timeError << " ("
295 0 : << dLong.getValue("arcsec") << " arcsec"
296 : << (longErrOnGreatCircle ? " along great circle" : "")
297 0 : << ")" << endl;
298 : }
299 0 : position << " --- " << labels.second << " " << latString;
300 0 : if (dLat.getValue() == 0) {
301 0 : position << " (fixed)" << endl;
302 : }
303 : else {
304 0 : position << " +/- " << dLat << endl;
305 : }
306 0 : if (dc) {
307 0 : const Vector<String> units = dc->worldAxisUnits();
308 0 : Vector<Double> world(dc->nWorldAxes(), 0);
309 0 : world[0] = longitude.getValue(units[0]);
310 0 : world[1] = lat.getValue(units[1]);
311 : // TODO do the pixel computations in another method
312 0 : pixelCoords.reset(new Vector<Double>());
313 0 : if (dc->toPixel(*pixelCoords, world)) {
314 0 : Vector<Double> increment = dc->increment();
315 0 : Double longPixErr = dLong.getValue() == 0
316 0 : ? 0 : abs(dLong.getValue(units[0])/increment[0]);
317 0 : Double latPixErr = dLat.getValue() == 0
318 0 : ? 0 : abs(dLat.getValue(units[1])/increment[1]);
319 0 : Vector<Double> longPix(2), latPix(2);
320 0 : longPix.set(roundDouble(longPixErr));
321 0 : latPix.set(roundDouble(latPixErr));
322 0 : precision = precisionForValueErrorPairs(longPix, latPix);
323 0 : position << std::fixed << std::setprecision(precision);
324 0 : position << " --- " << labels.first << " " << (*pixelCoords)[0];
325 0 : if (dLong.getValue() == 0) {
326 0 : position << " (fixed)" << endl;
327 : }
328 : else {
329 0 : position << " +/- " << longPixErr << " pixels" << endl;
330 : }
331 0 : position << " --- " << labels.second << " " << (*pixelCoords)[1];
332 0 : if (dLat.getValue() == 0) {
333 0 : position << " (fixed)" << endl;
334 : }
335 : else {
336 0 : position << " +/- " << latPixErr << " pixels" << endl;
337 : }
338 : }
339 : else {
340 0 : position << "unable to determine position in pixels:" << dc->errorMessage() << endl;
341 : }
342 : }
343 0 : return position.str();
344 : }
345 :
346 0 : std::pair<String, String> SkyComponent::_axisLabels(
347 : const DirectionCoordinate *const &dc
348 : ) {
349 0 : std::pair<String, String> labels;
350 0 : if (dc) {
351 0 : Vector<String> names = dc->worldAxisNames();
352 0 : for (uInt i=0; i<2; i++) {
353 0 : String label;
354 0 : names[i].downcase();
355 0 : if (names[i] == "right ascension") {
356 0 : label = "ra:";
357 : }
358 0 : else if (names[i] == "declination") {
359 0 : label = "dec:";
360 : }
361 0 : else if (names[i] == "longitude") {
362 0 : label = "long:";
363 : }
364 0 : else if (names[i] == "latitude") {
365 0 : label = "lat:";
366 : }
367 : else {
368 0 : label = names[i] + ":";
369 : }
370 0 : if (i == 0) {
371 0 : labels.first = label;
372 : }
373 : else {
374 0 : labels.second = label;
375 : }
376 : }
377 : }
378 : else {
379 0 : labels.first = "long:";
380 0 : labels.second = "lat: ";
381 : }
382 : typedef std::string::size_type size_type;
383 0 : size_type f = labels.first.size();
384 0 : size_type s = labels.second.size();
385 0 : if (f > s) {
386 0 : size_type d = f - s;
387 0 : for (size_type i=0; i<d ;i++) {
388 0 : labels.second += " ";
389 : }
390 : }
391 0 : else if (f < s) {
392 0 : size_type d = s - f;
393 0 : for (size_type i=0; i<d ;i++) {
394 0 : labels.first += " ";
395 : }
396 : }
397 0 : return labels;
398 : }
399 :
400 : }
401 :
|
