Line data Source code
1 : //# BeamCalc.cc: Implementation for BeamCalc
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 : //#include <stdio.h>
30 : //#include <complex.h>
31 : #include <cmath>
32 : #include <math.h>
33 : //#include <stdlib.h>
34 : //#include <string.h>
35 : #include <casacore/images/Images/TempImage.h>
36 : #include <synthesis/MeasurementEquations/AntennaResponses.h>
37 : #include <casacore/tables/Tables/TableProxy.h>
38 : #include <casacore/casa/Exceptions.h>
39 : #include <casacore/casa/Containers/ValueHolder.h>
40 : #include <casacore/casa/Arrays/Array.h>
41 : #include <synthesis/TransformMachines/SynthesisError.h>
42 : #include <synthesis/TransformMachines/BeamCalc.h>
43 : #include <casacore/casa/OS/Timer.h>
44 : #include <casacore/casa/System/AppState.h>
45 : #include <casacore/casa/OS/Directory.h>
46 : #include <casatools/Config/State.h>
47 : #ifdef _OPENMP
48 : #include <omp.h>
49 : #endif
50 : #if ((__GNUC__ >= 4) && (__GNUC_MINOR__ >= 4))
51 : #define GCC44x 1
52 : #else
53 : #define GCC44x 0
54 : #endif
55 :
56 :
57 : using namespace std;
58 : using namespace casacore;
59 : namespace casa{
60 :
61 : const Double BeamCalc::METER_INCH = 39.37008;
62 : const Double BeamCalc::INCH_METER = (1.0/BeamCalc::METER_INCH);
63 : const Double BeamCalc::NS_METER = 0.299792458; // Exact
64 : const Double BeamCalc::METER_NS = (1.0/BeamCalc::NS_METER);
65 : const Double BeamCalc::DEG_RAD = M_PI/180.0;
66 : const Double BeamCalc::RAD_DEG = 180.0/M_PI;
67 :
68 : BeamCalc* BeamCalc::instance_p = 0;
69 :
70 1 : BeamCalc::BeamCalc(const String repopath):
71 1 : obsName_p(""),
72 1 : antType_p(""),
73 1 : obsTime_p(),
74 1 : BeamCalc_NumBandCodes_p(0),
75 1 : BeamCalcGeometries_p(0),
76 1 : bandMinFreq_p(0),
77 1 : bandMaxFreq_p(0),
78 1 : antRespPath_p(""), repoPath_p(repopath){
79 1 : }
80 :
81 3424 : BeamCalc* BeamCalc::Instance(const String repopath){
82 3424 : if(instance_p==0){
83 1 : instance_p = new BeamCalc(repopath);
84 : }
85 3424 : return instance_p;
86 : }
87 :
88 : // initialise the beam calculation parameters
89 1733 : void BeamCalc::setBeamCalcGeometries(const String& obsName,
90 : const String& antType,
91 : const MEpoch& obsTime,
92 : const String& otherAntRayPath){
93 :
94 1733 : Unit uS("s");
95 1733 : Bool verbose = false;
96 :
97 :
98 5199 : if(obsName==obsName_p
99 1729 : && antType==antType_p
100 3462 : && obsTime.get(uS).getValue()==obsTime_p.get(uS).getValue()
101 3462 : && otherAntRayPath.empty()
102 : ){
103 1729 : return; // nothing to do (assuming the databases haven't changed)
104 : }
105 :
106 4 : cout << "Processing request for geometries from observatory " << obsName << ", antenna type " << antType << endl;
107 :
108 4 : LogIO os;
109 4 : os << LogOrigin("BeamCalc", "setBeamCalcGeometries()");
110 :
111 4 : if(obsName!=""){
112 4 : obsName_p = obsName;
113 : }
114 4 : if(antType!=""){
115 4 : antType_p = antType;
116 : }
117 4 : obsTime_p = obsTime;
118 :
119 :
120 4 : BeamCalcGeometries_p.resize(0);
121 :
122 4 : AntennaResponses aR;
123 4 : String antRespPath;
124 4 : String antRayPath = otherAntRayPath;
125 :
126 4 : Bool useInternal = false;
127 :
128 4 : os << LogIO::NORMAL << "Initialisation of geometries for observatory " << obsName_p
129 4 : << ", antenna type " << antType_p << LogIO::POST;
130 :
131 4 : if(otherAntRayPath.empty()){
132 4 : if(!MeasTable::AntennaResponsesPath(antRespPath, obsName_p)) {
133 4 : useInternal = true;
134 : }
135 : else{
136 0 : if(!aR.init(antRespPath)){
137 : // init failed
138 : String mesg="Initialisation of antenna response parameters for observatory "
139 0 : +obsName_p+" failed using path "+antRespPath;
140 0 : SynthesisError err(mesg);
141 0 : throw(err);
142 0 : }
143 : uInt respImageChannel;
144 0 : MFrequency respImageNomFreq;
145 : AntennaResponses::FuncTypes respImageFType;
146 0 : MVAngle respImageRotOffset;
147 :
148 0 : if(!aR.getImageName(antRayPath,
149 : respImageChannel,
150 : respImageNomFreq,
151 : respImageFType,
152 : respImageRotOffset,
153 : //
154 0 : obsName_p,
155 0 : obsTime_p,
156 0 : MFrequency(Quantity(0.,Unit("Hz")), MFrequency::TOPO), // any frequency
157 0 : AntennaResponses::INTERNAL,
158 0 : antType_p
159 : )
160 : ){ // no matching response found
161 : os << LogIO::NORMAL << "No matching antenna response found for observatory "
162 0 : << obsName_p << LogIO::POST;
163 0 : useInternal = true;
164 : }
165 0 : }
166 :
167 4 : if(useInternal){
168 :
169 4 : Bool found = False;
170 4 : String fullFileName;
171 4 : const std::list<std::string> &data_path = AppStateSource::fetch( ).dataPath( );
172 4 : const std::string distrodata_path = repoPath_p.size()> 0 ? repoPath_p : casatools::get_state().distroDataPath( );
173 : // The data path search need to be rewritten to adopt the recommanded
174 : // setting via python file for CASA6. For now, only the first path that
175 : // actually exist will be set to the data path (TT 2018.12.10)
176 4 : if (data_path.size() > 0 || distrodata_path.size() > 0) {
177 4 : for (std::list<std::string>::const_iterator it = data_path.begin();
178 4 : !found && it != data_path.end(); ++it) {
179 4 : Path lpath = Path(*it);
180 : // os<<"Here to datapath="<<lpath<<LogIO::POST;
181 : // Path lpath = Path(data_path);
182 4 : String slpath = lpath.absoluteName();
183 4 : String subdirname;
184 4 : if (obsName_p == "VLA" || obsName_p == "EVLA") {
185 4 : subdirname = "/nrao/VLA";
186 0 : } else if (obsName_p == "ALMA") {
187 0 : subdirname = "/alma/response";
188 : }
189 : // Directory ddir(slpath+subdirname);
190 : try {
191 :
192 4 : Directory ddir(slpath+subdirname);
193 4 : ddir.exists();
194 4 : found = True;
195 4 : fullFileName=slpath;
196 4 : break;
197 4 : }
198 0 : catch (...) {
199 0 : }
200 :
201 : //if (ddir.exists()) {
202 : // cerr<<" ddir exists:"<<slpath<<subdirname<<endl;
203 : // found = True;
204 : // fullFileName=slpath;
205 : // break;
206 : //}
207 12 : }
208 4 : if (!found && distrodata_path!="") {
209 0 : fullFileName = distrodata_path;
210 0 : found = True;
211 : }
212 : }
213 0 : else if(!found) {
214 0 : const char *sep=" ";
215 0 : char *aipsPath = strtok(getenv("CASAPATH"),sep);
216 0 : if (aipsPath == NULL)
217 0 : throw(SynthesisError("CASAPATH not found."));
218 0 : fullFileName=aipsPath;
219 0 : fullFileName+="/data";
220 : }
221 :
222 :
223 4 : if(obsName_p=="VLA" && antType_p=="STANDARD"){
224 2 : os << LogIO::NORMAL << "Will use default geometries for VLA STANDARD." << LogIO::POST;
225 2 : BeamCalc_NumBandCodes_p = VLABeamCalc_NumBandCodes;
226 2 : BeamCalcGeometries_p.resize(BeamCalc_NumBandCodes_p);
227 2 : bandMinFreq_p.resize(BeamCalc_NumBandCodes_p);
228 2 : bandMaxFreq_p.resize(BeamCalc_NumBandCodes_p);
229 16 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
230 14 : copyBeamCalcGeometry(&BeamCalcGeometries_p[i], &VLABeamCalcGeometryDefaults[i]);
231 14 : bandMinFreq_p[i] = VLABandMinFreqDefaults[i];
232 14 : bandMaxFreq_p[i] = VLABandMaxFreqDefaults[i];
233 : }
234 : //antRespPath_p = fullFileName + "/data/nrao/VLA";
235 2 : antRespPath_p = fullFileName + "/nrao/VLA";
236 : }
237 2 : else if(obsName_p=="EVLA" && antType_p=="STANDARD"){
238 2 : os << LogIO::NORMAL << "Will use default geometries for EVLA STANDARD." << LogIO::POST;
239 2 : BeamCalc_NumBandCodes_p = EVLABeamCalc_NumBandCodes;
240 2 : BeamCalcGeometries_p.resize(BeamCalc_NumBandCodes_p);
241 2 : bandMinFreq_p.resize(BeamCalc_NumBandCodes_p);
242 2 : bandMaxFreq_p.resize(BeamCalc_NumBandCodes_p);
243 20 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
244 18 : copyBeamCalcGeometry(&BeamCalcGeometries_p[i], &EVLABeamCalcGeometryDefaults[i]);
245 18 : bandMinFreq_p[i] = EVLABandMinFreqDefaults[i];
246 18 : bandMaxFreq_p[i] = EVLABandMaxFreqDefaults[i];
247 : }
248 : //antRespPath_p = fullFileName + "/data/nrao/VLA";
249 2 : antRespPath_p = fullFileName + "/nrao/VLA";
250 : }
251 0 : else if(obsName_p=="ALMA" && (antType_p=="DA" || antType_p=="DV" || antType_p=="PM")){
252 0 : os << LogIO::NORMAL << "Will use default geometries for ALMA DA, DV, and PM." << LogIO::POST;
253 0 : BeamCalc_NumBandCodes_p = ALMABeamCalc_NumBandCodes;
254 0 : BeamCalcGeometries_p.resize(BeamCalc_NumBandCodes_p);
255 0 : bandMinFreq_p.resize(BeamCalc_NumBandCodes_p);
256 0 : bandMaxFreq_p.resize(BeamCalc_NumBandCodes_p);
257 0 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
258 0 : copyBeamCalcGeometry(&BeamCalcGeometries_p[i], &ALMABeamCalcGeometryDefaults[i]);
259 0 : if(antType_p=="DA"){
260 0 : BeamCalcGeometries_p[i].legwidth *= -1.; // change from + to x shape
261 : }
262 0 : bandMinFreq_p[i] = ALMABandMinFreqDefaults[i];
263 0 : bandMaxFreq_p[i] = ALMABandMaxFreqDefaults[i];
264 : }
265 : //antRespPath_p = fullFileName + "/data/alma/responses";
266 0 : antRespPath_p = fullFileName + "/alma/responses";
267 : }
268 : else{
269 : String mesg="We don't have any antenna ray tracing parameters available for observatory "
270 0 : +obsName_p+", antenna type "+antType_p;
271 0 : SynthesisError err(mesg);
272 0 : throw(err);
273 0 : }
274 4 : return;
275 4 : } // end if(useInternal)
276 : }
277 :
278 :
279 0 : os << LogIO::NORMAL << "from file " << antRayPath << endl;
280 :
281 : try {
282 : // read temp table from ASCII file
283 0 : TableProxy antParTab = TableProxy(antRayPath, String(""), String("tempRayTraceTab.tab"),
284 0 : false, IPosition(), // autoheader, autoshape
285 0 : String(" "), // separator
286 0 : String("#"), // comment marker
287 : 0,-1, // first and last line
288 0 : Vector<String>(), Vector<String>());
289 :
290 0 : antParTab.deleteTable(true); // table will be deleted when it goes out of scope
291 :
292 : // read the table
293 0 : uInt nRows = antParTab.nrows();
294 0 : BeamCalc_NumBandCodes_p = nRows;
295 :
296 0 : BeamCalcGeometries_p.resize(BeamCalc_NumBandCodes_p);
297 0 : bandMinFreq_p.resize(BeamCalc_NumBandCodes_p);
298 0 : bandMaxFreq_p.resize(BeamCalc_NumBandCodes_p);
299 :
300 0 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
301 0 : sprintf(BeamCalcGeometries_p[i].name, "%s", antParTab.getCell("NAME", i).asString().c_str());
302 0 : bandMinFreq_p[i] = antParTab.getCell("MINFREQ", i).asDouble() * 1E9; // expect GHz
303 0 : bandMaxFreq_p[i] = antParTab.getCell("MAXFREQ", i).asDouble() * 1E9;
304 0 : BeamCalcGeometries_p[i].sub_h = antParTab.getCell("SUB_H", i).asDouble();
305 0 : Array<Double> ta1;
306 0 : ta1.assign(antParTab.getCell("FEEDPOS", i).asArrayDouble());
307 0 : for(uInt j=0; j<3;j++){
308 0 : BeamCalcGeometries_p[i].feedpos[j] = ta1(IPosition(1,j));
309 : }
310 0 : BeamCalcGeometries_p[i].subangle = antParTab.getCell("SUBANGLE", i).asDouble();
311 0 : BeamCalcGeometries_p[i].legwidth = antParTab.getCell("LEGWIDTH", i).asDouble();
312 0 : BeamCalcGeometries_p[i].legfoot = antParTab.getCell("LEGFOOT", i).asDouble();
313 0 : BeamCalcGeometries_p[i].legapex = antParTab.getCell("LEGAPEX", i).asDouble();
314 0 : BeamCalcGeometries_p[i].Rhole = antParTab.getCell("RHOLE", i).asDouble();
315 0 : BeamCalcGeometries_p[i].Rant = antParTab.getCell("RANT", i).asDouble();
316 0 : BeamCalcGeometries_p[i].reffreq = antParTab.getCell("REFFREQ", i).asDouble(); // stay in GHz
317 0 : Array<Double> ta2;
318 0 : ta2.assign(antParTab.getCell("TAPERPOLY", i).asArrayDouble());
319 0 : for(uInt j=0; j<5;j++){
320 0 : BeamCalcGeometries_p[i].taperpoly[j] = ta2(IPosition(1,j));
321 : }
322 0 : BeamCalcGeometries_p[i].ntaperpoly = antParTab.getCell("NTAPERPOLY", i).asInt();
323 0 : BeamCalcGeometries_p[i].astigm_0 = antParTab.getCell("ASTIGM_0", i).asDouble();
324 0 : BeamCalcGeometries_p[i].astigm_45 = antParTab.getCell("ASTIGM_45", i).asDouble();
325 0 : if(verbose){
326 : cout << "i name bandMinFreq_p bandMaxFreq_p sub_h feedpos feedpos feedpos subangle legwidth legfoot legapex"
327 0 : << " Rhole Rant reffreq taperpoly taperpoly taperpoly taperpoly taperpoly ntaperpoly astigm0 astigm45" << endl;
328 0 : cout << i << " " << BeamCalcGeometries_p[i].name << " " << bandMinFreq_p[i] << " " << bandMaxFreq_p[i]
329 0 : << " " << BeamCalcGeometries_p[i].sub_h
330 0 : << " " << BeamCalcGeometries_p[i].feedpos[0] << " " << BeamCalcGeometries_p[i].feedpos[1]
331 0 : << " " << BeamCalcGeometries_p[i].feedpos[2]
332 0 : << " " << BeamCalcGeometries_p[i].subangle << " " << BeamCalcGeometries_p[i].legwidth
333 0 : << " " << BeamCalcGeometries_p[i].legfoot << " " << BeamCalcGeometries_p[i].legapex
334 0 : << " " << BeamCalcGeometries_p[i].Rhole << " " << BeamCalcGeometries_p[i].Rant << " " << BeamCalcGeometries_p[i].reffreq
335 0 : << " " << BeamCalcGeometries_p[i].taperpoly[0] << " " << BeamCalcGeometries_p[i].taperpoly[1]
336 0 : << " " << BeamCalcGeometries_p[i].taperpoly[2] << " " << BeamCalcGeometries_p[i].taperpoly[3]
337 0 : << " " << BeamCalcGeometries_p[i].taperpoly[4] << " " << BeamCalcGeometries_p[i].ntaperpoly
338 0 : << " " << BeamCalcGeometries_p[i].astigm_0 << " " << BeamCalcGeometries_p[i].astigm_45 << endl;
339 : }
340 0 : }
341 :
342 0 : } catch (AipsError x) {
343 0 : String mesg="Initialisation of antenna ray tracing parameters for observatory "+obsName_p
344 0 : +" failed using path "+antRayPath+"\n with message "+x.getMesg();
345 0 : BeamCalcGeometries_p.resize(0);
346 0 : SynthesisError err(mesg);
347 0 : throw(err);
348 0 : }
349 :
350 0 : if(antRespPath.empty()){ // use containing directory of the antRayPath
351 0 : antRespPath_p = Path(antRayPath).dirName();
352 : }
353 : else{
354 0 : antRespPath_p = Path(antRespPath).dirName();
355 : }
356 :
357 0 : os << LogIO::NORMAL << "... successful." << LogIO::POST;
358 :
359 0 : return;
360 :
361 1733 : }
362 :
363 1733 : Int BeamCalc::getBandID(Double freq, // in Hz
364 : const String& obsName,
365 : const String& bandName,
366 : const String& antType,
367 : const MEpoch& obsTime,
368 : const String& otherAntRayPath){
369 :
370 1733 : setBeamCalcGeometries(obsName, antType, obsTime, otherAntRayPath);
371 :
372 : // Check against bandName if it is a non-NULL string. Otherwise
373 : // check against frequency range. The latter method is only for
374 : // backward compatibility reasons and for older MSes which did not
375 : // have correct band names in the SPW sub-table.
376 1733 : if (bandName != "")
377 2167 : for(uInt i=0; i<BeamCalcGeometries_p.nelements(); i++)
378 2125 : if(String(BeamCalcGeometries_p[i].bandName)==bandName) return i;
379 :
380 : // If the control flow gets here, the given bandName did not match
381 : // in SPW names in the MS. So use the fall-back option of using
382 : // the reference frequency to get the band ID (this will lead to
383 : // incorrect ID for the edge SPWs which might overlap in frequecy
384 : // with the adjecent band).
385 62 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
386 62 : if((bandMinFreq_p[i]<=freq)&&(freq<=bandMaxFreq_p[i])){
387 42 : return i;
388 : }
389 : }
390 0 : ostringstream mesg;
391 0 : mesg << obsName << "/" << bandName << "/" << antType << "/" << freq << "(Hz) combination not recognized.";
392 0 : throw(SynthesisError(mesg.str()));
393 :
394 0 : }
395 :
396 :
397 :
398 1691 : calcAntenna* BeamCalc::newAntenna(Double sub_h, Double feed_x, Double feed_y, Double feed_z,
399 : Double ftaper, Double thmax, const char *geomfile)
400 : {
401 : calcAntenna *a;
402 : Int i;
403 : Double d, r, m, z;
404 : FILE *in;
405 1691 : String fullFileName(antRespPath_p);
406 1691 : fullFileName = fullFileName + String("/") + geomfile;
407 :
408 1691 : in = fopen(fullFileName.c_str(), "r");
409 :
410 1691 : if(!in)
411 : {
412 0 : String msg = "File " + fullFileName
413 0 : + " not found.\n Did you forget to install package data repository?\n";
414 0 : throw(SynthesisError(msg));
415 0 : }
416 :
417 1691 : a = (calcAntenna *)malloc(sizeof(calcAntenna));
418 :
419 2117132 : for(i = 0; i < MAXGEOM; i++)
420 : {
421 2117132 : if(fscanf(in, "%lf%lf%lf", &r, &z, &m) != 3) break;
422 2115441 : a->z[i] = z;
423 2115441 : a->m[i] = m;
424 2115441 : a->radius = r;
425 : }
426 1691 : fclose(in);
427 1691 : a->ngeom = i;
428 1691 : a->zedge = z;
429 1691 : a->deltar = a->radius/(float)(a->ngeom-1.0);
430 1691 : a->bestparabola = a->zedge/(a->radius*a->radius);
431 1691 : if(i < 3)
432 : {
433 0 : fprintf(stderr, "geom file not valid\n");
434 0 : free(a);
435 0 : return 0;
436 : }
437 :
438 1691 : z = sub_h-feed_z;
439 :
440 1691 : a->sub_h = sub_h;
441 1691 : a->feed[0] = feed_x;
442 1691 : a->feed[1] = feed_y;
443 1691 : a->feed[2] = feed_z;
444 1691 : d = std::sqrt((double)(feed_x*feed_x + feed_y*feed_y + z*z));
445 1691 : if(z > 0.0)
446 : {
447 1691 : a->K = sub_h + d;
448 1691 : a->feeddir[0] = -feed_x/d;
449 1691 : a->feeddir[1] = -feed_y/d;
450 1691 : a->feeddir[2] = (sub_h-feed_z)/d;
451 : }
452 : else
453 : {
454 0 : a->K = std::sqrt((double(feed_x*feed_x + feed_y*feed_y + feed_z*feed_z)));
455 0 : a->feeddir[0] = -feed_x/d;
456 0 : a->feeddir[1] = -feed_y/d;
457 0 : a->feeddir[2] = (sub_h-feed_z)/d;
458 : }
459 6764 : for(i = 0; i < 3; i++) a->pfeeddir[i] = a->feeddir[i];
460 1691 : a->ftaper = fabs(ftaper);
461 1691 : a->thmax = thmax;
462 1691 : a->fa2pi = 2.0*M_PI*std::sqrt((double)ftaper)*0.1874/sin(thmax*M_PI/180.0);
463 1691 : a->legwidth = 0.0;
464 1691 : a->legfoot = a->radius/2.0;
465 1691 : a->legapex = sub_h*1.2;
466 1691 : a->legthick = 0.0;
467 1691 : a->hole_radius = 0.0;
468 1691 : a->dir[0] = a->dir[1] = 0.0;
469 1691 : a->dir[2] = 1.0;
470 1691 : strcpy(a->name, "unnamed");
471 1691 : a->k[0] = a->k[1] = a->k[2] = 0.0;
472 : /* default to no polarization state */
473 1691 : Antennasetfreq(a, 1.0);
474 1691 : Antennasetdir(a, 0); /* compute hhat and vhat */
475 1691 : a->gridsize = 0;
476 1691 : dishvalue(a, a->legfoot, &a->legfootz, 0);
477 :
478 1691 : return a;
479 1691 : }
480 :
481 1691 : void BeamCalc::deleteAntenna(calcAntenna *a)
482 : {
483 1691 : if(!a) return;
484 :
485 1691 : free(a);
486 : }
487 :
488 3382 : void BeamCalc::Antennasetfreq(calcAntenna *a, Double freq)
489 : {
490 : Int i;
491 :
492 3382 : a->freq = freq;
493 3382 : a->lambda = NS_METER/freq;
494 13528 : for(i = 0; i < 3; i++) a->k[i] = -2.0*M_PI*a->dir[i]/a->lambda;
495 3382 : }
496 :
497 3382 : void BeamCalc::Antennasetdir(calcAntenna *a, const Double *dir)
498 : {
499 : Double hmag;
500 : Int i;
501 :
502 3382 : if(dir)
503 : {
504 6764 : for(i = 0; i < 3; i++) a->dir[i] = dir[i];
505 1691 : if(a->dir[0] == 0.0 && a->dir[1] == 0.0)
506 : {
507 1691 : a->hhat[0] = 1.0;
508 1691 : a->hhat[1] = a->hhat[2] = 0.0;
509 1691 : a->vhat[1] = 1.0;
510 1691 : a->vhat[0] = a->vhat[2] = 0.0;
511 : }
512 : else
513 : {
514 0 : a->hhat[0] = a->dir[1];
515 0 : a->hhat[1] = -a->dir[0];
516 0 : a->hhat[2] = 0.0;
517 0 : hmag = sqrt(a->hhat[0]*a->hhat[0]
518 0 : + a->hhat[1]*a->hhat[1]);
519 0 : a->hhat[0] /= hmag;
520 0 : a->hhat[1] /= hmag;
521 :
522 0 : a->vhat[0] = a->hhat[1]*a->dir[2]
523 0 : - a->hhat[2]*a->dir[1];
524 0 : a->vhat[1] = a->hhat[2]*a->dir[0]
525 0 : - a->hhat[0]*a->dir[2];
526 0 : a->vhat[2] = a->hhat[0]*a->dir[1]
527 0 : - a->hhat[1]*a->dir[0];
528 : }
529 : }
530 13528 : for(i = 0; i < 3; i++) a->k[i] = -2.0*M_PI*a->dir[i]/a->lambda;
531 3382 : }
532 :
533 : /* sets feeddir after pathology is considered */
534 0 : void BeamCalc::alignfeed(calcAntenna *a, const Pathology *p)
535 : {
536 : Int i, j;
537 0 : Double f[3], s0[3], s[3], d[3], m=0.0;
538 :
539 0 : for(i = 0; i < 3; i++) f[i] = a->feed[i] + p->feedshift[i];
540 0 : for(i = 0; i < 3; i++) s0[i] = -p->subrotpoint[i];
541 0 : s0[2] += a->sub_h;
542 0 : for(i = 0; i < 3; i++)
543 : {
544 0 : s[i] = 0.0;
545 0 : for(j = 0; j < 3; j++)
546 0 : s[i] += p->subrot[i][j]*s0[j];
547 0 : s[i] += p->subrotpoint[i] + p->subshift[i];
548 0 : d[i] = s[i]-f[i];
549 0 : m += d[i]*d[i];
550 : }
551 0 : m = sqrt(m);
552 0 : for(i = 0; i < 3; i++) a->feeddir[i] = d[i]/m;
553 0 : }
554 :
555 5073 : void BeamCalc::getfeedbasis(const calcAntenna *a, Double B[3][3])
556 : {
557 : Int i;
558 : Double *dir, *vhat, *hhat;
559 :
560 5073 : hhat = B[0];
561 5073 : vhat = B[1];
562 5073 : dir = B[2];
563 :
564 20292 : for(i = 0; i < 3; i++) dir[i] = a->pfeeddir[i];
565 :
566 5073 : if(dir[0] == 0.0 && dir[1] == 0.0)
567 : {
568 0 : vhat[0] = 1.0;
569 0 : vhat[1] = vhat[2] = 0.0;
570 0 : hhat[1] = 1.0;
571 0 : hhat[0] = hhat[2] = 0.0;
572 : }
573 : else
574 : {
575 5073 : vhat[0] = dir[1];
576 5073 : vhat[1] = -dir[0];
577 5073 : vhat[2] = 0.0;
578 5073 : norm3(vhat);
579 :
580 5073 : hhat[0] = vhat[1]*dir[2] - vhat[2]*dir[1];
581 5073 : hhat[1] = vhat[2]*dir[0] - vhat[0]*dir[2];
582 5073 : hhat[2] = vhat[0]*dir[1] - vhat[1]*dir[0];
583 : }
584 5073 : }
585 :
586 3382 : void BeamCalc::Efield(const calcAntenna *a, const Complex *pol, Complex *E)
587 : {
588 : Double B[3][3];
589 : Double *hhat, *vhat;
590 :
591 3382 : getfeedbasis(a, B);
592 3382 : hhat = B[0];
593 3382 : vhat = B[1];
594 :
595 13528 : for(Int i = 0; i < 3; i++)
596 10146 : E[i] = Complex(hhat[i],0) * pol[0] + Complex(vhat[i],0) * pol[1];
597 3382 : }
598 :
599 0 : Int BeamCalc::Antennasetfeedpattern(calcAntenna* /*a*/,
600 : const char* /*filename*/,
601 : Double /*scale*/)
602 : {
603 : #if 0
604 : Int i, N, Nmax;
605 : Double x, delta;
606 : VecArray pat;
607 :
608 : a->feedpatterndelta = 0.0;
609 : if(a->feedpattern) deleteVector(a->feedpattern);
610 :
611 : if(filename == 0) return 1;
612 :
613 : pat = VecArrayfromfile(filename, 2);
614 :
615 : if(!pat) return 0;
616 : N = VectorSize(pat[0]);
617 : g_assert(N > 2);
618 : g_assert(pat[0][0] == 0.0);
619 :
620 : delta = pat[0][1];
621 : g_assert(delta > 0.0);
622 : for(i = 2; i < N; i++)
623 : {
624 : x = pat[0][i]-pat[0][i-1]-delta;
625 : g_assert(fabs(x) < delta/10000.0);
626 : }
627 :
628 : /* convert to radians */
629 : delta *= M_PI/180.0;
630 :
631 : /* and scale it */
632 : if(scale > 0.0) delta *= scale;
633 :
634 : /* Do we need to truncate the pattern? */
635 : Nmax = M_PI/delta;
636 : if(N > Nmax)
637 : {
638 : a->feedpattern = newVector(Nmax);
639 : for(i = 0; i < Nmax; i++)
640 : a->feedpattern[i] = fabs(pat[1][i]);
641 : deleteVector(pat[1]);
642 : }
643 : else a->feedpattern = pat[1];
644 :
645 : a->feedpatterndelta = delta;
646 : deleteVector(pat[0]);
647 : deleteVecArray(pat);
648 : #endif
649 0 : return 1;
650 : }
651 :
652 1691 : calcAntenna* BeamCalc::newAntennafromApertureCalcParams(ApertureCalcParams *ap)
653 : {
654 : calcAntenna *a;
655 1691 : Double dir[3] = {0.0, 0.0, 1.0};
656 : Double sub_h, feed_x, feed_y, feed_z, thmax, ftaper;
657 : char geomfile[128];//, *feedfile;
658 : BeamCalcGeometry *geom;
659 : Int i;
660 : Double x, freq, df;
661 :
662 : //LogIO os;
663 1691 : if((0<=ap->band) && (ap->band<(Int)BeamCalcGeometries_p.size())){
664 1691 : geom = &(BeamCalcGeometries_p[ap->band]);
665 : //os << "Using antenna parameters for " << geom->bandName << " band" << LogIO::POST;
666 : }
667 : else{
668 0 : SynthesisError err(String("Internal Error: attempt to access beam geometry for non-existing band."));
669 0 : throw(err);
670 0 : }
671 :
672 1691 : sub_h = geom->sub_h;
673 1691 : feed_x = geom->feedpos[0]; feed_x = -feed_x;
674 1691 : feed_y = geom->feedpos[1];
675 1691 : feed_z = geom->feedpos[2];
676 : //feedfile = 0;
677 1691 : thmax = geom->subangle;
678 :
679 1691 : freq = ap->freq;
680 1691 : if(freq <= 0.0) freq = geom->reffreq;
681 :
682 :
683 1691 : df = freq-geom->reffreq;
684 1691 : x = 1.0;
685 1691 : ftaper = 0.0;
686 3614 : for(i = 0; i < geom->ntaperpoly; i++)
687 : {
688 1923 : ftaper += geom->taperpoly[i]*x;
689 1923 : x *= df;
690 : }
691 1691 : sprintf(geomfile, "%s.surface", geom->name);
692 :
693 1691 : a = newAntenna(sub_h, feed_x, feed_y, feed_z, ftaper, thmax, geomfile);
694 1691 : if(!a) return 0;
695 :
696 1691 : strcpy(a->name, geom->name);
697 :
698 : /* feed pattern file is two column text file containing
699 : * angle (in degrees) and power (in dBi)
700 : */
701 :
702 : // if(feedfile != 0)
703 : // {
704 : // Double scale;
705 : // scale = getKeyValueDouble(kv, "feedpatternscale");
706 : // if(!Antennasetfeedpattern(a, feedfile, scale))
707 : // {
708 : // deleteAntenna(a);
709 : // fprintf(stderr, "Problem with feed file <%s>\n",
710 : // feedfile);
711 : // return 0;
712 : // }
713 : // }
714 :
715 1691 : Antennasetfreq(a, ap->freq);
716 :
717 1691 : a->legwidth = geom->legwidth;
718 1691 : a->legfoot = geom->legfoot;
719 1691 : a->legapex = geom->legapex;
720 :
721 1691 : a->hole_radius = geom->Rhole;
722 :
723 1691 : a->astigm_0 = geom->astigm_0;
724 1691 : a->astigm_45 = geom->astigm_45;
725 :
726 1691 : Antennasetdir(a, dir);
727 :
728 1691 : return a;
729 : }
730 :
731 269959130 : Int BeamCalc::dishvalue(const calcAntenna *a, Double r, Double *z, Double *m)
732 : {
733 : Double ma, mb, mc, zav, A, B, C, D;
734 : Double x, d, dd;
735 269959130 : Double s = 1.0;
736 : Int n;
737 :
738 269959130 : if(r == 0)
739 : {
740 0 : *z = a->z[0];
741 0 : *m = 0.0;
742 0 : return 1;
743 : }
744 :
745 269959130 : if(r < 0)
746 : {
747 0 : s = -1.0;
748 0 : r = -r;
749 : }
750 269959130 : d = a->deltar;
751 269959130 : dd = d*d;
752 :
753 269959130 : n = (Int)floor(r/d + 0.5); /* the middle point */
754 269959130 : if(n > a->ngeom-2) n = a->ngeom-2;
755 :
756 269959130 : x = r - n*d;
757 :
758 269959130 : if(n == 0)
759 : {
760 15219 : mc = a->m[1];
761 15219 : ma = -mc;
762 15219 : mb = 0.0;
763 15219 : zav = 2.0*a->z[1] + a->z[0];
764 : }
765 : else
766 : {
767 269943911 : ma = a->m[n-1];
768 269943911 : mb = a->m[n];
769 269943911 : mc = a->m[n+1];
770 269943911 : zav = a->z[n-1] + a->z[n] + a->z[n+1];
771 : }
772 :
773 269959130 : A = mb;
774 269959130 : B = 0.5*(mc - ma)/d;
775 269959130 : C = 0.5*(mc - 2.0*mb + ma)/dd;
776 :
777 269959130 : D = (zav - B*dd)/3.0;
778 :
779 269959130 : if(m) *m = s*(A + B*x + C*x*x);
780 269959130 : if(z) *z = s*(D + A*x + B*x*x/2.0 + C*x*x*x/3.0);
781 :
782 269959130 : return 1;
783 : }
784 :
785 269957439 : Int BeamCalc::astigdishvalue(const calcAntenna *a, Double x, Double y, Double *z, Double *m)
786 : {
787 : Double ma, mb, mc, zav, A, B, C, D;
788 : Double r, rr, theta, xp, d, dd, z5, z6, astigm, dastigm;
789 269957439 : Double s = 1.0;
790 : Int n;
791 :
792 269957439 : rr = x*x + y*y;
793 269957439 : r = sqrt(rr);
794 :
795 269957439 : if(r==0. || (a->astigm_0==0. && a->astigm_45==0.))
796 : {
797 269957439 : return dishvalue(a, r, z, m);
798 : }
799 :
800 : // the Zernike polynomials Z5 and Z6
801 : Double sin2th, cos2th, rho, rho2;
802 :
803 0 : theta = atan2(y,x);
804 0 : sin2th = sin(2.*theta);
805 0 : cos2th = cos(2.*theta);
806 0 : rho = r / a->radius;
807 0 : rho2 = rho*rho;
808 :
809 0 : z5 = sqrt(6.) * rho2 * sin2th;
810 0 : z6 = sqrt(6.) * rho2 * cos2th;
811 :
812 0 : astigm = 1. + a->astigm_45 * z5 + a->astigm_0 * z6;
813 0 : dastigm = 2.* rho2/r * sqrt(6.)*(a->astigm_45*sin2th + a->astigm_0*cos2th);
814 :
815 0 : d = a->deltar;
816 0 : dd = d*d;
817 :
818 0 : n = (Int)floor(r/d + 0.5); /* the middle point */
819 0 : if(n > a->ngeom-2) n = a->ngeom-2;
820 :
821 0 : xp = r - n*d;
822 :
823 0 : if(n == 0)
824 : {
825 0 : mc = a->m[1];
826 0 : ma = -mc;
827 0 : mb = 0.0;
828 0 : zav = 2.0*a->z[1] + a->z[0];
829 : }
830 : else
831 : {
832 0 : ma = a->m[n-1];
833 0 : mb = a->m[n];
834 0 : mc = a->m[n+1];
835 0 : zav = a->z[n-1] + a->z[n] + a->z[n+1];
836 : }
837 :
838 0 : A = mb;
839 0 : B = 0.5*(mc - ma)/d;
840 0 : C = 0.5*(mc - 2.0*mb + ma)/dd;
841 :
842 0 : D = (zav - B*dd)/3.0;
843 :
844 :
845 0 : Double zn = s*(D + A*xp + B*xp*xp/2.0 + C*xp*xp*xp/3.0);
846 0 : if(z) *z = zn * astigm;
847 0 : if(m) *m = s*(A + B*xp + C*xp*xp) * astigm + dastigm * zn;
848 :
849 0 : return 1;
850 : }
851 :
852 : /* Returns position of subreflector piece (x, y, z) and
853 : * its normal (u, v, w)
854 : */
855 29995271 : Int BeamCalc::subfromdish(const calcAntenna *a, Double x, Double y, Double *subpoint)
856 : {
857 : Double r, z, m, u, v, w;
858 : Double dx, dy, dz, dl, t;
859 : Double n[3], sf[3], sd[3];
860 : Int i;
861 :
862 29995271 : r = sqrt(x*x + y*y);
863 :
864 29995271 : if(r == 0)
865 : {
866 0 : subpoint[0] = 0.0;
867 0 : subpoint[1] = 0.0;
868 0 : subpoint[2] = a->sub_h;
869 : }
870 : else
871 : {
872 29995271 : astigdishvalue(a, x, y, &z, &m);
873 :
874 : /* Compute reflected unit vector direction */
875 29995271 : m = tan(2.0*atan(m));
876 29995271 : w = 1.0/sqrt(1.0+m*m);
877 29995271 : u = -m*(x/r)*w;
878 29995271 : v = -m*(y/r)*w;
879 :
880 29995271 : dx = a->feed[0]-x;
881 29995271 : dy = a->feed[1]-y;
882 29995271 : dz = a->feed[2]-z;
883 29995271 : dl = a->K + z;
884 :
885 29995271 : t = 0.5*(dx*dx + dy*dy + dz*dz - dl*dl)
886 29995271 : / (-dl + u*dx + v*dy + w*dz);
887 :
888 29995271 : subpoint[0] = x + u*t;
889 29995271 : subpoint[1] = y + v*t;
890 29995271 : subpoint[2] = z + w*t;
891 : }
892 :
893 119981084 : for(i = 0; i < 3; i++) sf[i] = a->feed[i] - subpoint[i];
894 29995271 : sd[0] = x - subpoint[0];
895 29995271 : sd[1] = y - subpoint[1];
896 29995271 : sd[2] = z - subpoint[2];
897 :
898 29995271 : norm3(sf);
899 29995271 : norm3(sd);
900 :
901 119981084 : for(i = 0; i < 3; i++) n[i] = sd[i]+sf[i];
902 :
903 29995271 : norm3(n);
904 :
905 119981084 : for(i = 0; i < 3; i++) subpoint[3+i] = n[i];
906 :
907 29995271 : return 1;
908 : }
909 :
910 0 : Int BeamCalc::dishfromsub(const calcAntenna *a, Double x, Double y, Double *dishpoint)
911 : {
912 :
913 : Double x1, y1, dx, dy, mx, my, r, d;
914 0 : const Double eps = 0.001;
915 0 : Int iter, niter=500;
916 : Double sub[5][6];
917 :
918 0 : x1 = x;
919 0 : y1 = y;
920 :
921 0 : for(iter = 0; iter < niter; iter++)
922 : {
923 0 : subfromdish(a, x1, y1, sub[0]);
924 0 : subfromdish(a, x1-eps, y1, sub[1]);
925 0 : subfromdish(a, x1+eps, y1, sub[2]);
926 0 : subfromdish(a, x1, y1-eps, sub[3]);
927 0 : subfromdish(a, x1, y1+eps, sub[4]);
928 0 : mx = 0.5*(sub[2][0]-sub[1][0])/eps;
929 0 : my = 0.5*(sub[4][1]-sub[3][1])/eps;
930 0 : dx = (x-sub[0][0])/mx;
931 0 : dy = (y-sub[0][1])/my;
932 0 : if(fabs(dx) > a->radius/7.0)
933 : {
934 0 : if(dx < 0) dx = -a->radius/7.0;
935 0 : else dx = a->radius/7.0;
936 : }
937 0 : if(fabs(dy) > a->radius/7.0)
938 : {
939 0 : if(dy < 0) dy = -a->radius/7.0;
940 0 : else dy = a->radius/7.0;
941 : }
942 0 : r = sqrt(x1*x1 + y1*y1);
943 0 : if(r >= a->radius)
944 0 : if(x1*dx + y1*dy > 0.0) return 0;
945 0 : x1 += 0.5*dx;
946 0 : y1 += 0.5*dy;
947 0 : if(fabs(dx) < 0.005*eps && fabs(dy) < 0.005*eps) break;
948 : }
949 0 : if(iter == niter) return 0;
950 :
951 0 : r = sqrt(x1*x1 + y1*y1);
952 0 : dishpoint[0] = x1;
953 0 : dishpoint[1] = y1;
954 : // dishpoint[2] = polyvalue(a->shape, r);
955 0 : dishpoint[3] = sub[0][0];
956 0 : dishpoint[4] = sub[0][1];
957 0 : dishpoint[5] = sub[0][2];
958 0 : d = sqrt(1.0+mx*mx+my*my);
959 0 : dishpoint[6] = mx/d;
960 0 : dishpoint[7] = my/d;
961 0 : dishpoint[8] = 1.0/d;
962 0 : dishpoint[9] = sub[0][3];
963 0 : dishpoint[10] = sub[0][4];
964 0 : dishpoint[11] = sub[0][5];
965 :
966 0 : if(r > a->radius) return 0;
967 0 : else return 1;
968 : }
969 :
970 0 : void BeamCalc::printAntenna(const calcAntenna *a)
971 : {
972 0 : printf("Antenna: %s %p\n", a->name, a);
973 0 : printf(" freq = %f GHz lambda = %f m\n", a->freq, a->lambda);
974 0 : printf(" feeddir = %f, %f, %f\n",
975 0 : a->feeddir[0], a->feeddir[1], a->feeddir[2]);
976 0 : printf(" pfeeddir = %f, %f, %f\n",
977 0 : a->pfeeddir[0], a->pfeeddir[1], a->pfeeddir[2]);
978 0 : }
979 :
980 29995271 : Ray * BeamCalc::newRay(const Double *sub)
981 : {
982 : Ray *ray;
983 : Int i;
984 :
985 29995271 : ray = (Ray *)malloc(sizeof(Ray));
986 209966897 : for(i = 0; i < 6; i++) ray->sub[i] = sub[i];
987 :
988 29995271 : return ray;
989 : }
990 :
991 29995271 : void BeamCalc::deleteRay(Ray *ray)
992 : {
993 29995271 : if(ray) free(ray);
994 29995271 : }
995 :
996 1691 : Pathology* BeamCalc::newPathology()
997 : {
998 : Pathology *P;
999 : Int i, j;
1000 :
1001 1691 : P = (Pathology *)malloc(sizeof(Pathology));
1002 :
1003 6764 : for(i = 0; i < 3; i++) P->subrotpoint[i] = P->subshift[i] = P->feedshift[i] = 0.0;
1004 21983 : for(i = 0; i < 3; i++) for(j = 0; j < 3; j++)
1005 15219 : P->feedrot[i][j] = P->subrot[i][j] = 0.0;
1006 6764 : for(i = 0; i < 3; i++) P->feedrot[i][i] = P->subrot[i][i] = 1.0;
1007 :
1008 1691 : P->az_offset = 0.0;
1009 1691 : P->el_offset = 0.0;
1010 1691 : P->phase_offset = 0.0;
1011 1691 : P->focus = 0.0;
1012 :
1013 1691 : return P;
1014 : }
1015 :
1016 1691 : Pathology* BeamCalc::newPathologyfromApertureCalcParams(ApertureCalcParams* /*ap*/)
1017 : {
1018 : Pathology *P;
1019 :
1020 1691 : P = newPathology();
1021 :
1022 1691 : return P;
1023 : }
1024 :
1025 1691 : void BeamCalc::deletePathology(Pathology *P)
1026 : {
1027 1691 : if(P) free(P);
1028 1691 : }
1029 :
1030 9547368 : void BeamCalc::normvec(const Double *a, const Double *b, Double *c)
1031 : {
1032 : Int i;
1033 : Double r;
1034 38189472 : for(i = 0; i < 3; i++) c[i] = a[i] - b[i];
1035 9547368 : r = sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
1036 38189472 : for(i = 0; i < 3; i++) c[i] /= r;
1037 9547368 : }
1038 :
1039 0 : Double BeamCalc::dAdOmega(const calcAntenna *a, const Ray *ray1, const Ray *ray2,
1040 : const Ray *ray3, const Pathology *p)
1041 : {
1042 : Double A, Omega;
1043 : Double n1[3], n2[3], n3[3], f[3], ci, cj, ck;
1044 : Int i;
1045 :
1046 : /* Area in aperture is in a plane z = const */
1047 0 : A = 0.5*fabs(
1048 0 : (ray1->aper[0]-ray2->aper[0])*(ray1->aper[1]-ray3->aper[1]) -
1049 0 : (ray1->aper[0]-ray3->aper[0])*(ray1->aper[1]-ray2->aper[1]) );
1050 :
1051 0 : for(i = 0; i < 3; i++) f[i] = a->feed[i] + p->feedshift[i];
1052 :
1053 0 : normvec(ray1->sub, f, n1);
1054 0 : normvec(ray2->sub, f, n2);
1055 0 : normvec(ray3->sub, f, n3);
1056 :
1057 0 : for(i = 0; i < 3; i++)
1058 : {
1059 0 : n1[i] -= n3[i];
1060 0 : n2[i] -= n3[i];
1061 : }
1062 :
1063 0 : ci = n1[1]*n2[2] - n1[2]*n2[1];
1064 0 : cj = n1[2]*n2[0] - n1[0]*n2[2];
1065 0 : ck = n1[0]*n2[1] - n1[1]*n2[0];
1066 :
1067 0 : Omega = 0.5*sqrt(ci*ci + cj*cj + ck*ck);
1068 :
1069 0 : return A/Omega;
1070 : }
1071 :
1072 3182456 : Double BeamCalc::dOmega(const calcAntenna *a, const Ray *ray1, const Ray *ray2,
1073 : const Ray *ray3, const Pathology *p)
1074 : {
1075 : Double Omega;
1076 : Double n1[3], n2[3], n3[3], f[3], ci, cj, ck;
1077 : Int i;
1078 :
1079 12729824 : for(i = 0; i < 3; i++) f[i] = a->feed[i] + p->feedshift[i];
1080 :
1081 3182456 : normvec(ray1->sub, f, n1);
1082 3182456 : normvec(ray2->sub, f, n2);
1083 3182456 : normvec(ray3->sub, f, n3);
1084 :
1085 12729824 : for(i = 0; i < 3; i++)
1086 : {
1087 9547368 : n1[i] -= n3[i];
1088 9547368 : n2[i] -= n3[i];
1089 : }
1090 :
1091 3182456 : ci = n1[1]*n2[2] - n1[2]*n2[1];
1092 3182456 : cj = n1[2]*n2[0] - n1[0]*n2[2];
1093 3182456 : ck = n1[0]*n2[1] - n1[1]*n2[0];
1094 :
1095 3182456 : Omega = 0.5*sqrt(ci*ci + cj*cj + ck*ck);
1096 :
1097 3182456 : return Omega;
1098 : }
1099 :
1100 3184147 : Double BeamCalc::Raylen(const Ray *ray)
1101 : {
1102 : Double len, d[3];
1103 : Int i;
1104 :
1105 : /* feed to subreflector */
1106 12736588 : for(i = 0; i < 3; i++) d[i] = ray->feed[i] - ray->sub[i];
1107 3184147 : len = sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
1108 :
1109 : /* subreflector to dish */
1110 12736588 : for(i = 0; i < 3; i++) d[i] = ray->sub[i] - ray->dish[i];
1111 3184147 : len += sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
1112 :
1113 : /* dish to aperture */
1114 12736588 : for(i = 0; i < 3; i++) d[i] = ray->dish[i] - ray->aper[i];
1115 3184147 : len += sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
1116 :
1117 3184147 : return len;
1118 : }
1119 :
1120 29995271 : void BeamCalc::Pathologize(Double *sub, const Pathology *p)
1121 : {
1122 : Int i;
1123 : Int j;
1124 : Double tmp[6];
1125 :
1126 119981084 : for(i = 0; i < 3; i++) sub[i] -= p->subrotpoint[i];
1127 :
1128 119981084 : for(i = 0; i < 3; i++)
1129 : {
1130 89985813 : tmp[i] = 0.0;
1131 89985813 : tmp[i+3] = 0.0;
1132 359943252 : for(j = 0; j < 3; j++) tmp[i] += p->subrot[i][j]*sub[j];
1133 359943252 : for(j = 0; j < 3; j++) tmp[i+3] += p->subrot[i][j]*sub[j+3];
1134 : }
1135 :
1136 119981084 : for(i = 0; i < 3; i++)
1137 89985813 : sub[i] = tmp[i] + p->subrotpoint[i] + p->subshift[i];
1138 89985813 : for(i = 4; i < 6; i++)
1139 59990542 : sub[i] = tmp[i];
1140 29995271 : }
1141 :
1142 0 : void BeamCalc::applyPathology(Pathology *P, calcAntenna *a)
1143 : {
1144 : Double dx[3];
1145 : Int i, j;
1146 :
1147 0 : if(P->focus != 0.0)
1148 : {
1149 0 : dx[0] = -a->feed[0];
1150 0 : dx[1] = -a->feed[1];
1151 0 : dx[2] = a->sub_h-a->feed[2];
1152 0 : norm3(dx);
1153 0 : for(i = 0; i < 3; i++) P->feedshift[i] += P->focus*dx[i];
1154 :
1155 0 : P->focus = 0.0;
1156 : }
1157 :
1158 0 : for(i = 0; i < 3; i++) a->pfeeddir[i] = 0.0;
1159 0 : for(j = 0; j < 3; j++) for(i = 0; i < 3; i++)
1160 0 : a->pfeeddir[j] += P->feedrot[j][i]*a->feeddir[i];
1161 0 : }
1162 :
1163 :
1164 29995271 : void BeamCalc::intersectdish(const calcAntenna *a, const Double *sub, const Double *unitdir,
1165 : Double *dish, Int niter)
1166 : {
1167 : Double A, B, C, t, m, r;
1168 : Double x[3], n[3];
1169 : Int i, iter;
1170 :
1171 : /* First intersect with ideal paraboloid */
1172 29995271 : A = a->bestparabola*(unitdir[0]*unitdir[0] + unitdir[1]*unitdir[1]);
1173 29995271 : B = 2.0*a->bestparabola*(unitdir[0]*sub[0] + unitdir[1]*sub[1])
1174 29995271 : -unitdir[2];
1175 29995271 : C = a->bestparabola*(sub[0]*sub[0] + sub[1]*sub[1]) - sub[2];
1176 29995271 : t = 0.5*(sqrt(B*B-4.0*A*C) - B)/A; /* take greater root */
1177 :
1178 29995271 : for(iter = 0; ; iter++)
1179 : {
1180 : /* get position (x) and normal (n) on the real dish */
1181 719886504 : for(i = 0; i < 2; i++) x[i] = sub[i] + t*unitdir[i];
1182 239962168 : r = sqrt(x[0]*x[0] + x[1]*x[1]);
1183 239962168 : astigdishvalue(a, x[0], x[1], &(x[2]), &m);
1184 239962168 : n[2] = 1.0/sqrt(1.0+m*m);
1185 239962168 : n[0] = -m*(x[0]/r)*n[2];
1186 239962168 : n[1] = -m*(x[1]/r)*n[2];
1187 :
1188 239962168 : if(iter >= niter) break;
1189 :
1190 209966897 : A = B = 0;
1191 839867588 : for(i = 0; i < 3; i++)
1192 : {
1193 629900691 : A += n[i]*(x[i]-sub[i]); /* n dot (x-sub) */
1194 629900691 : B += n[i]*unitdir[i]; /* n dot unitdir */
1195 : }
1196 209966897 : t = A/B;
1197 : }
1198 :
1199 119981084 : for(i = 0; i < 3; i++)
1200 : {
1201 89985813 : dish[i] = x[i];
1202 89985813 : dish[i+3] = n[i];
1203 : }
1204 29995271 : }
1205 :
1206 29995271 : void BeamCalc::intersectaperture(const calcAntenna *a, const Double *dish,
1207 : const Double *unitdir, Double *aper)
1208 : {
1209 : Double t;
1210 : Int i;
1211 :
1212 29995271 : t = (a->zedge-dish[2])/unitdir[2];
1213 119981084 : for(i = 0; i < 3; i++) aper[i] = dish[i] + t*unitdir[i];
1214 :
1215 29995271 : aper[3] = aper[4] = 0.0;
1216 29995271 : aper[5] = 1.0;
1217 29995271 : }
1218 :
1219 : /* gain in power */
1220 0 : Double BeamCalc::feedfunc(const calcAntenna *a, Double theta)
1221 : {
1222 : Double stheta;
1223 :
1224 0 : stheta = sin(theta);
1225 0 : return exp(2.0*(-0.083)*a->fa2pi*a->fa2pi*stheta*stheta);
1226 : }
1227 :
1228 : /* gain in power */
1229 3182456 : Double BeamCalc::feedgain(const calcAntenna *a, const Ray *ray, const Pathology */*p*/)
1230 : {
1231 3182456 : Double costheta = 0.0;
1232 : Double v[3];
1233 : Int i;
1234 :
1235 12729824 : for(i = 0; i < 3; i++) v[i] = ray->sub[i] - ray->feed[i];
1236 3182456 : norm3(v);
1237 :
1238 12729824 : for(i = 0; i < 3; i++)
1239 : {
1240 9547368 : costheta += a->pfeeddir[i]*v[i];
1241 : }
1242 :
1243 :
1244 3182456 : return exp(2.0*(-0.083)*a->fa2pi*a->fa2pi*(1.0-costheta*costheta));
1245 : }
1246 :
1247 29995271 : Ray* BeamCalc::trace(const calcAntenna *a, Double x, Double y, const Pathology *p)
1248 : {
1249 : Ray *ray;
1250 : Double idealsub[12];
1251 29995271 : Double fu[3], du[3], au[3], ndotf=0.0, ndotd=0.0;
1252 : Int i;
1253 29995271 : const Int niter = 7;
1254 :
1255 29995271 : subfromdish(a, x, y, idealsub);
1256 :
1257 29995271 : ray = newRay(idealsub);
1258 :
1259 29995271 : Pathologize(ray->sub, p);
1260 :
1261 29995271 : if(ray->sub[5] < -1.0 || ray->sub[5] > -0.0)
1262 : {
1263 0 : deleteRay(ray);
1264 0 : return 0;
1265 : }
1266 :
1267 119981084 : for(i = 0; i < 3; i++) ray->feed[i] = a->feed[i] + p->feedshift[i];
1268 :
1269 : /* now determine unit vector pointing to dish */
1270 :
1271 : /* unit toward feed */
1272 119981084 : for(i = 0; i < 3; i++) fu[i] = ray->feed[i] - ray->sub[i];
1273 29995271 : norm3(fu);
1274 :
1275 : /* unit toward dish */
1276 119981084 : for(i = 0; i < 3; i++) ndotf += ray->sub[i+3]*fu[i];
1277 119981084 : for(i = 0; i < 3; i++) du[i] = 2.0*ray->sub[i+3]*ndotf - fu[i];
1278 :
1279 : /* dish point */
1280 29995271 : intersectdish(a, ray->sub, du, ray->dish, niter);
1281 :
1282 : /* unit toward aperture */
1283 119981084 : for(i = 0; i < 3; i++) ndotd += ray->dish[i+3]*du[i];
1284 119981084 : for(i = 0; i < 3; i++) au[i] = du[i] - 2.0*ray->dish[i+3]*ndotd;
1285 :
1286 29995271 : intersectaperture(a, ray->dish, au, ray->aper);
1287 :
1288 29995271 : return ray;
1289 : }
1290 :
1291 6364912 : void BeamCalc::tracepol(Complex *E0, const Ray *ray, Complex *E1)
1292 : {
1293 6364912 : Complex fac;
1294 : Double v1[3], v2[3], v3[3];
1295 : Double r[3];
1296 : Int i;
1297 :
1298 25459648 : for(i = 0; i < 3; i++)
1299 : {
1300 19094736 : v1[i] = ray->sub[i] - ray->feed[i];
1301 19094736 : v2[i] = ray->dish[i] - ray->sub[i];
1302 19094736 : v3[i] = ray->aper[i] - ray->dish[i];
1303 19094736 : E1[i] = E0[i];
1304 : }
1305 6364912 : norm3(v1);
1306 6364912 : norm3(v2);
1307 6364912 : norm3(v3);
1308 :
1309 25459648 : for(i = 0; i < 3; i++) r[i] = v1[i] - v2[i];
1310 6364912 : norm3(r);
1311 6364912 : fac = Complex(r[0],0)*E1[0] + Complex(r[1],0)*E1[1] + Complex(r[2],0)*E1[2];
1312 25459648 : for(i = 0; i < 3; i++) E1[i] = Complex(r[i],0)*fac*2.0 - E1[i];
1313 :
1314 25459648 : for(i = 0; i < 3; i++) r[i] = v2[i] - v3[i];
1315 6364912 : norm3(r);
1316 6364912 : fac = Complex(r[0],0)*E1[0] + Complex(r[1],0)*E1[1] + Complex(r[2],0)*E1[2];
1317 25459648 : for(i = 0; i < 3; i++) E1[i] = Complex(r[i],0)*fac*2.0 - E1[i];
1318 6364912 : }
1319 :
1320 0 : Int BeamCalc::legplanewaveblock(const calcAntenna *a, Double x, Double y)
1321 : {
1322 : /* outside the leg foot area, the blockage is spherical wave */
1323 0 : if(x*x + y*y > a->legfoot*a->legfoot) return 0;
1324 :
1325 0 : if(a->legwidth == 0.0) return 0;
1326 :
1327 0 : if(strcmp(a->name, "VLBA") == 0)
1328 : {
1329 0 : const Double s=1.457937;
1330 0 : const Double c=1.369094;
1331 0 : if(fabs(s*x+c*y) < -a->legwidth) return 1;
1332 0 : if(fabs(s*x-c*y) < -a->legwidth) return 1;
1333 : }
1334 0 : else if(a->legwidth < 0.0) /* "x shaped" legs */
1335 : {
1336 0 : if(fabs(x-y)*M_SQRT2 < -a->legwidth) return 1;
1337 0 : if(fabs(x+y)*M_SQRT2 < -a->legwidth) return 1;
1338 : }
1339 0 : else if(a->legwidth > 0.0) /* "+ shaped" legs */
1340 : {
1341 0 : if(fabs(x)*2.0 < a->legwidth) return 1;
1342 0 : if(fabs(y)*2.0 < a->legwidth) return 1;
1343 : }
1344 :
1345 0 : return 0;
1346 : }
1347 :
1348 3375524 : Int BeamCalc::legplanewaveblock2(const calcAntenna *a, const Ray *ray)
1349 : {
1350 : Int i, n;
1351 : Double dr2;
1352 : // phi set but not used
1353 : Double theta /*, phi*/;
1354 : Double r0[3], dr[3], l0[3], l1[3], dl[3], D[3];
1355 : Double D2, N[3], ll, rr;
1356 3375524 : const Double thetaplus[4] =
1357 : {0, M_PI/2.0, M_PI, 3.0*M_PI/2.0};
1358 3375524 : const Double thetacross[4] =
1359 : {0.25*M_PI, 0.75*M_PI, 1.25*M_PI, 1.75*M_PI};
1360 3375524 : const Double thetavlba[4] =
1361 : {0.816817, 2.3247756, 3.9584096, 5.466368};
1362 : const Double *thetalut;
1363 :
1364 3375524 : if(a->legwidth == 0.0) return 0;
1365 :
1366 3375524 : if(strcmp(a->name, "VLBA") == 0) thetalut = thetavlba;
1367 3375524 : else if(a->legwidth < 0.0) thetalut = thetacross;
1368 3375524 : else thetalut = thetaplus;
1369 :
1370 : /* inside the leg feet is plane wave blockage */
1371 3375524 : dr2 = ray->dish[0]*ray->dish[0] + ray->dish[1]*ray->dish[1];
1372 3375524 : if(dr2 >= a->legfoot*a->legfoot) return 0;
1373 :
1374 4469424 : for(i = 0; i < 3; i++)
1375 : {
1376 3352068 : r0[i] = ray->dish[i];
1377 3352068 : dr[i] = ray->aper[i] - r0[i];
1378 : }
1379 1117356 : rr = r0[0]*r0[0] + r0[1]*r0[1];
1380 :
1381 1117356 : l0[2] = a->legfootz;
1382 1117356 : l1[0] = l1[1] = 0.0;
1383 1117356 : l1[2] = a->legapex;
1384 : // phi set but not used
1385 : // phi = atan2(r0[1], r0[0]);
1386 5483080 : for(n = 0; n < 4; n++)
1387 : {
1388 4407204 : theta = thetalut[n];
1389 4407204 : l0[0] = a->legfoot*cos(theta);
1390 4407204 : l0[1] = a->legfoot*sin(theta);
1391 4407204 : ll = l0[0]*l0[0] + l0[1]*l0[1];
1392 4407204 : if((l0[0]*r0[0] + l0[1]*r0[1])/sqrt(ll*rr) < 0.7) continue;
1393 4543680 : for(i = 0; i < 3; i++) dl[i] = l1[i] - l0[i];
1394 4543680 : for(i = 0; i < 3; i++) D[i] = r0[i] - l0[i];
1395 :
1396 1135920 : N[0] = dr[1]*dl[2] - dr[2]*dl[1];
1397 1135920 : N[1] = dr[2]*dl[0] - dr[0]*dl[2];
1398 1135920 : N[2] = dr[0]*dl[1] - dr[1]*dl[0];
1399 1135920 : norm3(N);
1400 :
1401 1135920 : D2 = D[0]*N[0] + D[1]*N[1] + D[2]*N[2];
1402 :
1403 1135920 : if(fabs(D2) <= 0.5*fabs(a->legwidth)) return 1;
1404 : }
1405 :
1406 1075876 : return 0;
1407 : }
1408 :
1409 3334044 : Int BeamCalc::legsphericalwaveblock(const calcAntenna *a, const Ray *ray)
1410 : {
1411 : Int i, n;
1412 : Double dr2;
1413 : // phi set but not used
1414 : Double theta /*, phi*/;
1415 : Double r0[3], dr[3], l0[3], l1[3], dl[3], D[3];
1416 : Double D2, N[3], ll, rr;
1417 3334044 : const Double thetaplus[4] =
1418 : {0, M_PI/2.0, M_PI, 3.0*M_PI/2.0};
1419 3334044 : const Double thetacross[4] =
1420 : {0.25*M_PI, 0.75*M_PI, 1.25*M_PI, 1.75*M_PI};
1421 3334044 : const Double thetavlba[4] =
1422 : {0.816817, 2.3247756, 3.9584096, 5.466368};
1423 : const Double *thetalut;
1424 :
1425 3334044 : if(a->legwidth == 0.0) return 0;
1426 :
1427 3334044 : if(strcmp(a->name, "VLBA") == 0) thetalut = thetavlba;
1428 3334044 : else if(a->legwidth < 0.0) thetalut = thetacross;
1429 3334044 : else thetalut = thetaplus;
1430 :
1431 : /* inside the leg feet is plane wave blockage */
1432 3334044 : dr2 = ray->dish[0]*ray->dish[0] + ray->dish[1]*ray->dish[1];
1433 3334044 : if(dr2 < a->legfoot*a->legfoot) return 0;
1434 :
1435 9032672 : for(i = 0; i < 3; i++)
1436 : {
1437 6774504 : r0[i] = ray->dish[i];
1438 6774504 : dr[i] = ray->sub[i] - r0[i];
1439 : }
1440 2258168 : rr = r0[0]*r0[0] + r0[1]*r0[1];
1441 :
1442 2258168 : l0[2] = a->legfootz;
1443 2258168 : l1[0] = l1[1] = 0.0;
1444 2258168 : l1[2] = a->legapex;
1445 : // phi set but not used
1446 : // phi = atan2(r0[1], r0[0]);
1447 10911870 : for(n = 0; n < 4; n++)
1448 : {
1449 8805290 : theta = thetalut[n];
1450 8805290 : l0[0] = a->legfoot*cos(theta);
1451 8805290 : l0[1] = a->legfoot*sin(theta);
1452 8805290 : ll = l0[0]*l0[0] + l0[1]*l0[1];
1453 8805290 : if((l0[0]*r0[0] + l0[1]*r0[1])/sqrt(ll*rr) < 0.7) continue;
1454 9159536 : for(i = 0; i < 3; i++) dl[i] = l1[i] - l0[i];
1455 9159536 : for(i = 0; i < 3; i++) D[i] = r0[i] - l0[i];
1456 :
1457 2289884 : N[0] = dr[1]*dl[2] - dr[2]*dl[1];
1458 2289884 : N[1] = dr[2]*dl[0] - dr[0]*dl[2];
1459 2289884 : N[2] = dr[0]*dl[1] - dr[1]*dl[0];
1460 2289884 : norm3(N);
1461 :
1462 2289884 : D2 = D[0]*N[0] + D[1]*N[1] + D[2]*N[2];
1463 :
1464 2289884 : if(fabs(D2) <= 0.5*fabs(a->legwidth)) return 1;
1465 : }
1466 :
1467 2106580 : return 0;
1468 : }
1469 :
1470 :
1471 32 : void BeamCalc::copyBeamCalcGeometry(BeamCalcGeometry* to, BeamCalcGeometry* from){
1472 32 : sprintf(to->name, "%s", from->name);
1473 32 : sprintf(to->bandName, "%s", from->bandName);
1474 32 : to->sub_h = from->sub_h;
1475 128 : for(uInt j=0; j<3;j++){
1476 96 : to->feedpos[j] = from->feedpos[j];
1477 : }
1478 32 : to->subangle = from->subangle;
1479 32 : to->legwidth = from->legwidth;
1480 32 : to->legfoot = from->legfoot;
1481 32 : to->legapex = from->legapex;
1482 32 : to->Rhole = from->Rhole;
1483 32 : to->Rant = from->Rant;
1484 32 : to->reffreq = from->reffreq;
1485 192 : for(uInt j=0; j<5;j++){
1486 160 : to->taperpoly[j] = from->taperpoly[j];
1487 : }
1488 32 : to->ntaperpoly = from->ntaperpoly;
1489 32 : to->astigm_0 = from->astigm_0;
1490 32 : to->astigm_45 = from->astigm_45;
1491 :
1492 32 : }
1493 :
1494 :
1495 : /* The meat of the calculation */
1496 :
1497 :
1498 1691 : Int BeamCalc::calculateAperture(ApertureCalcParams *ap)
1499 : {
1500 : // not used
1501 : // Complex fp, Exr, Eyr, Exl, Eyl;
1502 1691 : Complex Er[3], El[3];
1503 1691 : Complex Pr[2], Pl[2];
1504 1691 : Complex q[2];
1505 : //Double dx, dy, Rhole, Rant, x0, y0, R2, H2, eps;
1506 : //Complex rr, rl, lr, ll, tmp;
1507 : Double L0, phase;
1508 : Double sp, cp;
1509 : Double B[3][3];
1510 : calcAntenna *a;
1511 : Pathology *p;
1512 : Int nx, ny, os;
1513 : Int i, j;
1514 : //Double pac, pas; /* parallactic angle cosine / sine */
1515 1691 : Complex Iota; Iota=Complex(0,1);
1516 :
1517 : //UNUSED: Complex E1[3];
1518 : //UNUSED: Double x,y, r2, L, amp, dP, dA, d0, x1, y1, dx1, dy1, dx2, dy2, dO;
1519 : //UNUSED: Ray *ray, *rayx, *rayy;
1520 : //UNUSED: Int iter;
1521 : //UNUSED: Int niter=6;
1522 :
1523 1691 : a = newAntennafromApertureCalcParams(ap);
1524 1691 : p = newPathologyfromApertureCalcParams(ap);
1525 :
1526 : /* compute central ray pathlength */
1527 : {
1528 : Ray *tmpRay;
1529 1691 : tmpRay = trace(a, 0.0, 0.00001, p);
1530 1691 : L0 = Raylen(tmpRay);
1531 1691 : deleteRay(tmpRay);
1532 : }
1533 :
1534 : //pac = cos(ap->pa+M_PI/2);
1535 : //pas = sin(ap->pa+M_PI/2);
1536 :
1537 1691 : if(obsName_p=="EVLA" || obsName_p=="VLA"){
1538 : /* compute polarization vectors in circular basis */
1539 1691 : Pr[0] = 1.0/M_SQRT2; Pr[1] = Iota/M_SQRT2;
1540 1691 : Pl[0] = 1.0/M_SQRT2; Pl[1] = -Iota/M_SQRT2;
1541 :
1542 : /* compensate for feed orientation */
1543 1691 : getfeedbasis(a, B);
1544 1691 : phase = atan2(B[0][1], B[0][0]);
1545 1691 : cp = cos(phase);
1546 1691 : sp = sin(phase);
1547 :
1548 1691 : q[0] = Pr[0];
1549 1691 : q[1] = Pr[1];
1550 1691 : Pr[0] = Complex(cp,0)*q[0] + Complex(sp,0)*q[1];
1551 1691 : Pr[1] = -Complex(sp,0)*q[0] + Complex(cp,0)*q[1];
1552 1691 : q[0] = Pl[0];
1553 1691 : q[1] = Pl[1];
1554 1691 : Pl[0] = Complex(cp,0)*q[0] + Complex(sp,0)*q[1];
1555 1691 : Pl[1] = -Complex(sp,0)*q[0] + Complex(cp,0)*q[1];
1556 : }
1557 : else{
1558 : /* in linear basis */
1559 0 : Pr[0] = 1.0; Pr[1] = 0.0;
1560 0 : Pl[0] = 0.0; Pl[1] = 1.0;
1561 : }
1562 :
1563 :
1564 :
1565 : /* compute 3-vector feed efields for the two polarizations */
1566 1691 : Efield(a, Pr, Er);
1567 1691 : Efield(a, Pl, El);
1568 1691 : ap->aperture->set(Complex(0.0));
1569 :
1570 1691 : os = ap->oversamp;
1571 1691 : nx = ap->nx*os;
1572 1691 : ny = ap->ny*os;
1573 : //dx = ap->dx/os;
1574 : //dy = ap->dy/os;
1575 : //x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
1576 : //y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
1577 : //Rant = a->radius;
1578 : //Rhole = a->hole_radius;
1579 : //R2 = Rant*Rant;
1580 : //H2 = Rhole*Rhole;
1581 :
1582 : //eps = dx/4.0;
1583 :
1584 1691 : IPosition pos(4);
1585 : // shape = ap->aperture->shape();
1586 :
1587 :
1588 : // cerr << "max threads " << omp_get_max_threads()
1589 : // << " threads available " << omp_get_num_threads()
1590 : // << endl;
1591 : //Int Nth=1;
1592 : //#ifdef _OPENMP
1593 : //Nth=max(omp_get_max_threads()-2,1);
1594 : //#endif
1595 : // Timer tim;
1596 : // tim.mark();
1597 : //#pragma omp parallel default(none) firstprivate(Er, El, nx, ny) private(i,j) shared(ap, a, p, L0) num_threads(Nth)
1598 : {
1599 : //#pragma omp for
1600 2980187 : for(j = 0; j < ny; j++)
1601 : {
1602 6630526656 : for(i = 0; i < nx; i++)
1603 : {
1604 6627548160 : computePixelValues(ap, a, p, L0, Er, El, i,j);
1605 : }
1606 : }
1607 : }
1608 : // tim.show("BeamCalc:");
1609 :
1610 1691 : deletePathology(p);
1611 1691 : deleteAntenna(a);
1612 :
1613 1691 : return 1;
1614 1691 : }
1615 :
1616 6627548160 : void BeamCalc::computePixelValues(const ApertureCalcParams *ap,
1617 : const calcAntenna *a, const Pathology *p,
1618 : const Double &L0,
1619 : Complex *Er, Complex *El,
1620 : const Int &i, const Int &j)
1621 : {
1622 6627548160 : Complex fp, Exr, Eyr, Exl, Eyl;
1623 : // Complex Er[3], El[3];
1624 6627548160 : Complex E1[3];
1625 : Double dx, dy, Rant, x0, y0, x, y, r2, R2, H2, eps, Rhole;
1626 6627548160 : Complex rr, rl, lr, ll, tmp;
1627 : Double L, amp, dP, dA, dO, x1, y1, dx1, dy1, dx2, dy2, phase;
1628 : //Int nx, ny;
1629 : Int os;
1630 6627548160 : Int niter=6;
1631 : Double pac, pas;
1632 : Double cp,sp; /* parallactic angle cosine / sine */
1633 6627548160 : Complex Iota; Iota=Complex(0,1);
1634 6627548160 : IPosition pos(4);pos=0;
1635 :
1636 6627548160 : Ray *ray=0, *rayx=0, *rayy=0;
1637 : /* determine parallactic angle rotated coordinates */
1638 :
1639 6627548160 : os = ap->oversamp;
1640 : //nx = ap->nx*os;
1641 : //ny = ap->ny*os;
1642 6627548160 : dx = ap->dx/os;
1643 6627548160 : dy = ap->dy/os;
1644 6627548160 : x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
1645 6627548160 : y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
1646 6627548160 : Rant = a->radius;
1647 6627548160 : Rhole = a->hole_radius;
1648 6627548160 : R2 = Rant*Rant;
1649 6627548160 : H2 = Rhole*Rhole;
1650 : // for(Int i=0; i < nx; ++i)
1651 : {
1652 6627548160 : eps = dx/4.0;
1653 6627548160 : pac = cos(ap->pa+M_PI/2);
1654 6627548160 : pas = sin(ap->pa+M_PI/2);
1655 :
1656 6627548160 : x = pac*(x0 + i*dx) - pas*(y0 + j*dy);
1657 6627548160 : y = pas*(x0 + i*dx) + pac*(y0 + j*dy);
1658 6627548160 : x = -x;
1659 :
1660 6627548160 : if(fabs(x) > Rant) goto nextpoint;
1661 122325652 : if(fabs(y) > Rant) goto nextpoint;
1662 4528707 : r2 = x*x + y*y;
1663 4528707 : if(r2 > R2) goto nextpoint;
1664 3555387 : if(r2 < H2) goto nextpoint;
1665 :
1666 3375524 : ray = rayx = rayy = 0;
1667 :
1668 3375524 : x1 = x;
1669 3375524 : y1 = y;
1670 :
1671 23628668 : for(Int iter = 0; iter < niter; iter++)
1672 : {
1673 20253144 : ray = trace(a, x1, y1, p);
1674 20253144 : if(!ray) goto nextpoint;
1675 20253144 : x1 += (x - ray->aper[0]);
1676 20253144 : y1 += (y - ray->aper[1]);
1677 20253144 : deleteRay(ray);
1678 20253144 : ray = 0;
1679 : }
1680 :
1681 3375524 : ray = trace(a, x1, y1, p);
1682 :
1683 : /* check for leg blockage */
1684 3375524 : if(legplanewaveblock2(a, ray))
1685 41480 : goto nextpoint;
1686 3334044 : if(legsphericalwaveblock(a, ray))
1687 151588 : goto nextpoint;
1688 :
1689 3182456 : if(y < 0) rayy = trace(a, x1, y1+eps, p);
1690 1591228 : else rayy = trace(a, x1, y1-eps, p);
1691 :
1692 3182456 : if(x < 0) rayx = trace(a, x1+eps, y1, p);
1693 1591228 : else rayx = trace(a, x1-eps, y1, p);
1694 :
1695 3182456 : if(ray == 0 || rayx == 0 || rayy == 0)
1696 0 : goto nextpoint;
1697 :
1698 : /* compute solid angle subtended at the feed */
1699 3182456 : dx1 = rayx->aper[0]-ray->aper[0];
1700 3182456 : dy1 = rayx->aper[1]-ray->aper[1];
1701 3182456 : dx2 = rayy->aper[0]-ray->aper[0];
1702 3182456 : dy2 = rayy->aper[1]-ray->aper[1];
1703 :
1704 3182456 : dA = 0.5*fabs(dx1*dy2 - dx2*dy1);
1705 3182456 : dO = (dOmega(a, rayx, rayy, ray, p)/dA)*dx*dx;
1706 3182456 : dP = dO*feedgain(a, ray, p);
1707 3182456 : amp = sqrt(dP);
1708 :
1709 3182456 : L = Raylen(ray);
1710 :
1711 3182456 : phase = 2.0*M_PI*(L-L0)/a->lambda;
1712 :
1713 : /* phase retard the wave */
1714 3182456 : cp = cos(phase);
1715 3182456 : sp = sin(phase);
1716 : // fp = cp + sp*1.0i;
1717 :
1718 3182456 : fp = Complex(cp,sp);
1719 :
1720 :
1721 3182456 : tracepol(Er, ray, E1);
1722 3182456 : Exr = fp*amp*E1[0];
1723 3182456 : Eyr = fp*amp*E1[1];
1724 :
1725 : // rr = Exr - 1.0i*Eyr;
1726 : // rl = Exr + 1.0i*Eyr;
1727 3182456 : rr = Exr - Iota*Eyr;
1728 3182456 : rl = Exr + Iota*Eyr;
1729 :
1730 3182456 : tracepol(El, ray, E1);
1731 3182456 : Exl = fp*amp*E1[0];
1732 3182456 : Eyl = fp*amp*E1[1];
1733 : // lr = Exl - 1.0i*Eyl;
1734 : // ll = Exl + 1.0i*Eyl;
1735 3182456 : lr = Exl - Iota*Eyl;
1736 3182456 : ll = Exl + Iota*Eyl;
1737 : // pos(0)=(Int)((j/os) - (25.0/dy/os)/2 + shape(0)/2 - 0.5);
1738 : // pos(1)=(Int)((i/os) - (25.0/dx/os)/2 + shape(1)/2 - 0.5);
1739 : // Following 3 lines go with ANT tag in VLACalc.....
1740 : // Double antRadius=BeamCalcGeometryDefaults[ap->band].Rant;
1741 : // pos(0)=(Int)((j/os) - (antRadius/dy/os) + shape(0)/2 - 0.5);
1742 : // pos(1)=(Int)((i/os) - (antRadius/dx/os) + shape(1)/2 - 0.5);
1743 : // Following 2 lines go with the PIX tag in VLACalc...
1744 3182456 : pos(0)=(Int)((j/os));
1745 3182456 : pos(1)=(Int)((i/os));
1746 3182456 : pos(3)=0;
1747 :
1748 3182456 : pos(2)=0;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+rr,pos);
1749 3182456 : pos(2)=1;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+rl,pos);
1750 3182456 : pos(2)=2;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+lr,pos);
1751 3182456 : pos(2)=3;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+ll,pos);
1752 6627548160 : nextpoint:
1753 6627548160 : if(ray) deleteRay(ray);
1754 6627548160 : if(rayx) deleteRay(rayx);
1755 6627548160 : if(rayy) deleteRay(rayy);
1756 : }
1757 6627548160 : }
1758 : //
1759 : //----------------------------------------------------------------------
1760 : // Compute only the required polarizations.
1761 : //
1762 0 : Int BeamCalc::calculateAperture(ApertureCalcParams *ap, const Int& whichPoln)
1763 : {
1764 : //Complex fp, Exr, Eyr, Exl, Eyl;
1765 0 : Complex Er[3], El[3];
1766 0 : Complex Pr[2], Pl[2];
1767 0 : Complex q[2];
1768 : //Double dx, dy, Rhole, Rant;//x0, y0, R2, H2, eps;
1769 : //Complex rr, rl, lr, ll, tmp;
1770 : Double L0, phase;
1771 : Double sp, cp;
1772 : Double B[3][3];
1773 : calcAntenna *a;
1774 : Pathology *p;
1775 : Int nx, ny, os;
1776 : Int i, j;
1777 : //Double pac, pas; /* parallactic angle cosine / sine */
1778 0 : Complex Iota; Iota=Complex(0,1);
1779 :
1780 : //UNUSED: Complex E1[3];
1781 : //UNUSED: Double x,y, r2, L, amp, dP, dA, d0, x1, y1, dx1, dy1, dx2, dy2, dO;
1782 : //UNUSED: Ray *ray, *rayx, *rayy;
1783 : //UNUSED: Int iter;
1784 : //UNUSED: Int niter=6;
1785 :
1786 0 : a = newAntennafromApertureCalcParams(ap);
1787 0 : p = newPathologyfromApertureCalcParams(ap);
1788 :
1789 : /* compute central ray pathlength */
1790 : {
1791 : Ray *tmpRay;
1792 0 : tmpRay = trace(a, 0.0, 0.00001, p);
1793 0 : L0 = Raylen(tmpRay);
1794 0 : deleteRay(tmpRay);
1795 : }
1796 :
1797 : //pac = cos(ap->pa+M_PI/2);
1798 : //pas = sin(ap->pa+M_PI/2);
1799 :
1800 0 : if(obsName_p=="EVLA" || obsName_p=="VLA"){
1801 : /* compute polarization vectors in circular basis */
1802 0 : Pr[0] = 1.0/M_SQRT2; Pr[1] = Iota/M_SQRT2;
1803 0 : Pl[0] = 1.0/M_SQRT2; Pl[1] = -Iota/M_SQRT2;
1804 :
1805 : /* compensate for feed orientation */
1806 0 : getfeedbasis(a, B);
1807 0 : phase = atan2(B[0][1], B[0][0]);
1808 0 : cp = cos(phase);
1809 0 : sp = sin(phase);
1810 :
1811 0 : q[0] = Pr[0];
1812 0 : q[1] = Pr[1];
1813 0 : Pr[0] = Complex(cp,0)*q[0] + Complex(sp,0)*q[1];
1814 0 : Pr[1] = -Complex(sp,0)*q[0] + Complex(cp,0)*q[1];
1815 0 : q[0] = Pl[0];
1816 0 : q[1] = Pl[1];
1817 0 : Pl[0] = Complex(cp,0)*q[0] + Complex(sp,0)*q[1];
1818 0 : Pl[1] = -Complex(sp,0)*q[0] + Complex(cp,0)*q[1];
1819 : }
1820 : else{
1821 : /* in linear basis */
1822 0 : Pr[0] = 1.0; Pr[1] = 0.0;
1823 0 : Pl[0] = 0.0; Pl[1] = 1.0;
1824 : }
1825 :
1826 :
1827 :
1828 : /* compute 3-vector feed efields for the two polarizations */
1829 0 : if ((whichPoln == Stokes::RR) || (whichPoln == Stokes::XX))
1830 0 : Efield(a, Pr, Er);
1831 0 : else if ((whichPoln == Stokes::LL) || (whichPoln == Stokes::YY))
1832 0 : Efield(a, Pl, El);
1833 : else
1834 0 : {Efield(a, Pr, Er); Efield(a, Pl, El);}
1835 :
1836 0 : ap->aperture->set(Complex(0.0));
1837 :
1838 0 : os = ap->oversamp;
1839 0 : nx = ap->nx*os;
1840 0 : ny = ap->ny*os;
1841 : // dx = ap->dx/os;
1842 : // dy = ap->dy/os;
1843 : //x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
1844 : //y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
1845 : // Rant = a->radius;
1846 : // Rhole = a->hole_radius;
1847 : //R2 = Rant*Rant;
1848 : //H2 = Rhole*Rhole;
1849 :
1850 : //eps = dx/4.0;
1851 :
1852 0 : IPosition pos(4);
1853 : // shape = ap->aperture->shape();
1854 :
1855 :
1856 : // cerr << "max threads " << omp_get_max_threads()
1857 : // << " threads available " << omp_get_num_threads()
1858 : // << endl;
1859 : // Int Nth=1, localWhichPoln=whichPoln;
1860 0 : Int localWhichPoln=whichPoln;
1861 : //#ifdef _OPENMP
1862 : //Nth=max(omp_get_max_threads()-2,1);
1863 : //#endif
1864 : // Timer tim;
1865 : // tim.mark();
1866 : //#if GCC44x
1867 : //#pragma omp parallel default(none) firstprivate(Er, El, nx, ny, localWhichPoln) private(i,j) shared(ap, a, p, L0) num_threads(Nth)
1868 : //#else
1869 : //#pragma omp parallel default(none) firstprivate(Er, El, nx, ny, localWhichPoln) private(i,j) shared(ap, a, p, L0) num_threads(Nth)
1870 : //#endif
1871 : {
1872 : //#pragma omp for
1873 0 : for(j = 0; j < ny; j++)
1874 : {
1875 0 : for(i = 0; i < nx; i++)
1876 : {
1877 0 : computePixelValues(ap, a, p, L0, Er, El, i,j,localWhichPoln);
1878 : }
1879 : }
1880 : }
1881 : // tim.show("BeamCalc:");
1882 :
1883 0 : deletePathology(p);
1884 0 : deleteAntenna(a);
1885 :
1886 0 : return 1;
1887 0 : }
1888 :
1889 0 : void BeamCalc::computePixelValues(const ApertureCalcParams *ap,
1890 : const calcAntenna *a, const Pathology *p,
1891 : const Double &L0,
1892 : Complex *Er, Complex *El,
1893 : const Int &i, const Int &j,
1894 : const Int& whichPoln)
1895 : {
1896 0 : Complex fp, Exr, Eyr, Exl, Eyl;
1897 : // Complex Er[3], El[3];
1898 0 : Complex E1[3];
1899 : Double dx, dy, x0, y0, x, y, r2, Rhole, Rant, R2, H2, eps;
1900 0 : Complex rr, rl, lr, ll, tmp;
1901 : Double L, amp, dP, dA, dO, x1, y1, dx1, dy1, dx2, dy2, phase;
1902 : //Int nx, ny;
1903 : Int os;
1904 0 : Int niter=6;
1905 : Double pac, pas, cp,sp; /* parallactic angle cosine / sine */
1906 0 : Complex Iota; Iota=Complex(0,1);
1907 0 : IPosition pos(4);pos=0;
1908 :
1909 0 : Ray *ray=0, *rayx=0, *rayy=0;
1910 : /* determine parallactic angle rotated coordinates */
1911 :
1912 0 : os = ap->oversamp;
1913 : //nx = ap->nx*os;
1914 : //ny = ap->ny*os;
1915 0 : dx = ap->dx/os;
1916 0 : dy = ap->dy/os;
1917 0 : x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
1918 0 : y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
1919 0 : Rant = a->radius;
1920 0 : Rhole = a->hole_radius;
1921 0 : R2 = Rant*Rant;
1922 0 : H2 = Rhole*Rhole;
1923 : // for(Int i=0; i < nx; ++i)
1924 : {
1925 0 : eps = dx/4.0;
1926 0 : pac = cos(ap->pa+M_PI/2);
1927 0 : pas = sin(ap->pa+M_PI/2);
1928 :
1929 0 : x = pac*(x0 + i*dx) - pas*(y0 + j*dy);
1930 0 : y = pas*(x0 + i*dx) + pac*(y0 + j*dy);
1931 0 : x = -x;
1932 :
1933 0 : if(fabs(x) > Rant) goto nextpoint;
1934 0 : if(fabs(y) > Rant) goto nextpoint;
1935 0 : r2 = x*x + y*y;
1936 0 : if(r2 > R2) goto nextpoint;
1937 0 : if(r2 < H2) goto nextpoint;
1938 :
1939 0 : ray = rayx = rayy = 0;
1940 :
1941 0 : x1 = x;
1942 0 : y1 = y;
1943 :
1944 0 : for(Int iter = 0; iter < niter; iter++)
1945 : {
1946 0 : ray = trace(a, x1, y1, p);
1947 0 : if(!ray) goto nextpoint;
1948 0 : x1 += (x - ray->aper[0]);
1949 0 : y1 += (y - ray->aper[1]);
1950 0 : deleteRay(ray);
1951 0 : ray = 0;
1952 : }
1953 :
1954 0 : ray = trace(a, x1, y1, p);
1955 :
1956 : /* check for leg blockage */
1957 0 : if(legplanewaveblock2(a, ray))
1958 0 : goto nextpoint;
1959 0 : if(legsphericalwaveblock(a, ray))
1960 0 : goto nextpoint;
1961 :
1962 0 : if(y < 0) rayy = trace(a, x1, y1+eps, p);
1963 0 : else rayy = trace(a, x1, y1-eps, p);
1964 :
1965 0 : if(x < 0) rayx = trace(a, x1+eps, y1, p);
1966 0 : else rayx = trace(a, x1-eps, y1, p);
1967 :
1968 0 : if(ray == 0 || rayx == 0 || rayy == 0)
1969 0 : goto nextpoint;
1970 :
1971 : /* compute solid angle subtended at the feed */
1972 0 : dx1 = rayx->aper[0]-ray->aper[0];
1973 0 : dy1 = rayx->aper[1]-ray->aper[1];
1974 0 : dx2 = rayy->aper[0]-ray->aper[0];
1975 0 : dy2 = rayy->aper[1]-ray->aper[1];
1976 :
1977 0 : dA = 0.5*fabs(dx1*dy2 - dx2*dy1);
1978 0 : dO = (dOmega(a, rayx, rayy, ray, p)/dA)*dx*dx;
1979 0 : dP = dO*feedgain(a, ray, p);
1980 0 : amp = sqrt(dP);
1981 :
1982 0 : L = Raylen(ray);
1983 :
1984 0 : phase = 2.0*M_PI*(L-L0)/a->lambda;
1985 :
1986 : /* phase retard the wave */
1987 0 : cp = cos(phase);
1988 0 : sp = sin(phase);
1989 : // fp = cp + sp*1.0i;
1990 :
1991 0 : fp = Complex(cp,sp);
1992 :
1993 :
1994 0 : tracepol(Er, ray, E1);
1995 0 : Exr = fp*amp*E1[0];
1996 0 : Eyr = fp*amp*E1[1];
1997 : // rr = Exr - 1.0i*Eyr;
1998 : // rl = Exr + 1.0i*Eyr;
1999 0 : rr = Exr - Iota*Eyr;
2000 0 : rl = Exr + Iota*Eyr;
2001 :
2002 0 : tracepol(El, ray, E1);
2003 0 : Exl = fp*amp*E1[0];
2004 0 : Eyl = fp*amp*E1[1];
2005 : // lr = Exl - 1.0i*Eyl;
2006 : // ll = Exl + 1.0i*Eyl;
2007 0 : lr = Exl - Iota*Eyl;
2008 0 : ll = Exl + Iota*Eyl;
2009 :
2010 : // pos(0)=(Int)((j/os) - (25.0/dy/os)/2 + shape(0)/2 - 0.5);
2011 : // pos(1)=(Int)((i/os) - (25.0/dx/os)/2 + shape(1)/2 - 0.5);
2012 : // Following 3 lines go with ANT tag in VLACalc.....
2013 : // Double antRadius=BeamCalcGeometryDefaults[ap->band].Rant;
2014 : // pos(0)=(Int)((j/os) - (antRadius/dy/os) + shape(0)/2 - 0.5);
2015 : // pos(1)=(Int)((i/os) - (antRadius/dx/os) + shape(1)/2 - 0.5);
2016 : // Following 2 lines go with the PIX tag in VLACalc...
2017 0 : pos(0)=(Int)((j/os));
2018 0 : pos(1)=(Int)((i/os));
2019 0 : pos(2)=0;
2020 0 : pos(3)=0;
2021 :
2022 0 : if (whichPoln==Stokes::RR)
2023 0 : {tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+rr,pos);}
2024 0 : else if (whichPoln==Stokes::RL)
2025 0 : {tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+rl,pos);}
2026 0 : else if (whichPoln==Stokes::LR)
2027 0 : {tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+lr,pos);}
2028 0 : else if (whichPoln==Stokes::LL)
2029 0 : {tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+ll,pos);}
2030 : else {
2031 0 : SynthesisError err(String("BeamCalc::computePixelValues: Cannot handle Stokes ")+String(whichPoln));
2032 0 : throw(err);
2033 0 : }
2034 :
2035 0 : nextpoint:
2036 0 : if(ray) deleteRay(ray);
2037 0 : if(rayx) deleteRay(rayx);
2038 0 : if(rayy) deleteRay(rayy);
2039 : }
2040 0 : }
2041 :
2042 : //
2043 : //----------------------------------------------------------------------
2044 : // Compute only the required polarizations.for linear polarization
2045 : //
2046 0 : Int BeamCalc::calculateApertureLinPol(ApertureCalcParams *ap, const Int& whichPoln)
2047 : {
2048 0 : Complex Ex[3], Ey[3];
2049 0 : Complex Px[2], Py[2];
2050 : //Double dx, dy, Rhole, Rant;//x0, y0, R2, H2, eps;
2051 : Double L0;
2052 : calcAntenna *a;
2053 : Pathology *p;
2054 : Int nx, ny, os;
2055 : Int i, j;
2056 : //Double pac, pas; /* parallactic angle cosine / sine */
2057 : // Complex Iota=Complex(0,1);
2058 :
2059 :
2060 0 : a = newAntennafromApertureCalcParams(ap);
2061 0 : p = newPathologyfromApertureCalcParams(ap);
2062 :
2063 : /* compute central ray pathlength */
2064 : {
2065 : Ray *tmpRay;
2066 0 : tmpRay = trace(a, 0.0, 0.00001, p);
2067 0 : L0 = Raylen(tmpRay);
2068 0 : deleteRay(tmpRay);
2069 : }
2070 :
2071 : //pac = cos(ap->pa+M_PI/2);
2072 : //pas = sin(ap->pa+M_PI/2);
2073 :
2074 : /* in linear basis */
2075 0 : Px[0] = 0.0; Px[1] = 1.0;
2076 0 : Py[0] = 1.0; Py[1] = 0.0;
2077 :
2078 0 : IPosition pos(4); pos=0;
2079 :
2080 : /* compute 3-vector feed efields for the two polarizations */
2081 0 : Efield(a, Py, Ey);
2082 0 : Efield(a, Px, Ex);
2083 :
2084 0 : if (whichPoln == Stokes::XX){
2085 0 : pos(2)=0;
2086 : }
2087 0 : else if (whichPoln == Stokes::YY){
2088 0 : pos(2)=3;
2089 : }
2090 0 : else if (whichPoln == Stokes::XY){
2091 0 : pos(2)=1;
2092 : }
2093 0 : else if (whichPoln == Stokes::YX){
2094 0 : pos(2)=2;
2095 : }
2096 : else {
2097 0 : SynthesisError err(String("BeamCalc::calculateApertureLinPol: Cannot handle Stokes ")+String(whichPoln));
2098 0 : throw(err);
2099 0 : }
2100 :
2101 : // set only the affected plane to zero
2102 0 : for(j = 0; j < ap->nx; j++){
2103 0 : pos(0)= j;
2104 0 : for(i = 0; i < ap->ny; i++){
2105 0 : pos(1)= i;
2106 0 : ap->aperture->putAt(Complex(0.),pos);
2107 : }
2108 : }
2109 :
2110 0 : os = ap->oversamp;
2111 0 : nx = ap->nx*os;
2112 0 : ny = ap->ny*os;
2113 : // dx = ap->dx/os;
2114 : // dy = ap->dy/os;
2115 : //x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
2116 : //y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
2117 : // Rant = a->radius;
2118 : // Rhole = a->hole_radius;
2119 : //R2 = Rant*Rant;
2120 : //H2 = Rhole*Rhole;
2121 :
2122 : //eps = dx/4.0;
2123 :
2124 : // cerr << "max threads " << omp_get_max_threads()
2125 : // << " threads available " << omp_get_num_threads()
2126 : // << endl;
2127 0 : Int Nth=1, localWhichPoln=whichPoln;
2128 : #ifdef _OPENMP
2129 0 : Nth=max(omp_get_max_threads()-2,1);
2130 : #endif
2131 : // Timer tim;
2132 : // tim.mark();
2133 : #if GCC44x
2134 : #pragma omp parallel default(none) firstprivate(Ex, Ey, nx, ny, localWhichPoln) private(i,j) shared(ap, a, p, L0) num_threads(Nth)
2135 : #else
2136 0 : #pragma omp parallel default(none) firstprivate(Ex, Ey, nx, ny, localWhichPoln) private(i,j) shared(ap, a, p, L0) num_threads(Nth)
2137 : #endif
2138 : {
2139 : #pragma omp for
2140 : for(j = 0; j < ny; j++)
2141 : {
2142 : for(i = 0; i < nx; i++)
2143 : {
2144 : computePixelValuesLinPol(ap, a, p, L0, Ex, Ey, i,j,localWhichPoln);
2145 : }
2146 : }
2147 : }
2148 : // tim.show("BeamCalc:");
2149 :
2150 0 : deletePathology(p);
2151 0 : deleteAntenna(a);
2152 :
2153 0 : return 1;
2154 0 : }
2155 :
2156 0 : void BeamCalc::computePixelValuesLinPol(const ApertureCalcParams *ap,
2157 : const calcAntenna *a, const Pathology *p,
2158 : const Double &L0,
2159 : Complex *Ex, Complex *Ey,
2160 : const Int &i, const Int &j,
2161 : const Int& whichPoln)
2162 : {
2163 0 : Complex fp, Exx, Eyx, Exy, Eyy;
2164 :
2165 0 : Complex E1[3];
2166 : Double dx, dy, x0, y0, x, y, r2, Rhole, Rant, R2, H2, eps;
2167 0 : Complex xx, xy, yx, yy, tmp;
2168 : Double L, amp, dP, dA, dO, x1, y1, dx1, dy1, dx2, dy2, phase;
2169 : //Int nx, ny;
2170 : Int os;
2171 0 : Int niter=6;
2172 : Double pac, pas, cp,sp; /* parallactic angle cosine / sine */
2173 0 : Complex Iota; Iota=Complex(0,1);
2174 0 : IPosition pos(4);pos=0;
2175 :
2176 0 : Ray *ray=0, *rayx=0, *rayy=0;
2177 : /* determine parallactic angle rotated coordinates */
2178 :
2179 0 : os = ap->oversamp;
2180 : //nx = ap->nx*os;
2181 : //ny = ap->ny*os;
2182 0 : dx = ap->dx/os;
2183 0 : dy = ap->dy/os;
2184 0 : x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
2185 0 : y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
2186 0 : Rant = a->radius;
2187 0 : Rhole = a->hole_radius;
2188 0 : R2 = Rant*Rant;
2189 0 : H2 = Rhole*Rhole;
2190 : // for(Int i=0; i < nx; ++i)
2191 : {
2192 0 : eps = dx/4.0;
2193 0 : pac = cos(ap->pa+M_PI/2);
2194 0 : pas = sin(ap->pa+M_PI/2);
2195 :
2196 0 : x = pac*(x0 + i*dx) - pas*(y0 + j*dy);
2197 0 : y = pas*(x0 + i*dx) + pac*(y0 + j*dy);
2198 0 : x = -x;
2199 :
2200 0 : if(fabs(x) > Rant) goto nextpoint;
2201 0 : if(fabs(y) > Rant) goto nextpoint;
2202 0 : r2 = x*x + y*y;
2203 0 : if(r2 > R2) goto nextpoint;
2204 0 : if(r2 < H2) goto nextpoint;
2205 :
2206 0 : ray = rayx = rayy = 0;
2207 :
2208 0 : x1 = x;
2209 0 : y1 = y;
2210 :
2211 0 : for(Int iter = 0; iter < niter; iter++)
2212 : {
2213 0 : ray = trace(a, x1, y1, p);
2214 0 : if(!ray) goto nextpoint;
2215 0 : x1 += (x - ray->aper[0]);
2216 0 : y1 += (y - ray->aper[1]);
2217 0 : deleteRay(ray);
2218 0 : ray = 0;
2219 : }
2220 :
2221 0 : ray = trace(a, x1, y1, p);
2222 :
2223 : /* check for leg blockage */
2224 0 : if(legplanewaveblock2(a, ray))
2225 0 : goto nextpoint;
2226 0 : if(legsphericalwaveblock(a, ray))
2227 0 : goto nextpoint;
2228 :
2229 0 : if(y < 0) rayy = trace(a, x1, y1+eps, p);
2230 0 : else rayy = trace(a, x1, y1-eps, p);
2231 :
2232 0 : if(x < 0) rayx = trace(a, x1+eps, y1, p);
2233 0 : else rayx = trace(a, x1-eps, y1, p);
2234 :
2235 0 : if(ray == 0 || rayx == 0 || rayy == 0)
2236 0 : goto nextpoint;
2237 :
2238 : /* compute solid angle subtended at the feed */
2239 0 : dx1 = rayx->aper[0]-ray->aper[0];
2240 0 : dy1 = rayx->aper[1]-ray->aper[1];
2241 0 : dx2 = rayy->aper[0]-ray->aper[0];
2242 0 : dy2 = rayy->aper[1]-ray->aper[1];
2243 :
2244 0 : dA = 0.5*fabs(dx1*dy2 - dx2*dy1);
2245 0 : dO = (dOmega(a, rayx, rayy, ray, p)/dA)*dx*dx;
2246 0 : dP = dO*feedgain(a, ray, p);
2247 0 : amp = sqrt(dP);
2248 :
2249 0 : L = Raylen(ray);
2250 :
2251 0 : phase = 2.0*M_PI*(L-L0)/a->lambda;
2252 :
2253 : /* phase retard the wave */
2254 0 : cp = cos(phase);
2255 0 : sp = sin(phase);
2256 :
2257 0 : fp = Complex(cp,sp);
2258 :
2259 :
2260 0 : tracepol(Ex, ray, E1);
2261 0 : Exx = fp*amp*E1[0];
2262 0 : Eyx = fp*amp*E1[1];
2263 :
2264 0 : tracepol(Ey, ray, E1);
2265 0 : Exy = fp*amp*E1[0];
2266 0 : Eyy = fp*amp*E1[1];
2267 :
2268 :
2269 0 : xx = Exx;
2270 0 : xy = Complex(0.);
2271 0 : yx = Complex(0.);
2272 0 : yy = Eyy;
2273 :
2274 : // pos(0)=(Int)((j/os) - (25.0/dy/os)/2 + shape(0)/2 - 0.5);
2275 : // pos(1)=(Int)((i/os) - (25.0/dx/os)/2 + shape(1)/2 - 0.5);
2276 : // Following 3 lines go with ANT tag in VLACalc.....
2277 : // Double antRadius=BeamCalcGeometryDefaults[ap->band].Rant;
2278 : // pos(0)=(Int)((j/os) - (antRadius/dy/os) + shape(0)/2 - 0.5);
2279 : // pos(1)=(Int)((i/os) - (antRadius/dx/os) + shape(1)/2 - 0.5);
2280 : // Following 2 lines go with the PIX tag in VLACalc...
2281 0 : pos(0)=(Int)((j/os));
2282 0 : pos(1)=(Int)((i/os));
2283 0 : pos(3)=0;
2284 :
2285 0 : if (whichPoln==Stokes::XX){
2286 0 : pos(2)=0;
2287 0 : tmp=ap->aperture->getAt(pos);
2288 0 : ap->aperture->putAt(tmp+xx,pos);
2289 : }
2290 :
2291 0 : else if (whichPoln==Stokes::XY){
2292 0 : pos(2)=1;
2293 0 : tmp=ap->aperture->getAt(pos);
2294 0 : ap->aperture->putAt(tmp+xy,pos);
2295 : }
2296 :
2297 0 : else if (whichPoln==Stokes::YX){
2298 0 : pos(2)=2;
2299 0 : tmp=ap->aperture->getAt(pos);
2300 0 : ap->aperture->putAt(tmp+yx,pos);
2301 : }
2302 :
2303 0 : else if (whichPoln==Stokes::YY){
2304 0 : pos(2)=3;
2305 0 : tmp=ap->aperture->getAt(pos);
2306 0 : ap->aperture->putAt(tmp+yy,pos);
2307 : }
2308 :
2309 0 : nextpoint:
2310 0 : if(ray) deleteRay(ray);
2311 0 : if(rayx) deleteRay(rayx);
2312 0 : if(rayy) deleteRay(rayy);
2313 : }
2314 0 : }
2315 : };
2316 :
|
