Line data Source code
1 : //# RIorAParray.cc: Implementation of RI/AP on-demand converter
2 : //# Copyright (C) 2012
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 adressed 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 : //#
27 :
28 : #include <synthesis/CalTables/RIorAParray.h>
29 :
30 : #include <casacore/casa/aips.h>
31 : #include <casacore/casa/BasicSL/Constants.h>
32 : #include <casacore/casa/Arrays/Array.h>
33 : #include <casacore/casa/IO/ArrayIO.h>
34 : #include <casacore/casa/Arrays/ArrayIter.h>
35 : #include <casacore/casa/Arrays/Matrix.h>
36 : #include <casacore/casa/Arrays/ArrayMath.h>
37 : #include <casacore/casa/Logging/LogMessage.h>
38 : #include <casacore/casa/Logging/LogSink.h>
39 :
40 : #define RIORAPVERB false
41 :
42 : using namespace casacore;
43 : namespace casa { //# NAMESPACE CASA - BEGIN
44 :
45 :
46 : // Construct empty
47 0 : RIorAPArray::RIorAPArray() :
48 0 : c_ok_(false),
49 0 : f_ok_(false),
50 0 : phaseTracked_(false),
51 0 : c_(),
52 0 : f_()
53 0 : {}
54 :
55 : // Construct from external Complex Array
56 604763 : RIorAPArray::RIorAPArray(const Array<Complex>& c) :
57 604763 : c_ok_(false),
58 604763 : f_ok_(false),
59 604763 : phaseTracked_(false),
60 604763 : c_(),
61 604763 : f_()
62 : {
63 : if (RIORAPVERB) cout << "ctor(A<Complex>))" << endl;
64 :
65 : // Delegate to setData
66 604763 : this->setData(c);
67 604763 : }
68 :
69 : // Construct from external Float Array
70 88851 : RIorAPArray::RIorAPArray(const Array<Float>& f) :
71 88851 : c_ok_(false),
72 88851 : f_ok_(false),
73 88851 : phaseTracked_(false),
74 88851 : c_(),
75 88851 : f_()
76 : {
77 : if (RIORAPVERB) cout << "ctor(A<Float>))" << endl;
78 :
79 : // Delegate to setData
80 88851 : setData(f);
81 88851 : }
82 :
83 : // Destructor
84 693614 : RIorAPArray::~RIorAPArray() {};
85 :
86 :
87 604763 : void RIorAPArray::setData(const Array<Complex>& c) {
88 :
89 : // Discard any existing float part; will be created on-demand, if nec.
90 604763 : f_.resize();
91 604763 : f_ok_=false;
92 : // Reference incoming data array
93 604763 : if (c.ndim()==1) {
94 0 : IPosition ip=c.shape();
95 0 : ip.prepend(IPosition(1,1));
96 0 : c_.reference(c.reform(ip));
97 0 : }
98 : else
99 604763 : c_.reference(c);
100 : // Complex version now ok
101 604763 : c_ok_=true;
102 604763 : }
103 :
104 88851 : void RIorAPArray::setData(const Array<Float>& f) {
105 :
106 : // Discard any existing complex part; will be created on-demand, if nec.
107 88851 : c_.resize();
108 88851 : c_ok_=false;
109 :
110 : // Insist that incoming Float Array has ndim>1
111 88851 : if (f.ndim()<2)
112 0 : throw(AipsError("RIorAPArray: Input Float Array must be at least 2D."));
113 :
114 : // Reference incoming data array
115 88851 : f_.reference(f);
116 88851 : f_ok_=true;
117 88851 : }
118 :
119 : // State
120 0 : void RIorAPArray::state(Bool verbose) {
121 :
122 0 : cout << boolalpha;
123 0 : cout << "--state--" << endl;
124 0 : cout << "f_: ok=" << f_ok_ << " &=" << f_.data() << " sh=" << f_.shape() << " nrefs=" << f_.nrefs() << endl;
125 0 : cout << "c_: ok=" << c_ok_ << " &=" << c_.data() << " sh=" << c_.shape() << " nrefs=" << c_.nrefs() << endl;
126 :
127 0 : if (verbose) {
128 0 : cout.precision(10);
129 0 : cout << "f_ = " << f_ << endl;
130 0 : cout << "c_ = " << c_ << endl;
131 : }
132 0 : cout << "---------" << endl;
133 :
134 0 : }
135 :
136 :
137 : // Render Complex version (calc from Float, if necessary)
138 88851 : Array<Complex> RIorAPArray::c() {
139 88851 : if (!c_ok_) { // not already calculated
140 88851 : resizec_();
141 88851 : calc_c(); // calc internally
142 : }
143 88851 : return c_; // return the array
144 : }
145 :
146 :
147 : // Render Float version (calc from Complex, if necessary)
148 604763 : Array<Float> RIorAPArray::f(Bool trackphase) {
149 : // form it, if not already calculated or phase needs tracking
150 : // TBD optimize already-calc'd/needs phasetracked case
151 604763 : if (!f_ok_ || (trackphase && !phaseTracked_)) {
152 604763 : resizef_();
153 604763 : calc_f(trackphase);
154 : }
155 604763 : return f_; // return the array
156 : }
157 :
158 88851 : void RIorAPArray::resizec_() {
159 : if (RIORAPVERB) cout << "resizec_()" << endl;
160 88851 : IPosition cip(f_.shape());
161 88851 : cip(0)/=2; // First axis float->complex (half as long)
162 88851 : c_.resize(cip);
163 88851 : }
164 604763 : void RIorAPArray::resizef_() {
165 : if (RIORAPVERB) cout << "resizef_()" << endl;
166 604763 : IPosition fip(c_.shape());
167 604763 : fip(0)*=2; // First axis complex->float (twice as long)
168 604763 : f_.resize(fip);
169 :
170 : // Ensure that at least 2 axes...
171 604763 : if (fip.size()<2)
172 0 : throw(AipsError("RIorAPArray: Internal Float array, f_, must have ndim>1"));
173 604763 : }
174 :
175 : // Calc Complex version from Float info
176 : // (assumes c_ already correct size)
177 88851 : void RIorAPArray::calc_c() {
178 :
179 : if (RIORAPVERB) cout << "calc_c()" << endl;
180 :
181 88851 : if (!f_ok_)
182 0 : throw(AipsError("RIorAParray::f(): Can't calculate complex version from absent float version."));
183 :
184 88851 : Int ndim=f_.ndim();
185 88851 : Array<Float> amp;
186 88851 : Array<Float> ph;
187 88851 : IPosition blc(ndim,0),trc(f_.endPosition()),stp(ndim,1);
188 88851 : stp(0)=2; // by 2 in first axis
189 88851 : amp=f_(blc,trc,stp);
190 88851 : blc(0)=1; // phase is 2nd value on first axis
191 88851 : ph=f_(blc,trc,stp);
192 :
193 : // Form Float array with real/imag parts (not amp/phase)
194 88851 : Array<Float> ftmp(f_.shape());
195 88851 : blc(0)=0;
196 88851 : ftmp(blc,trc,stp)=amp*cos(ph);
197 88851 : blc(0)=1;
198 88851 : ftmp(blc,trc,stp)=amp*sin(ph);
199 :
200 : // Convert Float R/I array to Complex array
201 : // c_=RealToComplex(ftmp);
202 88851 : RealToComplex(c_,ftmp);
203 88851 : c_ok_=true;
204 :
205 88851 : }
206 :
207 : // Calc Float version from Complex info
208 : // (assumes f_ already correct size)
209 604763 : void RIorAPArray::calc_f(Bool trackphase) {
210 :
211 : if (RIORAPVERB) cout << "calc_f(" << boolalpha << trackphase << ")" << endl;
212 :
213 604763 : if (!c_ok_)
214 0 : throw(AipsError("RIorAParray::f(): Can't calculate float version from absent complex version."));
215 :
216 604763 : Int ndim=f_.ndim();
217 604763 : IPosition blc(ndim,0),trc(f_.endPosition()),stp(ndim,1);
218 604763 : stp(0)=2; // by 2 in first axis
219 604763 : Array<Float> amp(f_(blc,trc,stp));
220 604763 : blc(0)=1; // phase is 2nd value on first axis
221 604763 : Array<Float> ph(f_(blc,trc,stp));
222 :
223 604763 : amp=amplitude(c_).reform(amp.shape());
224 604763 : ph=phase(c_).reform(amp.shape());
225 604763 : f_ok_=true;
226 :
227 604763 : phaseTracked_=false;
228 :
229 604763 : if (trackphase)
230 604763 : trackPhase(ph);
231 :
232 604763 : }
233 :
234 : // Track phase _on_last_axis_
235 604763 : void RIorAPArray::trackPhase(Array<Float>& ph) {
236 :
237 : if (RIORAPVERB) cout << "trackPhase()" << endl;
238 :
239 604763 : IPosition phsh(ph.shape());
240 604763 : ArrayIterator<Float> phiter(ph,IPosition(1,ph.ndim()-1));
241 604763 : Vector<Float> ph1;
242 1293202 : while (!phiter.pastEnd()) {
243 688439 : ph1.reference(phiter.array());
244 837173 : for (uInt j=1;j<ph1.nelements();++j) {
245 148734 : if (ph1(j)>ph1(j-1))
246 77700 : while (ph1(j)>(ph1(j-1)+C::pi)) ph1(j)-=C::_2pi;
247 148734 : if (ph1(j)<ph1(j-1))
248 50614 : while (ph1(j)<(ph1(j-1)-C::pi)) ph1(j)+=C::_2pi;
249 : }
250 688439 : phiter.next();
251 : }
252 :
253 604763 : phaseTracked_=true;
254 :
255 604763 : }
256 :
257 :
258 :
259 : /*
260 :
261 : ****TBD: Handle filling supplied arrays? (to avoid copies)
262 :
263 :
264 : // Render Complex version onto supplied Array (calc from Float, if necessary)
265 : void RIorAPArray::c(Array<Complex>& c) {
266 :
267 :
268 : if (c_ok_) { // Already calculated internally
269 : if (c.conform(c_))
270 : c=c_;
271 : }
272 : else { // Do the calculation
273 :
274 : // f_ *must* be ok if we are going to calculate c_
275 : if (!f_ok_)
276 : throw(AipsError("RIorAParray::c(): Can't serve complex version from absent float version."));
277 :
278 : // Resize supplied storage
279 : IPosition cip(f_.shape());
280 : cip.getLast(cip.nelements()-1);
281 : c.resize(cip); // the _supplied_ array (could be external)
282 : if (&c_!=&c)
283 : c_.reference(c); // no-op if c_ and c are same array?
284 : // Do the calculation
285 : calc_c();
286 : }
287 :
288 : }
289 :
290 : // Render Float version onto supplied Array(calc from Complex, if necessary)
291 : void RIorAPArray::f(Array<Float>& f,Bool trackphase) {
292 : // c_ must be ok if we are going to calculate f_
293 : if (!c_ok_)
294 : throw(AipsError("RIorAParray::c(): Can't serve complex version from absent float version."));
295 :
296 : // Resize storage
297 : IPosition fip(c_.shape());
298 : fip.prepend(IPosition(1,2));
299 : f.resize(fip); // the supplied array (could be external)
300 : f_.reference(f); // no-op if f_ and f are same array?
301 :
302 : // Do the calculation
303 : calc_f(trackphase);
304 : }
305 : */
306 :
307 : } //# NAMESPACE CASA - END
|
