Line data Source code
1 : // -*- C++ -*-
2 : //# AWConvFuncEPJones.cc: Implementation of the AWConvFuncEPJones class
3 : //# Copyright (C) 1997,1998,1999,2000,2001,2002,2003
4 : //# Associated Universities, Inc. Washington DC, USA.
5 : //#
6 : //# This library is free software; you can redistribute it and/or modify it
7 : //# under the terms of the GNU Library General Public License as published by
8 : //# the Free Software Foundation; either version 2 of the License, or (at your
9 : //# option) any later version.
10 : //#
11 : //# This library is distributed in the hope that it will be useful, but WITHOUT
12 : //# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
13 : //# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
14 : //# License for more details.
15 : //#
16 : //# You should have received a copy of the GNU Library General Public License
17 : //# along with this library; if not, write to the Free Software Foundation,
18 : //# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
19 : //#
20 : //# Correspondence concerning AIPS++ should be addressed as follows:
21 : //# Internet email: aips2-request@nrao.edu.
22 : //# Postal address: AIPS++ Project Office
23 : //# National Radio Astronomy Observatory
24 : //# 520 Edgemont Road
25 : //# Charlottesville, VA 22903-2475 USA
26 : //#
27 : //# $Id$
28 : //
29 : #include <synthesis/TransformMachines/AWConvFuncEPJones.h>
30 : #include <synthesis/TransformMachines/SynthesisError.h>
31 : #include <casacore/images/Images/ImageInterface.h>
32 : #include <synthesis/TransformMachines/Utils.h>
33 : #include <synthesis/TransformMachines/BeamCalc.h>
34 : #include <synthesis/TransformMachines/CFStore.h>
35 : #include <synthesis/TransformMachines/CFStore2.h>
36 : #include <synthesis/TransformMachines/PSTerm.h>
37 : #include <synthesis/TransformMachines/WTerm.h>
38 : #include <synthesis/TransformMachines/ATerm.h>
39 : #include <synthesis/TransformMachines/VLACalcIlluminationConvFunc.h>
40 : #include <synthesis/TransformMachines/ConvolutionFunction.h>
41 : #include <synthesis/TransformMachines/PolOuterProduct.h>
42 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
43 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
44 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
45 : #include <casacore/lattices/LatticeMath/LatticeFFT.h>
46 : #include <casacore/casa/Utilities/CompositeNumber.h>
47 : #include <casacore/measures/Measures/MeasTable.h>
48 : #include <ostream>
49 :
50 : #define MAX_FREQ 1e30
51 :
52 : using namespace casacore;
53 : namespace casa{
54 : //
55 : //----------------------------------------------------------------------
56 : //
57 0 : AWConvFuncEPJones& AWConvFuncEPJones::operator=(const AWConvFuncEPJones& other)
58 : {
59 0 : if(this!=&other)
60 : {
61 0 : AWConvFunc::operator=(other);
62 0 : imageDC_p = other.imageDC_p;
63 0 : imageObsInfo_p = other.imageObsInfo_p;
64 : }
65 0 : return *this;
66 : }
67 : //
68 : //----------------------------------------------------------------------
69 : // Find the offset between the VB and the image phase center
70 : //
71 0 : Vector<Double> AWConvFuncEPJones::findPointingOffset(const ImageInterface<Complex>& image,
72 : const VisBuffer& vb)
73 : {
74 0 : storeImageParams(image,vb);
75 : //where in the image in pixels is this pointing
76 0 : pixFieldGrad_p.resize(2);
77 0 : toPix(vb);
78 0 : pixFieldGrad_p=thePix_p;
79 :
80 0 : MDirection fieldDir=direction1_p;
81 : //shift from center
82 0 : pixFieldGrad_p(0) = pixFieldGrad_p(0) - Double(nx_p / 2);
83 0 : pixFieldGrad_p(1) = pixFieldGrad_p(1) - Double(ny_p / 2);
84 :
85 : // Int convSampling=aTerm_p->getOversampling();
86 0 : Int convSampling=getOversampling(*psTerm_p,*wTerm_p,*aTerm_p);
87 :
88 : //phase gradient per pixel to apply
89 0 : pixFieldGrad_p(0) = -pixFieldGrad_p(0)*2.0*C::pi/Double(nx_p)/Double(convSampling);
90 0 : pixFieldGrad_p(1) = -pixFieldGrad_p(1)*2.0*C::pi/Double(ny_p)/Double(convSampling);
91 :
92 0 : return pixFieldGrad_p;
93 : }
94 : //
95 : //----------------------------------------------------------------------
96 : //
97 0 : void AWConvFuncEPJones::makeConvFunction(const ImageInterface<Complex>& image,
98 : const VisBuffer& vb,
99 : const Int wConvSize,
100 : const CountedPtr<PolOuterProduct>& pop,
101 : const Float pa,
102 : const Float dpa,
103 : const Vector<Double>& uvScale, const Vector<Double>& uvOffset,
104 : const Matrix<Double>& spwFreqSel,
105 : CFStore2& cfs,
106 : CFStore2& cfwts,
107 : Bool fillCF)
108 : {
109 0 : findPointingOffset(image,vb);
110 0 : AWConvFunc::makeConvFunction(image,vb,wConvSize,pop,pa,dpa,uvScale,uvOffset,spwFreqSel,cfs,cfwts,fillCF);
111 0 : }
112 : //
113 : //----------------------------------------------------------------------
114 : //
115 0 : void AWConvFuncEPJones::toPix(const VisBuffer& vb)
116 : {
117 0 : thePix_p.resize(2);
118 :
119 0 : if(dc_p.directionType() != MDirection::castType(vb.direction1()(0).getRef().getType())){
120 : //pointToPix_p.setModel(theDir);
121 :
122 0 : MEpoch timenow(Quantity(vb.time()(0), timeUnit_p), timeMType_p);
123 : //cout << "Ref " << vb.direction1()(0).getRefString() << " ep "
124 : //<< timenow.getRefString() << " time " <<
125 : //MVTime(timenow.getValue().getTime()).string(MVTime::YMD) <<
126 : //endl;
127 0 : pointFrame_p.resetEpoch(timenow);
128 : //////////////////////////
129 : //pointToPix holds pointFrame_p by reference...
130 : //thus good to go for conversion
131 0 : direction1_p=pointToPix_p(vb.direction1()(0));
132 0 : direction2_p=pointToPix_p(vb.direction2()(0));
133 0 : dc_p.toPixel(thePix_p, direction1_p);
134 :
135 : }
136 : else{
137 0 : direction1_p=vb.direction1()(0);
138 0 : direction2_p=vb.direction2()(0);
139 0 : dc_p.toPixel(thePix_p, vb.direction1()(0));
140 : }
141 0 : }
142 : //
143 : //----------------------------------------------------------------------
144 : //
145 0 : void AWConvFuncEPJones::storeImageParams(const ImageInterface<Complex>& iimage,
146 : const VisBuffer& vb){
147 : //image signature changed...rather simplistic for now
148 0 : if((iimage.shape().product() != nx_p*ny_p*nchan_p*npol_p) || nchan_p < 1){
149 0 : csys_p=iimage.coordinates();
150 0 : Int coordIndex=csys_p.findCoordinate(Coordinate::DIRECTION);
151 0 : AlwaysAssert(coordIndex>=0, AipsError);
152 0 : directionIndex_p=coordIndex;
153 0 : dc_p=csys_p.directionCoordinate(directionIndex_p);
154 0 : ObsInfo imInfo=csys_p.obsInfo();
155 0 : String tel= imInfo.telescope();
156 0 : MPosition pos;
157 0 : if (vb.msColumns().observation().nrow() > 0) {
158 0 : tel = vb.msColumns().observation().telescopeName()(vb.msColumns().observationId()(0));
159 : }
160 0 : if (tel.length() == 0 || !tel.contains("VLA") ||
161 0 : !MeasTable::Observatory(pos,tel)) {
162 : // unknown observatory, use first antenna
163 0 : Int ant1=vb.antenna1()(0);
164 0 : pos=vb.msColumns().antenna().positionMeas()(ant1);
165 : }
166 : //cout << "TELESCOPE " << tel << endl;
167 : //Store this to build epochs via the time access of visbuffer later
168 0 : timeMType_p=MEpoch::castType(vb.msColumns().timeMeas()(0).getRef().getType());
169 0 : timeUnit_p=Unit(vb.msColumns().timeMeas().measDesc().getUnits()(0).getName());
170 : // timeUnit_p=Unit("s");
171 : //cout << "UNIT " << timeUnit_p.getValue() << " name " << timeUnit_p.getName() << endl;
172 0 : pointFrame_p=MeasFrame(imInfo.obsDate(), pos);
173 0 : MDirection::Ref elRef(dc_p.directionType(), pointFrame_p);
174 : //For now we set the conversion from this direction
175 0 : pointToPix_p=MDirection::Convert( MDirection(), elRef);
176 0 : nx_p=iimage.shape()(coordIndex);
177 0 : ny_p=iimage.shape()(coordIndex+1);
178 : // pointingPix_p.resize(nx_p, ny_p);
179 : // pointingPix_p.set(false);
180 0 : coordIndex=csys_p.findCoordinate(Coordinate::SPECTRAL);
181 0 : Int pixAxis=csys_p.pixelAxes(coordIndex)[0];
182 0 : nchan_p=iimage.shape()(pixAxis);
183 0 : coordIndex=csys_p.findCoordinate(Coordinate::STOKES);
184 0 : pixAxis=csys_p.pixelAxes(coordIndex)[0];
185 0 : npol_p=iimage.shape()(pixAxis);
186 : }
187 :
188 0 : }
189 : }
|
