Line data Source code
1 : //# VLAIlluminationConvFunc.h: Definition for VLAIlluminationConvFunc
2 : //# Copyright (C) 1996,1997,1998,1999,2000,2002
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 adressed 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 : //# $Id$
28 :
29 : #ifndef SYNTHESIS_BEAMCALC_H
30 : #define SYNTHESIS_BEAMCALC_H
31 :
32 : //#include <casa/complex.h>
33 : #include <casacore/images/Images/TempImage.h>
34 : #include <casacore/casa/Exceptions.h>
35 : #include <casacore/casa/Logging/LogIO.h>
36 :
37 : namespace casa
38 : {
39 :
40 : #define MAXGEOM 2000
41 :
42 : typedef struct /* all dimensions in meters, GHz */
43 : {
44 : char name[16]; /* name of antenna, e.g., VLA */
45 : char bandName[16];
46 : casacore::Double sub_h; /* subreflector vertex height above primary vertex */
47 : casacore::Double feedpos[3]; /* position of feed */
48 : casacore::Double subangle; /* angle subtended by the subreflector */
49 : casacore::Double legwidth; /* strut width */
50 : casacore::Double legfoot; /* distance from optic axis of leg foot */
51 : casacore::Double legapex; /* hight of leg intersection */
52 : casacore::Double Rhole; /* radius of central hole */
53 : casacore::Double Rant; /* antenna radius */
54 : casacore::Double reffreq; /* a reference frequency */
55 : casacore::Double taperpoly[5]; /* polynomial expanded about reffreq */
56 : casacore::Int ntaperpoly; /* number of terms in polynomial */
57 :
58 : casacore::Double astigm_0; /* astigmatism: coefficient of Zernike Polyn. Z6 a.k.a. 0-90 */
59 : casacore::Double astigm_45; /* astigmatism: coefficient of Zernike Polyn. Z5 a.k.a. 45-135 */
60 :
61 : /* to be added later
62 : casacore::Double focus;
63 : casacore::Double dfeedpos[3];
64 : casacore::Double dsub_z;
65 : */
66 : } BeamCalcGeometry;
67 :
68 : typedef struct
69 : {
70 : casacore::Int oversamp; /* average this many points per cell */
71 : casacore::ImageInterface<casacore::Complex> *aperture=nullptr; /* Jones planes [Nx,Ny,NStokes,NFreq]*/
72 : std::shared_ptr<casacore::ImageInterface<casacore::Complex> > apertureptr;
73 : casacore::Double x0, y0; /* center of cell 0, 0, meters */
74 : casacore::Double dx, dy; /* increment in meters */
75 : casacore::Int nx, ny; /* calculation plane size in cells */
76 : /* last cell is at coordinates:
77 : X = x0 + (nx-1)*dx
78 : Y = y0 + (ny-1)*dy
79 : */
80 : casacore::Double pa; /* Parallactic angle, radians */
81 : casacore::Double freq; /* GHz */
82 : casacore::Int band;
83 : } ApertureCalcParams;
84 :
85 : /*
86 : * calcAntenna parameters
87 : */
88 : typedef struct
89 : {
90 : casacore::Double sub_h; /* height of subreflector (on axis) */
91 : casacore::Double feed[3]; /* position of the feed */
92 : casacore::Double feeddir[3]; /* unit vector pointing along feed */
93 : casacore::Double pfeeddir[3]; /* same as above after pathology applied */
94 : casacore::Double radius; /* antenna radius (m) */
95 : casacore::Double K;
96 : casacore::Double deltar;
97 : casacore::Double zedge; /* height at the edge of the dish */
98 : casacore::Double bestparabola; /* best fit parabola quadratic coef */
99 : casacore::Double ftaper; /* taper of feed */
100 : casacore::Double thmax; /* maximum angle of feed */
101 : casacore::Double fa2pi;
102 : casacore::Double legwidth;
103 : casacore::Double legfoot, legfootz;
104 : casacore::Double legapex;
105 : casacore::Double legthick;
106 : casacore::Double hole_radius;
107 : casacore::Double freq, lambda;
108 : casacore::Double astigm_0; /* astigmatism: coefficient of Zernike Polyn. Z6 a.k.a. 0-90 */
109 : casacore::Double astigm_45; /* astigmatism: coefficient of Zernike Polyn. Z5 a.k.a. 45-135 */
110 : casacore::Double dir[3];
111 : casacore::Double hhat[3], vhat[3]; /* unit vectors orthogonal to dir */
112 : casacore::Double z[MAXGEOM];
113 : casacore::Double m[MAXGEOM];
114 : casacore::Double k[3];
115 : casacore::Int ngeom;
116 : char name[16];
117 : casacore::Int gridsize;
118 : } calcAntenna;
119 :
120 : typedef struct
121 : {
122 : casacore::Double subrot[3][3]; /* 3x3 matrix rotating x,y,z or nx,ny,nz */
123 : casacore::Double feedrot[3][3]; /* 3x3 matrix rotating x,y,z or nx,ny,nz */
124 : casacore::Double subshift[3]; /* 3 length vector */
125 : casacore::Double feedshift[3]; /* 3 length vector */
126 : casacore::Double subrotpoint[3]; /* 3 vector describing point to rotate sub. */
127 : casacore::Double az_offset; /* azimuth pointing offset (radians) */
128 : casacore::Double el_offset; /* elevation pointing offset (radians) */
129 : casacore::Double phase_offset; /* DC offset in phase (radians) */
130 : casacore::Double focus; /* meters out of focus toward subreflector */
131 : } Pathology;
132 :
133 : typedef struct
134 : {
135 : casacore::Double aper[6]; /* aperture x, y, z, nx, ny, nz */
136 : casacore::Double dish[6]; /* dish x, y, z, nx, ny, nz */
137 : casacore::Double sub[6]; /* subreflector x, y, z, nx, ny, nz */
138 : casacore::Double feed[3]; /* feed x, y, z */
139 : } Ray;
140 :
141 : enum VLABeamCalcBandCode{
142 : BeamCalc_VLA_L = 0,
143 : BeamCalc_VLA_C,
144 : BeamCalc_VLA_X,
145 : BeamCalc_VLA_U,
146 : BeamCalc_VLA_K,
147 : BeamCalc_VLA_Q,
148 : BeamCalc_VLA_4,
149 :
150 : VLABeamCalc_NumBandCodes /* this line last */
151 : };
152 :
153 : enum EVLABeamCalcBandCode{
154 : BeamCalc_EVLA_L = 0,
155 : BeamCalc_EVLA_S,
156 : BeamCalc_EVLA_C,
157 : BeamCalc_EVLA_X,
158 : BeamCalc_EVLA_U,
159 : BeamCalc_EVLA_K,
160 : BeamCalc_EVLA_A,
161 : BeamCalc_EVLA_Q,
162 : BeamCalc_EVLA_4,
163 :
164 : EVLABeamCalc_NumBandCodes /* this line last */
165 : };
166 :
167 :
168 : enum ALMABeamCalcBandCode{
169 : BeamCalc_ALMA_1 = 0,
170 : BeamCalc_ALMA_2,
171 : BeamCalc_ALMA_3,
172 : BeamCalc_ALMA_4,
173 : BeamCalc_ALMA_5,
174 : BeamCalc_ALMA_6,
175 : BeamCalc_ALMA_7,
176 : BeamCalc_ALMA_8,
177 : BeamCalc_ALMA_9,
178 : BeamCalc_ALMA_10,
179 :
180 : ALMABeamCalc_NumBandCodes /* this line last */
181 : };
182 :
183 : extern casacore::Double VLABandMinFreqDefaults[VLABeamCalc_NumBandCodes];
184 : extern casacore::Double VLABandMaxFreqDefaults[VLABeamCalc_NumBandCodes];
185 : extern BeamCalcGeometry VLABeamCalcGeometryDefaults[VLABeamCalc_NumBandCodes];
186 : extern casacore::Double EVLABandMinFreqDefaults[EVLABeamCalc_NumBandCodes];
187 : extern casacore::Double EVLABandMaxFreqDefaults[EVLABeamCalc_NumBandCodes];
188 : extern BeamCalcGeometry EVLABeamCalcGeometryDefaults[EVLABeamCalc_NumBandCodes];
189 : extern casacore::Double ALMABandMinFreqDefaults[ALMABeamCalc_NumBandCodes];
190 : extern casacore::Double ALMABandMaxFreqDefaults[ALMABeamCalc_NumBandCodes];
191 : extern BeamCalcGeometry ALMABeamCalcGeometryDefaults[ALMABeamCalc_NumBandCodes];
192 :
193 : class BeamCalc
194 : {
195 : public:
196 :
197 : // This is a SINGLETON class
198 : static BeamCalc* Instance();
199 :
200 : void setBeamCalcGeometries(const casacore::String& obsName, // (the observatory name, e.g. "ALMA" or "ACA")
201 : const casacore::String& antennaType = "STANDARD",
202 : const casacore::MEpoch& obsTime = casacore::MEpoch(casacore::Quantity(50000., "d")),
203 : const casacore::String& otherAntRayPath=""); // override the AntennaResponses casacore::Table in Observatories
204 :
205 : casacore::Int getBandID(const casacore::Double freq, // in Hz
206 : const casacore::String& obsname,
207 : const casacore::String& bandName,
208 : const casacore::String& antennaType="STANDARD",
209 190 : const casacore::MEpoch& obsTime = casacore::MEpoch(casacore::Quantity(50000., "d")),
210 : const casacore::String& otherAntRayPath=""); // override the AntennaResponses casacore::Table in Observatories
211 :
212 : casacore::Int calculateAperture(ApertureCalcParams *ap);
213 : casacore::Int calculateAperture(ApertureCalcParams *ap, const casacore::Int& whichStokes);
214 : casacore::Int calculateApertureLinPol(ApertureCalcParams *ap, const casacore::Int& whichStokes);
215 :
216 : protected:
217 : BeamCalc();
218 :
219 : private:
220 :
221 : void computePixelValues(const ApertureCalcParams *ap, const calcAntenna *a, const Pathology *p,
222 : const casacore::Double &L0, casacore::Complex *Er, casacore::Complex *El,
223 : const casacore::Int &i, const casacore::Int &j);
224 : void computePixelValues(const ApertureCalcParams *ap,
225 : const calcAntenna *a, const Pathology *p,
226 : const casacore::Double &L0,
227 : casacore::Complex *Er, casacore::Complex *El,
228 : const casacore::Int &i, const casacore::Int &j,
229 : const casacore::Int& whichStokes);
230 : void computePixelValuesLinPol(const ApertureCalcParams *ap,
231 : const calcAntenna *a, const Pathology *p,
232 : const casacore::Double &L0,
233 : casacore::Complex *Ex, casacore::Complex *Ey,
234 : const casacore::Int &i, const casacore::Int &j,
235 : const casacore::Int& whichStokes);
236 : //normalizes a "vector" of 3 Doubles in the vector sense
237 152300604 : inline void norm3(casacore::Double *v)
238 : {
239 : casacore::Double s;
240 152300604 : s = sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
241 152300604 : v[0] /= s;
242 152300604 : v[1] /= s;
243 152300604 : v[2] /= s;
244 152300604 : }
245 :
246 : calcAntenna *newAntenna(casacore::Double sub_h, casacore::Double feed_x, casacore::Double feed_y, casacore::Double feed_z,
247 : casacore::Double ftaper, casacore::Double thmax, const char *geomfile);
248 :
249 : void deleteAntenna(calcAntenna *a);
250 :
251 : void Antennasetfreq(calcAntenna *a, casacore::Double freq);
252 :
253 : void Antennasetdir(calcAntenna *a, const casacore::Double *dir);
254 :
255 : // sets feeddir after pathology is considered
256 : void alignfeed(calcAntenna *a, const Pathology *p);
257 :
258 : void getfeedbasis(const calcAntenna *a, casacore::Double B[3][3]);
259 :
260 : void Efield(const calcAntenna *a, const casacore::Complex *pol, casacore::Complex *E);
261 :
262 : casacore::Int Antennasetfeedpattern(calcAntenna *a, const char *filename, casacore::Double scale);
263 :
264 : calcAntenna *newAntennafromApertureCalcParams(ApertureCalcParams *ap);
265 :
266 : casacore::Int dishvalue(const calcAntenna *a, casacore::Double r, casacore::Double *z, casacore::Double *m);
267 :
268 : casacore::Int astigdishvalue(const calcAntenna *a, casacore::Double x, casacore::Double y, casacore::Double *z, casacore::Double *m);
269 :
270 : // Returns position of subreflector piece (x, y, z) and
271 : // its normal (u, v, w)
272 : casacore::Int subfromdish(const calcAntenna *a, casacore::Double x, casacore::Double y, casacore::Double *subpoint);
273 :
274 : casacore::Int dishfromsub(const calcAntenna *a, casacore::Double x, casacore::Double y, casacore::Double *dishpoint);
275 :
276 : void printAntenna(const calcAntenna *a);
277 :
278 : Ray* newRay(const casacore::Double *sub);
279 :
280 : void deleteRay(Ray *ray);
281 :
282 : Pathology* newPathology();
283 :
284 : Pathology* newPathologyfromApertureCalcParams(ApertureCalcParams *ap);
285 :
286 : void deletePathology(Pathology *P);
287 :
288 : void normvec(const casacore::Double *a, const casacore::Double *b, casacore::Double *c);
289 :
290 : casacore::Double dAdOmega(const calcAntenna *a, const Ray *ray1, const Ray *ray2,
291 : const Ray *ray3, const Pathology *p);
292 :
293 : casacore::Double dOmega(const calcAntenna *a, const Ray *ray1, const Ray *ray2,
294 : const Ray *ray3, const Pathology *p);
295 :
296 : casacore::Double Raylen(const Ray *ray);
297 :
298 : void Pathologize(casacore::Double *sub, const Pathology *p);
299 :
300 : void applyPathology(Pathology *P, calcAntenna *a);
301 :
302 : void intersectdish(const calcAntenna *a, const casacore::Double *sub, const casacore::Double *unitdir,
303 : casacore::Double *dish, casacore::Int niter);
304 :
305 : void intersectaperture(const calcAntenna *a, const casacore::Double *dish,
306 : const casacore::Double *unitdir, casacore::Double *aper);
307 :
308 : casacore::Double feedfunc(const calcAntenna *a, casacore::Double theta);
309 :
310 : casacore::Double feedgain(const calcAntenna *a, const Ray *ray, const Pathology *p);
311 :
312 : Ray* trace(const calcAntenna *a, casacore::Double x, casacore::Double y, const Pathology *p);
313 :
314 : void tracepol(casacore::Complex *E0, const Ray *ray, casacore::Complex *E1);
315 :
316 : casacore::Int legplanewaveblock(const calcAntenna *a, casacore::Double x, casacore::Double y);
317 :
318 : casacore::Int legplanewaveblock2(const calcAntenna *a, const Ray *ray);
319 :
320 : casacore::Int legsphericalwaveblock(const calcAntenna *a, const Ray *ray);
321 :
322 : void copyBeamCalcGeometry(BeamCalcGeometry* to, BeamCalcGeometry* from);
323 :
324 : static BeamCalc* instance_p;
325 :
326 : casacore::String obsName_p;
327 : casacore::String antType_p;
328 : casacore::MEpoch obsTime_p;
329 : casacore::uInt BeamCalc_NumBandCodes_p;
330 : casacore::Vector<BeamCalcGeometry> BeamCalcGeometries_p;
331 : casacore::Vector<casacore::Double> bandMinFreq_p; // in Hz
332 : casacore::Vector<casacore::Double> bandMaxFreq_p; // in Hz
333 : casacore::String antRespPath_p;
334 :
335 : static const casacore::Double METER_INCH;
336 : static const casacore::Double INCH_METER;
337 : static const casacore::Double NS_METER;
338 : static const casacore::Double METER_NS;
339 : static const casacore::Double DEG_RAD;
340 : static const casacore::Double RAD_DEG;
341 :
342 : };
343 :
344 : };
345 :
346 : #endif
|