Line data Source code
1 : //# tGridFT.cc: Tests the Synthesis model data serving
2 : //# Copyright (C) 2016
3 : //# Associated Universities, Inc. Washington DC, USA.
4 : //#
5 : //# This program is free software; you can redistribute it and/or modify it
6 : //# under the terms of the GNU General Public License as published by the Free
7 : //# Software Foundation; either version 2 of the License, or (at your option)
8 : //# any later version.
9 : //#
10 : //# This program 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 General Public License for
13 : //# more details.
14 : //#
15 : //# You should have received a copy of the GNU General Public License along
16 : //# with this program; if not, write to the Free Software Foundation, Inc.,
17 : //# 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 : //# $Id$
27 : #include <casacore/measures/Measures/Stokes.h>
28 : #include <casacore/coordinates/Coordinates/CoordinateSystem.h>
29 : #include <casacore/coordinates/Coordinates/DirectionCoordinate.h>
30 : #include <casacore/coordinates/Coordinates/SpectralCoordinate.h>
31 : #include <casacore/coordinates/Coordinates/StokesCoordinate.h>
32 : #include <casacore/coordinates/Coordinates/Projection.h>
33 : #include <casacore/casa/Arrays/ArrayMath.h>
34 : #include <casacore/images/Images/ImageInterface.h>
35 : #include <casacore/images/Images/PagedImage.h>
36 : #include <casacore/images/Images/TempImage.h>
37 : #include <components/ComponentModels/ComponentList.h>
38 : #include <components/ComponentModels/ComponentShape.h>
39 : #include <components/ComponentModels/Flux.h>
40 : #include <casacore/tables/TaQL/ExprNode.h>
41 : #include <casacore/measures/Measures/MeasTable.h>
42 : #include <casacore/ms/MSSel/MSSelection.h>
43 : #include <synthesis/TransformMachines2/FTMachine.h>
44 : #include <synthesis/TransformMachines2/GridFT.h>
45 : #include <synthesis/TransformMachines2/SetJyGridFT.h>
46 : #include <synthesis/TransformMachines2/WProjectFT.h>
47 : #include <synthesis/TransformMachines2/SimpleComponentFTMachine.h>
48 : #include <msvis/MSVis/VisImagingWeight.h>
49 : #include <msvis/MSVis/VisibilityIterator2.h>
50 : #include <msvis/MSVis/VisBuffer2.h>
51 : #include <casacore/casa/namespace.h>
52 : #include <casacore/casa/OS/Directory.h>
53 : #include <casacore/casa/Utilities/Regex.h>
54 : #include "MakeMS.h"
55 :
56 : using namespace std;
57 : using namespace casa;
58 : using namespace casacore;
59 : using namespace casa::refim;
60 : using namespace casa::test;
61 :
62 3 : refim::FTMachine* makeFTM(String ftmachine, MPosition loc, MDirection dir){
63 3 : if(ftmachine=="GridFT")
64 1 : return new refim::GridFT(1000000, 16, "SF", loc, 1.0, false, false);
65 2 : if(ftmachine=="SetJyGridFT")
66 1 : return new refim::SetJyGridFT(1000000, 16, "SF", loc, dir, 1.0, false, false);
67 1 : if(ftmachine=="WProjectFT")
68 1 : return new refim::WProjectFT(50, loc, 1000000, 16);
69 0 : return nullptr;
70 :
71 : }
72 :
73 1 : Int main(/*int argc, char **argv*/){
74 : try{
75 2 : MDirection thedir(Quantity(20.0, "deg"), Quantity(20.0, "deg"));
76 1 : String msname("Test.ms");
77 1 : MakeMS::makems(msname, thedir);
78 1 : MeasurementSet thems(msname, Table::Update);
79 1 : thems.markForDelete();
80 2 : vi::VisibilityIterator2 vi2(thems,vi::SortColumns(),true);
81 1 : vi::VisBuffer2 *vb=vi2.getVisBuffer();
82 1 : VisImagingWeight viw("natural");
83 1 : vi2.useImagingWeight(viw);
84 1 : ComponentList cl;
85 1 : SkyComponent otherPoint(ComponentType::POINT);
86 1 : otherPoint.flux() = Flux<Double>(6.66e-2, 0.0, 0.0, 0.00000);
87 1 : otherPoint.shape().setRefDirection(thedir);
88 1 : cl.add(otherPoint);
89 1 : Matrix<Double> xform(2,2);
90 1 : xform = 0.0;
91 1 : xform.diagonal() = 1.0;
92 2 : DirectionCoordinate dc(MDirection::J2000, Projection::SIN, Quantity(20.0,"deg"), Quantity(20.0, "deg"),
93 2 : Quantity(0.5, "arcsec"), Quantity(0.5,"arcsec"),
94 0 : xform, 50.0, 50.0, 999.0,
95 4 : 999.0);
96 1 : Vector<Int> whichStokes(1, Stokes::I);
97 1 : StokesCoordinate stc(whichStokes);
98 1 : SpectralCoordinate spc(MFrequency::LSRK, 1.5e9, 1e6, 0.0 , 1.420405752E9);
99 1 : CoordinateSystem cs;
100 1 : cs.addCoordinate(dc); cs.addCoordinate(stc); cs.addCoordinate(spc);
101 2 : TempImage<Complex> im(IPosition(4,100,100,1,1), cs);//, "gulu.image");
102 5 : std::vector<String> ftmachines={"GridFT", "WProjectFT", "SetJyGridFT"};
103 1 : refim::FTMachine * ftm=nullptr;
104 4 : for (uInt k=0; k < ftmachines.size(); ++k){
105 3 : cerr << "########### Test " << ftmachines[k] << endl;
106 3 : im.set(Complex(0.0));
107 3 : MPosition loc;
108 3 : MeasTable::Observatory(loc, MSColumns(thems).observation().telescopeName()(0));
109 3 : ftm=makeFTM(ftmachines[k], loc, thedir);
110 3 : Matrix<Float> weight;
111 3 : vi2.originChunks();
112 3 : vi2.origin();
113 3 : ftm->initializeToSky(im, weight, *vb);
114 3 : refim::SimpleComponentFTMachine cft;
115 :
116 6 : for (vi2.originChunks();vi2.moreChunks(); vi2.nextChunk()){
117 603 : for (vi2.origin(); vi2.more(); vi2.next()){
118 600 : cft.get(*vb, cl);
119 600 : vb->setVisCube(vb->visCubeModel());
120 600 : ftm->put(*vb);
121 : }
122 : }
123 3 : ftm->finalizeToSky();
124 :
125 3 : ftm->getImage(weight, true);
126 :
127 : // cerr << "val at center " << im.getAt(IPosition(4, 50, 50, 0, 0)) << endl;
128 3 : AlwaysAssertExit(near(6.66e-2, real( im.getAt(IPosition(4, 50, 50, 0, 0))), 1.0e-5));
129 : ///Let us degrid now
130 3 : im.set(Complex(0.0));
131 3 : im.putAt(Complex(999.0, 0.0),IPosition(4, 50, 50, 0, 0) );
132 3 : vi2.originChunks();
133 3 : vi2.origin();
134 3 : ftm->initializeToVis(im, *vb);
135 6 : for (vi2.originChunks();vi2.moreChunks(); vi2.nextChunk()){
136 603 : for (vi2.origin(); vi2.more(); vi2.next()){
137 600 : ftm->get(*vb);
138 : //cerr << "mod " << (vb->visCubeModel())(0,0,0) << endl;
139 600 : AlwaysAssertExit(near(999.0, abs((vb->visCubeModel())(0,0,0)), ftmachines[k]=="WProjectFT" ? 1e-3 : 1.0e-5));
140 : }
141 : }
142 3 : }
143 : //detach the ms for cleaning up
144 1 : thems=MeasurementSet();
145 1 : } catch (AipsError x) {
146 0 : cout << "Caught exception " << endl;
147 0 : cout << x.getMesg() << endl;
148 0 : return(1);
149 0 : }
150 1 : cerr <<"OK" << endl;
151 1 : exit(0);
152 : }
|
