% $Id: doMp.m,v 1.2 2005/04/18 19:46:46 vnadkar Exp $ 
% 
% Load in the data and sort through to find just the good data
% i.e. remove the NaN's. Need to have first run doMpMap.bat on
% a data set in the current directory.
%
% Note this has been done in several stages with clears to reduce the
% amount of memory needed (a 24 hour session produces very large
% files!).
%
% This is a work in progress.
%
% Stuart Riley 4/23/01
% Copyright Trimble Navigation
%


load elev.out
e = elev(:,2:33);
clear elev;
i = find( isfinite(e));
e_good = e(i);
clear e;

load azi.out
a = azi(:,2:33);
clear azi
a_good = a(i);
clear a;

load mp1.out
m1 = mp1(:,2:33);
clear mp1
m1_good = m1(i);
clear m1;

load mp2.out
m2 = mp2(:,2:33);
clear mp2
m2_good = m2(i);
clear m2;

load sn1.out
s1 = sn1(:,2:33);
clear sn1
sn1_good = s1(i);
clear s1;

load sn2.out
s2 = sn2(:,2:33);
clear sn2
sn2_good = s2(i);
clear s2;

% 
% Now form elevation azimuth bins
%

k=1;
all=[e_good,a_good,m1_good,m2_good,sn1_good,sn2_good];

% 
% There's got to be a better way to use find here to do everything in
% one step..... need to think about it some more....
for el=10:10:90
  m=1;
  for az=-180:10:180
    i = find(   (all(:,1) >= el ) ...
              & (all(:,1) < el+5) ...
              & (all(:,2) >= az ) ...
              & (all(:,2) < az+5) );
    if(length(i) > 0)

      % MP1
      tmp = all(i,3);
      % Because of missing epochs there may be some
      % NaN's which need removing otherwise the sigma
      % will be NaN
      good = find (isfinite(tmp));
      if(length(good) > 0)
        mp_out(k,m,1) = std(tmp(good));
        mp_out(k,m,4) = length(tmp(good));

        % L1 CNR
        tmp = all(i,5);
        good = find (isfinite(tmp));
        mp_out(k,m,7) = mean(tmp(good));
      else
        mp_out(k,m,1) = NaN;
        mp_out(k,m,4) = NaN;
        mp_out(k,m,7) = NaN;
      end
      
      % MP2
      tmp = all(i,4);
      % Because of missing epochs there may be some
      % NaN's which need removing otherwise the sigma
      % will be NaN
      good = find (isfinite(tmp));
      if(length(good) > 0)
        mp_out(k,m,5) = std(tmp(good));
        mp_out(k,m,6) = length(tmp(good));
        
        % L2 CNR
        tmp = all(i,6);
        good = find (isfinite(tmp));
        mp_out(k,m,8) = mean(tmp(good));
      else
        mp_out(k,m,5) = NaN;
        mp_out(k,m,6) = NaN;
        mp_out(k,m,8) = NaN;
      end

      mp_out(k,m,2) = el;
      mp_out(k,m,3) = az;
    else
      mp_out(k,m,1) = NaN;
      mp_out(k,m,2) = el;
      mp_out(k,m,3) = az;
      mp_out(k,m,4) = NaN;
      mp_out(k,m,5) = NaN;
      mp_out(k,m,6) = NaN;
      mp_out(k,m,7) = NaN;
      mp_out(k,m,8) = NaN;
    end
    m = m + 1;
  end
  k = k + 1;
end

i = find( ~isfinite(mp_out));
mp_out(i) = 0;

%
% Now plot the data in various ways, mesh, surf, polar (with a contour
% plot superimposed) etc...
%

