Line data Source code
1 : //# CalSet.cc: Implementation of Calibration parameter cache
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: 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 :
27 : #include <synthesis/CalTables/CalSet.h>
28 :
29 : #include <casacore/casa/Arrays.h>
30 : #include <casacore/casa/BasicSL/String.h>
31 : #include <casacore/casa/Utilities/GenSort.h>
32 : #include <casacore/casa/Exceptions/Error.h>
33 :
34 : #include <sstream>
35 :
36 : #include <casacore/casa/Logging/LogMessage.h>
37 : #include <casacore/casa/Logging/LogSink.h>
38 :
39 : using namespace casacore;
40 : namespace casa { //# NAMESPACE CASA - BEGIN
41 :
42 :
43 0 : void smooth(CalSet<Complex>& cs,
44 : const String& smtype,
45 : const Double& smtime,
46 : Vector<Int> selfields) {
47 :
48 : // half-width
49 0 : Double thw(smtime/2.0);
50 :
51 : // Workspace
52 0 : Vector<Double> times;
53 0 : Vector<Int> fields;
54 0 : Vector<Int> slotidx;
55 0 : Vector<Complex> p;
56 0 : Vector<Bool> pOK, newpOK;
57 0 : Vector<Float> amp;
58 0 : Vector<Float> pha;
59 0 : Vector<Bool> mask;
60 0 : Float newamp(0.0), newpha(0.0);
61 :
62 0 : IPosition blc(4,0,0,0,0);
63 0 : IPosition trc(4,0,0,0,0);
64 0 : IPosition vec(1,0);
65 :
66 : // For each spw
67 0 : for (Int ispw=0;ispw<cs.nSpw();++ispw) {
68 :
69 0 : Int nSlot=cs.nTime(ispw);
70 :
71 : // Only if more than one slot in this spw
72 0 : if (nSlot>1) {
73 :
74 0 : vec(0)=nSlot;
75 0 : trc(3)=nSlot-1;
76 :
77 0 : times.reference(cs.time(ispw));
78 0 : fields.reference(cs.fieldId(ispw));
79 0 : slotidx.resize(nSlot);
80 0 : indgen(slotidx);
81 0 : newpOK.resize(nSlot);
82 :
83 : // Discern how many fields we must do
84 0 : Int nFld=selfields.nelements();
85 : // Do all fields present, if none explicitly specificed
86 0 : if (nFld==0) {
87 0 : selfields=fields;
88 0 : nFld=genSort(selfields,Sort::Ascending,(Sort::QuickSort | Sort::NoDuplicates));
89 0 : selfields.resize(nFld,true);
90 : }
91 :
92 : // Arrange to mask/index each field
93 0 : PtrBlock<Vector<Bool>*> fldmask(nSlot,NULL);
94 0 : PtrBlock<Vector<Int>*> fldidx(nSlot,NULL);
95 0 : for (Int ifld=0;ifld<nFld;++ifld) {
96 0 : fldmask[ifld] = new Vector<Bool>;
97 0 : (*fldmask[ifld]) = (fields==selfields(ifld));
98 0 : fldidx[ifld] = new Vector<Int>;
99 0 : (*fldidx[ifld]) = slotidx((*fldmask[ifld])).getCompressedArray();
100 : }
101 :
102 : // For each elem (ant or baseline)
103 0 : for (Int ielem=0;ielem<cs.nElem();++ielem) {
104 0 : blc(2)=trc(2)=ielem;
105 : // For each channel
106 0 : for (int ichan=0;ichan<cs.nChan(ispw);++ichan) {
107 0 : blc(1)=trc(1)=ichan;
108 : // For each param (pol)
109 0 : for (Int ipar=0;ipar<cs.nPar();++ipar) {
110 0 : blc(0)=trc(0)=ipar;
111 :
112 : // Reference slices of par/parOK
113 0 : p.reference(cs.par(ispw)(blc,trc).reform(vec));
114 0 : pOK.reference(cs.parOK(ispw)(blc,trc).reform(vec));
115 0 : newpOK=pOK;
116 :
117 0 : IPosition psh(p.shape());
118 0 : amp.resize(psh);
119 0 : pha.resize(psh);
120 :
121 : /*
122 : cout << ispw << " "
123 : << ielem << " "
124 : << ichan << " "
125 : << ipar << " "
126 : << "p.shape() = " << p.shape() << " "
127 : << "pOK.shape() = " << pOK.shape() << " "
128 : << "amp.shape() = " << amp.shape() << " "
129 : << "pha.shape() = " << pha.shape() << " "
130 : << endl;
131 : */
132 : // Copy out amp and phase for processing
133 0 : amplitude(amp,p);
134 0 : phase(pha,p);
135 :
136 : // Filter each field separately
137 0 : for (Int ifld=0;ifld<nFld;++ifld) {
138 :
139 0 : Vector<Int> idx;
140 0 : idx.reference(*fldidx[ifld]);
141 :
142 : // If more than one slot for this field
143 0 : Int nidx=idx.nelements();
144 0 : if (nidx>1) {
145 :
146 : // Remove any phase cycles
147 : // (TBD: improve this algorithm?)
148 0 : Float phdif(0.0);
149 0 : for (Int i=1;i<nidx;++i) {
150 0 : phdif=pha(idx(i))-pha(idx(i-1));
151 0 : if (phdif > C::pi) {
152 0 : pha(idx(i)) -= (2*C::pi);
153 : // cout << " **************cycle+++++++++++" << endl;
154 : }
155 0 : else if (phdif < -C::pi) {
156 0 : pha(idx(i)) += (2*C::pi);
157 : // cout << " **************cycle-----------" << endl;
158 : }
159 : }
160 :
161 0 : Vector<Bool> mask;
162 0 : for (Int i=0;i<nidx;++i) {
163 : // Make mask
164 0 : mask = (*fldmask[ifld]);
165 0 : mask = (mask && pOK);
166 0 : mask = (mask && ( (times > (times(idx(i))-thw)) &&
167 0 : (times <= (times(idx(i))+thw)) ) );
168 :
169 : // Avoid explicit zeros, for now
170 0 : mask = (mask && amp>=FLT_MIN);
171 :
172 :
173 : //cout << " " << ifld << " " << i << " " << idx(i) << " ";
174 : //for (Int j=0;j<mask.nelements();++j)
175 : // cout << mask(j);
176 : //cout << endl;
177 :
178 0 : if (ntrue(mask)>0) {
179 0 : if (smtype=="mean") {
180 0 : newamp = mean(amp(mask));
181 0 : newpha = mean(pha(mask));
182 : }
183 0 : else if (smtype=="median") {
184 0 : newamp = median(amp(mask),false);
185 0 : newpha = median(pha(mask),false);
186 : }
187 0 : p(idx(i)) = Complex(cos(newpha),sin(newpha))*newamp;
188 0 : newpOK(idx(i))=true;
189 : }
190 : else
191 0 : newpOK(idx(i))=false;
192 :
193 : } // i
194 : } // nidx>1
195 : } // ifld
196 : // keep new ok info
197 0 : pOK=newpOK;
198 : } // ipar
199 : } // ichan
200 : } // ielem
201 :
202 : // Delete the PtrBlocks
203 0 : for (Int ifld=0;ifld<nFld;++ifld) {
204 0 : delete fldmask[ifld];
205 0 : delete fldidx[ifld];
206 : }
207 0 : fldmask=NULL;
208 0 : fldidx=NULL;
209 : } // nSlot>1
210 : } // ispw
211 :
212 0 : }
213 :
214 : } //# NAMESPACE CASA - END
|
