Line data Source code
1 : //# SkyCompRep.cc: this defines SkyCompRep
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002
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: SkyCompRep.cc 21292 2012-11-28 14:58:19Z gervandiepen $
27 :
28 : #include <components/ComponentModels/ComponentShape.h>
29 : #include <components/ComponentModels/ComponentType.h>
30 : #include <components/ComponentModels/ConstantSpectrum.h>
31 : #include <components/ComponentModels/PointShape.h>
32 : #include <casacore/scimath/Mathematics/GaussianBeam.h>
33 : #include <components/ComponentModels/GaussianShape.h>
34 : #include <components/ComponentModels/DiskShape.h>
35 : #include <components/ComponentModels/LimbDarkenedDiskShape.h>
36 : #include <components/ComponentModels/SkyCompRep.h>
37 : #include <casacore/casa/Arrays/ArrayMath.h>
38 : #include <casacore/casa/Arrays/Cube.h>
39 : #include <casacore/casa/Arrays/Matrix.h>
40 : #include <casacore/casa/Arrays/MatrixIter.h>
41 : #include <casacore/casa/Arrays/Vector.h>
42 : #include <casacore/casa/Containers/Record.h>
43 : #include <casacore/casa/Containers/RecordFieldId.h>
44 : #include <casacore/casa/Containers/RecordInterface.h>
45 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
46 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
47 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
48 : #include <casacore/casa/Exceptions/Error.h>
49 : #include <casacore/casa/Logging/LogIO.h>
50 : #include <casacore/casa/Logging/LogOrigin.h>
51 : #include <casacore/casa/BasicMath/Math.h>
52 : #include <casacore/casa/BasicSL/Complex.h>
53 : #include <casacore/measures/Measures/MDirection.h>
54 : #include <casacore/measures/Measures/MFrequency.h>
55 : #include <casacore/measures/Measures/Stokes.h>
56 : #include <casacore/casa/Quanta/MVAngle.h>
57 : #include <casacore/casa/Quanta/Quantum.h>
58 : #include <casacore/casa/Quanta/Unit.h>
59 : #include <casacore/casa/Quanta/UnitMap.h>
60 : #include <casacore/casa/Utilities/Assert.h>
61 : #include <casacore/casa/Utilities/DataType.h>
62 : #include <components/ComponentModels/SpectralModel.h>
63 : #include <iostream>
64 :
65 : using namespace casacore;
66 : namespace casa { //# NAMESPACE CASA - BEGIN
67 :
68 51202 : SkyCompRep::SkyCompRep()
69 51202 : :itsShapePtr(new PointShape),
70 51202 : itsSpectrumPtr(new ConstantSpectrum),
71 51202 : itsFlux(),
72 51202 : itsLabel(),
73 102404 : itsOptParms()
74 :
75 : {
76 51202 : AlwaysAssert(ok(), AipsError);
77 :
78 51202 : }
79 :
80 2 : SkyCompRep::SkyCompRep(const ComponentType::Shape& shape)
81 2 : :itsShapePtr(ComponentType::construct(shape)),
82 2 : itsSpectrumPtr(new ConstantSpectrum),
83 2 : itsFlux(),
84 2 : itsLabel(),
85 4 : itsOptParms()
86 : {
87 2 : AlwaysAssert(ok(), AipsError);
88 2 : }
89 :
90 1 : SkyCompRep::SkyCompRep(const ComponentType::Shape& shape,
91 1 : const ComponentType::SpectralShape& spectrum)
92 1 : :itsShapePtr(ComponentType::construct(shape)),
93 1 : itsSpectrumPtr(ComponentType::construct(spectrum)),
94 1 : itsFlux(),
95 1 : itsLabel(),
96 2 : itsOptParms()
97 : {
98 1 : AlwaysAssert(ok(), AipsError);
99 1 : }
100 :
101 1 : SkyCompRep::SkyCompRep(const Flux<Double>& flux,
102 : const ComponentShape& shape,
103 1 : const SpectralModel& spectrum)
104 1 : :itsShapePtr(shape.clone()),
105 1 : itsSpectrumPtr(spectrum.clone()),
106 1 : itsFlux(flux.copy()),
107 1 : itsLabel(),
108 2 : itsOptParms()
109 : {
110 1 : AlwaysAssert(ok(), AipsError);
111 1 : }
112 :
113 0 : SkyCompRep::SkyCompRep(const SkyCompRep& other)
114 : :SkyCompBase(),
115 0 : itsShapePtr(other.itsShapePtr->clone()),
116 0 : itsSpectrumPtr(other.itsSpectrumPtr->clone()),
117 0 : itsFlux(other.itsFlux.copy()),
118 0 : itsLabel(other.itsLabel),
119 0 : itsOptParms(other.itsOptParms)
120 : {
121 0 : AlwaysAssert(ok(), AipsError);
122 0 : }
123 :
124 102410 : SkyCompRep::~SkyCompRep() {}
125 :
126 0 : SkyCompRep& SkyCompRep::operator=(const SkyCompRep& other) {
127 0 : if (this != &other) {
128 0 : itsShapePtr = other.itsShapePtr->clone();
129 0 : itsSpectrumPtr = other.itsSpectrumPtr->clone();
130 0 : itsFlux = other.itsFlux.copy();
131 0 : itsLabel = other.itsLabel;
132 0 : itsOptParms = other.itsOptParms;
133 : }
134 0 : AlwaysAssert(ok(), AipsError);
135 0 : return *this;
136 : }
137 :
138 0 : const Flux<Double>& SkyCompRep::flux() const {
139 0 : return itsFlux;
140 : }
141 :
142 51203 : Flux<Double>& SkyCompRep::flux() {
143 51203 : return itsFlux;
144 : }
145 :
146 0 : const ComponentShape& SkyCompRep::shape() const {
147 0 : DebugAssert(ok(), AipsError);
148 0 : return *itsShapePtr;
149 : }
150 :
151 52007 : ComponentShape& SkyCompRep::shape() {
152 52007 : return *itsShapePtr;
153 : }
154 :
155 0 : void SkyCompRep::setShape(const ComponentShape& newShape) {
156 0 : itsShapePtr = newShape.clone();
157 0 : }
158 :
159 0 : SpectralModel& SkyCompRep::spectrum() {
160 0 : return *itsSpectrumPtr;
161 : }
162 :
163 0 : const SpectralModel& SkyCompRep::spectrum() const {
164 0 : return *itsSpectrumPtr;
165 : }
166 :
167 0 : void SkyCompRep::setSpectrum(const SpectralModel& newSpectrum) {
168 0 : itsSpectrumPtr = newSpectrum.clone();
169 0 : }
170 :
171 0 : String& SkyCompRep::label() {
172 0 : return itsLabel;
173 : }
174 :
175 0 : const String& SkyCompRep::label() const {
176 0 : return itsLabel;
177 : }
178 :
179 0 : Vector<Double>& SkyCompRep::optionalParameters() {
180 0 : return itsOptParms;
181 : }
182 :
183 0 : const Vector<Double>& SkyCompRep::optionalParameters() const {
184 0 : return itsOptParms;
185 : }
186 :
187 0 : Bool SkyCompRep::isPhysical() const {
188 0 : Flux<Double> compFlux = flux().copy();
189 0 : compFlux.convertPol(ComponentType::STOKES);
190 0 : const Vector<DComplex>& iquv = compFlux.value();
191 0 : const DComplex& i = iquv(0);
192 0 : const DComplex& q = iquv(1);
193 0 : const DComplex& u = iquv(2);
194 0 : const DComplex& v = iquv(3);
195 0 : if (!nearAbs(i.imag(), 0.0) ||
196 0 : !nearAbs(q.imag(), 0.0) ||
197 0 : !nearAbs(u.imag(), 0.0) ||
198 0 : !nearAbs(v.imag(), 0.0) ) {
199 0 : return false;
200 : }
201 0 : if (square(i.real()) <
202 0 : square(q.real()) + square(u.real()) + square(v.real()) ) {
203 0 : return false;
204 : }
205 0 : return true;
206 0 : }
207 :
208 0 : Flux<Double> SkyCompRep::sample(const MDirection& direction,
209 : const MVAngle& pixelLatSize,
210 : const MVAngle& pixelLongSize,
211 : const MFrequency& centerFrequency) const {
212 0 : Double scale = itsShapePtr->sample(direction, pixelLatSize, pixelLongSize);
213 : //scale *= itsSpectrumPtr->sample(centerFrequency);
214 0 : Flux<Double> flux = itsFlux.copy();
215 0 : Vector<Double> iquv(4);
216 0 : flux.value(iquv);
217 0 : itsSpectrumPtr->sampleStokes(centerFrequency, iquv);
218 0 : flux.setValue(iquv);
219 0 : flux.convertPol(itsFlux.pol());
220 0 : flux.scaleValue(scale, scale, scale, scale);
221 0 : return flux;
222 0 : }
223 :
224 0 : void SkyCompRep::sample(Cube<Double>& samples, const Unit& reqUnit,
225 : const Vector<MVDirection>& directions,
226 : const MeasRef<MDirection>& dirRef,
227 : const MVAngle& pixelLatSize,
228 : const MVAngle& pixelLongSize,
229 : const Vector<MVFrequency>& frequencies,
230 : const MeasRef<MFrequency>& freqRef) const {
231 0 : const uInt nDirSamples = directions.nelements();
232 0 : const uInt nFreqSamples = frequencies.nelements();
233 0 : Flux<Double> f = itsFlux.copy();
234 0 : f.convertUnit(reqUnit);
235 0 : Vector<Double> fluxVal(4);
236 0 : f.value(fluxVal);
237 : /*const Double i = fluxVal(0);
238 : const Double q = fluxVal(1);
239 : const Double u = fluxVal(2);
240 : const Double v = fluxVal(3);
241 : */
242 :
243 0 : Vector<Double> dirScales(nDirSamples);
244 0 : itsShapePtr->sample(dirScales, directions, dirRef,
245 : pixelLatSize, pixelLongSize);
246 :
247 0 : Matrix<Double> freqIQUV(nFreqSamples, 4);
248 0 : for (uInt i=0; i<nFreqSamples; ++i) {
249 0 : freqIQUV.row(i) = fluxVal.copy();
250 : }
251 :
252 : //itsSpectrumPtr->sample(freqScales, frequencies, freqRef);
253 0 : itsSpectrumPtr->sampleStokes(freqIQUV, frequencies, freqRef);
254 0 : for (uInt d=0; d<nDirSamples; ++d) {
255 0 : const Double thisDirScale = dirScales(d);
256 0 : if (thisDirScale != 0) {
257 0 : for (uInt f=0; f<nFreqSamples; ++f) {
258 0 : for (uInt stok=0; stok <4; ++stok){
259 0 : samples(stok, d, f) += thisDirScale * freqIQUV(f, stok);
260 : }
261 : }
262 : }
263 : }
264 0 : }
265 :
266 0 : Flux<Double> SkyCompRep::visibility(const Vector<Double>& uvw,
267 : const Double& frequency) const {
268 0 : Flux<Double> flux = itsFlux.copy();
269 0 : Double scale = itsShapePtr->visibility(uvw, frequency).real();
270 0 : MFrequency freq(Quantity(frequency, "Hz"));
271 : //scale *= itsSpectrumPtr->sample(freq);
272 0 : Vector<Double> iquv(4);
273 0 : flux.value(iquv);
274 0 : itsSpectrumPtr->sampleStokes(freq, iquv);
275 0 : flux.setValue(iquv);
276 0 : flux.convertPol(itsFlux.pol());
277 0 : flux.scaleValue(scale, scale, scale, scale);
278 0 : return flux;
279 0 : }
280 :
281 25600 : void SkyCompRep::visibility(Cube<DComplex>& visibilities,
282 : const Matrix<Double>& uvws,
283 : const Vector<Double>& frequencies) const {
284 25600 : const uInt nFreq = frequencies.nelements();
285 25600 : const uInt nVis = uvws.ncolumn();
286 :
287 25600 : Vector<Double> uvw(3);
288 : //Block<DComplex> flux(4);
289 25600 : Vector<Double> iquv(4);
290 25600 : Flux<Double> flux=itsFlux.copy();
291 25600 : flux.value(iquv);
292 : /*for (uInt p = 0; p < 4; p++) {
293 : flux[p] = itsFlux.value(p);
294 : }
295 : */
296 :
297 25600 : Matrix<Double> fIQUV(frequencies.size(), 4);
298 25600 : Vector<MVFrequency> mvFreq(frequencies.size());
299 2585600 : for (uInt f = 0; f < nFreq; f++) {
300 2560000 : fIQUV.row(f) = iquv.copy();
301 2560000 : mvFreq(f)=MVFrequency(frequencies(f));
302 : }
303 :
304 : // It's not clear how we would get the right information down here.
305 : // In any event, it's probably not a pressing concern for most cases.
306 :
307 : //Indeed ...write something complex and make believe that
308 : // that transformation from different frames can happen and let it bite when some
309 : // poor sucker try to use it
310 : //At least for now making it that the frequency is expected to be in the frame of
311 : // the component
312 25600 : MeasRef<MFrequency> measRef(itsSpectrumPtr->refFrequency().getRef());
313 : //itsSpectrumPtr->sample(fscale, mvFreq, measRef);
314 25600 : itsSpectrumPtr->sampleStokes(fIQUV, mvFreq, measRef);
315 25600 : Vector<Flux<Double> > stokesFlux(nFreq);
316 2585600 : for (uInt f=0; f < nFreq; ++f){
317 2560000 : stokesFlux(f) = Flux<Double>(fIQUV.row(f));
318 2560000 : stokesFlux(f).convertPol(itsFlux.pol());
319 : }
320 25600 : Matrix<DComplex> scales(nVis, nFreq);
321 25600 : itsShapePtr->visibility(scales, uvws, frequencies);
322 : /*Matrix<DComplex> scales2(nFreq, nVis);
323 : for(uInt k=0; k < nFreq; ++k ){
324 : for (uInt j=0; j < nVis; ++j){
325 : scales2(k,j)=scales(j,k)*fscale(k);
326 : }
327 : }
328 : for (uInt p = 0; p < 4; ++p) {
329 : visibilities.yzPlane(p) = flux[p]*scales2;
330 : }
331 : */
332 306400 : for (uInt v = 0; v < nVis; ++v) {
333 28360800 : for (uInt f=0; f < nFreq; ++f ){
334 140400000 : for (uInt p = 0; p < 4; ++p) {
335 112320000 : visibilities(p, f, v)=scales(v, f)*stokesFlux[f].value(p);
336 : }
337 : }
338 : }
339 :
340 :
341 : /*
342 : for (uInt v = 0; v < nVis; v++) {
343 : uvw=uvws.column(v);
344 : Matrix<Complex> plane;
345 : plane.reference(visibilities.xyPlane(v));
346 : // Scale by the specified frequency behaviour
347 : for (uInt f = 0; f < nFreq; f++) {
348 : Double scale = itsShapePtr->visibility(uvw, frequencies(f)).real();
349 : scale *= fscale[f];
350 : for (uInt p = 0; p < 4; p++) {
351 : visibilities(p, f, v) = flux[p] * scale;
352 : }
353 : }
354 : }
355 : */
356 :
357 :
358 25600 : }
359 :
360 25600 : Bool SkyCompRep::fromRecord(String& errorMessage,
361 : const RecordInterface& record) {
362 : {
363 25600 : const String fluxString("flux");
364 25600 : if (record.isDefined(fluxString)) {
365 25600 : const RecordFieldId flux(fluxString);
366 25600 : if (record.dataType(flux) != TpRecord) {
367 0 : errorMessage += "The 'flux' field must be a record\n";
368 0 : return false;
369 : }
370 25600 : const Record& fluxRec = record.asRecord(flux);
371 25600 : if (!itsFlux.fromRecord(errorMessage, fluxRec)) {
372 0 : errorMessage += "Problem parsing the 'flux' field\n";
373 0 : return false;
374 : }
375 25600 : } else {
376 0 : LogIO logErr(LogOrigin("SkyCompRep", "fromRecord()"));
377 : logErr << LogIO::WARN
378 : << "The component does not have a 'flux' field." << endl
379 : << "The default is 1.0 Jy in I and 0.0 in Q, U & V"
380 0 : << LogIO::POST;
381 0 : itsFlux = Flux<Double>(1);
382 0 : }
383 25600 : }
384 : {
385 25600 : const String shapeString("shape");
386 25600 : if (record.isDefined(shapeString)) {
387 25600 : const RecordFieldId shape(shapeString);
388 25600 : if (record.dataType(shape) != TpRecord) {
389 0 : errorMessage += "\nThe 'shape' field must be a record";
390 0 : return false;
391 : }
392 25600 : const Record& shapeRec = record.asRecord(shape);
393 : const ComponentType::Shape recType =
394 25600 : ComponentShape::getType(errorMessage, shapeRec);
395 25600 : if (recType >= ComponentType::UNKNOWN_SHAPE) {
396 0 : errorMessage += String("Cannot create a component with a '" +
397 0 : ComponentType::name(recType) + "' shape\n");
398 0 : return false;
399 : }
400 25600 : if (recType != itsShapePtr->type()) {
401 0 : ComponentShape* newShape = ComponentType::construct(recType);
402 0 : AlwaysAssert(newShape != 0, AipsError);
403 0 : setShape(*newShape);
404 0 : delete newShape;
405 : }
406 25600 : if (!itsShapePtr->fromRecord(errorMessage, shapeRec)) {
407 0 : errorMessage += "Problem parsing the 'shape' field\n";
408 0 : return false;
409 : }
410 25600 : } else {
411 0 : LogIO logErr(LogOrigin("SkyCompRep", "fromRecord()"));
412 : logErr << LogIO::WARN
413 : << "The component does not have a 'shape' field." << endl
414 : << "The default is a point component at the J2000 north pole"
415 0 : << LogIO::POST;
416 0 : const Unit deg("deg");
417 0 : itsShapePtr = new PointShape(MDirection(Quantum<Double>(0.0, deg),
418 0 : Quantum<Double>(90.0, deg),
419 0 : MDirection::J2000));
420 0 : }
421 25600 : }
422 : {
423 25600 : const String spectrumString("spectrum");
424 25600 : if (record.isDefined(spectrumString)) {
425 25600 : const RecordFieldId spectrum(spectrumString);
426 25600 : if (record.dataType(spectrum) != TpRecord) {
427 0 : errorMessage += "\nThe 'spectrum' field must be a record";
428 0 : return false;
429 : }
430 25600 : const Record& spectrumRec = record.asRecord(spectrum);
431 : const ComponentType::SpectralShape recType =
432 25600 : SpectralModel::getType(errorMessage, spectrumRec);
433 25600 : if (recType >= ComponentType::UNKNOWN_SPECTRAL_SHAPE) {
434 0 : errorMessage += String("Cannot create a component with a '" +
435 0 : ComponentType::name(recType) + "' spectrum\n");
436 0 : return false;
437 : }
438 25600 : if (recType != itsSpectrumPtr->type()) {
439 0 : SpectralModel* newSpectrum = ComponentType::construct(recType);
440 0 : AlwaysAssert(newSpectrum != 0, AipsError);
441 0 : setSpectrum(*newSpectrum);
442 0 : delete newSpectrum;
443 : }
444 25600 : if (!itsSpectrumPtr->fromRecord(errorMessage, spectrumRec)) {
445 0 : return false;
446 : }
447 25600 : } else {
448 0 : LogIO logErr(LogOrigin("SkyCompRep", "fromRecord()"));
449 : logErr << LogIO::WARN
450 : << "The component does not have a 'spectrum' field." << endl
451 : << "The default is a constant spectrum"
452 0 : << LogIO::POST;
453 0 : itsSpectrumPtr = new ConstantSpectrum;
454 0 : }
455 25600 : }
456 : {
457 25600 : const String labelString("label");
458 25600 : if (record.isDefined(labelString)) {
459 25600 : const RecordFieldId label(labelString);
460 25600 : if (record.dataType(label) != TpString) {
461 0 : errorMessage += "\nThe 'label' field must be a string";
462 0 : return false;
463 : }
464 25600 : if (record.shape(label) != IPosition(1,1)) {
465 0 : errorMessage += "\nThe 'label' field must have only 1 element";
466 0 : return false;
467 : }
468 25600 : itsLabel = record.asString(label);
469 25600 : }
470 25600 : }
471 : {
472 25600 : const String optParmString("optionalParameters");
473 25600 : if (record.isDefined(optParmString)) {
474 0 : const RecordFieldId optionalParameters(optParmString);
475 0 : if (record.dataType(optionalParameters) != TpArrayDouble) {
476 0 : errorMessage += "\nThe 'optionalParameters' field must be a Double Array";
477 0 : return false;
478 : }
479 0 : itsOptParms = record.asArrayDouble(optionalParameters);
480 0 : }
481 25600 : }
482 25600 : return true;
483 : }
484 :
485 800 : Bool SkyCompRep::toRecord(String& errorMessage,
486 : RecordInterface& record) const {
487 : {
488 800 : Record fluxRec;
489 800 : if (!itsFlux.toRecord(errorMessage, fluxRec)) {
490 0 : return false;
491 : }
492 800 : record.defineRecord(RecordFieldId("flux"), fluxRec);
493 800 : }
494 : {
495 800 : Record shapeRec;
496 800 : if (!itsShapePtr->toRecord(errorMessage, shapeRec)) {
497 0 : return false;
498 : }
499 800 : record.defineRecord(RecordFieldId("shape"), shapeRec);
500 800 : }
501 : {
502 800 : Record spectrumRec;
503 800 : if (!itsSpectrumPtr->toRecord(errorMessage, spectrumRec)) {
504 0 : return false;
505 : }
506 800 : record.defineRecord(RecordFieldId("spectrum"), spectrumRec);
507 800 : }
508 : {
509 800 : Record optParmRec;
510 800 : if (!itsOptParms.empty()) {
511 0 : record.define(RecordFieldId("optionalParameters"), itsOptParms);
512 : }
513 800 : }
514 800 : record.define(RecordFieldId("label"), itsLabel);
515 800 : return true;
516 : }
517 :
518 51206 : Bool SkyCompRep::ok() const {
519 51206 : if (itsShapePtr.null()) {
520 0 : LogIO logErr(LogOrigin("SkyCompRep", "ok()"));
521 : logErr << LogIO::SEVERE << "Shape pointer is null"
522 0 : << LogIO::POST;
523 0 : return false;
524 0 : }
525 51206 : if (itsShapePtr->ok() == false) {
526 0 : LogIO logErr(LogOrigin("SkyCompRep", "ok()"));
527 : logErr << LogIO::SEVERE << "The component shape is not ok"
528 0 : << LogIO::POST;
529 0 : return false;
530 0 : }
531 51206 : if (itsSpectrumPtr.null()) {
532 0 : LogIO logErr(LogOrigin("SkyCompRep", "ok()"));
533 : logErr << LogIO::SEVERE << "Spectrum pointer is null"
534 0 : << LogIO::POST;
535 0 : return false;
536 0 : }
537 51206 : if (itsSpectrumPtr->ok() == false) {
538 0 : LogIO logErr(LogOrigin("SkyCompRep", "ok()"));
539 : logErr << LogIO::SEVERE << "The component spectrum is not ok"
540 0 : << LogIO::POST;
541 0 : return false;
542 0 : }
543 51206 : return true;
544 : }
545 :
546 :
547 :
548 0 : void SkyCompRep::fromPixel (
549 : Double& facToJy, const Vector<Double>& parameters,
550 : const Unit& brightnessUnitIn,
551 : const GaussianBeam& restoringBeam,
552 : const CoordinateSystem& cSys,
553 : ComponentType::Shape componentShape,
554 : Stokes::StokesTypes stokes
555 : ) {
556 : //
557 : // pars(0) = Flux Jy
558 : // pars(1) = x cen abs pix
559 : // pars(2) = y cen abs pix
560 : // pars(3) = major pix
561 : // pars(4) = minor pix
562 : // pars(5) = pa radians
563 : //
564 0 : LogIO os(LogOrigin("SkyCompRep", "fromPixel()"));
565 :
566 0 : ThrowIf(
567 : ! cSys.hasDirectionCoordinate(),
568 : "CoordinateSystem does not contain a DirectionCoordinate"
569 : );
570 0 : const DirectionCoordinate& dirCoord = cSys.directionCoordinate();
571 : //
572 : // We need to find the ratio that converts the input peak brightness
573 : // from whatevers/per whatever to Jy per whatever. E.g. mJy/beam to Jy/beam.
574 : // This ratio is passed back (for scaling errors) and is needed regardless of
575 : // the component type.
576 :
577 0 : facToJy = SkyCompRep::convertToJy (brightnessUnitIn);
578 :
579 : // Now proceed with type dependent conversions
580 0 : if (componentShape==ComponentType::POINT) {
581 0 : ThrowIf(
582 : parameters.nelements() != 3,
583 : "Wrong number of parameters for Point shape"
584 : );
585 0 : Vector<Double> pars(2);
586 0 : pars(0) = parameters(1);
587 0 : pars(1) = parameters(2);
588 0 : PointShape pointShape;
589 0 : pointShape.fromPixel(pars, dirCoord);
590 0 : setShape(pointShape);
591 : //
592 0 : Quantum<Double> value(parameters(0)*facToJy, Unit("Jy"));
593 0 : itsFlux.setUnit(Unit("Jy"));
594 0 : itsFlux.setValue (value, stokes);
595 0 : } else if (componentShape==ComponentType::GAUSSIAN || componentShape==ComponentType::DISK) {
596 0 : ThrowIf(
597 : parameters.nelements() != 6,
598 : "Wrong number of parameters for Gaussian or Point shape"
599 : );
600 :
601 : // Do x,y,major,minor,pa
602 0 : Vector<Double> pars(5);
603 0 : for (uInt i=0; i<5; i++) {
604 0 : pars(i) = parameters(i+1);
605 : }
606 0 : Quantum<Double> majorAxis, minorAxis, pa;
607 0 : if (componentShape==ComponentType::GAUSSIAN) {
608 0 : GaussianShape shp;
609 0 : shp.fromPixel (pars, dirCoord);
610 0 : setShape(shp);
611 0 : majorAxis = shp.majorAxis();
612 0 : minorAxis = shp.minorAxis();
613 0 : pa = shp.positionAngle();
614 0 : } else {
615 0 : DiskShape shp;
616 0 : shp.fromPixel (pars, dirCoord);
617 0 : setShape(shp);
618 0 : majorAxis = shp.majorAxis();
619 0 : minorAxis = shp.minorAxis();
620 0 : pa = shp.positionAngle();
621 0 : }
622 0 : Quantum<Double> peakFlux(parameters(0), brightnessUnitIn);
623 : Quantum<Double> integralFlux =
624 : SkyCompRep::peakToIntegralFlux (dirCoord, componentShape, peakFlux,
625 0 : majorAxis, minorAxis, restoringBeam);
626 0 : itsFlux.setUnit(integralFlux.getFullUnit());
627 0 : itsFlux.setValue (integralFlux, stokes);
628 0 : }
629 :
630 :
631 : // Spectrum; assumed constant !
632 0 : ConstantSpectrum constSpec;
633 0 : if (cSys.hasSpectralAxis()) {
634 0 : SpectralCoordinate specCoord = cSys.spectralCoordinate();
635 0 : MFrequency mFreq;
636 0 : ThrowIf(
637 : ! specCoord.toWorld(mFreq, 0.0),
638 : "SpectralCoordinate conversion failed because "
639 : + specCoord.errorMessage()
640 : );
641 0 : constSpec.setRefFrequency(mFreq);
642 0 : }
643 0 : setSpectrum(constSpec);
644 0 : }
645 :
646 0 : Vector<Double> SkyCompRep::toPixel (const Unit& brightnessUnitOut,
647 : const GaussianBeam& restoringBeam,
648 : const CoordinateSystem& cSys,
649 : Stokes::StokesTypes stokes) const
650 : //
651 : // pars(0) = FLux Jy
652 : // pars(1) = x cen abs pix
653 : // pars(2) = y cen abs pix
654 : // pars(3) = major pix
655 : // pars(4) = minor pix
656 : // pars(5) = pa radians
657 : //
658 : {
659 0 : LogIO os(LogOrigin("SkyCompRep", "toPixel()"));
660 :
661 : // Find DirectionCoordinate
662 :
663 0 : Int dirCoordinate = cSys.findCoordinate(Coordinate::DIRECTION);
664 0 : if (dirCoordinate==-1) {
665 0 : os << "CoordinateSystem does not contain a DirectionCoordinate" << LogIO::EXCEPTION;
666 : }
667 0 : const DirectionCoordinate& dirCoord = cSys.directionCoordinate(dirCoordinate);
668 :
669 : // Do x,y, and possibly major,minor,pa (disk/gaussian)
670 :
671 0 : const ComponentShape& componentShape = shape();
672 0 : Vector<Double> pars = componentShape.toPixel(dirCoord);
673 0 : Vector<Double> parameters(1+pars.nelements());
674 0 : for (uInt i=0; i<pars.nelements(); i++) parameters[i+1] = pars[i];
675 :
676 : // Now do Flux
677 :
678 0 : ComponentType::Shape type = componentShape.type();
679 0 : if (type==ComponentType::POINT) {
680 0 : Flux<Double> f = flux();
681 0 : Quantum<Double> fluxPeak = f.value (stokes, true);
682 0 : parameters(0) = fluxPeak.getValue(); // Jy
683 0 : } else if (type==ComponentType::GAUSSIAN || type==ComponentType::DISK) {
684 :
685 : // Define /beam and /pixel units.
686 :
687 0 : Bool integralInJy = true;
688 : Unit brightnessUnits = defineBrightnessUnits(os, brightnessUnitOut, dirCoord,
689 0 : restoringBeam, integralInJy);
690 :
691 : // Get Flux (integral) for particular Stokes.
692 :
693 0 : Flux<Double> f = flux();
694 0 : Quantum<Double> fluxIntegral = f.value (stokes, true);
695 :
696 : // Find peak value. Because we have defined /beam and /pixel units
697 : // above we can use Quanta mathematics to get the answer we want
698 :
699 : Double fac;
700 0 : if (type==ComponentType::GAUSSIAN) {
701 0 : fac = C::pi / 4.0 / log(2.0);
702 0 : } else if (type==ComponentType::DISK) {
703 0 : fac = C::pi;
704 : } else {
705 0 : fac = 1.0;
706 : }
707 : //
708 0 : const TwoSidedShape& ts = dynamic_cast<const TwoSidedShape&>(componentShape);
709 0 : Quantum<Double> major2 = ts.majorAxis();
710 0 : major2.convert(Unit("rad"));
711 0 : Quantum<Double> minor2 = ts.minorAxis();
712 0 : minor2.convert(Unit("rad"));
713 : //
714 0 : Quantum<Double> tmp = major2 * minor2;
715 0 : tmp.scale(fac);
716 0 : Quantum<Double> fluxPeak = fluxIntegral / tmp; // /beam or /pixel units divided out here
717 0 : fluxPeak.convert(brightnessUnits);
718 0 : parameters(0) = fluxPeak.getValue();
719 :
720 : // Undefine /beam and /pixel units
721 :
722 0 : SkyCompRep::undefineBrightnessUnits();
723 0 : }
724 0 : return parameters;
725 0 : }
726 :
727 0 : Unit SkyCompRep::defineBrightnessUnits (
728 : LogIO& os, const Unit& brightnessUnitIn,
729 : const DirectionCoordinate& dirCoord,
730 : const GaussianBeam& restoringBeam,
731 : const Bool integralIsJy
732 : ) {
733 : // Define "pixel" units
734 0 : Vector<String> units(2);
735 0 : units.set("rad");
736 0 : DirectionCoordinate dirCoord2 = dirCoord;
737 0 : dirCoord2.setWorldAxisUnits(units);
738 0 : Vector<Double> inc = dirCoord2.increment();
739 :
740 : // Define "beam" units if needed
741 : // ugh this gets confusing because of the static nature of UnitMap. In some
742 : // cases elsewhere beam and pixel are dimensionless
743 0 : if (brightnessUnitIn.getName().contains("beam")) {
744 0 : if (! restoringBeam.isNull()) {
745 : // GaussianBeam rB = restoringBeam;
746 : // Double fac = C::pi / 4.0 / log(2.0) * rB.getMajor().getValue(Unit("rad")) * rB.getMinor().getValue(Unit("rad"));
747 0 : Double fac = restoringBeam.getArea("rad2");
748 0 : UnitMap::putUser("beam", UnitVal(fac,String("rad2")));
749 : }
750 : else {
751 0 : throw AipsError("No beam defined even though the image brightness units are " + brightnessUnitIn.getName());
752 : }
753 : }
754 0 : UnitMap::putUser("pixel", UnitVal(abs(inc(0)*inc(1)), String("rad2")));
755 :
756 : // We must tell the old unit that it has been redefined
757 :
758 : // The need to do things this way rather than a simple
759 : // Unit unitOut = brightnessUnit is unclear to me but it is necessary,
760 : // at least it keeps the unit tests passing.
761 : // dmehring 04may2012
762 0 : Unit unitOut(brightnessUnitIn.getName());
763 :
764 : // Check integral units
765 :
766 0 : if (integralIsJy) {
767 0 : if (unitOut.empty()) {
768 : os << LogIO::WARN << "There are no image brightness units, assuming Jy/pixel"
769 0 : << LogIO::POST;
770 0 : unitOut = Unit("Jy/pixel");
771 : }
772 : else {
773 0 : Quantity t0(1.0, unitOut);
774 0 : Quantity t1(1.0, Unit("rad2"));
775 0 : Quantity t2 = t0 * t1;
776 0 : if (!t2.isConform(Unit("Jy"))) {
777 0 : os << LogIO::WARN << "The image units '" << unitOut.getName() << "' are not consistent " << endl;
778 0 : os << "with Jy when integrated over the sky. Assuming Jy/pixel" << LogIO::POST;
779 0 : unitOut = Unit("Jy/pixel");
780 : }
781 0 : }
782 : }
783 0 : return unitOut;
784 0 : }
785 :
786 0 : void SkyCompRep::undefineBrightnessUnits()
787 : {
788 0 : UnitMap::removeUser("beam");
789 0 : UnitMap::removeUser("pixel");
790 0 : UnitMap::clearCache();
791 0 : }
792 :
793 0 : Quantity SkyCompRep::peakToIntegralFlux (
794 : const DirectionCoordinate& dirCoord,
795 : const ComponentType::Shape componentShape,
796 : const Quantum<Double>& peakFlux,
797 : const Quantum<Double>& majorAxis,
798 : const Quantum<Double>& minorAxis,
799 : const GaussianBeam& restoringBeam) {
800 0 : LogIO os(LogOrigin("SkyCompRep", "peakToIntegralFlux()"));
801 :
802 : // Define /beam and /pixel units.
803 :
804 0 : Unit unitIn = peakFlux.getFullUnit();
805 0 : String unitName = unitIn.getName();
806 0 : Bool integralIsJy = ! (unitName == "Jy/beam.km/s" || unitName == "K");
807 :
808 : Unit brightnessUnit = SkyCompRep::defineBrightnessUnits(
809 : os, unitIn, dirCoord, restoringBeam, integralIsJy
810 0 : );
811 : // Scale to integrated
812 :
813 0 : Quantity tmp(peakFlux.getValue(), brightnessUnit);
814 0 : if (componentShape==ComponentType::GAUSSIAN) {
815 0 : tmp.scale(C::pi / 4.0 / log(2.0));
816 : }
817 0 : else if (componentShape==ComponentType::DISK) {
818 0 : tmp.scale(C::pi);
819 : }
820 : else {
821 0 : SkyCompRep::undefineBrightnessUnits();
822 0 : os << "Unrecognized shape for flux density conversion" << LogIO::EXCEPTION;
823 : }
824 0 : Quantity fluxIntegral;
825 0 : Quantity majorAxis2 = majorAxis;
826 0 : Quantity minorAxis2 = minorAxis;
827 0 : majorAxis2.convert(Unit("rad"));
828 0 : minorAxis2.convert(Unit("rad"));
829 0 : fluxIntegral = tmp * majorAxis2 * minorAxis2;
830 0 : if (fluxIntegral.isConform(Unit("Jy"))) {
831 0 : fluxIntegral.convert("Jy");
832 : }
833 0 : else if (fluxIntegral.isConform(Unit("Jy.m/s"))) {
834 0 : fluxIntegral.convert("Jy.km/s");
835 : }
836 0 : else if (fluxIntegral.isConform(Unit("K.rad2"))) {
837 : // do nothing, units are OK as is
838 : }
839 : else {
840 : os << LogIO::WARN << "Cannot convert units of Flux integral to Jy - will assume Jy"
841 0 : << LogIO::POST;
842 0 : fluxIntegral.setUnit(Unit("Jy"));
843 : }
844 0 : SkyCompRep::undefineBrightnessUnits();
845 0 : return fluxIntegral;
846 0 : }
847 :
848 0 : Quantity SkyCompRep::integralToPeakFlux (
849 : const DirectionCoordinate& dirCoord,
850 : const ComponentType::Shape componentShape,
851 : const Quantity& integralFlux,
852 : const Unit& brightnessUnit,
853 : const Quantity& majorAxis,
854 : const Quantity& minorAxis,
855 : const GaussianBeam& restoringBeam
856 : ) {
857 0 : LogIO os(LogOrigin("SkyCompRep", "integralToPeakFlux()"));
858 0 : Quantity tmp(integralFlux);
859 0 : if (tmp.isConform(Unit("Jy"))) {
860 0 : tmp.convert("Jy");
861 : }
862 0 : else if (tmp.isConform(Unit("Jy.m/s"))) {
863 0 : tmp.convert("Jy.km/s");
864 : }
865 0 : else if (tmp.isConform(Unit("K.rad2"))) {
866 : // do nothing units are OK as is
867 : }
868 : else {
869 0 : os << "Cannot convert units of Flux integral (" + integralFlux.getUnit() + ") to Jy"
870 0 : << LogIO::EXCEPTION;
871 : }
872 : // Define /beam and /pixel units.
873 :
874 0 : String unitName = brightnessUnit.getName();
875 0 : Bool integralIsJy = ! (unitName == "Jy/beam.km/s" || unitName == "K");
876 : Unit unitIn = SkyCompRep::defineBrightnessUnits(
877 : os, brightnessUnit, dirCoord,
878 : restoringBeam, integralIsJy
879 0 : );
880 :
881 : // Convert integral to Jy
882 :
883 : // Now scale
884 :
885 0 : if (componentShape==ComponentType::GAUSSIAN) {
886 0 : tmp.scale(4.0 * log(2.0) / C::pi);
887 : }
888 0 : else if (componentShape==ComponentType::DISK) {
889 0 : tmp.scale(1.0 / C::pi);
890 : }
891 : else {
892 0 : SkyCompRep::undefineBrightnessUnits();
893 0 : os << "Unrecognized shape for flux density conversion" << LogIO::EXCEPTION;
894 : }
895 :
896 : // And divide out shape
897 :
898 0 : Quantity peakFlux;
899 0 : Quantity majorAxis2(majorAxis);
900 0 : Quantity minorAxis2(minorAxis);
901 0 : majorAxis2.convert(Unit("rad"));
902 0 : minorAxis2.convert(Unit("rad"));
903 0 : peakFlux = tmp / majorAxis2 / minorAxis2;
904 0 : peakFlux.convert(unitIn);
905 0 : SkyCompRep::undefineBrightnessUnits();
906 0 : return peakFlux;
907 0 : }
908 :
909 0 : Double SkyCompRep::convertToJy (const Unit& brightnessUnit)
910 : {
911 0 : LogIO os(LogOrigin("SkyCompRep", __FUNCTION__));
912 : // We need to find the ratio that converts the input peak brightness
913 : // from whatevers/per whatever to Jy per whatever. E.g. mJy/beam
914 : // to Jy/beam. This ratio is passed back (for scaling errors) and
915 : // is needed regardless of the component type. So we start by
916 : // Defining /beam and /pixel units to be dimensionless
917 0 : Unit unitIn = brightnessUnit;
918 : // This is not thread safe, but in general anything
919 : // that relies on UnitMap is not because of UnitMap's
920 : // static nature
921 0 : SkyCompRep::undefineBrightnessUnits();
922 :
923 0 : UnitMap::putUser("pixel", UnitVal(1.0,String("")));
924 0 : UnitMap::putUser("beam", UnitVal(1.0,String("")));
925 0 : unitIn = Unit(unitIn.getName()); // Tell system to update this unit
926 0 : Quantum<Double> tmp(1.0, unitIn);
927 0 : Double facToJy = 1.0;
928 0 : if (tmp.isConform(Unit("Jy"))) {
929 0 : Quantum<Double> tmp2(tmp);
930 0 : tmp2.convert("Jy");
931 0 : facToJy = tmp2.getValue() / tmp.getValue();
932 0 : }
933 0 : else if (tmp.isConform(Unit("Jy.m/s"))) {
934 0 : Quantum<Double> tmp2(tmp);
935 0 : tmp2.convert("Jy.km/s");
936 0 : facToJy = tmp2.getValue() / tmp.getValue();
937 0 : }
938 0 : else if (tmp.isConform(Unit("K"))) {
939 0 : Quantum<Double> tmp2(tmp);
940 0 : tmp2.convert("K");
941 0 : facToJy = tmp2.getValue() / tmp.getValue();
942 0 : }
943 : else {
944 : os << LogIO::WARN << "Cannot convert units of brightness to Jy - will assume Jy"
945 0 : << LogIO::POST;
946 0 : facToJy = 1.0;
947 : }
948 :
949 : // Undefine /beam and /pixel
950 :
951 0 : SkyCompRep::undefineBrightnessUnits();
952 0 : return facToJy;
953 0 : }
954 :
955 : } //# NAMESPACE CASA - END
956 :
|
