Line data Source code
1 : //# SimpleComponentFTMachine.cc: Implementation of SimpleComponentFTMachine class
2 : //# Copyright (C) 1997,1998,1999,2000,2001
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$
27 :
28 : #include <synthesis/TransformMachines/SimpleComponentFTMachine.h>
29 : #include <casacore/scimath/Mathematics/RigidVector.h>
30 : #include <components/ComponentModels/ComponentShape.h>
31 : #include <components/ComponentModels/ComponentList.h>
32 : #include <components/ComponentModels/ComponentType.h>
33 : #include <components/ComponentModels/Flux.h>
34 : #include <components/ComponentModels/SkyComponent.h>
35 : #include <components/ComponentModels/SpectralModel.h>
36 : #include <msvis/MSVis/VisBuffer.h>
37 : #include <casacore/ms/MeasurementSets/MSColumns.h>
38 : #include <casacore/casa/Arrays/ArrayMath.h>
39 : #include <casacore/casa/Arrays/ArrayLogical.h>
40 : #include <casacore/casa/Arrays/Cube.h>
41 : #include <casacore/casa/Arrays/Matrix.h>
42 : #include <casacore/casa/Arrays/Vector.h>
43 : #include <casacore/casa/Arrays/IPosition.h>
44 : #include <casacore/casa/BasicSL/Complex.h>
45 : #include <casacore/casa/BasicSL/Constants.h>
46 : #include <casacore/measures/Measures/UVWMachine.h>
47 : #ifdef _OPENMP
48 : #include <omp.h>
49 : #endif
50 :
51 :
52 : using namespace casacore;
53 : namespace casa { //# NAMESPACE CASA - BEGIN
54 :
55 0 : void SimpleComponentFTMachine::get(VisBuffer& vb, SkyComponent& component,
56 : Int row)
57 : {
58 : // If row is -1 then we pass through all rows
59 : uInt startRow, endRow;
60 0 : if (row < 0) {
61 0 : startRow = 0;
62 0 : endRow = vb.nRow() - 1;
63 : } else {
64 0 : startRow = endRow = row;
65 : }
66 0 : const uInt nRow = endRow - startRow + 1;
67 :
68 0 : component.flux().convertUnit(Unit("Jy"));
69 : // Rotate the uvw
70 0 : Matrix<Double> uvw(3, nRow);
71 0 : Vector<Double> dphase(nRow);
72 : {
73 0 : const Vector<RigidVector<Double,3> >& uvwBuff = vb.uvw();
74 0 : for (uInt i = startRow, n = 0; i <= endRow; i++, n++) {
75 0 : const RigidVector<Double,3>& uvwValue = uvwBuff(i);
76 : //NEGATING THE UV to be consitent with the GridFT::get
77 0 : for (Int idim = 0; idim < 2; idim++) {
78 0 : uvw(idim, n) = -uvwValue(idim);
79 : }
80 0 : uvw(2, n) = uvwValue(2);
81 : }
82 0 : rotateUVW(uvw, dphase, vb, component.shape().refDirection());
83 0 : dphase *= -C::_2pi;
84 : }
85 0 : uInt npol=vb.modelVisCube().shape()(0);
86 0 : uInt nChan=vb.modelVisCube().shape()(1);
87 0 : Cube<Complex> modelData;
88 0 : modelData.reference(vb.modelVisCube());
89 0 : Vector<Complex> visibility(4);
90 : //UNUSED: Double phase;
91 : //UNUSED: Double phaseMult;
92 0 : Vector<Double>& frequency = vb.frequency();
93 0 : Vector<Double> invLambda = frequency/C::c;
94 :
95 : // Find the offsets in polarization.
96 0 : Vector<Int> corrType = vb.corrType().copy();
97 : {
98 0 : Int startPol = corrType(0);
99 0 : if((startPol > 4) && (startPol < 9)) {
100 0 : component.flux().convertPol(ComponentType::CIRCULAR);
101 0 : corrType -= 5;
102 : }
103 0 : else if((startPol > 8) && (startPol < 13)) {
104 0 : component.flux().convertPol(ComponentType::LINEAR);
105 0 : corrType -= 9;
106 : }
107 : else {
108 0 : component.flux().convertPol(ComponentType::STOKES);
109 0 : corrType -= startPol;
110 : }
111 : }
112 :
113 :
114 0 : uInt npart=1;
115 : #ifdef _OPENMP
116 0 : npart= numthreads_p <0 ? omp_get_max_threads() : min(numthreads_p, omp_get_max_threads());
117 : #endif
118 0 : if((nRow/npart)==0) npart=1;
119 0 : Block<Cube<DComplex> > dVisp(npart);
120 0 : Vector<uInt> nRowp(npart);
121 0 : Block<Matrix<Double> > uvwp(npart);
122 0 : Block<SkyComponent> compp(npart);
123 0 : Block<Int> startrow(npart);
124 0 : startrow[0]=0;
125 :
126 0 : nRowp.set(nRow/npart);
127 0 : nRowp(npart-1) += nRow%npart;
128 0 : Int sumrow=0;
129 0 : Record lala;
130 0 : String err;
131 0 : component.toRecord(err, lala);
132 0 : for (uInt k=0; k < npart; ++k){
133 0 : dVisp[k].resize(4, nChan, nRowp(k));
134 0 : uvwp[k].resize(3, nRowp(k));
135 0 : uvwp[k]=uvw(IPosition(2,0, sumrow ), IPosition(2, 2, sumrow+nRowp(k)-1));
136 0 : compp[k].fromRecord(err, lala);
137 0 : sumrow+=nRowp(k);
138 0 : startrow[k]=0;
139 0 : for (uInt j=0; j < k; ++j) startrow[k]+=nRowp(j);
140 : }
141 :
142 :
143 : //#pragma omp parallel default(none) firstprivate(npart) shared(frequency,dVisp, uvwp, compp) num_threads(npart)
144 0 : #pragma omp parallel firstprivate(npart) shared(dVisp, uvwp, compp) num_threads(npart)
145 : {
146 :
147 : #pragma omp for
148 : for (Int k=0; k < Int(npart); ++k){
149 : //Cube<DComplex> dVis(4, nChan, nRow);
150 : compp[k].visibility(dVisp[k], uvwp[k], frequency);
151 : }
152 :
153 : }
154 :
155 : // Loop over all rows
156 0 : sumrow=0;
157 :
158 : Bool isCopy;
159 0 : Complex *modData=modelData.getStorage(isCopy);
160 :
161 0 : #pragma omp parallel default(none) firstprivate(npart, npol, nChan, modData, corrType, nRowp, dphase, invLambda) shared(startrow, dVisp) num_threads(npart)
162 : {
163 : #pragma omp for
164 : for (Int k = 0; k < Int(npart); ++k) {
165 :
166 :
167 :
168 : applyPhasor(k, startrow, nRowp, dphase, invLambda, npol, nChan, corrType, dVisp[k], modData);
169 :
170 : }
171 :
172 : }
173 :
174 0 : modelData.putStorage(modData, isCopy);
175 :
176 0 : }
177 :
178 0 : void SimpleComponentFTMachine::applyPhasor(Int part, const Block<Int>& startrow, const Vector<uInt>& nRowp, const Vector<Double>& dphase, const Vector<Double>& invLambda, const Int npol, const Int nchan, const Vector<Int>& corrType, const Cube<DComplex>& dVis, Complex*& modData){
179 :
180 : Int r;
181 : Int rowoff;
182 : Double phaseMult;
183 : Double phase;
184 0 : for (uInt j=0; j< nRowp[part]; ++j){
185 0 : r=startrow[part]+j;
186 0 : rowoff=r*nchan*npol;
187 0 : phaseMult = dphase(r);
188 0 : for (Int chn = 0; chn < nchan; chn++) {
189 0 : phase = phaseMult * invLambda(chn);
190 0 : Complex phasor(cos(phase), sin(phase));
191 0 : for (Int pol=0; pol < npol; pol++) {
192 0 : const DComplex& val = dVis(corrType(pol), chn, j);
193 0 : modData[rowoff+chn*npol+pol] = Complex(val.real(), val.imag()) * conj(phasor);
194 : }
195 : }
196 : }
197 :
198 :
199 :
200 0 : }
201 :
202 :
203 200 : void SimpleComponentFTMachine::get(VisBuffer& vb, const ComponentList& compList,
204 : Int row)
205 : {
206 : // If row is -1 then we pass through all rows
207 : uInt startRow, endRow;
208 200 : if (row < 0) {
209 200 : startRow = 0;
210 200 : endRow = vb.nRow() - 1;
211 : } else {
212 0 : startRow = endRow = row;
213 : }
214 200 : const uInt nRow = endRow - startRow + 1;
215 :
216 : // Rotate the uvw
217 200 : Matrix<Double> uvw0(3, nRow);
218 200 : Matrix<Double> uvw(3, nRow);
219 200 : Vector<Double> dphase(nRow);
220 :
221 200 : const Vector<RigidVector<Double,3> >& uvwBuff = vb.uvw();
222 70400 : for (uInt i = startRow, n = 0; i <= endRow; i++, n++) {
223 70200 : const RigidVector<Double,3>& uvwValue = uvwBuff(i);
224 : //NEGATING the UV to be consitent with GridFT::get
225 210600 : for (Int idim = 0; idim < 2; idim++) {
226 140400 : uvw(idim, n) = -uvwValue(idim);
227 : }
228 70200 : uvw(2, n) = uvwValue(2);
229 : }
230 200 : uvw0=uvw.copy();
231 :
232 200 : uInt ncomponents=compList.nelements();
233 200 : uInt npol=vb.modelVisCube().shape()(0);
234 200 : uInt nChan=vb.modelVisCube().shape()(1);
235 200 : Cube<Complex> modelData;
236 200 : modelData.reference(vb.modelVisCube());
237 200 : modelData=0.0;
238 :
239 200 : Vector<Complex> visibility(4);
240 : //UNUSED: Double phase;
241 : //UNUSED: Double phaseMult;
242 200 : Vector<Double> frequency;
243 200 : frequency= vb.frequency();
244 200 : Vector<Double> invLambda = frequency/C::c;
245 : // Find the offsets in polarization.
246 200 : Vector<Int> corrTypeL = vb.corrType().copy();
247 200 : Vector<Int> corrTypeC = vb.corrType().copy();
248 200 : Vector<Int> corrType = vb.corrType().copy();
249 200 : corrTypeL -= 9;
250 200 : corrTypeC -= 5;
251 200 : Cube<DComplex> dVis(4, nChan, nRow);
252 200 : ComponentType::Polarisation poltype=ComponentType::CIRCULAR;
253 200 : if(anyGT(Int(Stokes::RR), vb.corrType())){
254 0 : poltype=ComponentType::STOKES;
255 0 : corrType = vb.corrType()-1;
256 : }
257 200 : else if(vb.polFrame()==MSIter::Linear) {
258 0 : poltype=ComponentType::LINEAR;
259 0 : corrType = corrTypeL;
260 : } else {
261 200 : poltype=ComponentType::CIRCULAR;
262 200 : corrType = corrTypeC;
263 : }
264 :
265 200 : uInt npart=1;
266 : #ifdef _OPENMP
267 200 : npart= numthreads_p <0 ? omp_get_max_threads() : min(numthreads_p, omp_get_max_threads());
268 : #endif
269 :
270 :
271 :
272 200 : if((nRow/npart)==0) npart=1;
273 200 : Block<Cube<DComplex> > dVisp(npart);
274 200 : Vector<uInt> nRowp(npart);
275 200 : Block<Matrix<Double> > uvwp(npart);
276 200 : Block<ComponentList> compp(npart);
277 200 : Block<MeasFrame> mFramep(npart);
278 200 : Block<MDirection> mDir(npart);
279 200 : Block<uInt> startrow(npart);
280 :
281 200 : nRowp.set(nRow/npart);
282 200 : nRowp(npart-1) += nRow%npart;
283 :
284 :
285 200 : Block<Vector<Double> > dphasecomp(ncomponents);
286 200 : Block<Matrix<Double> > uvwcomp(ncomponents);
287 200 : Block<Block<Vector<Double> > > dphasecomps(npart);
288 200 : Block<Block<Matrix<Double> > > uvwcomps(npart);
289 6600 : for (uInt jj=0; jj < npart; ++jj){
290 6400 : dphasecomps[jj].resize(ncomponents);
291 6400 : uvwcomps[jj].resize(ncomponents);
292 : }
293 : ///Have to do this in one thread as MeasFrame and thus uvwmachine is not thread safe
294 400 : for (uInt k=0; k < ncomponents; ++k){
295 200 : uvwcomp[k]=uvw;
296 200 : dphasecomp[k].resize(nRow);
297 200 : rotateUVW(uvwcomp[k], dphasecomp[k], vb, compList.component(k).shape().refDirection());
298 200 : dphasecomp[k] *= -C::_2pi;
299 200 : Int sumrow=0;
300 6600 : for (uInt jj=0; jj < npart; ++jj){
301 6400 : dphasecomps[jj][k]=dphasecomp[k](IPosition(1,sumrow), IPosition(1, sumrow+nRowp(jj)-1));
302 6400 : uvwcomps[jj][k]=uvwcomp[k](IPosition(2,0, sumrow ), IPosition(2, 2, sumrow+nRowp(jj)-1));
303 6400 : sumrow+=nRowp(jj);
304 : }
305 :
306 : }
307 :
308 :
309 :
310 : //UNUSED: Int sumrow=0;
311 200 : Record lala;
312 200 : String err;
313 200 : compList.toRecord(err, lala);
314 6600 : for (uInt k=0; k < npart; ++k){
315 6400 : compp[k].fromRecord(err, lala);
316 6400 : startrow[k]=0;
317 105600 : for (uInt j=0; j < k; ++j) startrow[k]+=nRowp(j);
318 : }
319 : Bool isCopy;
320 200 : Complex *modData=modelData.getStorage(isCopy);
321 200 : MDirection dataMDir=vb.phaseCenter(phaseCenterTime_p);
322 :
323 200 : #pragma omp parallel default(none) firstprivate(npart, npol, nChan, modData, corrType, nRowp, invLambda, frequency, poltype) shared(startrow, compp, uvwcomps, dphasecomps) num_threads(npart)
324 : {
325 : #pragma omp for
326 :
327 : for (Int k=0; k < Int(npart); ++k){
328 :
329 : predictVis(modData, invLambda, frequency, compp[k], poltype, corrType,
330 : startrow[k], nRowp[k], nChan, npol, uvwcomps[k], dphasecomps[k]);
331 :
332 : }
333 :
334 :
335 : }
336 200 : modelData.putStorage(modData, isCopy);
337 :
338 200 : }
339 :
340 :
341 6400 : void SimpleComponentFTMachine::predictVis(Complex*& modData, const Vector<Double>& invLambda,
342 : const Vector<Double>& frequency, const ComponentList& compList,
343 : ComponentType::Polarisation poltype, const Vector<Int>& corrType,
344 : const uInt startrow, const uInt nrows, const uInt nchan, const uInt npol, const Block<Matrix<Double> > & uvwcomp, const Block<Vector<Double> > & dphasecomp){
345 6400 : Cube<DComplex> dVis(4, nchan, nrows);
346 6400 : uInt ncomponents=compList.nelements();
347 6400 : Vector<Double> dphase(nrows);
348 :
349 12800 : for (uInt icomp=0;icomp<ncomponents; ++icomp) {
350 6400 : SkyComponent component=compList.component(icomp);
351 6400 : component.flux().convertUnit(Unit("Jy"));
352 6400 : component.flux().convertPol(poltype);
353 6400 : MDirection compdir=component.shape().refDirection();
354 6400 : Vector<Double> thisRow(3);
355 6400 : thisRow=0.0;
356 : //UNUSED: uInt i;
357 6400 : Matrix<Double> uvw;
358 6400 : uvw=uvwcomp[icomp];
359 6400 : Vector<Double> dphase;
360 6400 : dphase=dphasecomp[icomp];
361 6400 : component.visibility(dVis, uvw, frequency);
362 : Double phaseMult;
363 : Double phase;
364 : uInt realrow;
365 : //// Loop over all rows
366 76600 : for (uInt r = 0; r < nrows ; r++) {
367 70200 : phaseMult = dphase(r);
368 70200 : realrow=r+startrow;
369 7090200 : for (uInt chn = 0; chn < nchan; chn++) {
370 7020000 : phase = phaseMult * invLambda(chn);
371 7020000 : Complex phasor(cos(phase), sin(phase));
372 35100000 : for (uInt pol=0; pol < npol; pol++) {
373 28080000 : const DComplex& val = dVis(corrType(pol), chn, r);
374 28080000 : modData[realrow*npol*nchan+chn*npol + pol] += Complex(val.real(), val.imag()) * conj(phasor);
375 : }
376 : }
377 : }
378 :
379 :
380 6400 : }
381 :
382 6400 : }
383 :
384 :
385 :
386 : } //# NAMESPACE CASA - END
387 :
|