% MP1
i = find( ~isfinite(mp_out));
mp_out(i) = 0;
figure('Visible','off')
contour(10:10:90, -180:10:180, mp_out(:,:,1)');
colorbar;
xlabel('Elevation Angle [Degrees]');
Ylabel('Azimuth Angle [Degrees]');
title('MP1 Standard deviation [Metres]');
grid on
colorbar
print -dpng fig1.png


figure('Visible','off')
% Note centre of Polar plor is set to Zenith (90 degrees)
% If only I could get the bloody axis labeled correctly!
[X Y] = pol2cart( mp_out(:,:,3).*pi/180 , 90 - mp_out(:,:,2));
h = polar([0 2*pi], [10 90]);
delete(h)
hold on
contour(X,Y,mp_out(:,:,1));
colorbar
[x y] = pol2cart(pi,80);
text(x,y,'South 10 deg elev \rightarrow','HorizontalAlignment','right')
title('MP1 in Metres versus Elevation/Azimuth (Zenith centre of plot)');
print -dpng fig2.png

figure('Visible','off')
h = polar([0 2*pi], [10 90]);
delete(h)
hold on
surf(X,Y,mp_out(:,:,1));
colorbar
shading interp
[x y] = pol2cart(pi,80);
text(x,y,'South 10 deg elev \rightarrow','HorizontalAlignment','right')
title('MP1 in Metres versus Elevation/Azimuth (Zenith centre of plot)');
print -dpng fig2a.png


figure('Visible','off')
surf(10:10:90, -180:10:180, mp_out(:,:,1)');
view(0,90)
shading interp
xlabel('Elevation Angle [Degrees]');
Ylabel('Azimuth Angle [Degrees]');
title('MP1 Standard deviation [Metres]');
colorbar
print -dpng fig3.png

figure('Visible','off')
surf(10:10:90, -180:10:180, mp_out(:,:,4)');
view(0,90)
shading interp
xlabel('Elevation Angle [Degrees]');
Ylabel('Azimuth Angle [Degrees]');
title('Number of epochs of MP1 data');
colorbar
print -dpng fig4.png

% MP2
% i = find( ~isfinite(mp_out));
% mp_out(i) = 0;
figure('Visible','off')
contour(10:10:90, -180:10:180, mp_out(:,:,5)');
colorbar;
xlabel('Elevation Angle [Degrees]');
Ylabel('Azimuth Angle [Degrees]');
title('MP2 Standard deviation [Metres]');
grid on
colorbar
print -dpng fig5.png

figure('Visible','off')
% Note centre of Polar plor is set to Zenith (90 degrees)
% If only I could get the bloody axis labeled correctly!
[X Y] = pol2cart(  mp_out(:,:,3).*pi/180 , 90 - mp_out(:,:,2));
h = polar([0 2*pi], [10 90]);
delete(h)
hold on
contour(X,Y,mp_out(:,:,5));
colorbar
[x y] = pol2cart(pi,80);
text(x,y,'South 10 deg elev \rightarrow','HorizontalAlignment','right')
title('MP2 in Metres versus Elevation/Azimuth (Zenith centre of plot)');
print -dpng fig6.png

figure('Visible','off')
h = polar([0 2*pi], [10 90]);
delete(h)
hold on
surf(X,Y,mp_out(:,:,5));
colorbar
shading interp 
[x y] = pol2cart(pi,80);
text(x,y,'South 10 deg elev \rightarrow','HorizontalAlignment','right')
title('MP2 in Metres versus Elevation/Azimuth (Zenith centre of plot)');
print -dpng fig6a.png

figure('Visible','off')
surf(10:10:90, -180:10:180, mp_out(:,:,5)');
view(0,90)
shading interp
xlabel('Elevation Angle [Degrees]');
Ylabel('Azimuth Angle [Degrees]');
title('MP2 Standard deviation [Metres]');
colorbar
print -dpng fig7.png

figure('Visible','off')
surf(10:10:90, -180:10:180, mp_out(:,:,6)');
view(0,90)
shading interp
xlabel('Elevation Angle [Degrees]');
Ylabel('Azimuth Angle [Degrees]');
title('Number of epochs of MP2 data');
colorbar
print -dpng fig8.png

% CNR 1
figure('Visible','off')
contour(10:10:90, -180:10:180, mp_out(:,:,7)');
colorbar;
xlabel('Elevation Angle [Degrees]');
Ylabel('Azimuth Angle [Degrees]');
title('L1 CNR Standard deviation [dBHz]');
grid on
colorbar
print -dpng fig9.png


figure('Visible','off')
% Note centre of Polar plor is set to Zenith (90 degrees)
% If only I could get the bloody axis labeled correctly!
[X Y] = pol2cart( mp_out(:,:,3).*pi/180 , 90 - mp_out(:,:,2));
h = polar([0 2*pi], [10 90]);
delete(h)
hold on
contour(X,Y,mp_out(:,:,7));
colorbar
[x y] = pol2cart(pi,80);
text(x,y,'South 10 deg elev \rightarrow','HorizontalAlignment','right')
title('L1 CNR in dBHz versus Elevation/Azimuth (Zenith centre of plot)');
print -dpng fig10.png


% CNR 2
figure('Visible','off')
contour(10:10:90, -180:10:180, mp_out(:,:,8)');
colorbar;
xlabel('Elevation Angle [Degrees]');
Ylabel('Azimuth Angle [Degrees]');
title('L2 CNR Standard deviation [dBHz]');
grid on
colorbar
print -dpng fig11.png


figure('Visible','off')
% Note centre of Polar plor is set to Zenith (90 degrees)
% If only I could get the bloody axis labeled correctly!
[X Y] = pol2cart( mp_out(:,:,3).*pi/180 , 90 - mp_out(:,:,2));
h = polar([0 2*pi], [10 90]);
delete(h)
hold on
contour(X,Y,mp_out(:,:,8));
colorbar
[x y] = pol2cart(pi,80);
text(x,y,'South 10 deg elev \rightarrow','HorizontalAlignment','right')
title('L2 CNR in dBHz versus Elevation/Azimuth (Zenith centre of plot)');
print -dpng fig12.png


save



