%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% matlab script to illustrate computing C-band azimuth modulation from
% BYU .SIR azimuth modulation images computed from ASCAT data
%
% written by D.G. Long Jan 2011
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% location of BYU .sir files files
dloc='182-212-2010_az4/';
%dloc='182-212-2010_az4_nofilt/';

% define box around continent
x1=400; y1=460; x2=1600; y2=1480;

% create a circular color bar to show directions
cmap=colormap('gray');
cmap2(1:size(cmap,1)/2,:)=cmap(1:2:size(cmap,1),:);
cmap2((size(cmap,1)/2+1):size(cmap,1),:)=flipud(cmap(1:2:size(cmap,1),:));

% read and plot BYU .sir image data files
aname='msfa-a-Ant10-182-212.ave';
[a ahead]=loadsir([dloc aname]);

myfigure(1)
imagesc(a(y1:y2,x1:x2),[-25 0]);colorbar
title(sprintf('ASCAT H-pol a (%s) in dB @ 40 inc',aname));
axis off;
axis equal;
colormap('gray')
print -dpng ASCAT_a.png


bname='msfa-b-Ant10-182-212.ave';
[b bhead]=loadsir([dloc bname]);

myfigure(2)
imagesc(b(y1:y2,x1:x2),[-0.35 0]);colorbar
title(sprintf('ASCAT H-pol b (%s) in dB/deg',bname));
axis off;
axis equal;
colormap('gray')
print -dpng ASCAT_b.png


qname='msfa-q-Ant10-182-212.sir';
[q qhead]=loadsir([dloc qname]);

myfigure(3)
imagesc(q(y1:y2,x1:x2),[0 2]);colorbar
title(sprintf('ASCAT H-pol q (1st Az Har) (%s) in dB',qname));
axis off;
axis equal;
colormap('gray')
print -dpng ASCAT_q.png

rname='msfa-r-Ant10-182-212.sir';
[r rhead]=loadsir([dloc rname]);

myfigure(13)
imagesc(r(y1:y2,x1:x2),[0 360]);colorbar
title(sprintf('ASCAT H-pol r (1st Az Har phase) (%s) in deg',qname));
axis off;
axis equal;
colormap(cmap2)
print -dpng ASCAT_r.png


sname='msfa-s-Ant10-182-212.sir';
[s shead]=loadsir([dloc sname]);

myfigure(4)
imagesc(s(y1:y2,x1:x2),[0 6]);colorbar
title(sprintf('ASCAT H-pol s (2nd Az Har) (%s) in dB',sname));
axis off;
axis equal;
colormap('gray');
print -dpng ASCAT_s.png

tname='msfa-t-Ant10-182-212.sir';
[t thead]=loadsir([dloc tname]);

myfigure(14)
imagesc(t(y1:y2,x1:x2),[0 360]);colorbar
title(sprintf('ASCAT H-pol t (2nd Az Har phase) (%s) in deg',tname));
axis off;
axis equal;
colormap(cmap2);
print -dpng ASCAT_t.png


Qname='msfa-Q-Ant10-182-212.sir';
[Q thead]=loadsir([dloc Qname]);

myfigure(5)
imagesc(Q(y1:y2,x1:x2),[0 1]);colorbar
title(sprintf('ASCAT H-pol Q (3rd Az Har) (%s) in dB',Qname));
axis off;
axis equal;
colormap('gray')
print -dpng ASCAT_Q.png

Rname='msfa-R-Ant10-182-212.sir';
[R Rhead]=loadsir([dloc Rname]);

myfigure(15)
imagesc(R(y1:y2,x1:x2),[0 360]);colorbar
title(sprintf('ASCAT H-pol R (3rd Az Har phase) (%s) in deg',Rname));
axis off;
axis equal;
colormap(cmap2);
%colormap('gray');
print -dpng ASCAT_R.png


Sname='msfa-S-Ant10-182-212.sir';
[S thead]=loadsir([dloc Sname]);

myfigure(6)
imagesc(S(y1:y2,x1:x2),[0 2]);colorbar
title(sprintf('ASCAT H-pol S (4th Az Har) (%s) in dB',Sname));
axis off;
axis equal;
colormap('gray')
print -dpng ASCAT_S.png


Tname='msfa-T-Ant10-182-212.sir';
[T Thead]=loadsir([dloc Tname]);

myfigure(16)
imagesc(T(y1:y2,x1:x2),[0 360]);colorbar
title(sprintf('ASCAT H-pol T (4th Az Har phase) (%s) deg',Tname));
axis off;
axis equal;
colormap(cmap2);
colormap('gray');
print -dpng ASCAT_T.png


%%%%%%%
% show azimuth modulation curve from a particular location
%%%%%%%

% specify location
lat=-70.25;lon=124.45; % lindsley
[x, y] = latlon2pix(lon,lat,ahead);

% define box size for plotting location on image and plot
d=ceil(12.5/4.45)*2+1;
xx1=x-d; xx2=x+d;  yy1=y-d;  yy2=y+d;
xx=round(x);yy=round(y);

myfigure(11)
imagesc(a(y1:y2,x1:x2),[-25 0]);colorbar
title(sprintf('ASCAT V-pol (%s) sigma-0 in dB @ 40 deg inc',aname));
hold on;plot([xx1 xx1 xx2 xx2 xx1]-x1+1,[yy1 yy2 yy2 yy1 yy1]-y1+1,'m'); hold off;
hold on;plot(x-x1+1,y-y1+1,'r+'); hold off;
axis off;
axis equal;
colormap('gray')
print -dpng ASCAT_box.png

% print azimuth modulation image coefficients
disp(sprintf('ASCAT at location %4.2fS,%4.2fE',-lat,lon));
disp(sprintf('      a=%7.3f dB  b=%7.4f dB/deg',a(yy,xx),b(yy,xx)));
disp(sprintf(' 1st: q=%7.3f dB  r=%6.2f deg',q(yy,xx),r(yy,xx)));
disp(sprintf(' 2nd: s=%7.3f dB  t=%6.2f deg',s(yy,xx),t(yy,xx)));
disp(sprintf(' 3rd: Q=%7.3f dB  R=%6.2f deg',Q(yy,xx),R(yy,xx)));
disp(sprintf(' 4th: S=%7.3f dB  T=%6.2f deg',S(yy,xx),T(yy,xx)));

% compute azimuth modulation curve
DTR=pi/180;
ang=-180:180;
response=q(yy,xx)*cos((ang-r(yy,xx))*DTR)+s(yy,xx)*cos((2*ang-t(yy,xx))*DTR)+Q(yy,xx)*cos((3*ang-R(yy,xx))*DTR)+S(yy,xx)*cos((4*ang-T(yy,xx))*DTR);

% plot azimuth modulation curve
myfigure(12)
h=plot(ang,response);
set(h,'LineWidth',2);
xlabel('Relative azimuth angle (deg)');
ylabel('sigma-0 variation with azimuth (dB)');
title(sprintf('ASCAT at location %4.1fS %4.1fE  (a=%6.2f dB b=%0.4f dB/deg)',-lat,lon,a(yy,xx),b(yy,xx)));
axis([-180 180 -6 8])
print -dpng ASCAT_az.png
