%% polar Azimuth/Range/Amplitude plot, or Plan Position Indicator (PPI)

function [logA,map]=ppi(nc,txpol,rpix,mode)

if nargin<2, txpol='vv'; end
if nargin<3, rpix=200; end
if nargin<4, mode=''; end
seaonly=strcmp(mode,'seaonly');

% averaging to match window size
nsweep=size(nc{'adc_data'},1);
range=nc{'range'}(:);
xavg=max(1,ceil(nsweep/(4*pi*rpix)));
nrange=length(range); 
maxrange=max(range);

% process radar returns
AA=zeros(nrange,floor(nsweep/xavg));
ampl=zeros(xavg,floor(nsweep/xavg)); lenAmpl=prod(size(ampl));
for r=1:nrange,
  [I,Q,mnIQ,sdIQ]=ipixload(nc,txpol,r,'ohgr');
  absiq=abs(I+j*Q)*sqrt(prod(sdIQ));
  ampl(:)=absiq(1:lenAmpl);
  AA(r,:)=mean(ampl.^2,1);
  if seaonly, 
    % correct for power drop off
    AA(r,:)=AA(r,:)*(range(r)/range(end))^3; 
  end
end

% convert to polar coordinates
az=zeros(xavg,floor(nsweep/xavg)); lenAz=prod(size(az));
az(:)=unwrap(nc{'azimuth_angle'}(1:lenAz)*pi/180);
az=mean(az,1);
% remove corrupt azimuth angles
daz=median(diff(az));
indx=find(sign(diff(az))~=sign(daz));
if ~isempty(indx),
  az=az(1:indx(1)); AA=AA(:,1:indx(1));
end
if seaonly,
  % remove land
  indx=find(az>70*pi/180 & az<215*pi/180);
  az=az(indx); AA=AA(:,indx);
end
ppr=rpix/maxrange;
xPolar=ppr*range*cos(0.5*pi-az);
yPolar=ppr*range*sin(0.5*pi-az);
[xI,yI]=meshgrid(min(xPolar(:))-8:max(xPolar(:))+8,...
                 min(yPolar(:))-8:max(yPolar(:))+8);
xI=round(xI); yI=round(yI);
rI=sqrt(xI.^2+yI.^2); 
azI=0.5*pi-angle(xI+j*yI); indx=find(azI<min(az)); azI(indx)=azI(indx)+2*pi;
ok=find(rI>=ppr*min(range) & rI<=ppr*maxrange & azI>=min(az) & azI<=max(az));
AAI=interp2(az,ppr*range,AA,azI(ok),rI(ok));
logA=zeros(size(xI)); logA(ok)=log(AAI);

% add rticks
indx=find(rI>=2 & rI<=4);
logA(indx)=1;
rtick=[0:10:350]*pi/180;
indx=find(rtick>min(az) & rtick<max(az));
rtick=rtick(indx);
for r=1:length(rtick)
  indx1=find((rI>rpix & rI<=rpix+8));
  indx2=find(azI(indx1)>rtick(r)-0.75./rI(indx1) & azI(indx1)<rtick(r)+0.75./rI(indx1));
  logA(indx1(indx2))=1;
end
rtick=[0:90:350]*pi/180;
indx=find(rtick>min(az) & rtick<max(az));
rtick=rtick(indx);
for r=1:length(rtick)
  indx1=find(((rI>rpix & rI<=rpix+8) | (rI>4 & rI<ppr*min(range))));
  indx2=find(azI(indx1)>rtick(r)-1.5./rI(indx1) & azI(indx1)<rtick(r)+1.5./rI(indx1));
  logA(indx1(indx2))=1;
end

if seaonly, 
  % scale for best visual contrast
  ea=360-mean(nc{'elevation_angle'}(:));
  rmin=30/((ea+0.45)*pi/180);
  rmax=30/(ea*pi/180);
  best=find(rI>=rpix*rmin/maxrange & rI<=rpix*rmax/maxrange & ...
            azI>=min(az) & azI<=max(az));
  if length(best)<0.01*length(ok), best=ok; end
  sortA=sort(logA(best));
  len=length(sortA);
  med=sortA(round(len/2));
  mad=median(abs(sortA-med));
  mn=med-3*mad; mx=med+3*mad;
else 
  mn=min(logA(ok)); mx=max(logA(ok));
end
indx=find(logA(ok)<mn); logA(ok(indx))=mn;
indx=find(logA(ok)>mx); logA(ok(indx))=mx;
logA(ok)=2+(logA(ok)-mn)*61/(mx-mn);

if seaonly, map=gray; map(:,3)=map(:,3)*0.9; else map=jet; end
map(1,:)=[1 1 1];
map(2,:)=[0 0 0];
if ~nargout,
  colormap(map);
  image(logA+1,'cdatamapping','direct'); set(gca,'ydir','normal');
  set(gcf,'units','pixels','position',[60 20 100+size(logA,2) 80+size(logA,1)]);
  set(gca,'units','pixels','position',[60 60 size(logA,2) size(logA,1)]);
end
