% Rodrigo Parra - July 2002
% Adaptated for ERR170 course on April 14 2005

% Fill the vector L and l with the L atitude and l ongitude of the
% stations.
%
% Define the Hour angle range on the vector H. Nt is the number
% of samples in the interval
%
% The matrices u and v have size Nt times Nb
% contain the uv tracks of each baseline.

% Parameters -----------------------------------------------
Nt = 100;               % Number of points for time vector
c = 3e8;		% m/s	( speed of light )
f = 1667e6;		% Hz	( frequency )
lambda = c / f ;	% m	( wavelength )
Re = 6378e3;		% m	( radius of the Earth )
% ----------------------------------------------------------


% if you want snapshot mode uncomment H = 0
% if you want uv tracks, specify the range
% for the hour angle (and uncomment of course!)
% The rendition of the dirty beam in tracking mode
% can take some seconds.

% Tracking mode -------------------
% H = linspace ( -6 , 6 , Nt ) * pi / 12;

% Snapshot mode -------------------
 H = 0;


% Declination of the source ------------
dec = 30;		% deg.
dec = dec * pi / 180;   % rad.
% --------------------------------------


% Station geographic coordinates
% uncomment the array you want to analyze
% The geometries correspond with those described
% in the Problem number 3. For the case C you can
% specify the number of stations N to use.

% Randomly distributed stations usually
% have very good uv coverage.

% NORTH latitude in degrees (L)
% EAST longitude in degrees (l)

% CASE A -------------------------------
% L = [0 0 0 0 0 0 0 1 2 3];	  
% l = [-1 -2 -3 0 1 2 3 0 0 0];   

% CASE B -------------------------------
 L = [ 0 1 2 3 ...			     	 
   -sin(pi/6) -2*sin(pi/6) -3*sin(pi/6)...   	 
   -sin(pi/6) -2*sin(pi/6) -3*sin(pi/6)]+30; 	 
 l = [ 0 0 0 0 ...			     	 
   -cos(pi/6) -2*cos(pi/6) -3*cos(pi/6)...   	 
   +cos(pi/6) +2*cos(pi/6) +3*cos(pi/6)];    	 

% CASE C -------------------------------
% N = 10;		    	
% theta = [0 : N-1]*(2*pi/N)	
% L = cos( theta );	    	
% l = sin( theta );	    	
	
% RANDOMLY DISTRIBUTED STATIONS --------
% N = 10;
% L = ( rand( 1, N ) - 0.5 ) * 180;
% l = ( rand( 1 ,N ) - 0.5 ) * 180;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DON'T TOUCH ANYTHING BELOW THIS LINE!!!
% UNLESS YOU WANT TO BE BURNED IN HELL!!!
% ------------------------------------------------

% Number of stations
N = length ( l );

% Conversion to cartesian coordinates
[ x , y , z ] = sph2cart( pi * l / 180 , pi * L / 180 , Re );

% Possible combinations of 2 stations ( baselines )
I = nchoosek ( [ 1 : N ] , 2 );

I = [I;fliplr(I)];


% Number of baselines
Nb = length ( I );

% Pairs of coordinates for all combinations
x1 = x( I ); y1 = y( I ); z1 = z( I );

% Cartesian components of the baselines
dx = diff ( x1' )';
dy = diff ( y1' )';
dz = diff ( z1' )';

% Length of the baselines
d = sqrt ( dx.^2 + dy.^2 +dz.^2 );

% uv coordinates calculation
u = 	+dx * sin ( H )...
	+dy * cos ( H );

v =	-dx * cos ( H ) * sin ( dec ) +...
	+dy * sin ( H ) * sin ( dec ) +...
	+dz * ones ( size ( H ) ) * cos ( dec );
	
if length( H ) == Nt
  u = reshape( u , [ Nb * Nt 1 ] );
  v = reshape( v , [ Nb * Nt 1 ] );
end


uv = round([ u v ]/12000);

[ ROW , i_unique ] = unique( uv , 'rows' );

i_redundant = setdiff( [1:Nb] , i_unique ); 

u_nr = u( i_unique ); 
v_nr = v( i_unique );

u_r = u( i_redundant ); 
v_r = v( i_redundant );


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
plot(l,L,'o')
text(l,L,num2str([1:N]'))
axis equal
xlabel('East longitude');
ylabel('North latitude');
title('Array configuration')
figure(2)
plot( u_nr'/1e6/lambda , v_nr'/1e6/lambda , 'x' );
if H == 0
  hold on
  plot( u_r'/1e6/lambda , v_r'/1e6/lambda , 'g+' )
  hold off
end
axis equal; grid
xlabel ('u - Mega-wavelength')
ylabel ('v - Mega-wavelength')
title('uv coverage')
figure(3)
dirty_beam(u_nr,v_nr);
box
camproj('perspective')
title('Dirty beam')
