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: aips2-request@nrao.edu.
21 : //# Postal address: AIPS++ Project Office
22 : //# National Radio Astronomy Observatory
23 : //# 520 Edgemont Road
24 : //# Charlottesville, VA 22903-2475 USA
25 : //#
26 : //#
27 : //# $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():
71 : obsName_p(""),
72 : antType_p(""),
73 : obsTime_p(),
74 : BeamCalc_NumBandCodes_p(0),
75 : BeamCalcGeometries_p(0),
76 : bandMinFreq_p(0),
77 : bandMaxFreq_p(0),
78 1 : antRespPath_p(""){
79 1 : }
80 :
81 308 : BeamCalc* BeamCalc::Instance(){
82 308 : if(instance_p==0){
83 1 : instance_p = new BeamCalc();
84 : }
85 308 : return instance_p;
86 : }
87 :
88 : // initialise the beam calculation parameters
89 172 : void BeamCalc::setBeamCalcGeometries(const String& obsName,
90 : const String& antType,
91 : const MEpoch& obsTime,
92 : const String& otherAntRayPath){
93 :
94 344 : Unit uS("s");
95 172 : Bool verbose = false;
96 :
97 :
98 516 : if(obsName==obsName_p
99 171 : && antType==antType_p
100 343 : && obsTime.get(uS).getValue()==obsTime_p.get(uS).getValue()
101 343 : && otherAntRayPath.empty()
102 : ){
103 171 : return; // nothing to do (assuming the databases haven't changed)
104 : }
105 :
106 1 : cout << "Processing request for geometries from observatory " << obsName << ", antenna type " << antType << endl;
107 :
108 2 : LogIO os;
109 1 : os << LogOrigin("BeamCalc", "setBeamCalcGeometries()");
110 :
111 1 : if(obsName!=""){
112 1 : obsName_p = obsName;
113 : }
114 1 : if(antType!=""){
115 1 : antType_p = antType;
116 : }
117 1 : obsTime_p = obsTime;
118 :
119 :
120 1 : BeamCalcGeometries_p.resize(0);
121 :
122 2 : AntennaResponses aR;
123 2 : String antRespPath;
124 2 : String antRayPath = otherAntRayPath;
125 :
126 1 : Bool useInternal = false;
127 :
128 1 : os << LogIO::NORMAL << "Initialisation of geometries for observatory " << obsName_p
129 1 : << ", antenna type " << antType_p << LogIO::POST;
130 :
131 1 : if(otherAntRayPath.empty()){
132 1 : if(!MeasTable::AntennaResponsesPath(antRespPath, obsName_p)) {
133 1 : 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 : }
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 : }
166 :
167 1 : if(useInternal){
168 :
169 1 : Bool found = False;
170 2 : String fullFileName;
171 2 : const std::list<std::string> &data_path = AppStateSource::fetch( ).dataPath( );
172 1 : const std::string distrodata_path =casatools::get_state().distroDataPath( );
173 : //cerr<<"distrodata_path="<<distrodata_path<<endl;
174 : //cerr<<"DATA PATH==="<< *data_path <<endl;
175 : // The data path search need to be rewritten to adopt the recommanded setting via python
176 : // file for CASA6.
177 : // For now, only the first path that actually exist will be set to the data path (TT 2018.12.10)
178 1 : if (data_path.size() > 0 ) {
179 2 : for ( std::list<std::string>::const_iterator it=data_path.begin(); ! found && it != data_path.end(); ++it ) {
180 2 : Path lpath = Path(*it);
181 : //os<<"Here to datapath="<<lpath<<LogIO::POST;
182 : //Path lpath = Path(data_path);
183 2 : String slpath = lpath.absoluteName();
184 2 : String subdirname;
185 2 : if(obsName_p=="VLA" || obsName_p=="EVLA") {
186 2 : subdirname="/nrao/VLA";
187 : }
188 0 : else if(obsName_p=="ALMA"){
189 0 : subdirname="/alma/response";
190 : }
191 : //Directory ddir(slpath+subdirname);
192 : try {
193 4 : Directory ddir(slpath+subdirname);
194 1 : ddir.exists();
195 1 : found = True;
196 1 : fullFileName=slpath;
197 1 : break;
198 : }
199 1 : catch (...) {
200 : }
201 :
202 : //if (ddir.exists()) {
203 : // cerr<<" ddir exists:"<<slpath<<subdirname<<endl;
204 : // found = True;
205 : // fullFileName=slpath;
206 : // break;
207 : //}
208 : }
209 1 : if (!found && distrodata_path!="") {
210 0 : fullFileName = distrodata_path;
211 0 : found = True;
212 : }
213 : }
214 0 : else if(!found) {
215 0 : const char *sep=" ";
216 0 : char *aipsPath = strtok(getenv("CASAPATH"),sep);
217 0 : if (aipsPath == NULL)
218 0 : throw(SynthesisError("CASAPATH not found."));
219 0 : fullFileName=aipsPath;
220 0 : fullFileName+="/data";
221 : }
222 :
223 :
224 1 : if(obsName_p=="VLA" && antType_p=="STANDARD"){
225 0 : os << LogIO::NORMAL << "Will use default geometries for VLA STANDARD." << LogIO::POST;
226 0 : BeamCalc_NumBandCodes_p = VLABeamCalc_NumBandCodes;
227 0 : BeamCalcGeometries_p.resize(BeamCalc_NumBandCodes_p);
228 0 : bandMinFreq_p.resize(BeamCalc_NumBandCodes_p);
229 0 : bandMaxFreq_p.resize(BeamCalc_NumBandCodes_p);
230 0 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
231 0 : copyBeamCalcGeometry(&BeamCalcGeometries_p[i], &VLABeamCalcGeometryDefaults[i]);
232 0 : bandMinFreq_p[i] = VLABandMinFreqDefaults[i];
233 0 : bandMaxFreq_p[i] = VLABandMaxFreqDefaults[i];
234 : }
235 : //antRespPath_p = fullFileName + "/data/nrao/VLA";
236 0 : antRespPath_p = fullFileName + "/nrao/VLA";
237 : }
238 1 : else if(obsName_p=="EVLA" && antType_p=="STANDARD"){
239 1 : os << LogIO::NORMAL << "Will use default geometries for EVLA STANDARD." << LogIO::POST;
240 1 : BeamCalc_NumBandCodes_p = EVLABeamCalc_NumBandCodes;
241 1 : BeamCalcGeometries_p.resize(BeamCalc_NumBandCodes_p);
242 1 : bandMinFreq_p.resize(BeamCalc_NumBandCodes_p);
243 1 : bandMaxFreq_p.resize(BeamCalc_NumBandCodes_p);
244 10 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
245 9 : copyBeamCalcGeometry(&BeamCalcGeometries_p[i], &EVLABeamCalcGeometryDefaults[i]);
246 9 : bandMinFreq_p[i] = EVLABandMinFreqDefaults[i];
247 9 : bandMaxFreq_p[i] = EVLABandMaxFreqDefaults[i];
248 : }
249 : //antRespPath_p = fullFileName + "/data/nrao/VLA";
250 1 : antRespPath_p = fullFileName + "/nrao/VLA";
251 : }
252 0 : else if(obsName_p=="ALMA" && (antType_p=="DA" || antType_p=="DV" || antType_p=="PM")){
253 0 : os << LogIO::NORMAL << "Will use default geometries for ALMA DA, DV, and PM." << LogIO::POST;
254 0 : BeamCalc_NumBandCodes_p = ALMABeamCalc_NumBandCodes;
255 0 : BeamCalcGeometries_p.resize(BeamCalc_NumBandCodes_p);
256 0 : bandMinFreq_p.resize(BeamCalc_NumBandCodes_p);
257 0 : bandMaxFreq_p.resize(BeamCalc_NumBandCodes_p);
258 0 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
259 0 : copyBeamCalcGeometry(&BeamCalcGeometries_p[i], &ALMABeamCalcGeometryDefaults[i]);
260 0 : if(antType_p=="DA"){
261 0 : BeamCalcGeometries_p[i].legwidth *= -1.; // change from + to x shape
262 : }
263 0 : bandMinFreq_p[i] = ALMABandMinFreqDefaults[i];
264 0 : bandMaxFreq_p[i] = ALMABandMaxFreqDefaults[i];
265 : }
266 : //antRespPath_p = fullFileName + "/data/alma/responses";
267 0 : antRespPath_p = fullFileName + "/alma/responses";
268 : }
269 : else{
270 : String mesg="We don't have any antenna ray tracing parameters available for observatory "
271 0 : +obsName_p+", antenna type "+antType_p;
272 0 : SynthesisError err(mesg);
273 0 : throw(err);
274 : }
275 1 : return;
276 : } // end if(useInternal)
277 : }
278 :
279 :
280 0 : os << LogIO::NORMAL << "from file " << antRayPath << endl;
281 :
282 : try {
283 : // read temp table from ASCII file
284 0 : TableProxy antParTab = TableProxy(antRayPath, String(""), String("tempRayTraceTab.tab"),
285 0 : false, IPosition(), // autoheader, autoshape
286 0 : String(" "), // separator
287 0 : String("#"), // comment marker
288 : 0,-1, // first and last line
289 0 : Vector<String>(), Vector<String>());
290 :
291 0 : antParTab.deleteTable(true); // table will be deleted when it goes out of scope
292 :
293 : // read the table
294 0 : uInt nRows = antParTab.nrows();
295 0 : BeamCalc_NumBandCodes_p = nRows;
296 :
297 0 : BeamCalcGeometries_p.resize(BeamCalc_NumBandCodes_p);
298 0 : bandMinFreq_p.resize(BeamCalc_NumBandCodes_p);
299 0 : bandMaxFreq_p.resize(BeamCalc_NumBandCodes_p);
300 :
301 0 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
302 0 : sprintf(BeamCalcGeometries_p[i].name, "%s", antParTab.getCell("NAME", i).asString().c_str());
303 0 : bandMinFreq_p[i] = antParTab.getCell("MINFREQ", i).asDouble() * 1E9; // expect GHz
304 0 : bandMaxFreq_p[i] = antParTab.getCell("MAXFREQ", i).asDouble() * 1E9;
305 0 : BeamCalcGeometries_p[i].sub_h = antParTab.getCell("SUB_H", i).asDouble();
306 0 : Array<Double> ta1;
307 0 : ta1.assign(antParTab.getCell("FEEDPOS", i).asArrayDouble());
308 0 : for(uInt j=0; j<3;j++){
309 0 : BeamCalcGeometries_p[i].feedpos[j] = ta1(IPosition(1,j));
310 : }
311 0 : BeamCalcGeometries_p[i].subangle = antParTab.getCell("SUBANGLE", i).asDouble();
312 0 : BeamCalcGeometries_p[i].legwidth = antParTab.getCell("LEGWIDTH", i).asDouble();
313 0 : BeamCalcGeometries_p[i].legfoot = antParTab.getCell("LEGFOOT", i).asDouble();
314 0 : BeamCalcGeometries_p[i].legapex = antParTab.getCell("LEGAPEX", i).asDouble();
315 0 : BeamCalcGeometries_p[i].Rhole = antParTab.getCell("RHOLE", i).asDouble();
316 0 : BeamCalcGeometries_p[i].Rant = antParTab.getCell("RANT", i).asDouble();
317 0 : BeamCalcGeometries_p[i].reffreq = antParTab.getCell("REFFREQ", i).asDouble(); // stay in GHz
318 0 : Array<Double> ta2;
319 0 : ta2.assign(antParTab.getCell("TAPERPOLY", i).asArrayDouble());
320 0 : for(uInt j=0; j<5;j++){
321 0 : BeamCalcGeometries_p[i].taperpoly[j] = ta2(IPosition(1,j));
322 : }
323 0 : BeamCalcGeometries_p[i].ntaperpoly = antParTab.getCell("NTAPERPOLY", i).asInt();
324 0 : BeamCalcGeometries_p[i].astigm_0 = antParTab.getCell("ASTIGM_0", i).asDouble();
325 0 : BeamCalcGeometries_p[i].astigm_45 = antParTab.getCell("ASTIGM_45", i).asDouble();
326 0 : if(verbose){
327 : cout << "i name bandMinFreq_p bandMaxFreq_p sub_h feedpos feedpos feedpos subangle legwidth legfoot legapex"
328 0 : << " Rhole Rant reffreq taperpoly taperpoly taperpoly taperpoly taperpoly ntaperpoly astigm0 astigm45" << endl;
329 0 : cout << i << " " << BeamCalcGeometries_p[i].name << " " << bandMinFreq_p[i] << " " << bandMaxFreq_p[i]
330 0 : << " " << BeamCalcGeometries_p[i].sub_h
331 0 : << " " << BeamCalcGeometries_p[i].feedpos[0] << " " << BeamCalcGeometries_p[i].feedpos[1]
332 0 : << " " << BeamCalcGeometries_p[i].feedpos[2]
333 0 : << " " << BeamCalcGeometries_p[i].subangle << " " << BeamCalcGeometries_p[i].legwidth
334 0 : << " " << BeamCalcGeometries_p[i].legfoot << " " << BeamCalcGeometries_p[i].legapex
335 0 : << " " << BeamCalcGeometries_p[i].Rhole << " " << BeamCalcGeometries_p[i].Rant << " " << BeamCalcGeometries_p[i].reffreq
336 0 : << " " << BeamCalcGeometries_p[i].taperpoly[0] << " " << BeamCalcGeometries_p[i].taperpoly[1]
337 0 : << " " << BeamCalcGeometries_p[i].taperpoly[2] << " " << BeamCalcGeometries_p[i].taperpoly[3]
338 0 : << " " << BeamCalcGeometries_p[i].taperpoly[4] << " " << BeamCalcGeometries_p[i].ntaperpoly
339 0 : << " " << BeamCalcGeometries_p[i].astigm_0 << " " << BeamCalcGeometries_p[i].astigm_45 << endl;
340 : }
341 : }
342 :
343 0 : } catch (AipsError x) {
344 0 : String mesg="Initialisation of antenna ray tracing parameters for observatory "+obsName_p
345 0 : +" failed using path "+antRayPath+"\n with message "+x.getMesg();
346 0 : BeamCalcGeometries_p.resize(0);
347 0 : SynthesisError err(mesg);
348 0 : throw(err);
349 : }
350 :
351 0 : if(antRespPath.empty()){ // use containing directory of the antRayPath
352 0 : antRespPath_p = Path(antRayPath).dirName();
353 : }
354 : else{
355 0 : antRespPath_p = Path(antRespPath).dirName();
356 : }
357 :
358 0 : os << LogIO::NORMAL << "... successful." << LogIO::POST;
359 :
360 0 : return;
361 :
362 : }
363 :
364 172 : Int BeamCalc::getBandID(Double freq, // in Hz
365 : const String& obsName,
366 : const String& bandName,
367 : const String& antType,
368 : const MEpoch& obsTime,
369 : const String& otherAntRayPath){
370 :
371 172 : setBeamCalcGeometries(obsName, antType, obsTime, otherAntRayPath);
372 :
373 : // Check against bandName if it is a non-NULL string. Otherwise
374 : // check against frequency range. The latter method is only for
375 : // backward compatibility reasons and for older MSes which did not
376 : // have correct band names in the SPW sub-table.
377 172 : if (bandName != "")
378 1720 : for(uInt i=0; i<BeamCalcGeometries_p.nelements(); i++)
379 1548 : if(String(BeamCalcGeometries_p[i].bandName)==bandName) return i;
380 :
381 : // If the control flow gets here, the given bandName did not match
382 : // in SPW names in the MS. So use the fall-back option of using
383 : // the reference frequency to get the band ID (this will lead to
384 : // incorrect ID for the edge SPWs which might overlap in frequecy
385 : // with the adjecent band).
386 248 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
387 248 : if((bandMinFreq_p[i]<=freq)&&(freq<=bandMaxFreq_p[i])){
388 172 : return i;
389 : }
390 : }
391 0 : ostringstream mesg;
392 0 : mesg << obsName << "/" << bandName << "/" << antType << "/" << freq << "(Hz) combination not recognized.";
393 0 : throw(SynthesisError(mesg.str()));
394 :
395 : }
396 :
397 :
398 :
399 136 : calcAntenna* BeamCalc::newAntenna(Double sub_h, Double feed_x, Double feed_y, Double feed_z,
400 : Double ftaper, Double thmax, const char *geomfile)
401 : {
402 : calcAntenna *a;
403 : Int i;
404 : Double d, r, m, z;
405 : FILE *in;
406 272 : String fullFileName(antRespPath_p);
407 136 : fullFileName = fullFileName + String("/") + geomfile;
408 :
409 136 : in = fopen(fullFileName.c_str(), "r");
410 :
411 136 : if(!in)
412 : {
413 0 : String msg = "File " + fullFileName
414 0 : + " not found.\n Did you forget to install package data repository?\n";
415 0 : throw(SynthesisError(msg));
416 : }
417 :
418 136 : a = (calcAntenna *)malloc(sizeof(calcAntenna));
419 :
420 170272 : for(i = 0; i < MAXGEOM; i++)
421 : {
422 170272 : if(fscanf(in, "%lf%lf%lf", &r, &z, &m) != 3) break;
423 170136 : a->z[i] = z;
424 170136 : a->m[i] = m;
425 170136 : a->radius = r;
426 : }
427 136 : fclose(in);
428 136 : a->ngeom = i;
429 136 : a->zedge = z;
430 136 : a->deltar = a->radius/(float)(a->ngeom-1.0);
431 136 : a->bestparabola = a->zedge/(a->radius*a->radius);
432 136 : if(i < 3)
433 : {
434 0 : fprintf(stderr, "geom file not valid\n");
435 0 : free(a);
436 0 : return 0;
437 : }
438 :
439 136 : z = sub_h-feed_z;
440 :
441 136 : a->sub_h = sub_h;
442 136 : a->feed[0] = feed_x;
443 136 : a->feed[1] = feed_y;
444 136 : a->feed[2] = feed_z;
445 136 : d = std::sqrt((double)(feed_x*feed_x + feed_y*feed_y + z*z));
446 136 : if(z > 0.0)
447 : {
448 136 : a->K = sub_h + d;
449 136 : a->feeddir[0] = -feed_x/d;
450 136 : a->feeddir[1] = -feed_y/d;
451 136 : a->feeddir[2] = (sub_h-feed_z)/d;
452 : }
453 : else
454 : {
455 0 : a->K = std::sqrt((double(feed_x*feed_x + feed_y*feed_y + feed_z*feed_z)));
456 0 : a->feeddir[0] = -feed_x/d;
457 0 : a->feeddir[1] = -feed_y/d;
458 0 : a->feeddir[2] = (sub_h-feed_z)/d;
459 : }
460 544 : for(i = 0; i < 3; i++) a->pfeeddir[i] = a->feeddir[i];
461 136 : a->ftaper = fabs(ftaper);
462 136 : a->thmax = thmax;
463 136 : a->fa2pi = 2.0*M_PI*std::sqrt((double)ftaper)*0.1874/sin(thmax*M_PI/180.0);
464 136 : a->legwidth = 0.0;
465 136 : a->legfoot = a->radius/2.0;
466 136 : a->legapex = sub_h*1.2;
467 136 : a->legthick = 0.0;
468 136 : a->hole_radius = 0.0;
469 136 : a->dir[0] = a->dir[1] = 0.0;
470 136 : a->dir[2] = 1.0;
471 136 : strcpy(a->name, "unnamed");
472 136 : a->k[0] = a->k[1] = a->k[2] = 0.0;
473 : /* default to no polarization state */
474 136 : Antennasetfreq(a, 1.0);
475 136 : Antennasetdir(a, 0); /* compute hhat and vhat */
476 136 : a->gridsize = 0;
477 136 : dishvalue(a, a->legfoot, &a->legfootz, 0);
478 :
479 136 : return a;
480 : }
481 :
482 136 : void BeamCalc::deleteAntenna(calcAntenna *a)
483 : {
484 136 : if(!a) return;
485 :
486 136 : free(a);
487 : }
488 :
489 272 : void BeamCalc::Antennasetfreq(calcAntenna *a, Double freq)
490 : {
491 : Int i;
492 :
493 272 : a->freq = freq;
494 272 : a->lambda = NS_METER/freq;
495 1088 : for(i = 0; i < 3; i++) a->k[i] = -2.0*M_PI*a->dir[i]/a->lambda;
496 272 : }
497 :
498 272 : void BeamCalc::Antennasetdir(calcAntenna *a, const Double *dir)
499 : {
500 : Double hmag;
501 : Int i;
502 :
503 272 : if(dir)
504 : {
505 544 : for(i = 0; i < 3; i++) a->dir[i] = dir[i];
506 136 : if(a->dir[0] == 0.0 && a->dir[1] == 0.0)
507 : {
508 136 : a->hhat[0] = 1.0;
509 136 : a->hhat[1] = a->hhat[2] = 0.0;
510 136 : a->vhat[1] = 1.0;
511 136 : a->vhat[0] = a->vhat[2] = 0.0;
512 : }
513 : else
514 : {
515 0 : a->hhat[0] = a->dir[1];
516 0 : a->hhat[1] = -a->dir[0];
517 0 : a->hhat[2] = 0.0;
518 0 : hmag = sqrt(a->hhat[0]*a->hhat[0]
519 0 : + a->hhat[1]*a->hhat[1]);
520 0 : a->hhat[0] /= hmag;
521 0 : a->hhat[1] /= hmag;
522 :
523 0 : a->vhat[0] = a->hhat[1]*a->dir[2]
524 0 : - a->hhat[2]*a->dir[1];
525 0 : a->vhat[1] = a->hhat[2]*a->dir[0]
526 0 : - a->hhat[0]*a->dir[2];
527 0 : a->vhat[2] = a->hhat[0]*a->dir[1]
528 0 : - a->hhat[1]*a->dir[0];
529 : }
530 : }
531 1088 : for(i = 0; i < 3; i++) a->k[i] = -2.0*M_PI*a->dir[i]/a->lambda;
532 272 : }
533 :
534 : /* sets feeddir after pathology is considered */
535 0 : void BeamCalc::alignfeed(calcAntenna *a, const Pathology *p)
536 : {
537 : Int i, j;
538 0 : Double f[3], s0[3], s[3], d[3], m=0.0;
539 :
540 0 : for(i = 0; i < 3; i++) f[i] = a->feed[i] + p->feedshift[i];
541 0 : for(i = 0; i < 3; i++) s0[i] = -p->subrotpoint[i];
542 0 : s0[2] += a->sub_h;
543 0 : for(i = 0; i < 3; i++)
544 : {
545 0 : s[i] = 0.0;
546 0 : for(j = 0; j < 3; j++)
547 0 : s[i] += p->subrot[i][j]*s0[j];
548 0 : s[i] += p->subrotpoint[i] + p->subshift[i];
549 0 : d[i] = s[i]-f[i];
550 0 : m += d[i]*d[i];
551 : }
552 0 : m = sqrt(m);
553 0 : for(i = 0; i < 3; i++) a->feeddir[i] = d[i]/m;
554 0 : }
555 :
556 408 : void BeamCalc::getfeedbasis(const calcAntenna *a, Double B[3][3])
557 : {
558 : Int i;
559 : Double *dir, *vhat, *hhat;
560 :
561 408 : hhat = B[0];
562 408 : vhat = B[1];
563 408 : dir = B[2];
564 :
565 1632 : for(i = 0; i < 3; i++) dir[i] = a->pfeeddir[i];
566 :
567 408 : if(dir[0] == 0.0 && dir[1] == 0.0)
568 : {
569 0 : vhat[0] = 1.0;
570 0 : vhat[1] = vhat[2] = 0.0;
571 0 : hhat[1] = 1.0;
572 0 : hhat[0] = hhat[2] = 0.0;
573 : }
574 : else
575 : {
576 408 : vhat[0] = dir[1];
577 408 : vhat[1] = -dir[0];
578 408 : vhat[2] = 0.0;
579 408 : norm3(vhat);
580 :
581 408 : hhat[0] = vhat[1]*dir[2] - vhat[2]*dir[1];
582 408 : hhat[1] = vhat[2]*dir[0] - vhat[0]*dir[2];
583 408 : hhat[2] = vhat[0]*dir[1] - vhat[1]*dir[0];
584 : }
585 408 : }
586 :
587 272 : void BeamCalc::Efield(const calcAntenna *a, const Complex *pol, Complex *E)
588 : {
589 : Double B[3][3];
590 : Double *hhat, *vhat;
591 :
592 272 : getfeedbasis(a, B);
593 272 : hhat = B[0];
594 272 : vhat = B[1];
595 :
596 1088 : for(Int i = 0; i < 3; i++)
597 816 : E[i] = Complex(hhat[i],0) * pol[0] + Complex(vhat[i],0) * pol[1];
598 272 : }
599 :
600 0 : Int BeamCalc::Antennasetfeedpattern(calcAntenna* /*a*/,
601 : const char* /*filename*/,
602 : Double /*scale*/)
603 : {
604 : #if 0
605 : Int i, N, Nmax;
606 : Double x, delta;
607 : VecArray pat;
608 :
609 : a->feedpatterndelta = 0.0;
610 : if(a->feedpattern) deleteVector(a->feedpattern);
611 :
612 : if(filename == 0) return 1;
613 :
614 : pat = VecArrayfromfile(filename, 2);
615 :
616 : if(!pat) return 0;
617 : N = VectorSize(pat[0]);
618 : g_assert(N > 2);
619 : g_assert(pat[0][0] == 0.0);
620 :
621 : delta = pat[0][1];
622 : g_assert(delta > 0.0);
623 : for(i = 2; i < N; i++)
624 : {
625 : x = pat[0][i]-pat[0][i-1]-delta;
626 : g_assert(fabs(x) < delta/10000.0);
627 : }
628 :
629 : /* convert to radians */
630 : delta *= M_PI/180.0;
631 :
632 : /* and scale it */
633 : if(scale > 0.0) delta *= scale;
634 :
635 : /* Do we need to truncate the pattern? */
636 : Nmax = M_PI/delta;
637 : if(N > Nmax)
638 : {
639 : a->feedpattern = newVector(Nmax);
640 : for(i = 0; i < Nmax; i++)
641 : a->feedpattern[i] = fabs(pat[1][i]);
642 : deleteVector(pat[1]);
643 : }
644 : else a->feedpattern = pat[1];
645 :
646 : a->feedpatterndelta = delta;
647 : deleteVector(pat[0]);
648 : deleteVecArray(pat);
649 : #endif
650 0 : return 1;
651 : }
652 :
653 136 : calcAntenna* BeamCalc::newAntennafromApertureCalcParams(ApertureCalcParams *ap)
654 : {
655 : calcAntenna *a;
656 136 : Double dir[3] = {0.0, 0.0, 1.0};
657 : Double sub_h, feed_x, feed_y, feed_z, thmax, ftaper;
658 : char geomfile[128];//, *feedfile;
659 : BeamCalcGeometry *geom;
660 : Int i;
661 : Double x, freq, df;
662 :
663 : //LogIO os;
664 136 : if((0<=ap->band) && (ap->band<(Int)BeamCalcGeometries_p.size())){
665 136 : geom = &(BeamCalcGeometries_p[ap->band]);
666 : //os << "Using antenna parameters for " << geom->bandName << " band" << LogIO::POST;
667 : }
668 : else{
669 0 : SynthesisError err(String("Internal Error: attempt to access beam geometry for non-existing band."));
670 0 : throw(err);
671 : }
672 :
673 136 : sub_h = geom->sub_h;
674 136 : feed_x = geom->feedpos[0]; feed_x = -feed_x;
675 136 : feed_y = geom->feedpos[1];
676 136 : feed_z = geom->feedpos[2];
677 : //feedfile = 0;
678 136 : thmax = geom->subangle;
679 :
680 136 : freq = ap->freq;
681 136 : if(freq <= 0.0) freq = geom->reffreq;
682 :
683 136 : df = freq-geom->reffreq;
684 136 : x = 1.0;
685 136 : ftaper = 0.0;
686 464 : for(i = 0; i < geom->ntaperpoly; i++)
687 : {
688 328 : ftaper += geom->taperpoly[i]*x;
689 328 : x *= df;
690 : }
691 136 : sprintf(geomfile, "%s.surface", geom->name);
692 :
693 136 : a = newAntenna(sub_h, feed_x, feed_y, feed_z, ftaper, thmax, geomfile);
694 136 : if(!a) return 0;
695 :
696 136 : 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 136 : Antennasetfreq(a, ap->freq);
716 :
717 136 : a->legwidth = geom->legwidth;
718 136 : a->legfoot = geom->legfoot;
719 136 : a->legapex = geom->legapex;
720 :
721 136 : a->hole_radius = geom->Rhole;
722 :
723 136 : a->astigm_0 = geom->astigm_0;
724 136 : a->astigm_45 = geom->astigm_45;
725 :
726 136 : Antennasetdir(a, dir);
727 :
728 136 : return a;
729 : }
730 :
731 233292880 : 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 233292880 : Double s = 1.0;
736 : Int n;
737 :
738 233292880 : if(r == 0)
739 : {
740 0 : *z = a->z[0];
741 0 : *m = 0.0;
742 0 : return 1;
743 : }
744 :
745 233292880 : if(r < 0)
746 : {
747 0 : s = -1.0;
748 0 : r = -r;
749 : }
750 233292880 : d = a->deltar;
751 233292880 : dd = d*d;
752 :
753 233292880 : n = (Int)floor(r/d + 0.5); /* the middle point */
754 233292880 : if(n > a->ngeom-2) n = a->ngeom-2;
755 :
756 233292880 : x = r - n*d;
757 :
758 233292880 : if(n == 0)
759 : {
760 1224 : mc = a->m[1];
761 1224 : ma = -mc;
762 1224 : mb = 0.0;
763 1224 : zav = 2.0*a->z[1] + a->z[0];
764 : }
765 : else
766 : {
767 233291656 : ma = a->m[n-1];
768 233291656 : mb = a->m[n];
769 233291656 : mc = a->m[n+1];
770 233291656 : zav = a->z[n-1] + a->z[n] + a->z[n+1];
771 : }
772 :
773 233292880 : A = mb;
774 233292880 : B = 0.5*(mc - ma)/d;
775 233292880 : C = 0.5*(mc - 2.0*mb + ma)/dd;
776 :
777 233292880 : D = (zav - B*dd)/3.0;
778 :
779 233292880 : if(m) *m = s*(A + B*x + C*x*x);
780 233292880 : if(z) *z = s*(D + A*x + B*x*x/2.0 + C*x*x*x/3.0);
781 :
782 233292880 : return 1;
783 : }
784 :
785 233292744 : 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 233292744 : Double s = 1.0;
790 : Int n;
791 :
792 233292744 : rr = x*x + y*y;
793 233292744 : r = sqrt(rr);
794 :
795 233292744 : if(r==0. || (a->astigm_0==0. && a->astigm_45==0.))
796 : {
797 233292744 : 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 25921416 : 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 25921416 : r = sqrt(x*x + y*y);
863 :
864 25921416 : 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 25921416 : astigdishvalue(a, x, y, &z, &m);
873 :
874 : /* Compute reflected unit vector direction */
875 25921416 : m = tan(2.0*atan(m));
876 25921416 : w = 1.0/sqrt(1.0+m*m);
877 25921416 : u = -m*(x/r)*w;
878 25921416 : v = -m*(y/r)*w;
879 :
880 25921416 : dx = a->feed[0]-x;
881 25921416 : dy = a->feed[1]-y;
882 25921416 : dz = a->feed[2]-z;
883 25921416 : dl = a->K + z;
884 :
885 51842832 : t = 0.5*(dx*dx + dy*dy + dz*dz - dl*dl)
886 25921416 : / (-dl + u*dx + v*dy + w*dz);
887 :
888 25921416 : subpoint[0] = x + u*t;
889 25921416 : subpoint[1] = y + v*t;
890 25921416 : subpoint[2] = z + w*t;
891 : }
892 :
893 103685664 : for(i = 0; i < 3; i++) sf[i] = a->feed[i] - subpoint[i];
894 25921416 : sd[0] = x - subpoint[0];
895 25921416 : sd[1] = y - subpoint[1];
896 25921416 : sd[2] = z - subpoint[2];
897 :
898 25921416 : norm3(sf);
899 25921416 : norm3(sd);
900 :
901 103685664 : for(i = 0; i < 3; i++) n[i] = sd[i]+sf[i];
902 :
903 25921416 : norm3(n);
904 :
905 103685664 : for(i = 0; i < 3; i++) subpoint[3+i] = n[i];
906 :
907 25921416 : 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 25921416 : Ray * BeamCalc::newRay(const Double *sub)
981 : {
982 : Ray *ray;
983 : Int i;
984 :
985 25921416 : ray = (Ray *)malloc(sizeof(Ray));
986 181449912 : for(i = 0; i < 6; i++) ray->sub[i] = sub[i];
987 :
988 25921416 : return ray;
989 : }
990 :
991 25921416 : void BeamCalc::deleteRay(Ray *ray)
992 : {
993 25921416 : if(ray) free(ray);
994 25921416 : }
995 :
996 136 : Pathology* BeamCalc::newPathology()
997 : {
998 : Pathology *P;
999 : Int i, j;
1000 :
1001 136 : P = (Pathology *)malloc(sizeof(Pathology));
1002 :
1003 544 : for(i = 0; i < 3; i++) P->subrotpoint[i] = P->subshift[i] = P->feedshift[i] = 0.0;
1004 1768 : for(i = 0; i < 3; i++) for(j = 0; j < 3; j++)
1005 1224 : P->feedrot[i][j] = P->subrot[i][j] = 0.0;
1006 544 : for(i = 0; i < 3; i++) P->feedrot[i][i] = P->subrot[i][i] = 1.0;
1007 :
1008 136 : P->az_offset = 0.0;
1009 136 : P->el_offset = 0.0;
1010 136 : P->phase_offset = 0.0;
1011 136 : P->focus = 0.0;
1012 :
1013 136 : return P;
1014 : }
1015 :
1016 136 : Pathology* BeamCalc::newPathologyfromApertureCalcParams(ApertureCalcParams* /*ap*/)
1017 : {
1018 : Pathology *P;
1019 :
1020 136 : P = newPathology();
1021 :
1022 136 : return P;
1023 : }
1024 :
1025 136 : void BeamCalc::deletePathology(Pathology *P)
1026 : {
1027 136 : if(P) free(P);
1028 136 : }
1029 :
1030 8262576 : void BeamCalc::normvec(const Double *a, const Double *b, Double *c)
1031 : {
1032 : Int i;
1033 : Double r;
1034 33050304 : for(i = 0; i < 3; i++) c[i] = a[i] - b[i];
1035 8262576 : r = sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
1036 33050304 : for(i = 0; i < 3; i++) c[i] /= r;
1037 8262576 : }
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 2754192 : 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 11016768 : for(i = 0; i < 3; i++) f[i] = a->feed[i] + p->feedshift[i];
1080 :
1081 2754192 : normvec(ray1->sub, f, n1);
1082 2754192 : normvec(ray2->sub, f, n2);
1083 2754192 : normvec(ray3->sub, f, n3);
1084 :
1085 11016768 : for(i = 0; i < 3; i++)
1086 : {
1087 8262576 : n1[i] -= n3[i];
1088 8262576 : n2[i] -= n3[i];
1089 : }
1090 :
1091 2754192 : ci = n1[1]*n2[2] - n1[2]*n2[1];
1092 2754192 : cj = n1[2]*n2[0] - n1[0]*n2[2];
1093 2754192 : ck = n1[0]*n2[1] - n1[1]*n2[0];
1094 :
1095 2754192 : Omega = 0.5*sqrt(ci*ci + cj*cj + ck*ck);
1096 :
1097 2754192 : return Omega;
1098 : }
1099 :
1100 2754328 : Double BeamCalc::Raylen(const Ray *ray)
1101 : {
1102 : Double len, d[3];
1103 : Int i;
1104 :
1105 : /* feed to subreflector */
1106 11017312 : for(i = 0; i < 3; i++) d[i] = ray->feed[i] - ray->sub[i];
1107 2754328 : len = sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
1108 :
1109 : /* subreflector to dish */
1110 11017312 : for(i = 0; i < 3; i++) d[i] = ray->sub[i] - ray->dish[i];
1111 2754328 : len += sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
1112 :
1113 : /* dish to aperture */
1114 11017312 : for(i = 0; i < 3; i++) d[i] = ray->dish[i] - ray->aper[i];
1115 2754328 : len += sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
1116 :
1117 2754328 : return len;
1118 : }
1119 :
1120 25921416 : void BeamCalc::Pathologize(Double *sub, const Pathology *p)
1121 : {
1122 : Int i;
1123 : Int j;
1124 : Double tmp[6];
1125 :
1126 103685664 : for(i = 0; i < 3; i++) sub[i] -= p->subrotpoint[i];
1127 :
1128 103685664 : for(i = 0; i < 3; i++)
1129 : {
1130 77764248 : tmp[i] = 0.0;
1131 77764248 : tmp[i+3] = 0.0;
1132 311056992 : for(j = 0; j < 3; j++) tmp[i] += p->subrot[i][j]*sub[j];
1133 311056992 : for(j = 0; j < 3; j++) tmp[i+3] += p->subrot[i][j]*sub[j+3];
1134 : }
1135 :
1136 103685664 : for(i = 0; i < 3; i++)
1137 77764248 : sub[i] = tmp[i] + p->subrotpoint[i] + p->subshift[i];
1138 77764248 : for(i = 4; i < 6; i++)
1139 51842832 : sub[i] = tmp[i];
1140 25921416 : }
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 25921416 : 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 25921416 : A = a->bestparabola*(unitdir[0]*unitdir[0] + unitdir[1]*unitdir[1]);
1173 51842832 : B = 2.0*a->bestparabola*(unitdir[0]*sub[0] + unitdir[1]*sub[1])
1174 25921416 : -unitdir[2];
1175 25921416 : C = a->bestparabola*(sub[0]*sub[0] + sub[1]*sub[1]) - sub[2];
1176 25921416 : t = 0.5*(sqrt(B*B-4.0*A*C) - B)/A; /* take greater root */
1177 :
1178 25921416 : for(iter = 0; ; iter++)
1179 : {
1180 : /* get position (x) and normal (n) on the real dish */
1181 622113984 : for(i = 0; i < 2; i++) x[i] = sub[i] + t*unitdir[i];
1182 207371328 : r = sqrt(x[0]*x[0] + x[1]*x[1]);
1183 207371328 : astigdishvalue(a, x[0], x[1], &(x[2]), &m);
1184 207371328 : n[2] = 1.0/sqrt(1.0+m*m);
1185 207371328 : n[0] = -m*(x[0]/r)*n[2];
1186 207371328 : n[1] = -m*(x[1]/r)*n[2];
1187 :
1188 207371328 : if(iter >= niter) break;
1189 :
1190 181449912 : A = B = 0;
1191 725799648 : for(i = 0; i < 3; i++)
1192 : {
1193 544349736 : A += n[i]*(x[i]-sub[i]); /* n dot (x-sub) */
1194 544349736 : B += n[i]*unitdir[i]; /* n dot unitdir */
1195 : }
1196 181449912 : t = A/B;
1197 : }
1198 :
1199 103685664 : for(i = 0; i < 3; i++)
1200 : {
1201 77764248 : dish[i] = x[i];
1202 77764248 : dish[i+3] = n[i];
1203 : }
1204 25921416 : }
1205 :
1206 25921416 : void BeamCalc::intersectaperture(const calcAntenna *a, const Double *dish,
1207 : const Double *unitdir, Double *aper)
1208 : {
1209 : Double t;
1210 : Int i;
1211 :
1212 25921416 : t = (a->zedge-dish[2])/unitdir[2];
1213 103685664 : for(i = 0; i < 3; i++) aper[i] = dish[i] + t*unitdir[i];
1214 :
1215 25921416 : aper[3] = aper[4] = 0.0;
1216 25921416 : aper[5] = 1.0;
1217 25921416 : }
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 2754192 : Double BeamCalc::feedgain(const calcAntenna *a, const Ray *ray, const Pathology */*p*/)
1230 : {
1231 2754192 : Double costheta = 0.0;
1232 : Double v[3];
1233 : Int i;
1234 :
1235 11016768 : for(i = 0; i < 3; i++) v[i] = ray->sub[i] - ray->feed[i];
1236 2754192 : norm3(v);
1237 :
1238 11016768 : for(i = 0; i < 3; i++)
1239 : {
1240 8262576 : costheta += a->pfeeddir[i]*v[i];
1241 : }
1242 :
1243 :
1244 2754192 : return exp(2.0*(-0.083)*a->fa2pi*a->fa2pi*(1.0-costheta*costheta));
1245 : }
1246 :
1247 25921416 : Ray* BeamCalc::trace(const calcAntenna *a, Double x, Double y, const Pathology *p)
1248 : {
1249 : Ray *ray;
1250 : Double idealsub[12];
1251 25921416 : Double fu[3], du[3], au[3], ndotf=0.0, ndotd=0.0;
1252 : Int i;
1253 25921416 : const Int niter = 7;
1254 :
1255 25921416 : subfromdish(a, x, y, idealsub);
1256 :
1257 25921416 : ray = newRay(idealsub);
1258 :
1259 25921416 : Pathologize(ray->sub, p);
1260 :
1261 25921416 : if(ray->sub[5] < -1.0 || ray->sub[5] > -0.0)
1262 : {
1263 0 : deleteRay(ray);
1264 0 : return 0;
1265 : }
1266 :
1267 103685664 : 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 103685664 : for(i = 0; i < 3; i++) fu[i] = ray->feed[i] - ray->sub[i];
1273 25921416 : norm3(fu);
1274 :
1275 : /* unit toward dish */
1276 103685664 : for(i = 0; i < 3; i++) ndotf += ray->sub[i+3]*fu[i];
1277 103685664 : for(i = 0; i < 3; i++) du[i] = 2.0*ray->sub[i+3]*ndotf - fu[i];
1278 :
1279 : /* dish point */
1280 25921416 : intersectdish(a, ray->sub, du, ray->dish, niter);
1281 :
1282 : /* unit toward aperture */
1283 103685664 : for(i = 0; i < 3; i++) ndotd += ray->dish[i+3]*du[i];
1284 103685664 : for(i = 0; i < 3; i++) au[i] = du[i] - 2.0*ray->dish[i+3]*ndotd;
1285 :
1286 25921416 : intersectaperture(a, ray->dish, au, ray->aper);
1287 :
1288 25921416 : return ray;
1289 : }
1290 :
1291 5508384 : void BeamCalc::tracepol(Complex *E0, const Ray *ray, Complex *E1)
1292 : {
1293 5508384 : Complex fac;
1294 : Double v1[3], v2[3], v3[3];
1295 : Double r[3];
1296 : Int i;
1297 :
1298 22033536 : for(i = 0; i < 3; i++)
1299 : {
1300 16525152 : v1[i] = ray->sub[i] - ray->feed[i];
1301 16525152 : v2[i] = ray->dish[i] - ray->sub[i];
1302 16525152 : v3[i] = ray->aper[i] - ray->dish[i];
1303 16525152 : E1[i] = E0[i];
1304 : }
1305 5508384 : norm3(v1);
1306 5508384 : norm3(v2);
1307 5508384 : norm3(v3);
1308 :
1309 22033536 : for(i = 0; i < 3; i++) r[i] = v1[i] - v2[i];
1310 5508384 : norm3(r);
1311 5508384 : fac = Complex(r[0],0)*E1[0] + Complex(r[1],0)*E1[1] + Complex(r[2],0)*E1[2];
1312 22033536 : for(i = 0; i < 3; i++) E1[i] = Complex(r[i],0)*fac*2.0 - E1[i];
1313 :
1314 22033536 : for(i = 0; i < 3; i++) r[i] = v2[i] - v3[i];
1315 5508384 : norm3(r);
1316 5508384 : fac = Complex(r[0],0)*E1[0] + Complex(r[1],0)*E1[1] + Complex(r[2],0)*E1[2];
1317 22033536 : for(i = 0; i < 3; i++) E1[i] = Complex(r[i],0)*fac*2.0 - E1[i];
1318 5508384 : }
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 2916128 : 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 2916128 : const Double thetaplus[4] =
1357 : {0, M_PI/2.0, M_PI, 3.0*M_PI/2.0};
1358 2916128 : const Double thetacross[4] =
1359 : {0.25*M_PI, 0.75*M_PI, 1.25*M_PI, 1.75*M_PI};
1360 2916128 : const Double thetavlba[4] =
1361 : {0.816817, 2.3247756, 3.9584096, 5.466368};
1362 : const Double *thetalut;
1363 :
1364 2916128 : if(a->legwidth == 0.0) return 0;
1365 :
1366 2916128 : if(strcmp(a->name, "VLBA") == 0) thetalut = thetavlba;
1367 2916128 : else if(a->legwidth < 0.0) thetalut = thetacross;
1368 2916128 : else thetalut = thetaplus;
1369 :
1370 : /* inside the leg feet is plane wave blockage */
1371 2916128 : dr2 = ray->dish[0]*ray->dish[0] + ray->dish[1]*ray->dish[1];
1372 2916128 : if(dr2 >= a->legfoot*a->legfoot) return 0;
1373 :
1374 3827072 : for(i = 0; i < 3; i++)
1375 : {
1376 2870304 : r0[i] = ray->dish[i];
1377 2870304 : dr[i] = ray->aper[i] - r0[i];
1378 : }
1379 956768 : rr = r0[0]*r0[0] + r0[1]*r0[1];
1380 :
1381 956768 : l0[2] = a->legfootz;
1382 956768 : l1[0] = l1[1] = 0.0;
1383 956768 : l1[2] = a->legapex;
1384 : // phi set but not used
1385 : // phi = atan2(r0[1], r0[0]);
1386 4703120 : for(n = 0; n < 4; n++)
1387 : {
1388 3778640 : theta = thetalut[n];
1389 3778640 : l0[0] = a->legfoot*cos(theta);
1390 3778640 : l0[1] = a->legfoot*sin(theta);
1391 3778640 : ll = l0[0]*l0[0] + l0[1]*l0[1];
1392 3778640 : if((l0[0]*r0[0] + l0[1]*r0[1])/sqrt(ll*rr) < 0.7) continue;
1393 3877312 : for(i = 0; i < 3; i++) dl[i] = l1[i] - l0[i];
1394 3877312 : for(i = 0; i < 3; i++) D[i] = r0[i] - l0[i];
1395 :
1396 969328 : N[0] = dr[1]*dl[2] - dr[2]*dl[1];
1397 969328 : N[1] = dr[2]*dl[0] - dr[0]*dl[2];
1398 969328 : N[2] = dr[0]*dl[1] - dr[1]*dl[0];
1399 969328 : norm3(N);
1400 :
1401 969328 : D2 = D[0]*N[0] + D[1]*N[1] + D[2]*N[2];
1402 :
1403 969328 : if(fabs(D2) <= 0.5*fabs(a->legwidth)) return 1;
1404 : }
1405 :
1406 924480 : return 0;
1407 : }
1408 :
1409 2883840 : 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 2883840 : const Double thetaplus[4] =
1418 : {0, M_PI/2.0, M_PI, 3.0*M_PI/2.0};
1419 2883840 : const Double thetacross[4] =
1420 : {0.25*M_PI, 0.75*M_PI, 1.25*M_PI, 1.75*M_PI};
1421 2883840 : const Double thetavlba[4] =
1422 : {0.816817, 2.3247756, 3.9584096, 5.466368};
1423 : const Double *thetalut;
1424 :
1425 2883840 : if(a->legwidth == 0.0) return 0;
1426 :
1427 2883840 : if(strcmp(a->name, "VLBA") == 0) thetalut = thetavlba;
1428 2883840 : else if(a->legwidth < 0.0) thetalut = thetacross;
1429 2883840 : else thetalut = thetaplus;
1430 :
1431 : /* inside the leg feet is plane wave blockage */
1432 2883840 : dr2 = ray->dish[0]*ray->dish[0] + ray->dish[1]*ray->dish[1];
1433 2883840 : if(dr2 < a->legfoot*a->legfoot) return 0;
1434 :
1435 7837440 : for(i = 0; i < 3; i++)
1436 : {
1437 5878080 : r0[i] = ray->dish[i];
1438 5878080 : dr[i] = ray->sub[i] - r0[i];
1439 : }
1440 1959360 : rr = r0[0]*r0[0] + r0[1]*r0[1];
1441 :
1442 1959360 : l0[2] = a->legfootz;
1443 1959360 : l1[0] = l1[1] = 0.0;
1444 1959360 : l1[2] = a->legapex;
1445 : // phi set but not used
1446 : // phi = atan2(r0[1], r0[0]);
1447 9472680 : for(n = 0; n < 4; n++)
1448 : {
1449 7642968 : theta = thetalut[n];
1450 7642968 : l0[0] = a->legfoot*cos(theta);
1451 7642968 : l0[1] = a->legfoot*sin(theta);
1452 7642968 : ll = l0[0]*l0[0] + l0[1]*l0[1];
1453 7642968 : if((l0[0]*r0[0] + l0[1]*r0[1])/sqrt(ll*rr) < 0.7) continue;
1454 7934528 : for(i = 0; i < 3; i++) dl[i] = l1[i] - l0[i];
1455 7934528 : for(i = 0; i < 3; i++) D[i] = r0[i] - l0[i];
1456 :
1457 1983632 : N[0] = dr[1]*dl[2] - dr[2]*dl[1];
1458 1983632 : N[1] = dr[2]*dl[0] - dr[0]*dl[2];
1459 1983632 : N[2] = dr[0]*dl[1] - dr[1]*dl[0];
1460 1983632 : norm3(N);
1461 :
1462 1983632 : D2 = D[0]*N[0] + D[1]*N[1] + D[2]*N[2];
1463 :
1464 1983632 : if(fabs(D2) <= 0.5*fabs(a->legwidth)) return 1;
1465 : }
1466 :
1467 1829712 : return 0;
1468 : }
1469 :
1470 :
1471 9 : void BeamCalc::copyBeamCalcGeometry(BeamCalcGeometry* to, BeamCalcGeometry* from){
1472 9 : sprintf(to->name, "%s", from->name);
1473 9 : sprintf(to->bandName, "%s", from->bandName);
1474 9 : to->sub_h = from->sub_h;
1475 36 : for(uInt j=0; j<3;j++){
1476 27 : to->feedpos[j] = from->feedpos[j];
1477 : }
1478 9 : to->subangle = from->subangle;
1479 9 : to->legwidth = from->legwidth;
1480 9 : to->legfoot = from->legfoot;
1481 9 : to->legapex = from->legapex;
1482 9 : to->Rhole = from->Rhole;
1483 9 : to->Rant = from->Rant;
1484 9 : to->reffreq = from->reffreq;
1485 54 : for(uInt j=0; j<5;j++){
1486 45 : to->taperpoly[j] = from->taperpoly[j];
1487 : }
1488 9 : to->ntaperpoly = from->ntaperpoly;
1489 9 : to->astigm_0 = from->astigm_0;
1490 9 : to->astigm_45 = from->astigm_45;
1491 :
1492 9 : }
1493 :
1494 :
1495 : /* The meat of the calculation */
1496 :
1497 :
1498 136 : Int BeamCalc::calculateAperture(ApertureCalcParams *ap)
1499 : {
1500 : // not used
1501 : // Complex fp, Exr, Eyr, Exl, Eyl;
1502 136 : Complex Er[3], El[3];
1503 136 : Complex Pr[2], Pl[2];
1504 136 : 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 136 : 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 136 : a = newAntennafromApertureCalcParams(ap);
1524 136 : p = newPathologyfromApertureCalcParams(ap);
1525 :
1526 : /* compute central ray pathlength */
1527 : {
1528 : Ray *tmpRay;
1529 136 : tmpRay = trace(a, 0.0, 0.00001, p);
1530 136 : L0 = Raylen(tmpRay);
1531 136 : deleteRay(tmpRay);
1532 : }
1533 :
1534 : //pac = cos(ap->pa+M_PI/2);
1535 : //pas = sin(ap->pa+M_PI/2);
1536 :
1537 136 : if(obsName_p=="EVLA" || obsName_p=="VLA"){
1538 : /* compute polarization vectors in circular basis */
1539 136 : Pr[0] = 1.0/M_SQRT2; Pr[1] = Iota/M_SQRT2;
1540 136 : Pl[0] = 1.0/M_SQRT2; Pl[1] = -Iota/M_SQRT2;
1541 :
1542 : /* compensate for feed orientation */
1543 136 : getfeedbasis(a, B);
1544 136 : phase = atan2(B[0][1], B[0][0]);
1545 136 : cp = cos(phase);
1546 136 : sp = sin(phase);
1547 :
1548 136 : q[0] = Pr[0];
1549 136 : q[1] = Pr[1];
1550 136 : Pr[0] = Complex(cp,0)*q[0] + Complex(sp,0)*q[1];
1551 136 : Pr[1] = -Complex(sp,0)*q[0] + Complex(cp,0)*q[1];
1552 136 : q[0] = Pl[0];
1553 136 : q[1] = Pl[1];
1554 136 : Pl[0] = Complex(cp,0)*q[0] + Complex(sp,0)*q[1];
1555 136 : 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 136 : Efield(a, Pr, Er);
1567 136 : Efield(a, Pl, El);
1568 136 : ap->aperture->set(Complex(0.0));
1569 :
1570 136 : os = ap->oversamp;
1571 136 : nx = ap->nx*os;
1572 136 : 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 136 : 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 448648 : for(j = 0; j < ny; j++)
1601 : {
1602 2161563648 : for(i = 0; i < nx; i++)
1603 : {
1604 2161115136 : computePixelValues(ap, a, p, L0, Er, El, i,j);
1605 : }
1606 : }
1607 : }
1608 : // tim.show("BeamCalc:");
1609 :
1610 136 : deletePathology(p);
1611 136 : deleteAntenna(a);
1612 :
1613 272 : return 1;
1614 : }
1615 :
1616 2161115136 : 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 2161115136 : Complex fp, Exr, Eyr, Exl, Eyl;
1623 : // Complex Er[3], El[3];
1624 2161115136 : Complex E1[3];
1625 : Double dx, dy, Rant, x0, y0, x, y, r2, R2, H2, eps, Rhole;
1626 2161115136 : 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 2161115136 : Int niter=6;
1631 : Double pac, pas;
1632 : Double cp,sp; /* parallactic angle cosine / sine */
1633 2161115136 : Complex Iota; Iota=Complex(0,1);
1634 4322230272 : IPosition pos(4);pos=0;
1635 :
1636 2161115136 : Ray *ray=0, *rayx=0, *rayy=0;
1637 : /* determine parallactic angle rotated coordinates */
1638 :
1639 2161115136 : os = ap->oversamp;
1640 : //nx = ap->nx*os;
1641 : //ny = ap->ny*os;
1642 2161115136 : dx = ap->dx/os;
1643 2161115136 : dy = ap->dy/os;
1644 2161115136 : x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
1645 2161115136 : y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
1646 2161115136 : Rant = a->radius;
1647 2161115136 : Rhole = a->hole_radius;
1648 2161115136 : R2 = Rant*Rant;
1649 2161115136 : H2 = Rhole*Rhole;
1650 : // for(Int i=0; i < nx; ++i)
1651 : {
1652 2161115136 : eps = dx/4.0;
1653 2161115136 : pac = cos(ap->pa+M_PI/2);
1654 2161115136 : pas = sin(ap->pa+M_PI/2);
1655 :
1656 2161115136 : x = pac*(x0 + i*dx) - pas*(y0 + j*dy);
1657 2161115136 : y = pas*(x0 + i*dx) + pac*(y0 + j*dy);
1658 2161115136 : x = -x;
1659 :
1660 2161115136 : if(fabs(x) > Rant) goto nextpoint;
1661 80735620 : if(fabs(y) > Rant) goto nextpoint;
1662 3927880 : r2 = x*x + y*y;
1663 3927880 : if(r2 > R2) goto nextpoint;
1664 3084808 : if(r2 < H2) goto nextpoint;
1665 :
1666 2916128 : ray = rayx = rayy = 0;
1667 :
1668 2916128 : x1 = x;
1669 2916128 : y1 = y;
1670 :
1671 20412896 : for(Int iter = 0; iter < niter; iter++)
1672 : {
1673 17496768 : ray = trace(a, x1, y1, p);
1674 17496768 : if(!ray) goto nextpoint;
1675 17496768 : x1 += (x - ray->aper[0]);
1676 17496768 : y1 += (y - ray->aper[1]);
1677 17496768 : deleteRay(ray);
1678 17496768 : ray = 0;
1679 : }
1680 :
1681 2916128 : ray = trace(a, x1, y1, p);
1682 :
1683 : /* check for leg blockage */
1684 2916128 : if(legplanewaveblock2(a, ray))
1685 32288 : goto nextpoint;
1686 2883840 : if(legsphericalwaveblock(a, ray))
1687 129648 : goto nextpoint;
1688 :
1689 2754192 : if(y < 0) rayy = trace(a, x1, y1+eps, p);
1690 1377096 : else rayy = trace(a, x1, y1-eps, p);
1691 :
1692 2754192 : if(x < 0) rayx = trace(a, x1+eps, y1, p);
1693 1377096 : else rayx = trace(a, x1-eps, y1, p);
1694 :
1695 2754192 : if(ray == 0 || rayx == 0 || rayy == 0)
1696 0 : goto nextpoint;
1697 :
1698 : /* compute solid angle subtended at the feed */
1699 2754192 : dx1 = rayx->aper[0]-ray->aper[0];
1700 2754192 : dy1 = rayx->aper[1]-ray->aper[1];
1701 2754192 : dx2 = rayy->aper[0]-ray->aper[0];
1702 2754192 : dy2 = rayy->aper[1]-ray->aper[1];
1703 :
1704 2754192 : dA = 0.5*fabs(dx1*dy2 - dx2*dy1);
1705 2754192 : dO = (dOmega(a, rayx, rayy, ray, p)/dA)*dx*dx;
1706 2754192 : dP = dO*feedgain(a, ray, p);
1707 2754192 : amp = sqrt(dP);
1708 :
1709 2754192 : L = Raylen(ray);
1710 :
1711 2754192 : phase = 2.0*M_PI*(L-L0)/a->lambda;
1712 :
1713 : /* phase retard the wave */
1714 2754192 : cp = cos(phase);
1715 2754192 : sp = sin(phase);
1716 : // fp = cp + sp*1.0i;
1717 :
1718 2754192 : fp = Complex(cp,sp);
1719 :
1720 :
1721 2754192 : tracepol(Er, ray, E1);
1722 2754192 : Exr = fp*amp*E1[0];
1723 2754192 : Eyr = fp*amp*E1[1];
1724 :
1725 : // rr = Exr - 1.0i*Eyr;
1726 : // rl = Exr + 1.0i*Eyr;
1727 2754192 : rr = Exr - Iota*Eyr;
1728 2754192 : rl = Exr + Iota*Eyr;
1729 :
1730 2754192 : tracepol(El, ray, E1);
1731 2754192 : Exl = fp*amp*E1[0];
1732 2754192 : Eyl = fp*amp*E1[1];
1733 : // lr = Exl - 1.0i*Eyl;
1734 : // ll = Exl + 1.0i*Eyl;
1735 2754192 : lr = Exl - Iota*Eyl;
1736 2754192 : 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 2754192 : pos(0)=(Int)((j/os));
1745 2754192 : pos(1)=(Int)((i/os));
1746 2754192 : pos(3)=0;
1747 :
1748 2754192 : pos(2)=0;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+rr,pos);
1749 2754192 : pos(2)=1;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+rl,pos);
1750 2754192 : pos(2)=2;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+lr,pos);
1751 2754192 : pos(2)=3;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+ll,pos);
1752 2161115136 : nextpoint:
1753 2161115136 : if(ray) deleteRay(ray);
1754 2161115136 : if(rayx) deleteRay(rayx);
1755 2161115136 : if(rayy) deleteRay(rayy);
1756 : }
1757 2161115136 : }
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 : }
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 : }
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 : }
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 0 : #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 : #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 : }
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 :
|
