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 0 : BeamCalc::BeamCalc():
71 0 : obsName_p(""),
72 0 : antType_p(""),
73 0 : obsTime_p(),
74 0 : BeamCalc_NumBandCodes_p(0),
75 0 : BeamCalcGeometries_p(0),
76 0 : bandMinFreq_p(0),
77 0 : bandMaxFreq_p(0),
78 0 : antRespPath_p(""){
79 0 : }
80 :
81 0 : BeamCalc* BeamCalc::Instance(){
82 0 : if(instance_p==0){
83 0 : instance_p = new BeamCalc();
84 : }
85 0 : return instance_p;
86 : }
87 :
88 : // initialise the beam calculation parameters
89 0 : void BeamCalc::setBeamCalcGeometries(const String& obsName,
90 : const String& antType,
91 : const MEpoch& obsTime,
92 : const String& otherAntRayPath){
93 :
94 0 : Unit uS("s");
95 0 : Bool verbose = false;
96 :
97 :
98 0 : if(obsName==obsName_p
99 0 : && antType==antType_p
100 0 : && obsTime.get(uS).getValue()==obsTime_p.get(uS).getValue()
101 0 : && otherAntRayPath.empty()
102 : ){
103 0 : return; // nothing to do (assuming the databases haven't changed)
104 : }
105 :
106 0 : cout << "Processing request for geometries from observatory " << obsName << ", antenna type " << antType << endl;
107 :
108 0 : LogIO os;
109 0 : os << LogOrigin("BeamCalc", "setBeamCalcGeometries()");
110 :
111 0 : if(obsName!=""){
112 0 : obsName_p = obsName;
113 : }
114 0 : if(antType!=""){
115 0 : antType_p = antType;
116 : }
117 0 : obsTime_p = obsTime;
118 :
119 :
120 0 : BeamCalcGeometries_p.resize(0);
121 :
122 0 : AntennaResponses aR;
123 0 : String antRespPath;
124 0 : String antRayPath = otherAntRayPath;
125 :
126 0 : Bool useInternal = false;
127 :
128 0 : os << LogIO::NORMAL << "Initialisation of geometries for observatory " << obsName_p
129 0 : << ", antenna type " << antType_p << LogIO::POST;
130 :
131 0 : if(otherAntRayPath.empty()){
132 0 : if(!MeasTable::AntennaResponsesPath(antRespPath, obsName_p)) {
133 0 : 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 0 : if(useInternal){
168 :
169 0 : Bool found = False;
170 0 : String fullFileName;
171 0 : const std::list<std::string> &data_path = AppStateSource::fetch( ).dataPath( );
172 0 : 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 0 : if (data_path.size() > 0 ) {
179 0 : for ( std::list<std::string>::const_iterator it=data_path.begin(); ! found && it != data_path.end(); ++it ) {
180 0 : Path lpath = Path(*it);
181 : //os<<"Here to datapath="<<lpath<<LogIO::POST;
182 : //Path lpath = Path(data_path);
183 0 : String slpath = lpath.absoluteName();
184 0 : String subdirname;
185 0 : if(obsName_p=="VLA" || obsName_p=="EVLA") {
186 0 : 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 0 : Directory ddir(slpath+subdirname);
194 0 : ddir.exists();
195 0 : found = True;
196 0 : fullFileName=slpath;
197 0 : break;
198 0 : }
199 0 : catch (...) {
200 0 : }
201 :
202 : //if (ddir.exists()) {
203 : // cerr<<" ddir exists:"<<slpath<<subdirname<<endl;
204 : // found = True;
205 : // fullFileName=slpath;
206 : // break;
207 : //}
208 0 : }
209 0 : 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 0 : 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 0 : else if(obsName_p=="EVLA" && antType_p=="STANDARD"){
239 0 : os << LogIO::NORMAL << "Will use default geometries for EVLA STANDARD." << LogIO::POST;
240 0 : BeamCalc_NumBandCodes_p = EVLABeamCalc_NumBandCodes;
241 0 : BeamCalcGeometries_p.resize(BeamCalc_NumBandCodes_p);
242 0 : bandMinFreq_p.resize(BeamCalc_NumBandCodes_p);
243 0 : bandMaxFreq_p.resize(BeamCalc_NumBandCodes_p);
244 0 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
245 0 : copyBeamCalcGeometry(&BeamCalcGeometries_p[i], &EVLABeamCalcGeometryDefaults[i]);
246 0 : bandMinFreq_p[i] = EVLABandMinFreqDefaults[i];
247 0 : bandMaxFreq_p[i] = EVLABandMaxFreqDefaults[i];
248 : }
249 : //antRespPath_p = fullFileName + "/data/nrao/VLA";
250 0 : 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 0 : }
275 0 : return;
276 0 : } // 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 0 : }
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 0 : }
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 0 : }
363 :
364 0 : 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 0 : 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 0 : if (bandName != "")
378 0 : for(uInt i=0; i<BeamCalcGeometries_p.nelements(); i++)
379 0 : 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 0 : for(uInt i=0; i<BeamCalc_NumBandCodes_p; i++){
387 0 : if((bandMinFreq_p[i]<=freq)&&(freq<=bandMaxFreq_p[i])){
388 0 : 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 0 : }
396 :
397 :
398 :
399 0 : 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 0 : String fullFileName(antRespPath_p);
407 0 : fullFileName = fullFileName + String("/") + geomfile;
408 :
409 0 : in = fopen(fullFileName.c_str(), "r");
410 :
411 0 : 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 0 : }
417 :
418 0 : a = (calcAntenna *)malloc(sizeof(calcAntenna));
419 :
420 0 : for(i = 0; i < MAXGEOM; i++)
421 : {
422 0 : if(fscanf(in, "%lf%lf%lf", &r, &z, &m) != 3) break;
423 0 : a->z[i] = z;
424 0 : a->m[i] = m;
425 0 : a->radius = r;
426 : }
427 0 : fclose(in);
428 0 : a->ngeom = i;
429 0 : a->zedge = z;
430 0 : a->deltar = a->radius/(float)(a->ngeom-1.0);
431 0 : a->bestparabola = a->zedge/(a->radius*a->radius);
432 0 : if(i < 3)
433 : {
434 0 : fprintf(stderr, "geom file not valid\n");
435 0 : free(a);
436 0 : return 0;
437 : }
438 :
439 0 : z = sub_h-feed_z;
440 :
441 0 : a->sub_h = sub_h;
442 0 : a->feed[0] = feed_x;
443 0 : a->feed[1] = feed_y;
444 0 : a->feed[2] = feed_z;
445 0 : d = std::sqrt((double)(feed_x*feed_x + feed_y*feed_y + z*z));
446 0 : if(z > 0.0)
447 : {
448 0 : a->K = sub_h + d;
449 0 : a->feeddir[0] = -feed_x/d;
450 0 : a->feeddir[1] = -feed_y/d;
451 0 : 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 0 : for(i = 0; i < 3; i++) a->pfeeddir[i] = a->feeddir[i];
461 0 : a->ftaper = fabs(ftaper);
462 0 : a->thmax = thmax;
463 0 : a->fa2pi = 2.0*M_PI*std::sqrt((double)ftaper)*0.1874/sin(thmax*M_PI/180.0);
464 0 : a->legwidth = 0.0;
465 0 : a->legfoot = a->radius/2.0;
466 0 : a->legapex = sub_h*1.2;
467 0 : a->legthick = 0.0;
468 0 : a->hole_radius = 0.0;
469 0 : a->dir[0] = a->dir[1] = 0.0;
470 0 : a->dir[2] = 1.0;
471 0 : strcpy(a->name, "unnamed");
472 0 : a->k[0] = a->k[1] = a->k[2] = 0.0;
473 : /* default to no polarization state */
474 0 : Antennasetfreq(a, 1.0);
475 0 : Antennasetdir(a, 0); /* compute hhat and vhat */
476 0 : a->gridsize = 0;
477 0 : dishvalue(a, a->legfoot, &a->legfootz, 0);
478 :
479 0 : return a;
480 0 : }
481 :
482 0 : void BeamCalc::deleteAntenna(calcAntenna *a)
483 : {
484 0 : if(!a) return;
485 :
486 0 : free(a);
487 : }
488 :
489 0 : void BeamCalc::Antennasetfreq(calcAntenna *a, Double freq)
490 : {
491 : Int i;
492 :
493 0 : a->freq = freq;
494 0 : a->lambda = NS_METER/freq;
495 0 : for(i = 0; i < 3; i++) a->k[i] = -2.0*M_PI*a->dir[i]/a->lambda;
496 0 : }
497 :
498 0 : void BeamCalc::Antennasetdir(calcAntenna *a, const Double *dir)
499 : {
500 : Double hmag;
501 : Int i;
502 :
503 0 : if(dir)
504 : {
505 0 : for(i = 0; i < 3; i++) a->dir[i] = dir[i];
506 0 : if(a->dir[0] == 0.0 && a->dir[1] == 0.0)
507 : {
508 0 : a->hhat[0] = 1.0;
509 0 : a->hhat[1] = a->hhat[2] = 0.0;
510 0 : a->vhat[1] = 1.0;
511 0 : 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 0 : for(i = 0; i < 3; i++) a->k[i] = -2.0*M_PI*a->dir[i]/a->lambda;
532 0 : }
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 0 : void BeamCalc::getfeedbasis(const calcAntenna *a, Double B[3][3])
557 : {
558 : Int i;
559 : Double *dir, *vhat, *hhat;
560 :
561 0 : hhat = B[0];
562 0 : vhat = B[1];
563 0 : dir = B[2];
564 :
565 0 : for(i = 0; i < 3; i++) dir[i] = a->pfeeddir[i];
566 :
567 0 : 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 0 : vhat[0] = dir[1];
577 0 : vhat[1] = -dir[0];
578 0 : vhat[2] = 0.0;
579 0 : norm3(vhat);
580 :
581 0 : hhat[0] = vhat[1]*dir[2] - vhat[2]*dir[1];
582 0 : hhat[1] = vhat[2]*dir[0] - vhat[0]*dir[2];
583 0 : hhat[2] = vhat[0]*dir[1] - vhat[1]*dir[0];
584 : }
585 0 : }
586 :
587 0 : void BeamCalc::Efield(const calcAntenna *a, const Complex *pol, Complex *E)
588 : {
589 : Double B[3][3];
590 : Double *hhat, *vhat;
591 :
592 0 : getfeedbasis(a, B);
593 0 : hhat = B[0];
594 0 : vhat = B[1];
595 :
596 0 : for(Int i = 0; i < 3; i++)
597 0 : E[i] = Complex(hhat[i],0) * pol[0] + Complex(vhat[i],0) * pol[1];
598 0 : }
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 0 : calcAntenna* BeamCalc::newAntennafromApertureCalcParams(ApertureCalcParams *ap)
654 : {
655 : calcAntenna *a;
656 0 : 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 0 : if((0<=ap->band) && (ap->band<(Int)BeamCalcGeometries_p.size())){
665 0 : 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 0 : }
672 :
673 0 : sub_h = geom->sub_h;
674 0 : feed_x = geom->feedpos[0]; feed_x = -feed_x;
675 0 : feed_y = geom->feedpos[1];
676 0 : feed_z = geom->feedpos[2];
677 : //feedfile = 0;
678 0 : thmax = geom->subangle;
679 :
680 0 : freq = ap->freq;
681 0 : if(freq <= 0.0) freq = geom->reffreq;
682 :
683 : //cerr << "BEam CAlc freq "<< freq << " reffreq " << geom->reffreq << endl;
684 :
685 :
686 0 : df = freq-geom->reffreq;
687 0 : x = 1.0;
688 0 : ftaper = 0.0;
689 0 : for(i = 0; i < geom->ntaperpoly; i++)
690 : {
691 0 : ftaper += geom->taperpoly[i]*x;
692 0 : x *= df;
693 : }
694 0 : sprintf(geomfile, "%s.surface", geom->name);
695 :
696 0 : a = newAntenna(sub_h, feed_x, feed_y, feed_z, ftaper, thmax, geomfile);
697 0 : if(!a) return 0;
698 :
699 0 : strcpy(a->name, geom->name);
700 :
701 : /* feed pattern file is two column text file containing
702 : * angle (in degrees) and power (in dBi)
703 : */
704 :
705 : // if(feedfile != 0)
706 : // {
707 : // Double scale;
708 : // scale = getKeyValueDouble(kv, "feedpatternscale");
709 : // if(!Antennasetfeedpattern(a, feedfile, scale))
710 : // {
711 : // deleteAntenna(a);
712 : // fprintf(stderr, "Problem with feed file <%s>\n",
713 : // feedfile);
714 : // return 0;
715 : // }
716 : // }
717 :
718 0 : Antennasetfreq(a, ap->freq);
719 :
720 0 : a->legwidth = geom->legwidth;
721 0 : a->legfoot = geom->legfoot;
722 0 : a->legapex = geom->legapex;
723 :
724 0 : a->hole_radius = geom->Rhole;
725 :
726 0 : a->astigm_0 = geom->astigm_0;
727 0 : a->astigm_45 = geom->astigm_45;
728 :
729 0 : Antennasetdir(a, dir);
730 :
731 0 : return a;
732 : }
733 :
734 0 : Int BeamCalc::dishvalue(const calcAntenna *a, Double r, Double *z, Double *m)
735 : {
736 : Double ma, mb, mc, zav, A, B, C, D;
737 : Double x, d, dd;
738 0 : Double s = 1.0;
739 : Int n;
740 :
741 0 : if(r == 0)
742 : {
743 0 : *z = a->z[0];
744 0 : *m = 0.0;
745 0 : return 1;
746 : }
747 :
748 0 : if(r < 0)
749 : {
750 0 : s = -1.0;
751 0 : r = -r;
752 : }
753 0 : d = a->deltar;
754 0 : dd = d*d;
755 :
756 0 : n = (Int)floor(r/d + 0.5); /* the middle point */
757 0 : if(n > a->ngeom-2) n = a->ngeom-2;
758 :
759 0 : x = r - n*d;
760 :
761 0 : if(n == 0)
762 : {
763 0 : mc = a->m[1];
764 0 : ma = -mc;
765 0 : mb = 0.0;
766 0 : zav = 2.0*a->z[1] + a->z[0];
767 : }
768 : else
769 : {
770 0 : ma = a->m[n-1];
771 0 : mb = a->m[n];
772 0 : mc = a->m[n+1];
773 0 : zav = a->z[n-1] + a->z[n] + a->z[n+1];
774 : }
775 :
776 0 : A = mb;
777 0 : B = 0.5*(mc - ma)/d;
778 0 : C = 0.5*(mc - 2.0*mb + ma)/dd;
779 :
780 0 : D = (zav - B*dd)/3.0;
781 :
782 0 : if(m) *m = s*(A + B*x + C*x*x);
783 0 : if(z) *z = s*(D + A*x + B*x*x/2.0 + C*x*x*x/3.0);
784 :
785 0 : return 1;
786 : }
787 :
788 0 : Int BeamCalc::astigdishvalue(const calcAntenna *a, Double x, Double y, Double *z, Double *m)
789 : {
790 : Double ma, mb, mc, zav, A, B, C, D;
791 : Double r, rr, theta, xp, d, dd, z5, z6, astigm, dastigm;
792 0 : Double s = 1.0;
793 : Int n;
794 :
795 0 : rr = x*x + y*y;
796 0 : r = sqrt(rr);
797 :
798 0 : if(r==0. || (a->astigm_0==0. && a->astigm_45==0.))
799 : {
800 0 : return dishvalue(a, r, z, m);
801 : }
802 :
803 : // the Zernike polynomials Z5 and Z6
804 : Double sin2th, cos2th, rho, rho2;
805 :
806 0 : theta = atan2(y,x);
807 0 : sin2th = sin(2.*theta);
808 0 : cos2th = cos(2.*theta);
809 0 : rho = r / a->radius;
810 0 : rho2 = rho*rho;
811 :
812 0 : z5 = sqrt(6.) * rho2 * sin2th;
813 0 : z6 = sqrt(6.) * rho2 * cos2th;
814 :
815 0 : astigm = 1. + a->astigm_45 * z5 + a->astigm_0 * z6;
816 0 : dastigm = 2.* rho2/r * sqrt(6.)*(a->astigm_45*sin2th + a->astigm_0*cos2th);
817 :
818 0 : d = a->deltar;
819 0 : dd = d*d;
820 :
821 0 : n = (Int)floor(r/d + 0.5); /* the middle point */
822 0 : if(n > a->ngeom-2) n = a->ngeom-2;
823 :
824 0 : xp = r - n*d;
825 :
826 0 : if(n == 0)
827 : {
828 0 : mc = a->m[1];
829 0 : ma = -mc;
830 0 : mb = 0.0;
831 0 : zav = 2.0*a->z[1] + a->z[0];
832 : }
833 : else
834 : {
835 0 : ma = a->m[n-1];
836 0 : mb = a->m[n];
837 0 : mc = a->m[n+1];
838 0 : zav = a->z[n-1] + a->z[n] + a->z[n+1];
839 : }
840 :
841 0 : A = mb;
842 0 : B = 0.5*(mc - ma)/d;
843 0 : C = 0.5*(mc - 2.0*mb + ma)/dd;
844 :
845 0 : D = (zav - B*dd)/3.0;
846 :
847 :
848 0 : Double zn = s*(D + A*xp + B*xp*xp/2.0 + C*xp*xp*xp/3.0);
849 0 : if(z) *z = zn * astigm;
850 0 : if(m) *m = s*(A + B*xp + C*xp*xp) * astigm + dastigm * zn;
851 :
852 0 : return 1;
853 : }
854 :
855 : /* Returns position of subreflector piece (x, y, z) and
856 : * its normal (u, v, w)
857 : */
858 0 : Int BeamCalc::subfromdish(const calcAntenna *a, Double x, Double y, Double *subpoint)
859 : {
860 : Double r, z, m, u, v, w;
861 : Double dx, dy, dz, dl, t;
862 : Double n[3], sf[3], sd[3];
863 : Int i;
864 :
865 0 : r = sqrt(x*x + y*y);
866 :
867 0 : if(r == 0)
868 : {
869 0 : subpoint[0] = 0.0;
870 0 : subpoint[1] = 0.0;
871 0 : subpoint[2] = a->sub_h;
872 : }
873 : else
874 : {
875 0 : astigdishvalue(a, x, y, &z, &m);
876 :
877 : /* Compute reflected unit vector direction */
878 0 : m = tan(2.0*atan(m));
879 0 : w = 1.0/sqrt(1.0+m*m);
880 0 : u = -m*(x/r)*w;
881 0 : v = -m*(y/r)*w;
882 :
883 0 : dx = a->feed[0]-x;
884 0 : dy = a->feed[1]-y;
885 0 : dz = a->feed[2]-z;
886 0 : dl = a->K + z;
887 :
888 0 : t = 0.5*(dx*dx + dy*dy + dz*dz - dl*dl)
889 0 : / (-dl + u*dx + v*dy + w*dz);
890 :
891 0 : subpoint[0] = x + u*t;
892 0 : subpoint[1] = y + v*t;
893 0 : subpoint[2] = z + w*t;
894 : }
895 :
896 0 : for(i = 0; i < 3; i++) sf[i] = a->feed[i] - subpoint[i];
897 0 : sd[0] = x - subpoint[0];
898 0 : sd[1] = y - subpoint[1];
899 0 : sd[2] = z - subpoint[2];
900 :
901 0 : norm3(sf);
902 0 : norm3(sd);
903 :
904 0 : for(i = 0; i < 3; i++) n[i] = sd[i]+sf[i];
905 :
906 0 : norm3(n);
907 :
908 0 : for(i = 0; i < 3; i++) subpoint[3+i] = n[i];
909 :
910 0 : return 1;
911 : }
912 :
913 0 : Int BeamCalc::dishfromsub(const calcAntenna *a, Double x, Double y, Double *dishpoint)
914 : {
915 :
916 : Double x1, y1, dx, dy, mx, my, r, d;
917 0 : const Double eps = 0.001;
918 0 : Int iter, niter=500;
919 : Double sub[5][6];
920 :
921 0 : x1 = x;
922 0 : y1 = y;
923 :
924 0 : for(iter = 0; iter < niter; iter++)
925 : {
926 0 : subfromdish(a, x1, y1, sub[0]);
927 0 : subfromdish(a, x1-eps, y1, sub[1]);
928 0 : subfromdish(a, x1+eps, y1, sub[2]);
929 0 : subfromdish(a, x1, y1-eps, sub[3]);
930 0 : subfromdish(a, x1, y1+eps, sub[4]);
931 0 : mx = 0.5*(sub[2][0]-sub[1][0])/eps;
932 0 : my = 0.5*(sub[4][1]-sub[3][1])/eps;
933 0 : dx = (x-sub[0][0])/mx;
934 0 : dy = (y-sub[0][1])/my;
935 0 : if(fabs(dx) > a->radius/7.0)
936 : {
937 0 : if(dx < 0) dx = -a->radius/7.0;
938 0 : else dx = a->radius/7.0;
939 : }
940 0 : if(fabs(dy) > a->radius/7.0)
941 : {
942 0 : if(dy < 0) dy = -a->radius/7.0;
943 0 : else dy = a->radius/7.0;
944 : }
945 0 : r = sqrt(x1*x1 + y1*y1);
946 0 : if(r >= a->radius)
947 0 : if(x1*dx + y1*dy > 0.0) return 0;
948 0 : x1 += 0.5*dx;
949 0 : y1 += 0.5*dy;
950 0 : if(fabs(dx) < 0.005*eps && fabs(dy) < 0.005*eps) break;
951 : }
952 0 : if(iter == niter) return 0;
953 :
954 0 : r = sqrt(x1*x1 + y1*y1);
955 0 : dishpoint[0] = x1;
956 0 : dishpoint[1] = y1;
957 : // dishpoint[2] = polyvalue(a->shape, r);
958 0 : dishpoint[3] = sub[0][0];
959 0 : dishpoint[4] = sub[0][1];
960 0 : dishpoint[5] = sub[0][2];
961 0 : d = sqrt(1.0+mx*mx+my*my);
962 0 : dishpoint[6] = mx/d;
963 0 : dishpoint[7] = my/d;
964 0 : dishpoint[8] = 1.0/d;
965 0 : dishpoint[9] = sub[0][3];
966 0 : dishpoint[10] = sub[0][4];
967 0 : dishpoint[11] = sub[0][5];
968 :
969 0 : if(r > a->radius) return 0;
970 0 : else return 1;
971 : }
972 :
973 0 : void BeamCalc::printAntenna(const calcAntenna *a)
974 : {
975 0 : printf("Antenna: %s %p\n", a->name, a);
976 0 : printf(" freq = %f GHz lambda = %f m\n", a->freq, a->lambda);
977 0 : printf(" feeddir = %f, %f, %f\n",
978 0 : a->feeddir[0], a->feeddir[1], a->feeddir[2]);
979 0 : printf(" pfeeddir = %f, %f, %f\n",
980 0 : a->pfeeddir[0], a->pfeeddir[1], a->pfeeddir[2]);
981 0 : }
982 :
983 0 : Ray * BeamCalc::newRay(const Double *sub)
984 : {
985 : Ray *ray;
986 : Int i;
987 :
988 0 : ray = (Ray *)malloc(sizeof(Ray));
989 0 : for(i = 0; i < 6; i++) ray->sub[i] = sub[i];
990 :
991 0 : return ray;
992 : }
993 :
994 0 : void BeamCalc::deleteRay(Ray *ray)
995 : {
996 0 : if(ray) free(ray);
997 0 : }
998 :
999 0 : Pathology* BeamCalc::newPathology()
1000 : {
1001 : Pathology *P;
1002 : Int i, j;
1003 :
1004 0 : P = (Pathology *)malloc(sizeof(Pathology));
1005 :
1006 0 : for(i = 0; i < 3; i++) P->subrotpoint[i] = P->subshift[i] = P->feedshift[i] = 0.0;
1007 0 : for(i = 0; i < 3; i++) for(j = 0; j < 3; j++)
1008 0 : P->feedrot[i][j] = P->subrot[i][j] = 0.0;
1009 0 : for(i = 0; i < 3; i++) P->feedrot[i][i] = P->subrot[i][i] = 1.0;
1010 :
1011 0 : P->az_offset = 0.0;
1012 0 : P->el_offset = 0.0;
1013 0 : P->phase_offset = 0.0;
1014 0 : P->focus = 0.0;
1015 :
1016 0 : return P;
1017 : }
1018 :
1019 0 : Pathology* BeamCalc::newPathologyfromApertureCalcParams(ApertureCalcParams* /*ap*/)
1020 : {
1021 : Pathology *P;
1022 :
1023 0 : P = newPathology();
1024 :
1025 0 : return P;
1026 : }
1027 :
1028 0 : void BeamCalc::deletePathology(Pathology *P)
1029 : {
1030 0 : if(P) free(P);
1031 0 : }
1032 :
1033 0 : void BeamCalc::normvec(const Double *a, const Double *b, Double *c)
1034 : {
1035 : Int i;
1036 : Double r;
1037 0 : for(i = 0; i < 3; i++) c[i] = a[i] - b[i];
1038 0 : r = sqrt(c[0]*c[0] + c[1]*c[1] + c[2]*c[2]);
1039 0 : for(i = 0; i < 3; i++) c[i] /= r;
1040 0 : }
1041 :
1042 0 : Double BeamCalc::dAdOmega(const calcAntenna *a, const Ray *ray1, const Ray *ray2,
1043 : const Ray *ray3, const Pathology *p)
1044 : {
1045 : Double A, Omega;
1046 : Double n1[3], n2[3], n3[3], f[3], ci, cj, ck;
1047 : Int i;
1048 :
1049 : /* Area in aperture is in a plane z = const */
1050 0 : A = 0.5*fabs(
1051 0 : (ray1->aper[0]-ray2->aper[0])*(ray1->aper[1]-ray3->aper[1]) -
1052 0 : (ray1->aper[0]-ray3->aper[0])*(ray1->aper[1]-ray2->aper[1]) );
1053 :
1054 0 : for(i = 0; i < 3; i++) f[i] = a->feed[i] + p->feedshift[i];
1055 :
1056 0 : normvec(ray1->sub, f, n1);
1057 0 : normvec(ray2->sub, f, n2);
1058 0 : normvec(ray3->sub, f, n3);
1059 :
1060 0 : for(i = 0; i < 3; i++)
1061 : {
1062 0 : n1[i] -= n3[i];
1063 0 : n2[i] -= n3[i];
1064 : }
1065 :
1066 0 : ci = n1[1]*n2[2] - n1[2]*n2[1];
1067 0 : cj = n1[2]*n2[0] - n1[0]*n2[2];
1068 0 : ck = n1[0]*n2[1] - n1[1]*n2[0];
1069 :
1070 0 : Omega = 0.5*sqrt(ci*ci + cj*cj + ck*ck);
1071 :
1072 0 : return A/Omega;
1073 : }
1074 :
1075 0 : Double BeamCalc::dOmega(const calcAntenna *a, const Ray *ray1, const Ray *ray2,
1076 : const Ray *ray3, const Pathology *p)
1077 : {
1078 : Double Omega;
1079 : Double n1[3], n2[3], n3[3], f[3], ci, cj, ck;
1080 : Int i;
1081 :
1082 0 : for(i = 0; i < 3; i++) f[i] = a->feed[i] + p->feedshift[i];
1083 :
1084 0 : normvec(ray1->sub, f, n1);
1085 0 : normvec(ray2->sub, f, n2);
1086 0 : normvec(ray3->sub, f, n3);
1087 :
1088 0 : for(i = 0; i < 3; i++)
1089 : {
1090 0 : n1[i] -= n3[i];
1091 0 : n2[i] -= n3[i];
1092 : }
1093 :
1094 0 : ci = n1[1]*n2[2] - n1[2]*n2[1];
1095 0 : cj = n1[2]*n2[0] - n1[0]*n2[2];
1096 0 : ck = n1[0]*n2[1] - n1[1]*n2[0];
1097 :
1098 0 : Omega = 0.5*sqrt(ci*ci + cj*cj + ck*ck);
1099 :
1100 0 : return Omega;
1101 : }
1102 :
1103 0 : Double BeamCalc::Raylen(const Ray *ray)
1104 : {
1105 : Double len, d[3];
1106 : Int i;
1107 :
1108 : /* feed to subreflector */
1109 0 : for(i = 0; i < 3; i++) d[i] = ray->feed[i] - ray->sub[i];
1110 0 : len = sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
1111 :
1112 : /* subreflector to dish */
1113 0 : for(i = 0; i < 3; i++) d[i] = ray->sub[i] - ray->dish[i];
1114 0 : len += sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
1115 :
1116 : /* dish to aperture */
1117 0 : for(i = 0; i < 3; i++) d[i] = ray->dish[i] - ray->aper[i];
1118 0 : len += sqrt(d[0]*d[0] + d[1]*d[1] + d[2]*d[2]);
1119 :
1120 0 : return len;
1121 : }
1122 :
1123 0 : void BeamCalc::Pathologize(Double *sub, const Pathology *p)
1124 : {
1125 : Int i;
1126 : Int j;
1127 : Double tmp[6];
1128 :
1129 0 : for(i = 0; i < 3; i++) sub[i] -= p->subrotpoint[i];
1130 :
1131 0 : for(i = 0; i < 3; i++)
1132 : {
1133 0 : tmp[i] = 0.0;
1134 0 : tmp[i+3] = 0.0;
1135 0 : for(j = 0; j < 3; j++) tmp[i] += p->subrot[i][j]*sub[j];
1136 0 : for(j = 0; j < 3; j++) tmp[i+3] += p->subrot[i][j]*sub[j+3];
1137 : }
1138 :
1139 0 : for(i = 0; i < 3; i++)
1140 0 : sub[i] = tmp[i] + p->subrotpoint[i] + p->subshift[i];
1141 0 : for(i = 4; i < 6; i++)
1142 0 : sub[i] = tmp[i];
1143 0 : }
1144 :
1145 0 : void BeamCalc::applyPathology(Pathology *P, calcAntenna *a)
1146 : {
1147 : Double dx[3];
1148 : Int i, j;
1149 :
1150 0 : if(P->focus != 0.0)
1151 : {
1152 0 : dx[0] = -a->feed[0];
1153 0 : dx[1] = -a->feed[1];
1154 0 : dx[2] = a->sub_h-a->feed[2];
1155 0 : norm3(dx);
1156 0 : for(i = 0; i < 3; i++) P->feedshift[i] += P->focus*dx[i];
1157 :
1158 0 : P->focus = 0.0;
1159 : }
1160 :
1161 0 : for(i = 0; i < 3; i++) a->pfeeddir[i] = 0.0;
1162 0 : for(j = 0; j < 3; j++) for(i = 0; i < 3; i++)
1163 0 : a->pfeeddir[j] += P->feedrot[j][i]*a->feeddir[i];
1164 0 : }
1165 :
1166 :
1167 0 : void BeamCalc::intersectdish(const calcAntenna *a, const Double *sub, const Double *unitdir,
1168 : Double *dish, Int niter)
1169 : {
1170 : Double A, B, C, t, m, r;
1171 : Double x[3], n[3];
1172 : Int i, iter;
1173 :
1174 : /* First intersect with ideal paraboloid */
1175 0 : A = a->bestparabola*(unitdir[0]*unitdir[0] + unitdir[1]*unitdir[1]);
1176 0 : B = 2.0*a->bestparabola*(unitdir[0]*sub[0] + unitdir[1]*sub[1])
1177 0 : -unitdir[2];
1178 0 : C = a->bestparabola*(sub[0]*sub[0] + sub[1]*sub[1]) - sub[2];
1179 0 : t = 0.5*(sqrt(B*B-4.0*A*C) - B)/A; /* take greater root */
1180 :
1181 0 : for(iter = 0; ; iter++)
1182 : {
1183 : /* get position (x) and normal (n) on the real dish */
1184 0 : for(i = 0; i < 2; i++) x[i] = sub[i] + t*unitdir[i];
1185 0 : r = sqrt(x[0]*x[0] + x[1]*x[1]);
1186 0 : astigdishvalue(a, x[0], x[1], &(x[2]), &m);
1187 0 : n[2] = 1.0/sqrt(1.0+m*m);
1188 0 : n[0] = -m*(x[0]/r)*n[2];
1189 0 : n[1] = -m*(x[1]/r)*n[2];
1190 :
1191 0 : if(iter >= niter) break;
1192 :
1193 0 : A = B = 0;
1194 0 : for(i = 0; i < 3; i++)
1195 : {
1196 0 : A += n[i]*(x[i]-sub[i]); /* n dot (x-sub) */
1197 0 : B += n[i]*unitdir[i]; /* n dot unitdir */
1198 : }
1199 0 : t = A/B;
1200 : }
1201 :
1202 0 : for(i = 0; i < 3; i++)
1203 : {
1204 0 : dish[i] = x[i];
1205 0 : dish[i+3] = n[i];
1206 : }
1207 0 : }
1208 :
1209 0 : void BeamCalc::intersectaperture(const calcAntenna *a, const Double *dish,
1210 : const Double *unitdir, Double *aper)
1211 : {
1212 : Double t;
1213 : Int i;
1214 :
1215 0 : t = (a->zedge-dish[2])/unitdir[2];
1216 0 : for(i = 0; i < 3; i++) aper[i] = dish[i] + t*unitdir[i];
1217 :
1218 0 : aper[3] = aper[4] = 0.0;
1219 0 : aper[5] = 1.0;
1220 0 : }
1221 :
1222 : /* gain in power */
1223 0 : Double BeamCalc::feedfunc(const calcAntenna *a, Double theta)
1224 : {
1225 : Double stheta;
1226 :
1227 0 : stheta = sin(theta);
1228 0 : return exp(2.0*(-0.083)*a->fa2pi*a->fa2pi*stheta*stheta);
1229 : }
1230 :
1231 : /* gain in power */
1232 0 : Double BeamCalc::feedgain(const calcAntenna *a, const Ray *ray, const Pathology */*p*/)
1233 : {
1234 0 : Double costheta = 0.0;
1235 : Double v[3];
1236 : Int i;
1237 :
1238 0 : for(i = 0; i < 3; i++) v[i] = ray->sub[i] - ray->feed[i];
1239 0 : norm3(v);
1240 :
1241 0 : for(i = 0; i < 3; i++)
1242 : {
1243 0 : costheta += a->pfeeddir[i]*v[i];
1244 : }
1245 :
1246 :
1247 0 : return exp(2.0*(-0.083)*a->fa2pi*a->fa2pi*(1.0-costheta*costheta));
1248 : }
1249 :
1250 0 : Ray* BeamCalc::trace(const calcAntenna *a, Double x, Double y, const Pathology *p)
1251 : {
1252 : Ray *ray;
1253 : Double idealsub[12];
1254 0 : Double fu[3], du[3], au[3], ndotf=0.0, ndotd=0.0;
1255 : Int i;
1256 0 : const Int niter = 7;
1257 :
1258 0 : subfromdish(a, x, y, idealsub);
1259 :
1260 0 : ray = newRay(idealsub);
1261 :
1262 0 : Pathologize(ray->sub, p);
1263 :
1264 0 : if(ray->sub[5] < -1.0 || ray->sub[5] > -0.0)
1265 : {
1266 0 : deleteRay(ray);
1267 0 : return 0;
1268 : }
1269 :
1270 0 : for(i = 0; i < 3; i++) ray->feed[i] = a->feed[i] + p->feedshift[i];
1271 :
1272 : /* now determine unit vector pointing to dish */
1273 :
1274 : /* unit toward feed */
1275 0 : for(i = 0; i < 3; i++) fu[i] = ray->feed[i] - ray->sub[i];
1276 0 : norm3(fu);
1277 :
1278 : /* unit toward dish */
1279 0 : for(i = 0; i < 3; i++) ndotf += ray->sub[i+3]*fu[i];
1280 0 : for(i = 0; i < 3; i++) du[i] = 2.0*ray->sub[i+3]*ndotf - fu[i];
1281 :
1282 : /* dish point */
1283 0 : intersectdish(a, ray->sub, du, ray->dish, niter);
1284 :
1285 : /* unit toward aperture */
1286 0 : for(i = 0; i < 3; i++) ndotd += ray->dish[i+3]*du[i];
1287 0 : for(i = 0; i < 3; i++) au[i] = du[i] - 2.0*ray->dish[i+3]*ndotd;
1288 :
1289 0 : intersectaperture(a, ray->dish, au, ray->aper);
1290 :
1291 0 : return ray;
1292 : }
1293 :
1294 0 : void BeamCalc::tracepol(Complex *E0, const Ray *ray, Complex *E1)
1295 : {
1296 0 : Complex fac;
1297 : Double v1[3], v2[3], v3[3];
1298 : Double r[3];
1299 : Int i;
1300 :
1301 0 : for(i = 0; i < 3; i++)
1302 : {
1303 0 : v1[i] = ray->sub[i] - ray->feed[i];
1304 0 : v2[i] = ray->dish[i] - ray->sub[i];
1305 0 : v3[i] = ray->aper[i] - ray->dish[i];
1306 0 : E1[i] = E0[i];
1307 : }
1308 0 : norm3(v1);
1309 0 : norm3(v2);
1310 0 : norm3(v3);
1311 :
1312 0 : for(i = 0; i < 3; i++) r[i] = v1[i] - v2[i];
1313 0 : norm3(r);
1314 0 : fac = Complex(r[0],0)*E1[0] + Complex(r[1],0)*E1[1] + Complex(r[2],0)*E1[2];
1315 0 : for(i = 0; i < 3; i++) E1[i] = Complex(r[i],0)*fac*2.0 - E1[i];
1316 :
1317 0 : for(i = 0; i < 3; i++) r[i] = v2[i] - v3[i];
1318 0 : norm3(r);
1319 0 : fac = Complex(r[0],0)*E1[0] + Complex(r[1],0)*E1[1] + Complex(r[2],0)*E1[2];
1320 0 : for(i = 0; i < 3; i++) E1[i] = Complex(r[i],0)*fac*2.0 - E1[i];
1321 0 : }
1322 :
1323 0 : Int BeamCalc::legplanewaveblock(const calcAntenna *a, Double x, Double y)
1324 : {
1325 : /* outside the leg foot area, the blockage is spherical wave */
1326 0 : if(x*x + y*y > a->legfoot*a->legfoot) return 0;
1327 :
1328 0 : if(a->legwidth == 0.0) return 0;
1329 :
1330 0 : if(strcmp(a->name, "VLBA") == 0)
1331 : {
1332 0 : const Double s=1.457937;
1333 0 : const Double c=1.369094;
1334 0 : if(fabs(s*x+c*y) < -a->legwidth) return 1;
1335 0 : if(fabs(s*x-c*y) < -a->legwidth) return 1;
1336 : }
1337 0 : else if(a->legwidth < 0.0) /* "x shaped" legs */
1338 : {
1339 0 : if(fabs(x-y)*M_SQRT2 < -a->legwidth) return 1;
1340 0 : if(fabs(x+y)*M_SQRT2 < -a->legwidth) return 1;
1341 : }
1342 0 : else if(a->legwidth > 0.0) /* "+ shaped" legs */
1343 : {
1344 0 : if(fabs(x)*2.0 < a->legwidth) return 1;
1345 0 : if(fabs(y)*2.0 < a->legwidth) return 1;
1346 : }
1347 :
1348 0 : return 0;
1349 : }
1350 :
1351 0 : Int BeamCalc::legplanewaveblock2(const calcAntenna *a, const Ray *ray)
1352 : {
1353 : Int i, n;
1354 : Double dr2;
1355 : // phi set but not used
1356 : Double theta /*, phi*/;
1357 : Double r0[3], dr[3], l0[3], l1[3], dl[3], D[3];
1358 : Double D2, N[3], ll, rr;
1359 0 : const Double thetaplus[4] =
1360 : {0, M_PI/2.0, M_PI, 3.0*M_PI/2.0};
1361 0 : const Double thetacross[4] =
1362 : {0.25*M_PI, 0.75*M_PI, 1.25*M_PI, 1.75*M_PI};
1363 0 : const Double thetavlba[4] =
1364 : {0.816817, 2.3247756, 3.9584096, 5.466368};
1365 : const Double *thetalut;
1366 :
1367 0 : if(a->legwidth == 0.0) return 0;
1368 :
1369 0 : if(strcmp(a->name, "VLBA") == 0) thetalut = thetavlba;
1370 0 : else if(a->legwidth < 0.0) thetalut = thetacross;
1371 0 : else thetalut = thetaplus;
1372 :
1373 : /* inside the leg feet is plane wave blockage */
1374 0 : dr2 = ray->dish[0]*ray->dish[0] + ray->dish[1]*ray->dish[1];
1375 0 : if(dr2 >= a->legfoot*a->legfoot) return 0;
1376 :
1377 0 : for(i = 0; i < 3; i++)
1378 : {
1379 0 : r0[i] = ray->dish[i];
1380 0 : dr[i] = ray->aper[i] - r0[i];
1381 : }
1382 0 : rr = r0[0]*r0[0] + r0[1]*r0[1];
1383 :
1384 0 : l0[2] = a->legfootz;
1385 0 : l1[0] = l1[1] = 0.0;
1386 0 : l1[2] = a->legapex;
1387 : // phi set but not used
1388 : // phi = atan2(r0[1], r0[0]);
1389 0 : for(n = 0; n < 4; n++)
1390 : {
1391 0 : theta = thetalut[n];
1392 0 : l0[0] = a->legfoot*cos(theta);
1393 0 : l0[1] = a->legfoot*sin(theta);
1394 0 : ll = l0[0]*l0[0] + l0[1]*l0[1];
1395 0 : if((l0[0]*r0[0] + l0[1]*r0[1])/sqrt(ll*rr) < 0.7) continue;
1396 0 : for(i = 0; i < 3; i++) dl[i] = l1[i] - l0[i];
1397 0 : for(i = 0; i < 3; i++) D[i] = r0[i] - l0[i];
1398 :
1399 0 : N[0] = dr[1]*dl[2] - dr[2]*dl[1];
1400 0 : N[1] = dr[2]*dl[0] - dr[0]*dl[2];
1401 0 : N[2] = dr[0]*dl[1] - dr[1]*dl[0];
1402 0 : norm3(N);
1403 :
1404 0 : D2 = D[0]*N[0] + D[1]*N[1] + D[2]*N[2];
1405 :
1406 0 : if(fabs(D2) <= 0.5*fabs(a->legwidth)) return 1;
1407 : }
1408 :
1409 0 : return 0;
1410 : }
1411 :
1412 0 : Int BeamCalc::legsphericalwaveblock(const calcAntenna *a, const Ray *ray)
1413 : {
1414 : Int i, n;
1415 : Double dr2;
1416 : // phi set but not used
1417 : Double theta /*, phi*/;
1418 : Double r0[3], dr[3], l0[3], l1[3], dl[3], D[3];
1419 : Double D2, N[3], ll, rr;
1420 0 : const Double thetaplus[4] =
1421 : {0, M_PI/2.0, M_PI, 3.0*M_PI/2.0};
1422 0 : const Double thetacross[4] =
1423 : {0.25*M_PI, 0.75*M_PI, 1.25*M_PI, 1.75*M_PI};
1424 0 : const Double thetavlba[4] =
1425 : {0.816817, 2.3247756, 3.9584096, 5.466368};
1426 : const Double *thetalut;
1427 :
1428 0 : if(a->legwidth == 0.0) return 0;
1429 :
1430 0 : if(strcmp(a->name, "VLBA") == 0) thetalut = thetavlba;
1431 0 : else if(a->legwidth < 0.0) thetalut = thetacross;
1432 0 : else thetalut = thetaplus;
1433 :
1434 : /* inside the leg feet is plane wave blockage */
1435 0 : dr2 = ray->dish[0]*ray->dish[0] + ray->dish[1]*ray->dish[1];
1436 0 : if(dr2 < a->legfoot*a->legfoot) return 0;
1437 :
1438 0 : for(i = 0; i < 3; i++)
1439 : {
1440 0 : r0[i] = ray->dish[i];
1441 0 : dr[i] = ray->sub[i] - r0[i];
1442 : }
1443 0 : rr = r0[0]*r0[0] + r0[1]*r0[1];
1444 :
1445 0 : l0[2] = a->legfootz;
1446 0 : l1[0] = l1[1] = 0.0;
1447 0 : l1[2] = a->legapex;
1448 : // phi set but not used
1449 : // phi = atan2(r0[1], r0[0]);
1450 0 : for(n = 0; n < 4; n++)
1451 : {
1452 0 : theta = thetalut[n];
1453 0 : l0[0] = a->legfoot*cos(theta);
1454 0 : l0[1] = a->legfoot*sin(theta);
1455 0 : ll = l0[0]*l0[0] + l0[1]*l0[1];
1456 0 : if((l0[0]*r0[0] + l0[1]*r0[1])/sqrt(ll*rr) < 0.7) continue;
1457 0 : for(i = 0; i < 3; i++) dl[i] = l1[i] - l0[i];
1458 0 : for(i = 0; i < 3; i++) D[i] = r0[i] - l0[i];
1459 :
1460 0 : N[0] = dr[1]*dl[2] - dr[2]*dl[1];
1461 0 : N[1] = dr[2]*dl[0] - dr[0]*dl[2];
1462 0 : N[2] = dr[0]*dl[1] - dr[1]*dl[0];
1463 0 : norm3(N);
1464 :
1465 0 : D2 = D[0]*N[0] + D[1]*N[1] + D[2]*N[2];
1466 :
1467 0 : if(fabs(D2) <= 0.5*fabs(a->legwidth)) return 1;
1468 : }
1469 :
1470 0 : return 0;
1471 : }
1472 :
1473 :
1474 0 : void BeamCalc::copyBeamCalcGeometry(BeamCalcGeometry* to, BeamCalcGeometry* from){
1475 0 : sprintf(to->name, "%s", from->name);
1476 0 : sprintf(to->bandName, "%s", from->bandName);
1477 0 : to->sub_h = from->sub_h;
1478 0 : for(uInt j=0; j<3;j++){
1479 0 : to->feedpos[j] = from->feedpos[j];
1480 : }
1481 0 : to->subangle = from->subangle;
1482 0 : to->legwidth = from->legwidth;
1483 0 : to->legfoot = from->legfoot;
1484 0 : to->legapex = from->legapex;
1485 0 : to->Rhole = from->Rhole;
1486 0 : to->Rant = from->Rant;
1487 0 : to->reffreq = from->reffreq;
1488 0 : for(uInt j=0; j<5;j++){
1489 0 : to->taperpoly[j] = from->taperpoly[j];
1490 : }
1491 0 : to->ntaperpoly = from->ntaperpoly;
1492 0 : to->astigm_0 = from->astigm_0;
1493 0 : to->astigm_45 = from->astigm_45;
1494 :
1495 0 : }
1496 :
1497 :
1498 : /* The meat of the calculation */
1499 :
1500 :
1501 0 : Int BeamCalc::calculateAperture(ApertureCalcParams *ap)
1502 : {
1503 : // not used
1504 : // Complex fp, Exr, Eyr, Exl, Eyl;
1505 0 : Complex Er[3], El[3];
1506 0 : Complex Pr[2], Pl[2];
1507 0 : Complex q[2];
1508 : //Double dx, dy, Rhole, Rant, x0, y0, R2, H2, eps;
1509 : //Complex rr, rl, lr, ll, tmp;
1510 : Double L0, phase;
1511 : Double sp, cp;
1512 : Double B[3][3];
1513 : calcAntenna *a;
1514 : Pathology *p;
1515 : Int nx, ny, os;
1516 : Int i, j;
1517 : //Double pac, pas; /* parallactic angle cosine / sine */
1518 0 : Complex Iota; Iota=Complex(0,1);
1519 :
1520 : //UNUSED: Complex E1[3];
1521 : //UNUSED: Double x,y, r2, L, amp, dP, dA, d0, x1, y1, dx1, dy1, dx2, dy2, dO;
1522 : //UNUSED: Ray *ray, *rayx, *rayy;
1523 : //UNUSED: Int iter;
1524 : //UNUSED: Int niter=6;
1525 :
1526 0 : a = newAntennafromApertureCalcParams(ap);
1527 0 : p = newPathologyfromApertureCalcParams(ap);
1528 :
1529 : /* compute central ray pathlength */
1530 : {
1531 : Ray *tmpRay;
1532 0 : tmpRay = trace(a, 0.0, 0.00001, p);
1533 0 : L0 = Raylen(tmpRay);
1534 0 : deleteRay(tmpRay);
1535 : }
1536 :
1537 : //pac = cos(ap->pa+M_PI/2);
1538 : //pas = sin(ap->pa+M_PI/2);
1539 :
1540 0 : if(obsName_p=="EVLA" || obsName_p=="VLA"){
1541 : /* compute polarization vectors in circular basis */
1542 0 : Pr[0] = 1.0/M_SQRT2; Pr[1] = Iota/M_SQRT2;
1543 0 : Pl[0] = 1.0/M_SQRT2; Pl[1] = -Iota/M_SQRT2;
1544 :
1545 : /* compensate for feed orientation */
1546 0 : getfeedbasis(a, B);
1547 0 : phase = atan2(B[0][1], B[0][0]);
1548 0 : cp = cos(phase);
1549 0 : sp = sin(phase);
1550 :
1551 0 : q[0] = Pr[0];
1552 0 : q[1] = Pr[1];
1553 0 : Pr[0] = Complex(cp,0)*q[0] + Complex(sp,0)*q[1];
1554 0 : Pr[1] = -Complex(sp,0)*q[0] + Complex(cp,0)*q[1];
1555 0 : q[0] = Pl[0];
1556 0 : q[1] = Pl[1];
1557 0 : Pl[0] = Complex(cp,0)*q[0] + Complex(sp,0)*q[1];
1558 0 : Pl[1] = -Complex(sp,0)*q[0] + Complex(cp,0)*q[1];
1559 : }
1560 : else{
1561 : /* in linear basis */
1562 0 : Pr[0] = 1.0; Pr[1] = 0.0;
1563 0 : Pl[0] = 0.0; Pl[1] = 1.0;
1564 : }
1565 :
1566 :
1567 :
1568 : /* compute 3-vector feed efields for the two polarizations */
1569 0 : Efield(a, Pr, Er);
1570 0 : Efield(a, Pl, El);
1571 0 : ap->aperture->set(Complex(0.0));
1572 :
1573 0 : os = ap->oversamp;
1574 0 : nx = ap->nx*os;
1575 0 : ny = ap->ny*os;
1576 : //dx = ap->dx/os;
1577 : //dy = ap->dy/os;
1578 : //x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
1579 : //y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
1580 : //Rant = a->radius;
1581 : //Rhole = a->hole_radius;
1582 : //R2 = Rant*Rant;
1583 : //H2 = Rhole*Rhole;
1584 :
1585 : //eps = dx/4.0;
1586 :
1587 0 : IPosition pos(4);
1588 : // shape = ap->aperture->shape();
1589 :
1590 :
1591 : // cerr << "max threads " << omp_get_max_threads()
1592 : // << " threads available " << omp_get_num_threads()
1593 : // << endl;
1594 : //Int Nth=1;
1595 : //#ifdef _OPENMP
1596 : //Nth=max(omp_get_max_threads()-2,1);
1597 : //#endif
1598 : // Timer tim;
1599 : // tim.mark();
1600 : //#pragma omp parallel default(none) firstprivate(Er, El, nx, ny) private(i,j) shared(ap, a, p, L0) num_threads(Nth)
1601 : {
1602 : //#pragma omp for
1603 0 : for(j = 0; j < ny; j++)
1604 : {
1605 0 : for(i = 0; i < nx; i++)
1606 : {
1607 0 : computePixelValues(ap, a, p, L0, Er, El, i,j);
1608 : }
1609 : }
1610 : }
1611 : // tim.show("BeamCalc:");
1612 :
1613 0 : deletePathology(p);
1614 0 : deleteAntenna(a);
1615 :
1616 0 : return 1;
1617 0 : }
1618 :
1619 0 : void BeamCalc::computePixelValues(const ApertureCalcParams *ap,
1620 : const calcAntenna *a, const Pathology *p,
1621 : const Double &L0,
1622 : Complex *Er, Complex *El,
1623 : const Int &i, const Int &j)
1624 : {
1625 0 : Complex fp, Exr, Eyr, Exl, Eyl;
1626 : // Complex Er[3], El[3];
1627 0 : Complex E1[3];
1628 : Double dx, dy, Rant, x0, y0, x, y, r2, R2, H2, eps, Rhole;
1629 0 : Complex rr, rl, lr, ll, tmp;
1630 : Double L, amp, dP, dA, dO, x1, y1, dx1, dy1, dx2, dy2, phase;
1631 : //Int nx, ny;
1632 : Int os;
1633 0 : Int niter=6;
1634 : Double pac, pas;
1635 : Double cp,sp; /* parallactic angle cosine / sine */
1636 0 : Complex Iota; Iota=Complex(0,1);
1637 0 : IPosition pos(4);pos=0;
1638 :
1639 0 : Ray *ray=0, *rayx=0, *rayy=0;
1640 : /* determine parallactic angle rotated coordinates */
1641 :
1642 0 : os = ap->oversamp;
1643 : //nx = ap->nx*os;
1644 : //ny = ap->ny*os;
1645 0 : dx = ap->dx/os;
1646 0 : dy = ap->dy/os;
1647 0 : x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
1648 0 : y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
1649 0 : Rant = a->radius;
1650 0 : Rhole = a->hole_radius;
1651 0 : R2 = Rant*Rant;
1652 0 : H2 = Rhole*Rhole;
1653 : // for(Int i=0; i < nx; ++i)
1654 : {
1655 0 : eps = dx/4.0;
1656 0 : pac = cos(ap->pa+M_PI/2);
1657 0 : pas = sin(ap->pa+M_PI/2);
1658 :
1659 0 : x = pac*(x0 + i*dx) - pas*(y0 + j*dy);
1660 0 : y = pas*(x0 + i*dx) + pac*(y0 + j*dy);
1661 0 : x = -x;
1662 :
1663 0 : if(fabs(x) > Rant) goto nextpoint;
1664 0 : if(fabs(y) > Rant) goto nextpoint;
1665 0 : r2 = x*x + y*y;
1666 0 : if(r2 > R2) goto nextpoint;
1667 0 : if(r2 < H2) goto nextpoint;
1668 :
1669 0 : ray = rayx = rayy = 0;
1670 :
1671 0 : x1 = x;
1672 0 : y1 = y;
1673 :
1674 0 : for(Int iter = 0; iter < niter; iter++)
1675 : {
1676 0 : ray = trace(a, x1, y1, p);
1677 0 : if(!ray) goto nextpoint;
1678 0 : x1 += (x - ray->aper[0]);
1679 0 : y1 += (y - ray->aper[1]);
1680 0 : deleteRay(ray);
1681 0 : ray = 0;
1682 : }
1683 :
1684 0 : ray = trace(a, x1, y1, p);
1685 :
1686 : /* check for leg blockage */
1687 0 : if(legplanewaveblock2(a, ray))
1688 0 : goto nextpoint;
1689 0 : if(legsphericalwaveblock(a, ray))
1690 0 : goto nextpoint;
1691 :
1692 0 : if(y < 0) rayy = trace(a, x1, y1+eps, p);
1693 0 : else rayy = trace(a, x1, y1-eps, p);
1694 :
1695 0 : if(x < 0) rayx = trace(a, x1+eps, y1, p);
1696 0 : else rayx = trace(a, x1-eps, y1, p);
1697 :
1698 0 : if(ray == 0 || rayx == 0 || rayy == 0)
1699 0 : goto nextpoint;
1700 :
1701 : /* compute solid angle subtended at the feed */
1702 0 : dx1 = rayx->aper[0]-ray->aper[0];
1703 0 : dy1 = rayx->aper[1]-ray->aper[1];
1704 0 : dx2 = rayy->aper[0]-ray->aper[0];
1705 0 : dy2 = rayy->aper[1]-ray->aper[1];
1706 :
1707 0 : dA = 0.5*fabs(dx1*dy2 - dx2*dy1);
1708 0 : dO = (dOmega(a, rayx, rayy, ray, p)/dA)*dx*dx;
1709 0 : dP = dO*feedgain(a, ray, p);
1710 0 : amp = sqrt(dP);
1711 :
1712 0 : L = Raylen(ray);
1713 :
1714 0 : phase = 2.0*M_PI*(L-L0)/a->lambda;
1715 :
1716 : /* phase retard the wave */
1717 0 : cp = cos(phase);
1718 0 : sp = sin(phase);
1719 : // fp = cp + sp*1.0i;
1720 :
1721 0 : fp = Complex(cp,sp);
1722 :
1723 :
1724 0 : tracepol(Er, ray, E1);
1725 0 : Exr = fp*amp*E1[0];
1726 0 : Eyr = fp*amp*E1[1];
1727 :
1728 : // rr = Exr - 1.0i*Eyr;
1729 : // rl = Exr + 1.0i*Eyr;
1730 0 : rr = Exr - Iota*Eyr;
1731 0 : rl = Exr + Iota*Eyr;
1732 :
1733 0 : tracepol(El, ray, E1);
1734 0 : Exl = fp*amp*E1[0];
1735 0 : Eyl = fp*amp*E1[1];
1736 : // lr = Exl - 1.0i*Eyl;
1737 : // ll = Exl + 1.0i*Eyl;
1738 0 : lr = Exl - Iota*Eyl;
1739 0 : ll = Exl + Iota*Eyl;
1740 : // pos(0)=(Int)((j/os) - (25.0/dy/os)/2 + shape(0)/2 - 0.5);
1741 : // pos(1)=(Int)((i/os) - (25.0/dx/os)/2 + shape(1)/2 - 0.5);
1742 : // Following 3 lines go with ANT tag in VLACalc.....
1743 : // Double antRadius=BeamCalcGeometryDefaults[ap->band].Rant;
1744 : // pos(0)=(Int)((j/os) - (antRadius/dy/os) + shape(0)/2 - 0.5);
1745 : // pos(1)=(Int)((i/os) - (antRadius/dx/os) + shape(1)/2 - 0.5);
1746 : // Following 2 lines go with the PIX tag in VLACalc...
1747 0 : pos(0)=(Int)((j/os));
1748 0 : pos(1)=(Int)((i/os));
1749 0 : pos(3)=0;
1750 :
1751 0 : pos(2)=0;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+rr,pos);
1752 0 : pos(2)=1;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+rl,pos);
1753 0 : pos(2)=2;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+lr,pos);
1754 0 : pos(2)=3;tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+ll,pos);
1755 0 : nextpoint:
1756 0 : if(ray) deleteRay(ray);
1757 0 : if(rayx) deleteRay(rayx);
1758 0 : if(rayy) deleteRay(rayy);
1759 : }
1760 0 : }
1761 : //
1762 : //----------------------------------------------------------------------
1763 : // Compute only the required polarizations.
1764 : //
1765 0 : Int BeamCalc::calculateAperture(ApertureCalcParams *ap, const Int& whichPoln)
1766 : {
1767 : //Complex fp, Exr, Eyr, Exl, Eyl;
1768 0 : Complex Er[3], El[3];
1769 0 : Complex Pr[2], Pl[2];
1770 0 : Complex q[2];
1771 : //Double dx, dy, Rhole, Rant;//x0, y0, R2, H2, eps;
1772 : //Complex rr, rl, lr, ll, tmp;
1773 : Double L0, phase;
1774 : Double sp, cp;
1775 : Double B[3][3];
1776 : calcAntenna *a;
1777 : Pathology *p;
1778 : Int nx, ny, os;
1779 : Int i, j;
1780 : //Double pac, pas; /* parallactic angle cosine / sine */
1781 0 : Complex Iota; Iota=Complex(0,1);
1782 :
1783 : //UNUSED: Complex E1[3];
1784 : //UNUSED: Double x,y, r2, L, amp, dP, dA, d0, x1, y1, dx1, dy1, dx2, dy2, dO;
1785 : //UNUSED: Ray *ray, *rayx, *rayy;
1786 : //UNUSED: Int iter;
1787 : //UNUSED: Int niter=6;
1788 :
1789 0 : a = newAntennafromApertureCalcParams(ap);
1790 0 : p = newPathologyfromApertureCalcParams(ap);
1791 :
1792 : /* compute central ray pathlength */
1793 : {
1794 : Ray *tmpRay;
1795 0 : tmpRay = trace(a, 0.0, 0.00001, p);
1796 0 : L0 = Raylen(tmpRay);
1797 0 : deleteRay(tmpRay);
1798 : }
1799 :
1800 : //pac = cos(ap->pa+M_PI/2);
1801 : //pas = sin(ap->pa+M_PI/2);
1802 :
1803 0 : if(obsName_p=="EVLA" || obsName_p=="VLA"){
1804 : /* compute polarization vectors in circular basis */
1805 0 : Pr[0] = 1.0/M_SQRT2; Pr[1] = Iota/M_SQRT2;
1806 0 : Pl[0] = 1.0/M_SQRT2; Pl[1] = -Iota/M_SQRT2;
1807 :
1808 : /* compensate for feed orientation */
1809 0 : getfeedbasis(a, B);
1810 0 : phase = atan2(B[0][1], B[0][0]);
1811 0 : cp = cos(phase);
1812 0 : sp = sin(phase);
1813 :
1814 0 : q[0] = Pr[0];
1815 0 : q[1] = Pr[1];
1816 0 : Pr[0] = Complex(cp,0)*q[0] + Complex(sp,0)*q[1];
1817 0 : Pr[1] = -Complex(sp,0)*q[0] + Complex(cp,0)*q[1];
1818 0 : q[0] = Pl[0];
1819 0 : q[1] = Pl[1];
1820 0 : Pl[0] = Complex(cp,0)*q[0] + Complex(sp,0)*q[1];
1821 0 : Pl[1] = -Complex(sp,0)*q[0] + Complex(cp,0)*q[1];
1822 : }
1823 : else{
1824 : /* in linear basis */
1825 0 : Pr[0] = 1.0; Pr[1] = 0.0;
1826 0 : Pl[0] = 0.0; Pl[1] = 1.0;
1827 : }
1828 :
1829 :
1830 :
1831 : /* compute 3-vector feed efields for the two polarizations */
1832 0 : if ((whichPoln == Stokes::RR) || (whichPoln == Stokes::XX))
1833 0 : Efield(a, Pr, Er);
1834 0 : else if ((whichPoln == Stokes::LL) || (whichPoln == Stokes::YY))
1835 0 : Efield(a, Pl, El);
1836 : else
1837 0 : {Efield(a, Pr, Er); Efield(a, Pl, El);}
1838 :
1839 0 : ap->aperture->set(Complex(0.0));
1840 :
1841 0 : os = ap->oversamp;
1842 0 : nx = ap->nx*os;
1843 0 : ny = ap->ny*os;
1844 : // dx = ap->dx/os;
1845 : // dy = ap->dy/os;
1846 : //x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
1847 : //y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
1848 : // Rant = a->radius;
1849 : // Rhole = a->hole_radius;
1850 : //R2 = Rant*Rant;
1851 : //H2 = Rhole*Rhole;
1852 :
1853 : //eps = dx/4.0;
1854 :
1855 0 : IPosition pos(4);
1856 : // shape = ap->aperture->shape();
1857 :
1858 :
1859 : // cerr << "max threads " << omp_get_max_threads()
1860 : // << " threads available " << omp_get_num_threads()
1861 : // << endl;
1862 : // Int Nth=1, localWhichPoln=whichPoln;
1863 0 : Int localWhichPoln=whichPoln;
1864 : //#ifdef _OPENMP
1865 : //Nth=max(omp_get_max_threads()-2,1);
1866 : //#endif
1867 : // Timer tim;
1868 : // tim.mark();
1869 : //#if GCC44x
1870 : //#pragma omp parallel default(none) firstprivate(Er, El, nx, ny, localWhichPoln) private(i,j) shared(ap, a, p, L0) num_threads(Nth)
1871 : //#else
1872 : //#pragma omp parallel default(none) firstprivate(Er, El, nx, ny, localWhichPoln) private(i,j) shared(ap, a, p, L0) num_threads(Nth)
1873 : //#endif
1874 : {
1875 : //#pragma omp for
1876 0 : for(j = 0; j < ny; j++)
1877 : {
1878 0 : for(i = 0; i < nx; i++)
1879 : {
1880 0 : computePixelValues(ap, a, p, L0, Er, El, i,j,localWhichPoln);
1881 : }
1882 : }
1883 : }
1884 : // tim.show("BeamCalc:");
1885 :
1886 0 : deletePathology(p);
1887 0 : deleteAntenna(a);
1888 :
1889 0 : return 1;
1890 0 : }
1891 :
1892 0 : void BeamCalc::computePixelValues(const ApertureCalcParams *ap,
1893 : const calcAntenna *a, const Pathology *p,
1894 : const Double &L0,
1895 : Complex *Er, Complex *El,
1896 : const Int &i, const Int &j,
1897 : const Int& whichPoln)
1898 : {
1899 0 : Complex fp, Exr, Eyr, Exl, Eyl;
1900 : // Complex Er[3], El[3];
1901 0 : Complex E1[3];
1902 : Double dx, dy, x0, y0, x, y, r2, Rhole, Rant, R2, H2, eps;
1903 0 : Complex rr, rl, lr, ll, tmp;
1904 : Double L, amp, dP, dA, dO, x1, y1, dx1, dy1, dx2, dy2, phase;
1905 : //Int nx, ny;
1906 : Int os;
1907 0 : Int niter=6;
1908 : Double pac, pas, cp,sp; /* parallactic angle cosine / sine */
1909 0 : Complex Iota; Iota=Complex(0,1);
1910 0 : IPosition pos(4);pos=0;
1911 :
1912 0 : Ray *ray=0, *rayx=0, *rayy=0;
1913 : /* determine parallactic angle rotated coordinates */
1914 :
1915 0 : os = ap->oversamp;
1916 : //nx = ap->nx*os;
1917 : //ny = ap->ny*os;
1918 0 : dx = ap->dx/os;
1919 0 : dy = ap->dy/os;
1920 0 : x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
1921 0 : y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
1922 0 : Rant = a->radius;
1923 0 : Rhole = a->hole_radius;
1924 0 : R2 = Rant*Rant;
1925 0 : H2 = Rhole*Rhole;
1926 : // for(Int i=0; i < nx; ++i)
1927 : {
1928 0 : eps = dx/4.0;
1929 0 : pac = cos(ap->pa+M_PI/2);
1930 0 : pas = sin(ap->pa+M_PI/2);
1931 :
1932 0 : x = pac*(x0 + i*dx) - pas*(y0 + j*dy);
1933 0 : y = pas*(x0 + i*dx) + pac*(y0 + j*dy);
1934 0 : x = -x;
1935 :
1936 0 : if(fabs(x) > Rant) goto nextpoint;
1937 0 : if(fabs(y) > Rant) goto nextpoint;
1938 0 : r2 = x*x + y*y;
1939 0 : if(r2 > R2) goto nextpoint;
1940 0 : if(r2 < H2) goto nextpoint;
1941 :
1942 0 : ray = rayx = rayy = 0;
1943 :
1944 0 : x1 = x;
1945 0 : y1 = y;
1946 :
1947 0 : for(Int iter = 0; iter < niter; iter++)
1948 : {
1949 0 : ray = trace(a, x1, y1, p);
1950 0 : if(!ray) goto nextpoint;
1951 0 : x1 += (x - ray->aper[0]);
1952 0 : y1 += (y - ray->aper[1]);
1953 0 : deleteRay(ray);
1954 0 : ray = 0;
1955 : }
1956 :
1957 0 : ray = trace(a, x1, y1, p);
1958 :
1959 : /* check for leg blockage */
1960 0 : if(legplanewaveblock2(a, ray))
1961 0 : goto nextpoint;
1962 0 : if(legsphericalwaveblock(a, ray))
1963 0 : goto nextpoint;
1964 :
1965 0 : if(y < 0) rayy = trace(a, x1, y1+eps, p);
1966 0 : else rayy = trace(a, x1, y1-eps, p);
1967 :
1968 0 : if(x < 0) rayx = trace(a, x1+eps, y1, p);
1969 0 : else rayx = trace(a, x1-eps, y1, p);
1970 :
1971 0 : if(ray == 0 || rayx == 0 || rayy == 0)
1972 0 : goto nextpoint;
1973 :
1974 : /* compute solid angle subtended at the feed */
1975 0 : dx1 = rayx->aper[0]-ray->aper[0];
1976 0 : dy1 = rayx->aper[1]-ray->aper[1];
1977 0 : dx2 = rayy->aper[0]-ray->aper[0];
1978 0 : dy2 = rayy->aper[1]-ray->aper[1];
1979 :
1980 0 : dA = 0.5*fabs(dx1*dy2 - dx2*dy1);
1981 0 : dO = (dOmega(a, rayx, rayy, ray, p)/dA)*dx*dx;
1982 0 : dP = dO*feedgain(a, ray, p);
1983 0 : amp = sqrt(dP);
1984 :
1985 0 : L = Raylen(ray);
1986 :
1987 0 : phase = 2.0*M_PI*(L-L0)/a->lambda;
1988 :
1989 : /* phase retard the wave */
1990 0 : cp = cos(phase);
1991 0 : sp = sin(phase);
1992 : // fp = cp + sp*1.0i;
1993 :
1994 0 : fp = Complex(cp,sp);
1995 :
1996 :
1997 0 : tracepol(Er, ray, E1);
1998 0 : Exr = fp*amp*E1[0];
1999 0 : Eyr = fp*amp*E1[1];
2000 : // rr = Exr - 1.0i*Eyr;
2001 : // rl = Exr + 1.0i*Eyr;
2002 0 : rr = Exr - Iota*Eyr;
2003 0 : rl = Exr + Iota*Eyr;
2004 :
2005 0 : tracepol(El, ray, E1);
2006 0 : Exl = fp*amp*E1[0];
2007 0 : Eyl = fp*amp*E1[1];
2008 : // lr = Exl - 1.0i*Eyl;
2009 : // ll = Exl + 1.0i*Eyl;
2010 0 : lr = Exl - Iota*Eyl;
2011 0 : ll = Exl + Iota*Eyl;
2012 :
2013 : // pos(0)=(Int)((j/os) - (25.0/dy/os)/2 + shape(0)/2 - 0.5);
2014 : // pos(1)=(Int)((i/os) - (25.0/dx/os)/2 + shape(1)/2 - 0.5);
2015 : // Following 3 lines go with ANT tag in VLACalc.....
2016 : // Double antRadius=BeamCalcGeometryDefaults[ap->band].Rant;
2017 : // pos(0)=(Int)((j/os) - (antRadius/dy/os) + shape(0)/2 - 0.5);
2018 : // pos(1)=(Int)((i/os) - (antRadius/dx/os) + shape(1)/2 - 0.5);
2019 : // Following 2 lines go with the PIX tag in VLACalc...
2020 0 : pos(0)=(Int)((j/os));
2021 0 : pos(1)=(Int)((i/os));
2022 0 : pos(2)=0;
2023 0 : pos(3)=0;
2024 :
2025 0 : if (whichPoln==Stokes::RR)
2026 0 : {tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+rr,pos);}
2027 0 : else if (whichPoln==Stokes::RL)
2028 0 : {tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+rl,pos);}
2029 0 : else if (whichPoln==Stokes::LR)
2030 0 : {tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+lr,pos);}
2031 0 : else if (whichPoln==Stokes::LL)
2032 0 : {tmp=ap->aperture->getAt(pos);ap->aperture->putAt(tmp+ll,pos);}
2033 : else {
2034 0 : SynthesisError err(String("BeamCalc::computePixelValues: Cannot handle Stokes ")+String(whichPoln));
2035 0 : throw(err);
2036 0 : }
2037 :
2038 0 : nextpoint:
2039 0 : if(ray) deleteRay(ray);
2040 0 : if(rayx) deleteRay(rayx);
2041 0 : if(rayy) deleteRay(rayy);
2042 : }
2043 0 : }
2044 :
2045 : //
2046 : //----------------------------------------------------------------------
2047 : // Compute only the required polarizations.for linear polarization
2048 : //
2049 0 : Int BeamCalc::calculateApertureLinPol(ApertureCalcParams *ap, const Int& whichPoln)
2050 : {
2051 0 : Complex Ex[3], Ey[3];
2052 0 : Complex Px[2], Py[2];
2053 : //Double dx, dy, Rhole, Rant;//x0, y0, R2, H2, eps;
2054 : Double L0;
2055 : calcAntenna *a;
2056 : Pathology *p;
2057 : Int nx, ny, os;
2058 : Int i, j;
2059 : //Double pac, pas; /* parallactic angle cosine / sine */
2060 : // Complex Iota=Complex(0,1);
2061 :
2062 :
2063 0 : a = newAntennafromApertureCalcParams(ap);
2064 0 : p = newPathologyfromApertureCalcParams(ap);
2065 :
2066 : /* compute central ray pathlength */
2067 : {
2068 : Ray *tmpRay;
2069 0 : tmpRay = trace(a, 0.0, 0.00001, p);
2070 0 : L0 = Raylen(tmpRay);
2071 0 : deleteRay(tmpRay);
2072 : }
2073 :
2074 : //pac = cos(ap->pa+M_PI/2);
2075 : //pas = sin(ap->pa+M_PI/2);
2076 :
2077 : /* in linear basis */
2078 0 : Px[0] = 0.0; Px[1] = 1.0;
2079 0 : Py[0] = 1.0; Py[1] = 0.0;
2080 :
2081 0 : IPosition pos(4); pos=0;
2082 :
2083 : /* compute 3-vector feed efields for the two polarizations */
2084 0 : Efield(a, Py, Ey);
2085 0 : Efield(a, Px, Ex);
2086 :
2087 0 : if (whichPoln == Stokes::XX){
2088 0 : pos(2)=0;
2089 : }
2090 0 : else if (whichPoln == Stokes::YY){
2091 0 : pos(2)=3;
2092 : }
2093 0 : else if (whichPoln == Stokes::XY){
2094 0 : pos(2)=1;
2095 : }
2096 0 : else if (whichPoln == Stokes::YX){
2097 0 : pos(2)=2;
2098 : }
2099 : else {
2100 0 : SynthesisError err(String("BeamCalc::calculateApertureLinPol: Cannot handle Stokes ")+String(whichPoln));
2101 0 : throw(err);
2102 0 : }
2103 :
2104 : // set only the affected plane to zero
2105 0 : for(j = 0; j < ap->nx; j++){
2106 0 : pos(0)= j;
2107 0 : for(i = 0; i < ap->ny; i++){
2108 0 : pos(1)= i;
2109 0 : ap->aperture->putAt(Complex(0.),pos);
2110 : }
2111 : }
2112 :
2113 0 : os = ap->oversamp;
2114 0 : nx = ap->nx*os;
2115 0 : ny = ap->ny*os;
2116 : // dx = ap->dx/os;
2117 : // dy = ap->dy/os;
2118 : //x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
2119 : //y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
2120 : // Rant = a->radius;
2121 : // Rhole = a->hole_radius;
2122 : //R2 = Rant*Rant;
2123 : //H2 = Rhole*Rhole;
2124 :
2125 : //eps = dx/4.0;
2126 :
2127 : // cerr << "max threads " << omp_get_max_threads()
2128 : // << " threads available " << omp_get_num_threads()
2129 : // << endl;
2130 0 : Int Nth=1, localWhichPoln=whichPoln;
2131 : #ifdef _OPENMP
2132 0 : Nth=max(omp_get_max_threads()-2,1);
2133 : #endif
2134 : // Timer tim;
2135 : // tim.mark();
2136 : #if GCC44x
2137 : #pragma omp parallel default(none) firstprivate(Ex, Ey, nx, ny, localWhichPoln) private(i,j) shared(ap, a, p, L0) num_threads(Nth)
2138 : #else
2139 0 : #pragma omp parallel default(none) firstprivate(Ex, Ey, nx, ny, localWhichPoln) private(i,j) shared(ap, a, p, L0) num_threads(Nth)
2140 : #endif
2141 : {
2142 : #pragma omp for
2143 : for(j = 0; j < ny; j++)
2144 : {
2145 : for(i = 0; i < nx; i++)
2146 : {
2147 : computePixelValuesLinPol(ap, a, p, L0, Ex, Ey, i,j,localWhichPoln);
2148 : }
2149 : }
2150 : }
2151 : // tim.show("BeamCalc:");
2152 :
2153 0 : deletePathology(p);
2154 0 : deleteAntenna(a);
2155 :
2156 0 : return 1;
2157 0 : }
2158 :
2159 0 : void BeamCalc::computePixelValuesLinPol(const ApertureCalcParams *ap,
2160 : const calcAntenna *a, const Pathology *p,
2161 : const Double &L0,
2162 : Complex *Ex, Complex *Ey,
2163 : const Int &i, const Int &j,
2164 : const Int& whichPoln)
2165 : {
2166 0 : Complex fp, Exx, Eyx, Exy, Eyy;
2167 :
2168 0 : Complex E1[3];
2169 : Double dx, dy, x0, y0, x, y, r2, Rhole, Rant, R2, H2, eps;
2170 0 : Complex xx, xy, yx, yy, tmp;
2171 : Double L, amp, dP, dA, dO, x1, y1, dx1, dy1, dx2, dy2, phase;
2172 : //Int nx, ny;
2173 : Int os;
2174 0 : Int niter=6;
2175 : Double pac, pas, cp,sp; /* parallactic angle cosine / sine */
2176 0 : Complex Iota; Iota=Complex(0,1);
2177 0 : IPosition pos(4);pos=0;
2178 :
2179 0 : Ray *ray=0, *rayx=0, *rayy=0;
2180 : /* determine parallactic angle rotated coordinates */
2181 :
2182 0 : os = ap->oversamp;
2183 : //nx = ap->nx*os;
2184 : //ny = ap->ny*os;
2185 0 : dx = ap->dx/os;
2186 0 : dy = ap->dy/os;
2187 0 : x0 = ap->x0 - ap->dx/2.0 + dx/2.0;
2188 0 : y0 = ap->y0 - ap->dy/2.0 + dy/2.0;
2189 0 : Rant = a->radius;
2190 0 : Rhole = a->hole_radius;
2191 0 : R2 = Rant*Rant;
2192 0 : H2 = Rhole*Rhole;
2193 : // for(Int i=0; i < nx; ++i)
2194 : {
2195 0 : eps = dx/4.0;
2196 0 : pac = cos(ap->pa+M_PI/2);
2197 0 : pas = sin(ap->pa+M_PI/2);
2198 :
2199 0 : x = pac*(x0 + i*dx) - pas*(y0 + j*dy);
2200 0 : y = pas*(x0 + i*dx) + pac*(y0 + j*dy);
2201 0 : x = -x;
2202 :
2203 0 : if(fabs(x) > Rant) goto nextpoint;
2204 0 : if(fabs(y) > Rant) goto nextpoint;
2205 0 : r2 = x*x + y*y;
2206 0 : if(r2 > R2) goto nextpoint;
2207 0 : if(r2 < H2) goto nextpoint;
2208 :
2209 0 : ray = rayx = rayy = 0;
2210 :
2211 0 : x1 = x;
2212 0 : y1 = y;
2213 :
2214 0 : for(Int iter = 0; iter < niter; iter++)
2215 : {
2216 0 : ray = trace(a, x1, y1, p);
2217 0 : if(!ray) goto nextpoint;
2218 0 : x1 += (x - ray->aper[0]);
2219 0 : y1 += (y - ray->aper[1]);
2220 0 : deleteRay(ray);
2221 0 : ray = 0;
2222 : }
2223 :
2224 0 : ray = trace(a, x1, y1, p);
2225 :
2226 : /* check for leg blockage */
2227 0 : if(legplanewaveblock2(a, ray))
2228 0 : goto nextpoint;
2229 0 : if(legsphericalwaveblock(a, ray))
2230 0 : goto nextpoint;
2231 :
2232 0 : if(y < 0) rayy = trace(a, x1, y1+eps, p);
2233 0 : else rayy = trace(a, x1, y1-eps, p);
2234 :
2235 0 : if(x < 0) rayx = trace(a, x1+eps, y1, p);
2236 0 : else rayx = trace(a, x1-eps, y1, p);
2237 :
2238 0 : if(ray == 0 || rayx == 0 || rayy == 0)
2239 0 : goto nextpoint;
2240 :
2241 : /* compute solid angle subtended at the feed */
2242 0 : dx1 = rayx->aper[0]-ray->aper[0];
2243 0 : dy1 = rayx->aper[1]-ray->aper[1];
2244 0 : dx2 = rayy->aper[0]-ray->aper[0];
2245 0 : dy2 = rayy->aper[1]-ray->aper[1];
2246 :
2247 0 : dA = 0.5*fabs(dx1*dy2 - dx2*dy1);
2248 0 : dO = (dOmega(a, rayx, rayy, ray, p)/dA)*dx*dx;
2249 0 : dP = dO*feedgain(a, ray, p);
2250 0 : amp = sqrt(dP);
2251 :
2252 0 : L = Raylen(ray);
2253 :
2254 0 : phase = 2.0*M_PI*(L-L0)/a->lambda;
2255 :
2256 : /* phase retard the wave */
2257 0 : cp = cos(phase);
2258 0 : sp = sin(phase);
2259 :
2260 0 : fp = Complex(cp,sp);
2261 :
2262 :
2263 0 : tracepol(Ex, ray, E1);
2264 0 : Exx = fp*amp*E1[0];
2265 0 : Eyx = fp*amp*E1[1];
2266 :
2267 0 : tracepol(Ey, ray, E1);
2268 0 : Exy = fp*amp*E1[0];
2269 0 : Eyy = fp*amp*E1[1];
2270 :
2271 :
2272 0 : xx = Exx;
2273 0 : xy = Complex(0.);
2274 0 : yx = Complex(0.);
2275 0 : yy = Eyy;
2276 :
2277 : // pos(0)=(Int)((j/os) - (25.0/dy/os)/2 + shape(0)/2 - 0.5);
2278 : // pos(1)=(Int)((i/os) - (25.0/dx/os)/2 + shape(1)/2 - 0.5);
2279 : // Following 3 lines go with ANT tag in VLACalc.....
2280 : // Double antRadius=BeamCalcGeometryDefaults[ap->band].Rant;
2281 : // pos(0)=(Int)((j/os) - (antRadius/dy/os) + shape(0)/2 - 0.5);
2282 : // pos(1)=(Int)((i/os) - (antRadius/dx/os) + shape(1)/2 - 0.5);
2283 : // Following 2 lines go with the PIX tag in VLACalc...
2284 0 : pos(0)=(Int)((j/os));
2285 0 : pos(1)=(Int)((i/os));
2286 0 : pos(3)=0;
2287 :
2288 0 : if (whichPoln==Stokes::XX){
2289 0 : pos(2)=0;
2290 0 : tmp=ap->aperture->getAt(pos);
2291 0 : ap->aperture->putAt(tmp+xx,pos);
2292 : }
2293 :
2294 0 : else if (whichPoln==Stokes::XY){
2295 0 : pos(2)=1;
2296 0 : tmp=ap->aperture->getAt(pos);
2297 0 : ap->aperture->putAt(tmp+xy,pos);
2298 : }
2299 :
2300 0 : else if (whichPoln==Stokes::YX){
2301 0 : pos(2)=2;
2302 0 : tmp=ap->aperture->getAt(pos);
2303 0 : ap->aperture->putAt(tmp+yx,pos);
2304 : }
2305 :
2306 0 : else if (whichPoln==Stokes::YY){
2307 0 : pos(2)=3;
2308 0 : tmp=ap->aperture->getAt(pos);
2309 0 : ap->aperture->putAt(tmp+yy,pos);
2310 : }
2311 :
2312 0 : nextpoint:
2313 0 : if(ray) deleteRay(ray);
2314 0 : if(rayx) deleteRay(rayx);
2315 0 : if(rayy) deleteRay(rayy);
2316 : }
2317 0 : }
2318 : };
2319 :
|
