Line data Source code
1 : //# PBMath1DNumeric.cc: Implementation for PBMath1DNumeric
2 : //# Copyright (C) 1996,1997,1998,1999,2003
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 <casacore/casa/aips.h>
30 : #include <casacore/casa/BasicSL/Complex.h>
31 : #include <casacore/casa/BasicMath/Math.h>
32 : #include <casacore/scimath/Mathematics/FFTServer.h>
33 : #include <synthesis/TransformMachines/PBMath1DNumeric.h>
34 : #include <casacore/casa/Quanta.h>
35 : #include <casacore/measures/Measures.h>
36 :
37 :
38 :
39 : using namespace casacore;
40 : namespace casa { //# NAMESPACE CASA - BEGIN
41 :
42 0 : PBMath1DNumeric::PBMath1DNumeric(const Vector<Float>& numericArray,
43 : Quantity maxRad, Quantity refFreq,
44 : Bool isThisVP,
45 : BeamSquint squint,
46 0 : Bool useSymmetricBeam) :
47 : PBMath1D(maxRad, refFreq, isThisVP, squint, useSymmetricBeam),
48 0 : numericArray_p(numericArray)
49 : {
50 0 : fillPBArray();
51 :
52 0 : if (useSymmetricBeam) {
53 0 : symmetrizeSquintedBeam();
54 : }
55 0 : };
56 :
57 :
58 :
59 :
60 0 : PBMath1DNumeric::~PBMath1DNumeric()
61 : {
62 0 : };
63 :
64 :
65 0 : PBMath1DNumeric& PBMath1DNumeric::operator=(const PBMath1DNumeric& other)
66 : {
67 0 : if (this == &other)
68 0 : return *this;
69 :
70 0 : PBMath1D::operator=(other);
71 0 : numericArray_p = other.numericArray_p;
72 0 : return *this;
73 : };
74 :
75 :
76 :
77 :
78 0 : void PBMath1DNumeric::fillPBArray()
79 : {
80 :
81 0 : LogIO os(LogOrigin("PBMath1DNumeric", "fillPBArray"));
82 0 : Int nSamples=10000;
83 0 : vp_p.resize(nSamples);
84 0 : Vector<Float> vp_temp(2*nSamples);
85 :
86 0 : inverseIncrementRadius_p=Double(nSamples-1)/maximumRadius_p.getValue("'");
87 :
88 : // Sinc interpolation to fill up vp_p from numericArray_p:
89 : {
90 :
91 :
92 0 : Vector<Float> beam2(2*(numericArray_p.nelements()-1));
93 0 : beam2 = 0.0;
94 0 : uInt n2 = beam2.nelements();
95 : uInt i;
96 0 : for (i=0; i<numericArray_p.nelements()-1;i++) {
97 0 : beam2(n2/2+i) = numericArray_p(i);
98 0 : beam2(n2/2-i) = numericArray_p(i);
99 : }
100 0 : beam2(0) = numericArray_p(numericArray_p.nelements()-1);
101 :
102 0 : Vector<Complex> fft2;
103 :
104 :
105 0 : FFTServer<Float, Complex> server(IPosition(1, beam2.nelements()));
106 0 : server.fft0(fft2, beam2);
107 :
108 : // Now 0-pad fft2
109 0 : Vector<Complex> fft2pad(nSamples+1);
110 :
111 0 : for (i=0; i<fft2pad.nelements();i++) {
112 0 : fft2pad(i) = 0.0;
113 : }
114 :
115 0 : for (i=0; i<fft2.nelements();i++) {
116 0 : fft2pad(i) = fft2(i);
117 : }
118 :
119 : // and FT back into vp_temp
120 0 : server.fft0(vp_temp, fft2pad);
121 :
122 : uInt j;
123 : Float taper;
124 0 : for (i=vp_temp.nelements()/2,j=0; i< vp_temp.nelements();i++,j++) {
125 0 : taper = vp_temp(i)/vp_temp(vp_temp.nelements()/2);
126 0 : if (isThisVP_p) {
127 0 : vp_p(j) = taper;
128 : } else {
129 : // this will have sign problems! should do this via vp, not pb!
130 0 : if (taper <= 0.0) {
131 0 : vp_p(j) = 0.0;
132 : } else {
133 0 : vp_p(j) = sqrt(taper);
134 : }
135 : }
136 : }
137 :
138 0 : }
139 :
140 0 : };
141 :
142 :
143 : //Bool PBMath1DNumeric::flushToTable(Table& beamSubTable, Int iRow)
144 : //{};
145 :
146 :
147 : void
148 0 : PBMath1DNumeric::summary(Int nValues)
149 : {
150 0 : PBMath1D::summary(nValues);
151 0 : LogIO os(LogOrigin("PBMath1DNumeric", "summary"));
152 0 : for (uInt i=0; i < numericArray_p.nelements(); i++) {
153 0 : os << "Numeric values: " << i << " " << numericArray_p(i) << LogIO::POST;
154 : }
155 0 : };
156 :
157 0 : Int PBMath1DNumeric::support (const CoordinateSystem& cs){
158 :
159 0 : Int value=PBMath1D::support(cs);
160 : //As the maxRadius is determines where the numeric value stop...the support should be a bit larger
161 0 : value=Int(Double(value)*0.55)*2;
162 0 : return value;
163 :
164 :
165 : };
166 :
167 :
168 : } //# NAMESPACE CASA - END
169 :
|
