Line data Source code
1 : //# IlluminationConvFunc.cc: Implementation for IlluminationConvFunc
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 : #define USETABLES 0
30 : #include <synthesis/TransformMachines2/IlluminationConvFunc.h>
31 :
32 : using namespace casacore;
33 : namespace casa{
34 : namespace refim{
35 : //
36 : //------------------------------------------------------------------------
37 : //
38 : /*
39 : void IlluminationConvFunc::getValue(CF_TYPE *coord,
40 : CF_TYPE *raoff1, CF_TYPE *raoff2,
41 : CF_TYPE *decoff1, CF_TYPE *decoff2,
42 : CF_TYPE *area,
43 : std::complex<CF_TYPE>& weight,
44 : std::complex<CF_TYPE>& dweight1,
45 : std::complex<CF_TYPE>& dweight2
46 : )
47 : */
48 0 : CF_TYPE IlluminationConvFunc::getValue(Double *coord,
49 : Double *raoff1, Double *raoff2,
50 : Double *decoff1, Double *decoff2,
51 : Double *area,
52 : Int *doGrad,
53 : Complex& weight,
54 : Complex& dweight1,
55 : Complex& dweight2,
56 : Double& /*currentCFPA*/
57 : // ,Double lsigma
58 : )
59 : {
60 0 : DComplex Kij;
61 0 : Double Rij,u,v,dra,ddec,s, MY_PI=M_PI;
62 : // Double cpa,spa;
63 : //
64 : // If the PA of the convolution function changed, set up a new
65 : // rotation matrix
66 : //
67 : /*
68 : if (pa_p != currentCFPA)
69 : {
70 : pa_p = currentCFPA;
71 : cpa = cos(pa_p);
72 : spa = sin(pa_p);
73 : }
74 : */
75 : //
76 : // Compute the PA rotated (u,v)
77 : //
78 : // u = cpa*coord[0] - spa*coord[1];
79 : // v = spa*coord[0] + cpa*coord[1];
80 0 : u = coord[0];
81 0 : v = coord[1];
82 0 : s=4*sigma;
83 0 : dra=(*raoff1- *raoff2);
84 0 : ddec=(*decoff1 - *decoff2);
85 0 : Rij = (dra*dra + ddec*ddec)*s/2.0;
86 : // Kij = Complex(0,M_PI*(u*(*raoff1+*raoff2) + v*(*decoff1+*decoff2)));
87 0 : Kij = Complex(0,MY_PI*(u*(*raoff1+*raoff2) + v*(*decoff1+*decoff2)));
88 : // Eij0 = (u*u + v*v)*M_PI*M_PI/s;
89 :
90 : // Rij = (dra*dra + ddec*ddec);
91 0 : Rij = 0;
92 : #if(USETABLES)
93 : //
94 : // Compute it using lookup tables for exponentials.
95 : //
96 : Double tt=imag(Kij);
97 :
98 : // weight = ExpTable((-Eij0-Rij))*CExpTable(tt)/ *area;
99 : weight = ExpTable((-Rij))*CExpTable(tt)/ *area;
100 : #else
101 : //
102 : // Compute it exactly
103 : //
104 : // weight = (exp((-Eij0-Rij)+Kij)/ *area);
105 0 : weight = (exp((-Rij)+Kij)/ *area);
106 : // cout << exp(-Rij) << " " << exp(Kij) << " " << dra << " " << ddec << endl;
107 : #endif
108 :
109 0 : if (doGrad)
110 : {
111 : //
112 : // Following 2 lines when Rij is not zero
113 : //
114 : // dweight1 = (Complex(0,M_PI*u)-Complex(dra* s,0.0))*weight;
115 : // dweight2 = (Complex(0,M_PI*v)-Complex(ddec* s,0.0))*weight;
116 : //
117 : // dweight1 = -Complex(2*ddec,0) + (Complex(0,M_PI*u));//*weight;
118 : // dweight2 = -Complex(2*dra,0) + (Complex(0,M_PI*v));//*weight;
119 : //
120 : // Following 2 lines when Rij=0
121 : //
122 0 : dweight1 = (Complex(0,MY_PI*u));//*weight;
123 0 : dweight2 = (Complex(0,MY_PI*v));//*weight;
124 : }
125 : // cout << "Eij: " << weight << " " << dweight1 << " " << dweight2 << endl;
126 0 : return 1.0;
127 : }
128 : //
129 : //------------------------------------------------------------------------
130 : //
131 0 : CF_TYPE IlluminationConvFunc::area(Vector<Int>& convSupport, Vector<Double>& uvScale)
132 : {
133 0 : CF_TYPE u,v, area=0;
134 :
135 0 : for(Int i=-convSupport(0);i<convSupport(0);i++)
136 : {
137 0 : u = i*M_PI/uvScale(0);
138 0 : u *= u;
139 0 : for(Int j=-convSupport(0);j<convSupport(0);j++)
140 : {
141 0 : v = j*M_PI/uvScale(1);
142 0 : v *= v;
143 : #ifdef USETABLES
144 0 : area += ExpTable(-(u+v)/(4*sigma));
145 : #else
146 : area += exp(-(u+v)/(4*sigma));
147 : #endif
148 : }
149 : }
150 0 : return area;
151 : };
152 : //
153 : //------------------------------------------------------------------------
154 : //
155 0 : Vector<Int> IlluminationConvFunc::supportSize(Vector<Double>& uvScale)
156 : {
157 : CF_TYPE u;
158 0 : Vector<Int> supportSize(uvScale.shape());
159 :
160 0 : for(Int i=0;i<100;i++)
161 : {
162 0 : u=i*M_PI/uvScale(0);
163 0 : u=u*u/(4*sigma);
164 0 : u=exp(-u);
165 0 : if (u < 1e-5)
166 : {
167 0 : supportSize(0) = i;
168 0 : supportSize(0)++; // Increment it by one since this is
169 : // used in FORTRAN indexing
170 0 : break;
171 : }
172 : }
173 0 : return supportSize;
174 : };
175 : };
176 : };
|
