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