Line data Source code
1 : *=======================================================================
2 : * -*- FORTRAN -*-
3 : * Copyright (C) 1999,2001,2002
4 : * Associated Universities, Inc. Washington DC, USA.
5 : *
6 : * This library is free software; you can redistribute it and/or
7 : * modify it under the terms of the GNU Library General Public
8 : * License as published by the Free Software Foundation; either
9 : * version 2 of the License, or (at your option) any later version.
10 : *
11 : * This library is distributed in the hope that it will be useful,
12 : * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 : * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 : * GNU Library General Public License for more details.
15 : *
16 : * You should have received a copy of the GNU Library General Public
17 : * License along with this library; if not, write to the Free
18 : * Software Foundation, Inc., 675 Massachusetts Ave, Cambridge,
19 : * MA 02139, USA.
20 : *
21 : * Correspondence concerning AIPS++ should be addressed as follows:
22 : * Internet email: aips2-request@nrao.edu.
23 : * Postal address: AIPS++ Project Office
24 : * National Radio Astronomy Observatory
25 : * 520 Edgemont Road
26 : * Charlottesville, VA 22903-2475 USA
27 : *
28 : * $Id:$
29 : *-----------------------------------------------------------------------
30 : C The LeakyAntsol algorithm (A&A, v.375, p.344-350 (2001)).
31 : C
32 : C Modified calantso.FullArray.f file to include poln leakage terms.
33 : C
34 : C S.Bhatnagar 22 Sept. 2000
35 : C
36 : C Gain is set to a normal value when ByPass is 1 and to a smaller
37 : C value otherwise. This implies that for normal antsol, the gain
38 : C (and hence the search step size) is also normal and for
39 : C leakyantsol, the gain is much smaller (finer step size).
40 : C
41 : C S.Bhatnagar 30 Oct. 2001
42 : C
43 0 : subroutine leakyantso (corgain, corwt, nant, antgain, leakage,
44 : $ mode, oresid, nresid, scanwt,refant,stat,ByPass)
45 : C
46 : integer maxnant
47 : character*(*) routine
48 : parameter (routine = 'calantso')
49 : parameter (maxnant = 120)
50 :
51 : integer nant,mode,stat
52 : complex corgain(nant,*), antgain(*), leakage(*)
53 : real oresid, nresid, scanwt
54 : real corwt(nant, *)
55 : C
56 : integer iter, niter, iant, jant, ByPass
57 : integer refant
58 : C
59 : C Number of logical antennas that this routine can handle.
60 : C
61 : real eps, antwt, presid, gain
62 : complex gtop(maxnant),top(maxnant)
63 : real bottom(maxnant)
64 : data eps /1e-10/
65 : data niter /30000/
66 : data gain /0.001/
67 : C=====================================================================
68 : C
69 : C initialise gain only if the gain on entry is zero
70 : C
71 0 : oresid = 0.0
72 0 : nresid = 0.0
73 0 : scanwt = 0.0
74 :
75 : C
76 : C set initial estimates of gains
77 : C
78 : c$$$ do iant=1,nant
79 : c$$$ write(*,*) 'w=',(corwt(iant,jant),jant=1,nant)
80 : c$$$ enddo
81 0 : do iant = 1, nant
82 0 : antwt = 0.0
83 0 : do jant = 1, nant
84 0 : antwt = antwt + corwt(iant, jant)
85 : enddo
86 0 : if (antwt.gt.0.0) then
87 0 : scanwt = scanwt + antwt
88 0 : if (abs(antgain(iant)).eq.0.0) then
89 0 : do jant = 1, nant
90 : antgain(iant) = antgain(iant) + corgain(iant, jant) *
91 0 : 1 corwt(iant, jant)
92 : enddo
93 0 : antgain(iant) = antgain(iant) / antwt
94 : end if
95 : else
96 0 : antgain(iant) = cmplx(0.0,0.0)
97 : end if
98 : enddo
99 :
100 : C
101 : C Set initial leakage terms
102 : C
103 0 : if (ByPass .eq. 1) then
104 0 : gain=0.1
105 0 : goto 1212
106 : else
107 0 : gain=0.01
108 : endif
109 : c$$$ do iant = 1, nant
110 : c$$$ antwt=0.0
111 : c$$$ top(1)=0.0
112 : c$$$ top(2)=0.0
113 : c$$$ do jant = 1, nant
114 : c$$$ antwt = antwt + corwt(iant, jant)
115 : c$$$ enddo
116 : c$$$ if (abs(leakage(iant)) .eq. 0.0) then
117 : c$$$ do jant = 1, nant
118 : c$$$ if (iant .ne. jant) then
119 : c$$$ top(1)=top(1)+corgain(iant,jant)*corwt(iant,jant)
120 : c$$$ top(2)=top(2)+conjg(antgain(jant))*corwt(iant,jant)
121 : c$$$ endif
122 : c$$$ enddo
123 : c$$$ leakage(iant)=(top(1)-antgain(iant)*top(2))/antwt
124 : c$$$ endif
125 : c$$$ enddo
126 :
127 0 : do iant = 1, nant
128 0 : antwt = 0.0
129 0 : do jant = 1, nant
130 0 : antwt = antwt + corwt(iant, jant)
131 : enddo
132 0 : if (antwt .gt. 0.0) then
133 0 : if (abs(leakage(iant)) .eq. 0.0) then
134 0 : do jant = 1, nant
135 0 : if (iant .ne. jant) then
136 : leakage(iant) = leakage(iant) +
137 : 1 (corgain(iant,jant)-
138 : 1 antgain(iant)*conjg(antgain(jant)))*
139 0 : 1 corwt(iant,jant)
140 : endif
141 : enddo
142 0 : leakage(iant)=leakage(iant)/antwt
143 : endif
144 : else
145 0 : leakage(iant)=cmplx(0.0,0.0)
146 : endif
147 : enddo
148 : 1212 continue
149 : c
150 : c find antgains w.r.t. to the ref. ant.
151 : c
152 0 : if (scanwt.eq.0.0) then
153 0 : write (0,*)routine,': no valid data'
154 : end if
155 : C
156 : C Default mode is Amp-Phase solution.
157 : C
158 0 : if (mode.eq.0) then
159 : C
160 : C Phase only solution. Put antgain.Amp = 1.0
161 : C
162 0 : do iant = 1, nant
163 0 : if(abs(antgain(iant)).ne.0.0) then
164 0 : antgain(iant) = antgain(iant)/abs(antgain(iant))
165 : end if
166 0 : if(abs(leakage(iant)).ne.0.0) then
167 0 : leakage(iant) = leakage(iant)/abs(leakage(iant))
168 : end if
169 : enddo
170 0 : elseif (mode.eq.1) then
171 : C
172 : C Amp only solution. Put antgain.Phase = 0.0
173 : C
174 0 : do iant = 1, nant
175 0 : antgain(iant) = abs(antgain(iant))
176 0 : leakage(iant) = abs(leakage(iant))
177 : enddo
178 : end if
179 : c
180 : c find initial rms
181 : c
182 0 : presid = 0.0
183 0 : antwt=0.0
184 0 : do iant = 1, nant
185 0 : do jant = 1, nant
186 : oresid=oresid+abs(antgain(iant)*conjg(antgain(jant)) +
187 : 1 leakage(iant)*conjg(leakage(jant)) -
188 0 : 1 corgain(iant,jant))**2 * corwt(iant,jant)
189 0 : antwt = antwt + corwt(iant, jant)
190 : enddo
191 : enddo
192 0 : oresid=oresid/antwt
193 : c write(*,*)"Initial RMS=",oresid
194 : c
195 : c oresid = sqrt(oresid/scanwt)
196 : c oresid = 0.0001
197 0 : if (oresid .lt. eps) goto 110
198 : c
199 0 : nresid = 0.0
200 0 : do 100 iter = 1, niter
201 : c if (mod(iter,10) .eq. 0) write(*,*) iter
202 :
203 : C
204 : C
205 : C----------------------Gain terms--------------------------
206 : C
207 : C
208 0 : do iant = 1, nant
209 0 : top(iant) = cmplx (0.0, 0.0)
210 0 : gtop(iant)=cmplx(0,0)
211 0 : bottom(iant) = 0.0
212 0 : do jant=1,nant
213 0 : if (iant .ne. jant) then
214 : top(iant) = top(iant)
215 : 1 +corwt(iant, jant) * antgain(jant) *
216 0 : 1 corgain(iant, jant)
217 : gtop(iant)=gtop(iant)+conjg(leakage(jant))*
218 0 : 1 antgain(jant)*corwt(iant,jant)
219 : bottom(iant) = bottom(iant) + abs(antgain(jant))**2 *
220 0 : 1 corwt(iant, jant)
221 : endif
222 : enddo
223 0 : gtop(iant) = gtop(iant)*leakage(iant)
224 : enddo
225 : c
226 : c find new antenna gain
227 : c
228 0 : do iant = 1, nant
229 0 : if (bottom(iant).ne.0.0) then
230 : antgain(iant) = (1.0-gain) * antgain(iant) +
231 0 : 1 gain * (top(iant)-gtop(iant)) / bottom(iant)
232 : end if
233 : enddo
234 : C
235 : C
236 : C----------------------Leakage terms--------------------------
237 : C
238 : C
239 0 : if (ByPass .eq. 1) goto 2121
240 0 : do iant = 1, nant
241 0 : top(iant) = cmplx (0.0, 0.0)
242 0 : gtop(iant)=cmplx(0,0)
243 0 : bottom(iant) = 0.0
244 0 : do jant=1,nant
245 0 : if (iant .ne. jant) then
246 : top(iant) = top(iant)
247 : 1 +corwt(iant, jant) * leakage(jant) *
248 0 : 1 corgain(iant, jant)
249 : gtop(iant)=gtop(iant)+conjg(antgain(jant))*
250 0 : 1 leakage(jant)*corwt(iant,jant)
251 : bottom(iant) = bottom(iant) + abs(leakage(jant))**2 *
252 0 : 1 corwt(iant, jant)
253 : endif
254 : enddo
255 0 : gtop(iant) = gtop(iant)*antgain(iant)
256 : enddo
257 : c
258 : c find new leakage gain
259 : c
260 0 : do 301 iant = 1, nant
261 0 : if (bottom(iant).ne.0.0) then
262 : leakage(iant) = (1.0-gain) * leakage(iant) +
263 0 : 1 gain * (top(iant)-gtop(iant)) / bottom(iant)
264 : end if
265 0 : 301 continue
266 : 2121 continue
267 : C
268 : C
269 : C-------------------------------------------------------------
270 : C
271 : C
272 : c
273 : c is this a phase-only solution?
274 : c
275 0 : if (mode.eq.0) then
276 0 : do iant = 1, nant
277 0 : if(abs(antgain(iant)).ne.0.0) then
278 0 : antgain(iant) = antgain(iant)/abs(antgain(iant))
279 : end if
280 0 : if(abs(leakage(iant)).ne.0.0) then
281 0 : leakage(iant) = leakage(iant)/abs(leakage(iant))
282 : end if
283 : enddo
284 0 : elseif (mode.eq.1) then
285 0 : do iant = 1, nant
286 0 : antgain(iant) = abs(antgain(iant))
287 0 : leakage(iant) = abs(leakage(iant))
288 : enddo
289 : end if
290 : c
291 : c find residual
292 : c
293 0 : presid = nresid
294 0 : nresid = 0.0
295 0 : antwt = 0.0
296 0 : if (mode.eq.1)then
297 0 : do jant = 1, nant
298 0 : do iant = 1, nant
299 : nresid = nresid + corwt(iant, jant) *
300 : 1 abs(antgain(iant)*antgain(jant) +
301 : 2 leakage(iant)*leakage(jant) -
302 0 : 3 abs(corgain(iant, jant)))**2
303 0 : antwt = antwt + corwt(iant, jant)
304 : enddo
305 : enddo
306 : else
307 0 : do jant = 1, nant
308 0 : do iant = 1, nant
309 : nresid = nresid + corwt(iant, jant) *
310 : 1 abs(antgain(iant)*conjg(antgain(jant)) +
311 : 2 leakage(iant)*conjg(leakage(jant)) -
312 0 : 3 corgain(iant, jant))**2
313 0 : antwt = antwt + corwt(iant, jant)
314 : enddo
315 : enddo
316 : endif
317 0 : nresid = sqrt(nresid/antwt)
318 : c
319 0 : if (abs(nresid-presid).le.eps*oresid) goto 110
320 : c
321 0 : 100 continue
322 : c
323 : 110 continue
324 0 : if (iter.ge.niter) then
325 0 : stat = -iter
326 : else
327 0 : stat = iter
328 : endif
329 : c write(*,*)'NResid-PResid=',nresid,presid,abs(nresid-presid)
330 : C
331 : C Reference the phases with the reference antenna.
332 : C
333 0 : if (abs(antgain(refant)) .ne. 0.0) then
334 0 : top(1) = conjg(antgain(refant))/abs(antgain(refant))
335 : else
336 0 : write(0,*) "###Error: Reference antenna is dead for gains!"
337 0 : top(1)=cmplx(1,0)
338 : endif
339 :
340 0 : if (ByPass .eq. 0) then
341 0 : if (abs(leakage(refant)) .ne. 0.0) then
342 : c top(2) = conjg(leakage(refant))/abs(leakage(refant))
343 0 : top(2) = leakage(refant)
344 : else
345 0 : write(0,*) "###Error: Reference antenna is dead for leakages!"
346 0 : top(2)=cmplx(1,0)
347 : endif
348 : endif
349 :
350 0 : do iant=1,nant
351 0 : antgain(iant) = antgain(iant)*top(1)
352 0 : if (ByPass .eq. 0) leakage(iant) = leakage(iant)-top(2)
353 : enddo
354 : C
355 : C Make Antgain and Leakage orthogonal
356 : C
357 : c$$$ if (ByPass .ne. 1) then
358 : c$$$ top(1)=cmplx(0.0,0.0)
359 : c$$$ do iant = 1, nant
360 : c$$$ top(1)=top(1)+conjg(antgain(iant))*leakage(iant)
361 : c$$$ enddo
362 : c$$$
363 : c$$$ top(2)=(top(1))/abs(top(1))
364 : c$$$
365 : c$$$ bottom(1)=cmplx(0.0,0.0)
366 : c$$$ bottom(2)=cmplx(0.0,0.0)
367 : c$$$ do iant = 1, nant
368 : c$$$c bottom(1)=bottom(1)+abs(antgain(iant))**2
369 : c$$$c bottom(2)=bottom(1)+abs(leakage(iant))**2
370 : c$$$ bottom(1)=bottom(1)+antgain(iant)*conjg(antgain(iant))
371 : c$$$ bottom(2)=bottom(1)+leakage(iant)*conjg(leakage(iant))
372 : c$$$C antgain(iant)=antgain(iant)*top(2)
373 : c$$$ enddo
374 : c$$$
375 : c$$$ antwt=atan(2.0*abs(top(1))/(bottom(1)-bottom(2)))/2.0
376 : c$$$
377 : c$$$ do iant = 1, nant
378 : c$$$ top(1)= antgain(iant)*cos(antwt)+leakage(iant)*sin(antwt)
379 : c$$$ top(2)=-antgain(iant)*sin(antwt)+leakage(iant)*cos(antwt)
380 : c$$$ antgain(iant)=top(1)
381 : c$$$ leakage(iant)=top(2)
382 : c$$$ enddo
383 : c$$$
384 : c$$$ top(1)=cmplx(0,0)
385 : c$$$ do iant=1,nant
386 : c$$$ top(1)=conjg(antgain(iant))*leakage(iant)
387 : c$$$ enddo
388 : c$$$c$$$ write(*,*)'Gdag.Alpha=',abs(top(1)),
389 : c$$$c$$$ 1 57.295*atan2(imagpart(top(1)),realpart(top(1)))
390 : c$$$ endif
391 : 999 continue
392 0 : end
|
