-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathacm2skyimage.m
42 lines (37 loc) · 1.32 KB
/
acm2skyimage.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
function skymap = acm2skyimage(acm, xpos, ypos, zpos, freq, l, m)
% skymap = acm2skyimag(acm, xpos, ypos, freq, l, m)
%
% conversion of ACM to sky image
%
% arguments
% acm : nelem x nelem x nchannel array correlation matrix
% xpos : x-position of antennas (vector of length nelem)
% ypos : y-position of antennas (vector of length nelem)
% freq : frequencies in Hz (vector of length nchannel)
% l, m : l and m points to which to direct the array (vectors)
%
% return value
% skymap : length(l) x length(m) x nchannel matrix containing the
% resulting sky maps
%
% SJW, 2004
%
% modified April 19, 2005 by SJW: correct coordinate conventions
% modified July 20, 2006 by SJW: optimization for release in LOFAR package
[nelem, dummy, nchannel] = size(acm);
c = 2.9979245e8;
skymap = zeros(length(l), length(m), nchannel);
for nch = 1:nchannel
disp(['acm2skyimage: working on channel ' num2str(nch) ' of ' num2str(nchannel)]);
lambda = c / freq(nch);
k = 2 * pi / lambda;
wx = exp(-i * k * xpos(:) * l(:).');
wy = exp(-i * k * ypos(:) * m(:).');
wz = exp(-i * k * zpos(:) * (sqrt (1 - l(:).^2 - m(:).^2)).');
for lidx = 1:length(l)
for midx = 1:length(m)
weight = wx(:, lidx) .* wy(:, midx); % .* wz(:,midx);
skymap(lidx, midx, nch) = real(weight' * (acm(:, :, nch)) * weight);
end
end
end