Line data Source code
1 : //# FJones.cc: Implementation of Ionosphere
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2001,2002,2003,2011,2014
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/MeasurementComponents/FJones.h>
28 : //#include <synthesis/MeasurementComponents/CalCorruptor.h>
29 :
30 : #include <msvis/MSVis/VisibilityIterator2.h>
31 : #include <msvis/MSVis/VisBuffer2.h>
32 : #include <msvis/MSVis/VisBuffer.h> // still used in apply context
33 : #include <casacore/ms/MeasurementSets/MSColumns.h>
34 : #include <casacore/casa/BasicMath/Math.h>
35 :
36 : #include <casacore/casa/Arrays/ArrayMath.h>
37 : #include <casacore/casa/Arrays/MatrixMath.h>
38 : #include <casacore/casa/BasicSL/String.h>
39 : #include <casacore/casa/Utilities/Assert.h>
40 : #include <casacore/casa/Utilities/GenSort.h>
41 : #include <casacore/casa/Exceptions/Error.h>
42 : #include <casacore/casa/OS/Memory.h>
43 : #include <casacore/casa/System/Aipsrc.h>
44 : #include <casacore/images/Images/PagedImage.h>
45 :
46 : #include <sstream>
47 :
48 : #include <casacore/measures/Measures/MDirection.h>
49 : #include <casacore/measures/Measures/MEpoch.h>
50 : #include <casacore/measures/Measures/MeasTable.h>
51 :
52 : #include <casacore/casa/Logging/LogMessage.h>
53 : #include <casacore/casa/Logging/LogSink.h>
54 :
55 : using namespace casacore;
56 : namespace casa { //# NAMESPACE CASA - BEGIN
57 :
58 :
59 : // **********************************************************
60 : // FJones
61 : //
62 :
63 0 : FJones::FJones(VisSet& vs) :
64 : VisCal(vs),
65 : VisMueller(vs),
66 : SolvableVisJones(vs),
67 : tectype_(""),
68 : mframe_(),
69 : emm_(NULL),
70 : ionhgt_(450.,"km"),
71 : tecimname_("unspecified"),
72 : za_(),
73 0 : radper_(2.3649) // rad*Hz2/G
74 : {
75 0 : if (prtlev()>2) cout << "FJones::FJones(vs)" << endl;
76 :
77 : // Prepare zenith angle storage
78 0 : za().resize(nAnt());
79 0 : za().set(0.0);
80 :
81 : // Prepare line-of-sight B field (G) storage
82 0 : BlosG_.resize(nAnt());
83 0 : BlosG_.set(0.0);
84 :
85 0 : }
86 :
87 0 : FJones::FJones(String msname,Int MSnAnt,Int MSnSpw) :
88 : VisCal(msname,MSnAnt,MSnSpw),
89 : VisMueller(msname,MSnAnt,MSnSpw),
90 : SolvableVisJones(msname,MSnAnt,MSnSpw),
91 : tectype_(""),
92 : mframe_(),
93 : emm_(NULL),
94 : ionhgt_(450.,"km"),
95 : tecimname_("unspecified"),
96 : za_(),
97 0 : radper_(2.3649) // rad*Hz2/G
98 : {
99 0 : if (prtlev()>2) cout << "FJones::FJones(msname,MSnAnt,MSnSpw)" << endl;
100 :
101 : // Prepare zenith angle storage
102 0 : za().resize(MSnAnt);
103 0 : za().set(0.0);
104 :
105 : // Prepare line-of-sight B field (G) storage
106 0 : BlosG_.resize(MSnAnt);
107 0 : BlosG_.set(0.0);
108 :
109 0 : }
110 :
111 1 : FJones::FJones(const MSMetaInfoForCal& msmc) :
112 : VisCal(msmc),
113 : VisMueller(msmc),
114 : SolvableVisJones(msmc),
115 : tectype_(""),
116 : mframe_(),
117 : emm_(NULL),
118 : ionhgt_(450.,"km"),
119 : tecimname_("unspecified"),
120 : za_(),
121 1 : radper_(2.3649) // rad*Hz2/G
122 : {
123 1 : if (prtlev()>2) cout << "FJones::FJones(msmc)" << endl;
124 :
125 : // Prepare zenith angle storage
126 1 : za().resize(nAnt());
127 1 : za().set(0.0);
128 :
129 : // Prepare line-of-sight B field (G) storage
130 1 : BlosG_.resize(nAnt());
131 1 : BlosG_.set(0.0);
132 :
133 1 : }
134 :
135 2 : FJones::~FJones() {
136 1 : if (prtlev()>2) cout << "FJones::~FJones()" << endl;
137 :
138 1 : if (emm_)
139 0 : delete emm_;
140 :
141 2 : }
142 :
143 : // Enforce calwt=F and spwmap=[0]
144 0 : void FJones::setApply(const Record& apply) {
145 :
146 0 : cout << "FJones (iononsphere) apply now provisionally includes disp. delay." << endl;
147 :
148 : // Force spwmap to all 0 (FJones pars are not spw-dep
149 : // NB: this is required before calling parent, which sets
150 : // up spwmap-dep interpolation
151 0 : logSink() << " (" << this->typeName()
152 0 : << ": Overriding with spwmap=[0] since " << this->typeName()
153 : << " is not spw-dependent)"
154 0 : << LogIO::POST;
155 0 : spwMap().assign(Vector<Int>(1,0));
156 :
157 : // Remove spwmap from input Record, and pass to generic coe
158 0 : Record newapply;
159 0 : newapply=apply;
160 0 : if (newapply.isDefined("spwmap"))
161 0 : newapply.removeField("spwmap");
162 :
163 : // Do conventional things
164 0 : SolvableVisJones::setApply(newapply);
165 :
166 : // Extract ionhgt_ from table header, if available
167 0 : if (ct_) {
168 0 : const TableRecord& tr(ct_->keywordSet());
169 0 : if (tr.isDefined("IonosphereHeight(km)"))
170 0 : ionhgt_=Quantity(tr.asDouble("IonosphereHeight(km)"),"km");
171 : }
172 0 : logSink() << "Using ionosphere height = " << ionhgt_ << LogIO::POST;
173 :
174 :
175 : // Enforce no weight calibration (this is a phase-like correction)
176 0 : calWt()=false;
177 :
178 0 : }
179 :
180 :
181 0 : String FJones::applyinfo() {
182 :
183 0 : ostringstream o;
184 0 : o << typeName();
185 0 : o << ": table=" << calTableName();
186 :
187 0 : return String(o);
188 :
189 : }
190 :
191 1 : void FJones::setSpecify(const Record& specify) {
192 :
193 : // extract TEC retrieval mode... (madrigal, ionics, etc.)
194 1 : if (specify.isDefined("caltype")) {
195 1 : tectype_=upcase(specify.asString("caltype"));
196 : }
197 :
198 : // We are importing from a TEC image:
199 1 : if (tectype_!="TECIM")
200 0 : throw(AipsError("Unrecognized TEC mode: "+tectype_));
201 :
202 : // Extract TEC image name from infile
203 : // TBD: move this to private data and setSpecify
204 1 : if (specify.isDefined("infile"))
205 1 : tecimname_=specify.asString("infile");
206 :
207 1 : if (tecimname_=="")
208 0 : throw(AipsError("A valid TEC image file must be specified!"));
209 :
210 1 : if ( !Table::isReadable(tecimname_) )
211 0 : throw(AipsError("Cannot find expected TEC image: "+tecimname_));
212 : else
213 1 : logSink() << "Found required TEC image: "+tecimname_ << LogIO::POST;
214 :
215 : // extract intended cal table name...
216 1 : if (specify.isDefined("caltable")) {
217 1 : calTableName()=specify.asString("caltable");
218 :
219 1 : if (Table::isReadable(calTableName()))
220 : logSink() << "FYI: We are going to overwrite an existing CalTable: "
221 0 : << calTableName()
222 0 : << LogIO::POST;
223 : }
224 : // we are creating a table from scratch
225 1 : logSink() << "Creating " << typeName()
226 : << " table."
227 1 : << LogIO::POST;
228 :
229 : // Extract user's ionhgt_, if specified
230 1 : Vector<Double> parameters;
231 1 : if (specify.isDefined("parameter")) {
232 1 : parameters=specify.asArrayDouble("parameter");
233 1 : if (parameters.nelements()==1) {
234 0 : ionhgt_=Quantity(parameters[0],"km");
235 0 : logSink() << "Using user-specified ionosphere height = " << ionhgt_ << LogIO::POST;
236 : }
237 : else
238 1 : logSink() << "Using default ionosphere height ("+tectype_+") = " << ionhgt_ << LogIO::POST;
239 : }
240 : else
241 0 : logSink() << "Using default ionosphere height ("+tectype_+") = " << ionhgt_ << LogIO::POST;
242 :
243 : // Create a new caltable to fill up
244 1 : createMemCalTable();
245 :
246 : // Add ionhgt_ to table keyword
247 1 : if (ct_) {
248 :
249 1 : TableRecord& tr(ct_->rwKeywordSet());
250 1 : tr.define("IonosphereHeight(km)",ionhgt_.getValue("km"));
251 :
252 : }
253 :
254 :
255 : // Set up solveAllRPar
256 1 : initSolvePar();
257 :
258 : // Initialize parameter space
259 1 : solveAllRPar()=0.0;
260 1 : solveAllParOK()=true;
261 1 : solveAllParErr()=0.0; // what to use here?
262 1 : solveAllParSNR()=100.0; // what to use here?
263 :
264 1 : }
265 :
266 1 : void FJones::specify(const Record&) {
267 :
268 : //cout << "FJones-specific specify: " << tectype_ << endl;
269 :
270 1 : logSink() << "Deriving evolving line-of-sight TEC for whole dataset (this may take a little while)" << LogIO::POST;
271 :
272 2 : PagedImage<Float> tecim(tecimname_);
273 2 : CoordinateSystem teccs(tecim.coordinates());
274 :
275 : // Specific sort columns for FJones
276 2 : Block<Int> cols(5);
277 1 : cols[0]=MS::ARRAY_ID;
278 1 : cols[1]=MS::SCAN_NUMBER;
279 1 : cols[2]=MS::TIME;
280 1 : cols[3]=MS::FIELD_ID; // NB: assumed unique per timestamp!
281 1 : cols[4]=MS::DATA_DESC_ID; // required but irrelevant in this context
282 :
283 2 : MeasurementSet ms(msName());
284 2 : vi::SortColumns sc(cols);
285 2 : vi::VisibilityIterator2 vi2(ms,sc,true);
286 1 : vi::VisBuffer2 *vb = vi2.getVisBuffer();
287 1 : vi2.originChunks();
288 1 : vi2.origin();
289 :
290 2 : MSAntennaColumns antcol(ms.antenna());
291 :
292 : // Measures stuff for ionosphere coords
293 2 : MeasFrame mframe;
294 3 : MEpoch when(MVEpoch(Quantity(vb->time()(0),"s")),MEpoch::Ref(MEpoch::UTC));
295 1 : mframe.set(when);
296 2 : MPosition where;
297 1 : MeasTable::Observatory(where,"VLA");
298 2 : MVPosition vla=where.getValue();
299 1 : mframe.set(where);
300 2 : MDirection::Ref mdref=vb->phaseCenter().getRef();
301 1 : EarthMagneticMachine* emm = new EarthMagneticMachine(mdref,ionhgt_,mframe);
302 :
303 : // Workspace for geometry/TEC retrieval/interp
304 2 : MDirection md;
305 2 : MVPosition ionpos;
306 :
307 2 : Vector<Double> cf(3);
308 2 : Vector<Int> c(3);
309 2 : IPosition corner(3,0,0,0),lastcorner(3,0,0,0),twos(3,2,2,2);
310 :
311 1 : Double day0(-999.0);
312 2 : Vector<Double> refv(teccs.referenceValue());
313 1 : if (refv[2]>0.0)
314 1 : day0=0.0;
315 :
316 2 : Vector<Double> wc(3,0.0);
317 1 : Double &lng(wc(0)), &lat(wc(1));
318 :
319 : // wc0 evolves in iant loop below to place relevant factors in focus for calculations
320 2 : Vector<Double> wc0(3,0.0);
321 1 : wc0=teccs.toWorld(lastcorner);
322 : // Scalar referenes into wc0 for calculation below
323 : // Memory of wc0 can't be changed!
324 1 : Double &lng0(wc0(0)), &lat0(wc0(1)), &time0(wc0(2));
325 :
326 : // Fractional location in intervals
327 1 : Double dlng(0.0),dlat(0.0),dtime(0.0);
328 :
329 2 : Vector<Double> inc(teccs.increment());
330 1 : Double& dlng0(inc(0)), dlat0(inc(1)), dtime0(inc(2));
331 :
332 : // tecube evolves in iant loop below to place relevant factor in focus for calculations
333 2 : Cube<Float> tecube(2,2,2,0.0); // lng, lat, time
334 1 : tecube=tecim.getSlice(lastcorner,lastcorner+2,false);
335 : // Scalar referenes into tecube for calculation below
336 : // Memory of tecube can't be changed!
337 1 : Float &a0(tecube(0,0,0)), &a1(tecube(1,0,0)), &a2(tecube(1,1,0)), &a3(tecube(0,1,0));
338 1 : Float &b0(tecube(0,0,1)), &b1(tecube(1,0,1)), &b2(tecube(1,1,1)), &b3(tecube(0,1,1));
339 :
340 1 : Float t0(0.0), u0(0.0), v0(0.0);
341 1 : Float teca(0.0), tecb(0.0), ztec(4e18), tec(4e18);
342 :
343 1 : Double lasttime(-1);
344 1 : Int lastscan(-999);
345 1 : Int isol(0);
346 167 : for (vi2.originChunks(); vi2.moreChunks(); vi2.nextChunk()) {
347 :
348 332 : for (vi2.origin(); vi2.more(); vi2.next()) {
349 :
350 166 : refTime()=vb->time()(0);
351 166 : if (day0<0.0)
352 0 : day0=86400.0*floor(refTime()/86400.0);
353 :
354 166 : wc(2)=refTime()-day0;
355 :
356 166 : currScan()=vb->scan()(0);
357 :
358 : // only when the timestamp changes
359 166 : if (refTime()!=lasttime) {
360 :
361 : // Set meta-info (TBD: use syncMeta?)
362 83 : currSpw()=0; // only spw=0 for FJones (TEC is not freq-dep)
363 83 : currObs()=vb->observationId()(0);
364 83 : currScan()=vb->scan()(0);
365 83 : currField()=vb->fieldId()(0);
366 :
367 : // Calculate zenith angle for current time/direction
368 83 : za().resize(nAnt());
369 83 : Vector<MDirection> antazel(vb->azel(refTime()));
370 83 : Double* a=za().data();
371 1660 : for (Int iant=0;iant<nAnt();++iant,++a)
372 1577 : (*a)=C::pi_2 - antazel(iant).getAngle().getValue()(1);
373 :
374 : // This vb's direction measure...
375 83 : md=vb->phaseCenter();
376 :
377 : // Update current time in frame
378 83 : when.set(Quantity(refTime(),"s"));
379 83 : mframe.set(when);
380 :
381 1660 : for (Int iant=0;iant<nAnt();++iant) {
382 1577 : mframe.set(antcol.positionMeas()(iant));
383 :
384 : // Force measures calculation...
385 1577 : emm->set(mframe); // ensure frame is updating...
386 1577 : emm->calculate(md.getValue());
387 1577 : ionpos=emm->getPosition();
388 1577 : wc(0)=ionpos.getLong("deg").getValue();
389 1577 : wc(1)=ionpos.getLat("deg").getValue();
390 :
391 : // Find corner of 2x2x2 cube that we need
392 1577 : teccs.toPixel(cf,wc); // fractional pixel coord
393 1577 : convertArray(c,cf); // integer pixel coord
394 1577 : corner=IPosition(c); // as IPosition
395 1577 : if (corner!=lastcorner) { // If new
396 :
397 : // (2020/12/17 gmoellen)
398 : // This fails: (changes storage address in wc0!)
399 : //wc0=teccs.toWorld(corner);
400 :
401 : // (2020/12/17, gmoellen)
402 : // This works: (wc0's storage stays the same, values copied into it)
403 : // Scalar references used below still work!
404 422 : Vector<Double> toW(teccs.toWorld(corner));
405 211 : wc0=toW;
406 :
407 : // (2020/12/17 gmoellen)
408 : // This fails: (changes storage address of tecube)
409 : // tecube=tecim.getSlice(corner,twos,false);
410 :
411 : // (2020/12/17 gmoellen)
412 : // This works: (tecube's storage stays the same, values copied into it)
413 : // Scalar references used below still work!
414 422 : Cube<Float> newtecube(tecim.getSlice(corner,twos,false));
415 211 : tecube=newtecube;
416 :
417 211 : lastcorner=corner;
418 :
419 : /*
420 : cout << "***tecube = " << tecube << endl;
421 : cout << a0 << " " << a1 << " " << a2 << " " << a3
422 : << b0 << " " << b1 << " " << b2 << " " << b3 << endl;
423 :
424 :
425 : if (iant<4)
426 : cout << "**";
427 : }
428 : else
429 : if (iant<4)
430 : cout << " ";
431 : */
432 : }
433 :
434 1577 : dlng=lng-lng0;
435 1577 : dlat=lat-lat0;
436 1577 : dtime=wc(2)-time0;
437 :
438 1577 : t0=dlng/dlng0;
439 1577 : u0=dlat/dlat0;
440 1577 : v0=dtime/dtime0;
441 :
442 : if (false && iant==0) {
443 : cout << dlng << " " << dlng0 << " " << t0 << endl;
444 : cout << dlat << " " << dlat0 << " " << u0 << endl;
445 : cout << dtime << " " << dtime0 << " " << v0 << endl;
446 : }
447 :
448 : // Zenith TEC at pierce point
449 1577 : teca=(1.-t0)*(1.-u0)*a0 + t0*(1.-u0)*a1 + t0*u0*a2 + (1.-t0)*u0*a3;
450 1577 : tecb=(1.-t0)*(1.-u0)*b0 + t0*(1.-u0)*b1 + t0*u0*b2 + (1.-t0)*u0*b3;
451 1577 : ztec=(1.-v0)*teca + v0*tecb;
452 :
453 : // LOS tecat pierce point:
454 1577 : tec=ztec/cos(za()(iant));
455 :
456 : // retrieve time-, direction-, and antenna- specific vertical TEC...
457 1577 : solveAllRPar()(0,0,iant)=tec*1.e16; // in e-/m2
458 :
459 1577 : if (iant==0 && currScan()!=lastscan) {
460 5 : lastscan=currScan();
461 :
462 : /*
463 : MVDirection azel=vb->azel0(refTime()).getValue();
464 : cout.precision(12);
465 : cout << isol << " "
466 : << " obs=" << currObs()
467 : << " scan=" << currScan()
468 : << " fld=" << currField()
469 : << " time=" << wc(2)
470 : << " ant=" << antcol.name()(iant)+":"+antcol.station()(iant)
471 : << " long/lat=" << lng << "/" << lat
472 : //<< " pixel=" << cf << " " << c
473 : << " za=" << za()(iant)*180/C::pi
474 : << " ztec=" // << teca << " " << tecb << " "
475 : << ztec
476 : << " --> lostec=" << tec
477 : << endl;
478 : */
479 :
480 : }
481 : }
482 :
483 : // store in the caltable
484 83 : keepNCT();
485 :
486 83 : ++isol;
487 :
488 : }
489 :
490 166 : lasttime=refTime();
491 :
492 : }
493 : }
494 1 : }
495 :
496 : // FJones needs to know pol basis, and some geometry
497 0 : void FJones::syncMeta(const VisBuffer& vb) {
498 :
499 : // Call parent (sets currTime())
500 0 : SolvableVisJones::syncMeta(vb);
501 :
502 : // Basis
503 0 : if (vb.corrType()(0)==5) // Circulars
504 0 : pjonestype_=Jones::Diagonal;
505 0 : else if (vb.corrType()(0)==9) // Linears
506 0 : pjonestype_=Jones::General;
507 :
508 : // Geometry
509 0 : phasedir_p=vb.msColumns().field().phaseDirMeas(currField());
510 0 : antpos_p.reference(vb.msColumns().antenna().positionMeas());
511 :
512 0 : }
513 :
514 : // FJones needs to know pol basis, and some geometry
515 0 : void FJones::syncMeta2(const vi::VisBuffer2& vb) {
516 :
517 : // Call parent (sets currTime())
518 0 : SolvableVisJones::syncMeta2(vb);
519 :
520 : // Basis
521 0 : if (vb.correlationTypes()(0)==5) // Circulars
522 0 : pjonestype_=Jones::Diagonal;
523 0 : else if (vb.correlationTypes()(0)==9) // Linears
524 0 : pjonestype_=Jones::General;
525 :
526 0 : phasedir_p=vb.subtableColumns().field().phaseDirMeas(currField());
527 0 : antpos_p.reference(vb.subtableColumns().antenna().positionMeas());
528 :
529 0 : }
530 :
531 0 : void FJones::calcPar() {
532 :
533 0 : if (prtlev()>6) cout << " FJones::calcPar()" << endl;
534 :
535 : // set time in mframe_
536 0 : MEpoch epoch(Quantity(currTime(),"s"));
537 0 : epoch.setRefString("UTC");
538 0 : mframe_.set(epoch);
539 :
540 : // Set this antenna's position in mframe_
541 0 : const MPosition& antpos0 = antpos_p(0);
542 0 : mframe_.set(antpos0);
543 :
544 : // set direction ref in emm
545 0 : const MDirection::Ref phasedirR=phasedir_p.getRef();
546 0 : const MVDirection phasedirV=phasedir_p.getValue();
547 :
548 : // Construct workable EMM (if not yet done)
549 0 : if (!emm_)
550 0 : emm_ = new EarthMagneticMachine(phasedirR,ionhgt_,mframe_);
551 :
552 : // Calculate ant-dep mag field
553 0 : BlosG_.resize(nAnt());
554 :
555 : // loop over antennas
556 0 : for (Int iant=0;iant<nAnt();++iant) {
557 :
558 : // Set this antenna's position in mframe_
559 0 : const MPosition& antpos = antpos_p(iant);
560 0 : mframe_.resetPosition(antpos);
561 0 : emm_->set(mframe_); // seems to be needed to force new position in emm calculations
562 :
563 : // calculate this ant's field
564 0 : emm_->calculate(phasedirV);
565 :
566 0 : BlosG_(iant)=(emm_->getLOSField("G").getValue()); // Sign?
567 :
568 : }
569 :
570 : // cout.precision(16);
571 : // cout << "BlosG_ = " << BlosG_ << endl;
572 :
573 : // Get current zenith tec(t)
574 0 : SolvableVisCal::calcPar();
575 :
576 : // cout << "currRPar() = " << currRPar() << endl;
577 :
578 :
579 : // Pars now valid, matrices not yet
580 0 : validateP();
581 0 : invalidateJ(); // Force new calculation of matrix elements
582 :
583 0 : }
584 :
585 :
586 0 : void FJones::calcAllJones() {
587 :
588 : using casacore::operator*;
589 :
590 0 : if (prtlev()>6) cout << " FJones::calcAllJones()" << endl;
591 :
592 : // Nominally no ionosphere
593 0 : currJElem()=Complex(1.0);
594 0 : currJElemOK().set(true);
595 :
596 : // cout << currSpw() << " "
597 : // << currTime() << " "
598 : // << "currJElem().shape() = " << currJElem().shape() << endl;
599 :
600 0 : Complex* J=currJElem().data();
601 0 : Float* lostec=currRPar().data();
602 0 : Bool* lostecok=currParOK().data();
603 : Double f,rotpers2,tec,del,rot;
604 0 : Complex cdel;
605 :
606 :
607 0 : for (Int iant=0; iant<nAnt(); ++iant,++lostec,++lostecok) {
608 0 : if ((*lostecok)) {
609 0 : tec = Double(*lostec);
610 0 : rotpers2 = radper_*tec*BlosG_(iant);
611 :
612 0 : for (Int ich=0;ich<nChanMat();++ich) {
613 0 : f=currFreq()(ich)*1.0e9; // Hz
614 0 : del = -8.4483e-7*tec/f; // opposite sign cf rot!
615 0 : rot = rotpers2/f/f;
616 :
617 0 : switch (jonesType()) {
618 : // Circular version:
619 0 : case Jones::Diagonal: {
620 0 : J[0]=Complex(cos(del+rot),sin(del+rot));
621 0 : J[1]=Complex(cos(del-rot),sin(del-rot));
622 0 : J+=2;
623 0 : break;
624 : }
625 : // Linear version:
626 0 : case Jones::General: {
627 0 : cdel=Complex(cos(del),sin(del));
628 0 : J[0]=J[3]=cdel*cos(-rot);
629 0 : J[1]=cdel*sin(-rot);
630 0 : J[2]=-J[1];
631 0 : J+=4;
632 0 : break;
633 : }
634 0 : default:
635 0 : throw(AipsError("PJones doesn't know if it is Diag (Circ) or General (Lin)"));
636 : break;
637 :
638 : }
639 :
640 : }
641 : }
642 :
643 : }
644 : /*
645 : cout << "tec = " << tec << endl;
646 : cout << "rot = " << rot << " " << rotpers2 << endl;
647 : cout << "currJElem() = " << currJElem() << endl;
648 : */
649 0 : validateJ();
650 0 : }
651 :
652 :
653 :
654 : } //# NAMESPACE CASA - END
|
