Line data Source code
1 : //# CalInterp.cc: Implementation of Calibration Interpolation
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,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: 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 : #include <synthesis/CalTables/CalInterp.h>
28 :
29 : #include <casacore/casa/Arrays.h>
30 : #include <casacore/scimath/Mathematics/MatrixMathLA.h>
31 : #include <casacore/casa/BasicSL/String.h>
32 : #include <casacore/casa/Utilities/Assert.h>
33 : #include <casacore/casa/Utilities/BinarySearch.h>
34 : #include <casacore/casa/Exceptions/Error.h>
35 :
36 : #include <sstream>
37 :
38 : #include <casacore/casa/Logging/LogMessage.h>
39 : #include <casacore/casa/Logging/LogSink.h>
40 :
41 : using namespace casacore;
42 : namespace casa { //# NAMESPACE CASA - BEGIN
43 :
44 0 : CalInterp::CalInterp(CalSet<Complex>& cs,
45 : const String& timetype,
46 0 : const String& freqtype) :
47 0 : cs_(&cs),
48 0 : timeType_(timetype),
49 0 : freqType_(freqtype),
50 0 : spwMap_(cs.nSpw(),-1),
51 0 : spwOK_(cs.nSpw(),false),
52 0 : lastTime_(cs.nSpw(),-1.0e99),
53 0 : finit_(cs.nSpw(),false),
54 0 : nFreq_(cs.nSpw(),1),
55 0 : solFreq_(cs.nSpw(),NULL),
56 0 : datFreq_(cs.nSpw(),NULL),
57 0 : currSpw_(-1),
58 0 : currSlot_(cs.nSpw(),-1),
59 0 : exactTime_(false),
60 :
61 0 : ip4d_(cs.nSpw(),NULL),
62 0 : ip3d_(cs.nSpw(),NULL),
63 0 : ip2d_(cs.nSpw(),NULL),
64 :
65 0 : t0_(cs.nSpw(),0.0),
66 0 : tS_(cs.nSpw(),0.0),
67 0 : lastlo_(cs.nSpw(),-1),
68 :
69 0 : tAC_(cs.nSpw(),NULL),
70 0 : tPC_(cs.nSpw(),NULL),
71 0 : tCC_(cs.nSpw(),NULL),
72 0 : tOk_(cs.nSpw(),NULL),
73 :
74 0 : ch0_(cs.nSpw(),NULL),
75 0 : ef_(cs.nSpw(),NULL),
76 0 : df_(cs.nSpw(),NULL),
77 0 : fdf_(cs.nSpw(),NULL),
78 :
79 0 : fAC_(cs.nSpw(),NULL),
80 0 : fPC_(cs.nSpw(),NULL),
81 0 : fCC_(cs.nSpw(),NULL),
82 0 : fOk_(cs.nSpw(),NULL),
83 :
84 0 : verbose_(false)
85 :
86 : {
87 :
88 0 : if (verbose_) cout << "CalInterp::constructor" << endl;
89 :
90 0 : if (verbose_) cout << " timeType_ = " << timeType_ << endl;
91 0 : if (verbose_) cout << " freqType_ = " << freqType_ << endl;
92 :
93 : // Nominally, spwOK_ follows CalSet
94 0 : spwOK_ = cs_->spwOK();
95 :
96 : // Allocate (zero-size) working arrays and shapes
97 : // (will resize non-trivially as needed)
98 0 : for (currSpw_=0;currSpw_<cs.nSpw();currSpw_++) {
99 :
100 : // nFreq()=1 here, will set this in initFreqInterp
101 0 : ip4d_[currSpw()] = new IPosition(4,2,nPar(),nFreq(),nElem());
102 0 : ip3d_[currSpw()] = new IPosition(3,nPar(),nFreq(),nElem());
103 0 : ip2d_[currSpw()] = new IPosition(2,nFreq(),nElem());
104 :
105 0 : tAC_[currSpw()] = new Array<Float>();
106 0 : tOk_[currSpw()] = new Cube<Bool>();
107 0 : tPC_[currSpw()] = new Array<Float>();
108 0 : tCC_[currSpw()] = new Array<Complex>();
109 :
110 0 : fAC_[currSpw()] = new Array<Float>();
111 0 : fOk_[currSpw()] = new Cube<Bool>();
112 0 : fPC_[currSpw()] = new Array<Float>();
113 0 : fCC_[currSpw()] = new Array<Complex>();
114 :
115 0 : df_[currSpw()] = new Vector<Double>();
116 0 : fdf_[currSpw()] = new Vector<Double>();
117 0 : ch0_[currSpw()] = new Vector<Int>();
118 0 : ef_[currSpw()] = new Vector<Bool>();
119 : }
120 :
121 0 : currSpw_=-1;
122 :
123 0 : };
124 :
125 :
126 0 : CalInterp::~CalInterp() {
127 :
128 0 : if (verbose_) cout << "CalInterp::destructor" << endl;
129 0 : deflTimeC();
130 0 : deflFreqC();
131 0 : deflFreqA();
132 0 : }
133 :
134 :
135 0 : Bool CalInterp::interpolate(const Double& time,
136 : const Int& spw,
137 : const Vector<Double>& /*freq*/) {
138 :
139 0 : if (verbose_) cout << endl << "CalInterp::interpolate()" << endl;
140 :
141 : // TODO:
142 : // - catch case where requested spw has no solutions
143 :
144 : // Assume no change, for starters
145 0 : Bool newInterp(false);
146 :
147 : // If spw has changed, re-map spw
148 0 : currSpw()=spw;
149 :
150 0 : newInterp = interpTime(time);
151 :
152 : // Interpolate in Freq
153 : // if (newInterp) interpFreq(freq);
154 :
155 : // Finalize interpolation
156 0 : finalize();
157 :
158 0 : return newInterp;
159 :
160 : }
161 :
162 :
163 0 : Bool CalInterp::interpTime(const Double& time) {
164 :
165 : // Interpolate in Time
166 0 : if (verbose_) cout << "CalInterp::interpTime()" << endl;
167 :
168 0 : Bool newTime=false; // assume no new calculations needed
169 :
170 0 : Bool newSlot(false);
171 0 : if (time != lastTime() ) {
172 0 : lastTime()=time;
173 :
174 : // if not "nearest", must recalcuate time interp
175 0 : if ( !nearestT() ) newTime=true;
176 :
177 : // Find relevant time slot
178 0 : newSlot = findSlot(time);
179 :
180 0 : if (newSlot)
181 0 : newTime = true;
182 :
183 0 : if (!exactTime() && !nearestT())
184 0 : updTimeCoeff();
185 :
186 : }
187 :
188 0 : if (verbose_) cout << " " << boolalpha
189 0 : << "newTime = " << newTime << " "
190 0 : << "newSlot = " << newSlot << " "
191 0 : << "currSlot() = " << currSlot() << " "
192 0 : << "fieldId = " << cs_->fieldId(currSpwMap())(currSlot()) << " "
193 0 : << "exactTime() = " << exactTime() << " "
194 0 : << "nearestT() = " << nearestT() << " "
195 0 : << endl;
196 :
197 : // newTime=true here if we need to re-calc time interp (*any* mode)
198 :
199 0 : if (newTime) {
200 :
201 0 : if ( nearestT() || exactTime() ) {
202 0 : exactTime_=true; // Behave as exact
203 :
204 0 : if (verbose_) cout << " "
205 0 : << "FOUND EXACT TIME!" << endl;
206 : // Just reference CalSet parameter
207 0 : IPosition blc(4,0,0,0,currSlot());
208 0 : IPosition trc(4,nPar()-1,nChan()-1,nElem()-1,currSlot());
209 :
210 0 : Cube<Complex> t;
211 0 : t.reference(csPar()(blc,trc).reform(IPosition(3,nPar(),nChan(),nElem())));
212 0 : r_.reference(t);
213 0 : tOk().reference(csParOK()(blc,trc).reform(IPosition(3,nPar(),nChan(),nElem())));
214 :
215 0 : } else {
216 :
217 0 : if (verbose_) cout << " "
218 0 : << "NON-EXACT TIME." << endl;
219 :
220 : // Do non-trivial interpolation
221 : // r_.resize(0,0,0);
222 :
223 :
224 0 : interpTimeCalc(time);
225 : }
226 :
227 0 : if (verbose_)
228 0 : cout << " CalInterp addr: " << r_.data()
229 0 : << endl;
230 :
231 :
232 :
233 : }
234 :
235 0 : return newTime; // signals whether new interp calc required
236 :
237 : }
238 :
239 0 : void CalInterp::interpFreq(const Vector<Double>& freq) {
240 :
241 0 : if (verbose_) cout << "CalInterp::interpFreq()" << endl;
242 :
243 : // Only if more than one cal channel is freq interp (potentially) necessary
244 0 : if (nChan() > 1) {
245 :
246 : // setup freq interp info (~no-op if already done)
247 0 : initFreqInterp(freq);
248 :
249 : // if all exact freqs, just reference time interp result:
250 0 : if ( allEQ(ef(),true) ) {
251 :
252 : // un-reference r
253 0 : r.resize();
254 :
255 : } else {
256 : // do non-trivial freq interpolation
257 :
258 : // ensure correct info available from time interpolation
259 0 : calcAPC();
260 :
261 : // update coeffs
262 0 : updFreqCoeff();
263 :
264 : // do freq interp calc
265 0 : interpFreqCalc();
266 :
267 : }
268 :
269 : }
270 :
271 0 : }
272 :
273 0 : Bool CalInterp::findSlot(const Double& time) {
274 :
275 0 : if (verbose_) cout << "CalInterp::findSlot()" << endl;
276 :
277 0 : Int slot(-1);
278 :
279 0 : Bool newSlot(false); // Assume no change
280 :
281 : // If only one slot, we've found it
282 0 : if (nTime()==1) {
283 0 : slot=0;
284 0 : exactTime_=true;
285 :
286 : // More than one slot, find the right one:
287 : } else {
288 :
289 0 : Vector<Double> timelist(csTimes());
290 :
291 0 : if (exactTime_) newSlot=true;
292 :
293 : // Behave as nearest outside absolute range of available calibration
294 : // to avoid wild extrapolation, else search for the correct soln slot
295 0 : if (time<timelist(0)) {
296 : // Before first timestamp, use first slot as nearest one
297 0 : slot=0;
298 0 : exactTime_=true;
299 : }
300 0 : else if (time>timelist(nTime()-1)) {
301 : // After last timestamp, use last slot as nearest one
302 0 : slot=nTime()-1;
303 0 : exactTime_=true;
304 : }
305 : else
306 : // Find index in timelist where time would be:
307 0 : slot=binarySearch(exactTime_,timelist,time,nTime(),0);
308 :
309 : // cout << "time = " << time << " slot = " << slot << " nTime() = " << nTime() << endl;
310 : // If not already an exact match...
311 0 : if ( !exactTime_ ) {
312 :
313 : // identify this time via index just prior
314 0 : if (slot>0) slot--;
315 :
316 : // If nearest, fine-tune slot to actual nearest:
317 0 : if ( timeType()=="nearest" ) {
318 0 : exactTime_=true; // Nearest behaves like exact match
319 0 : if (slot!=nTime()-1 &&
320 0 : (timelist(slot+1)-time) < (time-timelist(slot)) )
321 0 : slot++;
322 : } else {
323 : // linear modes require one later slot
324 0 : if (slot==nTime()-1) slot--;
325 : }
326 : }
327 :
328 0 : }
329 0 : newSlot = (slot!=currSlot());
330 0 : if (newSlot) {
331 0 : currSlot_(currSpw_)=slot;
332 : }
333 :
334 0 : return newSlot;
335 :
336 : }
337 :
338 0 : void CalInterp::updTimeCoeff() {
339 :
340 0 : if (verbose_) cout << "CalInterp::updTimeCoeff()" << endl;
341 :
342 0 : if ( currSlot() != lastlo() ) {
343 0 : lastlo()=currSlot();
344 :
345 : // Resize Coefficient arrays
346 0 : IPosition ip4s(4,2,nPar(),nChan(),nElem());
347 0 : IPosition ip3s(3,nPar(),nChan(),nElem());
348 :
349 0 : tAC().resize(ip4s);
350 0 : tOk().resize(); // ensure not referencing csParOK!
351 0 : tOk().resize(ip3s);
352 :
353 0 : if ( timeType()=="linear")
354 0 : tPC().resize(ip4s);
355 0 : else if (timeType()=="aipslin")
356 0 : tCC().resize(ip4s);
357 :
358 : // Time ref/step for this interval
359 0 : t0()=csTimes()(currSlot());
360 0 : tS()=csTimes()(currSlot()+1)-t0();
361 :
362 : // For indexing parameter cache
363 0 : IPosition lo(4,0,0,0,currSlot()), hi(4,0,0,0,currSlot()+1);
364 0 : IPosition ref(4,0,0,0,0), slope(4,1,0,0,0);
365 :
366 0 : for (Int ielem=0;ielem<nElem();ielem++) {
367 0 : lo(2)=hi(2)=ref(3)=slope(3)=ielem;
368 0 : for (Int ichan=0;ichan<nChan();ichan++) {
369 0 : lo(1)=hi(1)=ref(2)=slope(2)=ichan;
370 0 : for (Int ipar=0;ipar<nPar();ipar++) {
371 0 : lo(0)=hi(0)=ref(1)=slope(1)=ipar;
372 :
373 0 : tOk()(ipar,ichan,ielem) = (csParOK()(lo) && csParOK()(hi));
374 :
375 0 : if (tOk()(ipar,ichan,ielem)) {
376 : // Intercept
377 0 : tAC()(ref) = abs(csPar()(lo));
378 0 : tAC()(slope) = abs(csPar()(hi)) - tAC()(ref);
379 :
380 0 : if (timeType()=="linear") {
381 0 : tPC()(ref) = arg(csPar()(lo));
382 :
383 : // Slope
384 0 : Float pslope = arg(csPar()(hi)) - tPC()(ref);
385 : // Catch simple phase wraps
386 0 : if (pslope > C::pi)
387 0 : pslope-=(2*C::pi);
388 0 : else if (pslope < -C::pi)
389 0 : pslope+=(2*C::pi);
390 0 : tPC()(slope) = pslope;
391 :
392 0 : } else if (timeType()=="aipslin") {
393 0 : tCC()(ref) = csPar()(lo);
394 0 : tCC()(slope) = csPar()(hi)-tCC()(ref);
395 :
396 : }
397 :
398 : } else {
399 0 : tAC()(ref)=1.0;
400 0 : tAC()(slope)=0.0;
401 0 : if (timeType()=="linear") {
402 0 : tPC()(ref)=0.0;
403 0 : tPC()(slope)=0.0;
404 0 : } else if (timeType()=="aipslin") {
405 0 : tCC()(ref)=Complex(1.0,0.0);
406 0 : tCC()(slope)=Complex(0.0,0.0);
407 : }
408 : }
409 : }
410 : }
411 : }
412 0 : }
413 :
414 0 : }
415 :
416 :
417 0 : void CalInterp::interpTimeCalc(const Double& time) {
418 :
419 0 : if (verbose_) cout << "CalInterp::interpTimeCalc()" << endl;
420 :
421 : // TODO:
422 : // a. Use matrix math instead of loops? (tOk() usage?)
423 :
424 : // Fractional time interval for this timestamp
425 0 : Float dt( Float( (time-t0())/tS() ) );
426 :
427 : // Ensure intermediate results cache is properly sized and ref'd
428 0 : if (tA_.nelements()==0) {
429 0 : tA_.resize(nPar(),nChan(),nElem());
430 0 : a.reference(tA_);
431 0 : ok.reference(tOk());
432 0 : if (timeType()=="linear") {
433 0 : tP_.resize(nPar(),nChan(),nElem());
434 0 : p.reference(tP_);
435 0 : c.resize();
436 : }
437 0 : else if (timeType()=="aipslin") {
438 0 : tC_.resize(nPar(),nChan(),nElem());
439 0 : c.reference(tC_);
440 0 : p.resize();
441 : }
442 : }
443 :
444 0 : IPosition ref(4,0,0,0,0),slope(4,1,0,0,0);
445 0 : for (Int ielem=ref(3)=slope(3)=0;
446 0 : ielem<nElem();
447 0 : ielem++,ref(3)++,slope(3)++) {
448 0 : for (Int ichan=ref(2)=slope(2)=0;
449 0 : ichan<nChan();
450 0 : ichan++,ref(2)++,slope(2)++) {
451 0 : for (Int ipar=ref(1)=slope(1)=0;
452 0 : ipar<nPar();
453 0 : ipar++,ref(1)++,slope(1)++) {
454 :
455 0 : if (tOk()(ipar,ichan,ielem)) {
456 :
457 0 : tA_(ipar,ichan,ielem) = tAC()(ref) + tAC()(slope)*dt;
458 0 : if (timeType()=="linear") {
459 0 : tP_(ipar,ichan,ielem) = tPC()(ref) + tPC()(slope)*dt;
460 0 : } else if (timeType()=="aipslin") {
461 0 : Complex tCtmp(tCC()(ref) + tCC()(slope)*dt);
462 0 : Float Amp(abs(tCtmp));
463 0 : if (Amp>0.0)
464 0 : tC_(ipar,ichan,ielem) = tCtmp/abs(tCtmp);
465 : else
466 0 : tC_(ipar,ichan,ielem) = Complex(1.0);
467 : }
468 : }
469 : else {
470 0 : tA_(ipar,ichan,ielem) = 1.0;
471 0 : if (timeType()=="linear") {
472 0 : tP_(ipar,ichan,ielem) = 0.0;
473 0 : } else if (timeType()=="aipslin") {
474 0 : tC_(ipar,ichan,ielem) = Complex(1.0,0.0);
475 : }
476 : } // tOk()
477 : } // ipar
478 : } // ichan
479 : } // ielem
480 :
481 0 : if (verbose_) {
482 0 : cout << "tA_ = " << tA_.nonDegenerate() << endl;
483 0 : if (timeType()=="linear")
484 0 : cout << "tP_ = " << tP_.nonDegenerate() << endl;
485 0 : else if (timeType()=="aipslin")
486 0 : cout << "tC_ = " << tC_.nonDegenerate() << endl;
487 : }
488 :
489 0 : }
490 :
491 :
492 0 : void CalInterp::initFreqInterp(const Vector<Double> freq) {
493 :
494 0 : if (verbose_) cout << "CalInterp::initFreqInterp()" << endl;
495 :
496 : // Initialize freq interpolation if
497 : // a. not yet initialized (first time thru), or,
498 : // b. # of requested frequencies not equal to previous, or,
499 : // (c. # of requested freqs same, but freqs different --- NYI)
500 :
501 :
502 0 : if (nFreq() != Int(freq.nelements()) ||
503 0 : !allEQ(freq,datFreq()) )
504 0 : finit()=false;
505 :
506 0 : if (!finit()) {
507 :
508 : // Remember number of requested frequencies
509 0 : nFreq() = freq.nelements();
510 :
511 : // Remember the frequencies:
512 0 : datFreq() = freq;
513 :
514 0 : ch0().resize(nFreq()); ch0()=-1;
515 0 : ef().resize(nFreq()); ef()=false;
516 :
517 0 : if (nChan()==1) {
518 : // one-to-many case (e.g., G)
519 : // need not support non-trivial freq interpolation
520 0 : ch0()=0;
521 0 : ef()=true;
522 :
523 : } else {
524 : // many-to-many case (e.g., B)
525 : // support non-trivial freq interpolation if some non-exact freqs
526 :
527 : // Fill ref ch indices for each input freq:
528 0 : Int ichan(0);
529 0 : for (Int i=0;i<nFreq();i++) {
530 :
531 0 : ichan = binarySearch(ef()(i),csFreq(),freq(i),nChan(),0);
532 :
533 : // interp-type adjustments:
534 0 : if (!ef()(i)) {
535 0 : if ( freqType()=="nearest" ) {
536 : // refer to true nearest, behave as exact match
537 0 : ef()(i)=true;
538 0 : if ( ichan==nChan() ||
539 0 : (ichan>0 &&
540 0 : abs(freq(i)-csFreq()(ichan-1))<abs(csFreq()(ichan)-freq(i))) )
541 0 : ichan--;
542 : } else {
543 : // refer to index at low and of pair that brackets this freq, within bounds
544 0 : if (ichan>0) ichan--;
545 0 : if (ichan==nChan()-1) ichan--;
546 : }
547 : }
548 0 : ch0()(i)=ichan;
549 :
550 : }
551 :
552 : // Fill in frac freq info if some non-exact freqs
553 0 : if ( !allEQ(ef(),true) ) {
554 :
555 0 : df().resize(nFreq()); df()=0.0;
556 0 : fdf().resize(nFreq()); fdf()=0.0;
557 :
558 : // Fill in frac freq info
559 0 : for (Int i=0;i<nFreq()-1;i++) {
560 0 : if (!ef()(i)) {
561 0 : df()(i) = freq(i)-csFreq()(ch0()(i));
562 0 : fdf()(i) = df()(i)/abs( csFreq()(ch0()(i)+1)-csFreq()(ch0()(i)) );
563 : }
564 : }
565 : }
566 : }
567 :
568 : }
569 :
570 0 : }
571 :
572 :
573 0 : void CalInterp::calcAPC() {
574 :
575 0 : if (verbose_) cout << "CalInterp::calcAPC()" << endl;
576 :
577 : // re-package interpTimeCalc result for updFreqCoeff, as necessary
578 :
579 0 : if (freqType()=="linear") {
580 0 : if (timeType()=="nearest") {
581 : // need a/p from r
582 0 : tA_.resize(nPar(),nChan(),nElem());
583 0 : tA_ = amplitude(r);
584 0 : tP_.resize(nPar(),nChan(),nElem());
585 0 : tP_ = phase(r);
586 0 : r.resize(); // drop reference to nearest cached solution
587 : }
588 0 : else if (timeType()=="aipslin") {
589 : // a ok, need p
590 0 : tP_.resize(nPar(),nChan(),nElem());
591 0 : tP_ = phase(tC_);
592 : }
593 0 : } else if (freqType()=="aipslin") {
594 0 : if (timeType()=="nearest") {
595 : // need a/c
596 0 : tA_.resize(nPar(),nChan(),nElem());
597 0 : tA_ = amplitude(r);
598 0 : tC_.resize(nPar(),nChan(),nElem());
599 0 : tC_ = r;
600 0 : Array<Float> rp,ip;
601 0 : rPart(tC_,rp);
602 0 : iPart(tC_,ip);
603 :
604 : // NEED TO CHECK for tA_=0 here!
605 :
606 0 : rp/=tA_;
607 0 : ip/=tA_;
608 :
609 0 : r.resize(); // drop reference to nearest cached solution
610 0 : }
611 0 : else if (timeType()=="linear") {
612 : // a ok, need c
613 0 : tC_.resize(nPar(),nChan(),nElem());
614 0 : Array<Float> rp,ip;
615 0 : rPart(tC_,rp);
616 0 : iPart(tC_,ip);
617 0 : rp = cos(tP_);
618 0 : ip = sin(tP_);
619 0 : }
620 : }
621 :
622 0 : }
623 :
624 0 : void CalInterp::updFreqCoeff() {
625 :
626 0 : if (verbose_) cout << "CalInterp::updFreqCoeff()" << endl;
627 :
628 : // Resize Coefficient arrays (no-op if already correct size)
629 0 : ip4d()(2)=ip2d()(0)=nFreq();
630 0 : fAC().resize(ip4d());
631 0 : fOk().resize(ip3d());
632 0 : if ( timeType()=="linear")
633 0 : fPC().resize(ip4d());
634 0 : else if (timeType()=="aipslin")
635 0 : fCC().resize(ip4d());
636 :
637 : // For indexing parameter cache
638 0 : IPosition lo(3,0,0,0), hi(3,0,0,0);
639 0 : IPosition ref(4,0,0,0,0), slope(4,1,0,0,0);
640 0 : for (Int ielem=0;ielem<nElem();ielem++) {
641 0 : lo(2)=hi(2)=ref(3)=slope(3)=ielem;
642 0 : for (Int ichan=0;ichan<nChan()-1;ichan++) {
643 0 : lo(1)=hi(1)=ref(2)=slope(2)=ichan;
644 :
645 0 : if (!ef()(ichan)) {
646 0 : for (Int ipar=0;ipar<nPar();ipar++) {
647 0 : lo(0)=hi(0)=ref(1)=slope(1)=ipar;
648 :
649 0 : fOk()(ipar,ichan,ielem) = tOk()(ipar,ichan,ielem) &&
650 0 : tOk()(ipar,ichan+1,ielem);
651 0 : if (fOk()(ipar,ichan,ielem)) {
652 :
653 : // Intercept
654 0 : fAC()(ref) = tA_(lo);
655 0 : fAC()(slope) = tA_(hi) - tA_(lo);
656 :
657 0 : if (freqType()=="linear") {
658 0 : fPC()(ref) = tP_(lo);
659 :
660 : // Slope
661 0 : Float pslope = tP_(hi) - tP_(lo);
662 : // Catch/remove simple phase wraps
663 0 : if (pslope > C::pi)
664 0 : pslope-=(2*C::pi);
665 0 : else if (pslope < -C::pi)
666 0 : pslope+=(2*C::pi);
667 0 : fPC()(slope) = pslope;
668 :
669 0 : } else if (freqType()=="aipslin") {
670 0 : fCC()(ref) = tC_(lo);
671 0 : fCC()(slope) = tC_(hi)-tC_(lo);
672 :
673 : }
674 : } else {
675 0 : fAC()(ref)=1.0;
676 0 : fAC()(slope)=0.0;
677 0 : if (timeType()=="linear") {
678 0 : fPC()(ref)=0.0;
679 0 : fPC()(slope)=0.0;
680 0 : } else if (timeType()=="aipslin") {
681 0 : fCC()(ref)=Complex(1.0,0.0);
682 0 : fCC()(slope)=Complex(0.0,0.0);
683 : }
684 : } // fOk()
685 : } // ipar
686 : } // !ef()(ichan)
687 : } // ichan
688 : } // ielem
689 :
690 0 : }
691 :
692 :
693 0 : void CalInterp::interpFreqCalc() {
694 :
695 0 : if (verbose_) cout << "CalInterp::interpFreqCalc()" << endl;
696 :
697 : // TODO:
698 : // a. Use matrix math instead of loops? (fOk(), ef() usage?)
699 :
700 0 : fA_.resize(nPar(),nFreq(),nElem());
701 0 : a.reference(fA_);
702 0 : ok.reference(fOk());
703 0 : if (freqType()=="linear") {
704 0 : fP_.resize(nPar(),nFreq(),nElem());
705 0 : p.reference(fP_);
706 0 : c.resize();
707 : }
708 0 : else if (freqType()=="aipslin") {
709 0 : fC_.resize(nPar(),nFreq(),nElem());
710 0 : c.reference(fC_);
711 0 : p.resize();
712 : }
713 :
714 0 : IPosition ref(4,0,0,0,0),slope(4,1,0,0,0);
715 0 : for (Int ielem=ref(3)=slope(3)=0;
716 0 : ielem<nElem();
717 0 : ielem++,ref(3)++,slope(3)++) {
718 0 : for (Int ichan=ref(2)=slope(2)=0;
719 0 : ichan<nFreq();
720 0 : ichan++,ref(2)++,slope(2)++) {
721 : // if non-exact channel and otherwise ok, calc interp per par
722 0 : if ( !ef()(ichan) ) {
723 0 : for (Int ipar=ref(1)=slope(1)=0;
724 0 : ipar<nPar();
725 0 : ipar++,ref(1)++,ref(2)++) {
726 0 : if (fOk()(ipar,ichan,ielem) ) {
727 0 : fA_(ipar,ichan,ielem) = fAC()(ref) + fAC()(slope)*fdf()(ichan);
728 0 : if ( freqType()=="linear") {
729 0 : fP_(ipar,ichan,ielem) = fPC()(ref) + fPC()(slope)*fdf()(ichan);
730 0 : } else if ( freqType()=="aipslin") {
731 0 : fC_(ipar,ichan,ielem) = fCC()(ref) + fCC()(slope)*fdf()(ichan);
732 : }
733 : }
734 : else {
735 : // copy from time interp
736 0 : fA_(ipar,ichan,ielem) = tA_(ipar,ch0()(ichan),ielem);
737 0 : if ( freqType()=="linear") {
738 0 : fP_(ipar,ichan,ielem) = tP_(ipar,ch0()(ichan),ielem);
739 0 : } else if ( freqType()=="aipslin") {
740 0 : fC_(ipar,ichan,ielem) = tC_(ipar,ch0()(ichan),ielem);
741 : }
742 : } // fOk()
743 : } // ipar
744 : } // !ef()
745 : } // ichan
746 : } // ielem
747 :
748 0 : }
749 :
750 0 : void CalInterp::finalize() {
751 :
752 0 : if (verbose_) cout << "CalInterp::finalize()" << endl;
753 :
754 0 : if (!exactTime()) {
755 0 : if (p.nelements() > 0 ) {
756 : // if we did phase interp, shape and set r/i parts of r_
757 : // cout << "CalInterp::finalize(): pppppp" << endl;
758 0 : c.resize(p.shape());
759 :
760 0 : Array<Float> rp,ip;
761 0 : rPart(c,rp);
762 0 : iPart(c,ip);
763 :
764 0 : rp = a*cos(p);
765 0 : ip = a*sin(p);
766 :
767 0 : }
768 0 : else if (c.nelements() > 0) {
769 :
770 : // just reference c
771 0 : Array<Float> rp,ip;
772 0 : rPart(c,rp);
773 0 : iPart(c,ip);
774 0 : rp*=a;
775 0 : ip*=a;
776 :
777 0 : }
778 0 : r_.reference(c);
779 : }
780 :
781 : // make public by reference
782 0 : r.reference(r_);
783 :
784 0 : ok.reference(tOk());
785 :
786 0 : }
787 :
788 :
789 : // Deflate in-focus time interpolation coeff cache
790 0 : void CalInterp::deflTimeC() {
791 :
792 0 : for (Int ispw=0; ispw<nSpw(); ispw++) {
793 0 : if (tAC_[ispw]) delete tAC_[ispw];
794 0 : if (tOk_[ispw]) delete tOk_[ispw];
795 0 : if (tPC_[ispw]) delete tPC_[ispw];
796 0 : if (tCC_[ispw]) delete tCC_[ispw];
797 :
798 0 : tAC_[ispw]=tPC_[ispw]= NULL;
799 0 : tOk_[ispw]=NULL;
800 0 : tCC_[ispw]=NULL;
801 : }
802 0 : }
803 :
804 : // Deflate in-focus freq interpolation coefficient cache
805 0 : void CalInterp::deflFreqC() {
806 :
807 0 : for (Int ispw=0;ispw<nSpw();ispw++) {
808 0 : if (fAC_[ispw]) delete fAC_[ispw];
809 0 : if (fOk_[ispw]) delete fOk_[ispw];
810 0 : if (fPC_[ispw]) delete fPC_[ispw];
811 0 : if (fCC_[ispw]) delete fCC_[ispw];
812 :
813 0 : fAC_[ispw]=fPC_[ispw]=NULL;
814 0 : fOk_[ispw]=NULL;
815 0 : fCC_[ispw]=NULL;
816 : }
817 0 : }
818 :
819 : // Deflate in-focus freq interpolation abscissa data
820 0 : void CalInterp::deflFreqA() {
821 :
822 0 : for (Int ispw=0;ispw<nSpw();ispw++) {
823 0 : if ( ch0_[ispw] ) delete ch0_[ispw];
824 0 : if ( ef_[ispw] ) delete ef_[ispw];
825 0 : if ( df_[ispw] ) delete df_[ispw];
826 0 : if ( fdf_[ispw] ) delete fdf_[ispw];
827 :
828 0 : ch0_[ispw] = NULL;
829 0 : ef_[ispw] = NULL;
830 0 : df_[ispw] = fdf_[ispw] = NULL;
831 : }
832 0 : }
833 :
834 0 : void CalInterp::asFloatArr(const Array<Complex>& in, Array<Float>& out) {
835 0 : IPosition ip1(in.shape());
836 0 : IPosition ip2(1,2);
837 0 : ip2.append(ip1);
838 0 : out.takeStorage(ip2,(Float*)(in.data()),SHARE);
839 0 : }
840 :
841 0 : void CalInterp::part(const Array<Complex>& c,
842 : const Int& which,
843 : Array<Float>& f) {
844 :
845 0 : IPosition ip1(c.shape());
846 0 : Array<Float> asfl;
847 0 : asFloatArr(c,asfl);
848 0 : IPosition ip2(asfl.shape());
849 0 : IPosition blc(ip2), trc(ip2);
850 0 : blc=0; trc-=1;
851 0 : blc(0)=trc(0)=which;
852 0 : f.reference(asfl(blc,trc).reform(ip1));
853 0 : }
854 :
855 0 : void CalInterp::setSpwOK() {
856 :
857 0 : for (Int i=0;i<nSpw();i++)
858 0 : if (spwMap(i)>-1)
859 0 : spwOK_(i) = cs_->spwOK()(spwMap(i));
860 : else
861 0 : spwOK_(i) = cs_->spwOK()(i);
862 :
863 : // cout << "CalInterp::spwOK() (spwmap) = " << boolalpha << spwOK() << endl;
864 :
865 0 : }
866 :
867 :
868 :
869 :
870 : } //# NAMESPACE CASA - END
871 :
|
