function radpat(ang1,R1,st1,ang2,R2,st2,ang3,R3,st3,ang4,R4,st4)
%RADPAT  Polar coordinate plot used for antenna radiation patterns.
%          RADPAT(ANG1,R1,ST1,ANG2,R2,ST2,ANG3,R3,ST3,ANG4,R4,ST4) 
%            plots up to four curves in dB format. 
%
%          ANGi are angles in degrees, 
%          Ri are radiation pattern values (radii for plot traces), &
%          STi are the linestyles.
%          See PLOT for a description of legal linestyles.
%
%          Ri can be in dB or not in dB (resulting plot is in dB).
%          Axis labels can be placed on horizontal or vertical axis.
%          Choice of normalized or unnormalized (show gains) patterns.
%          Minimum dB level at plot center can be specified.
%          Maximum dB level at outermost plot circle can be specified 
%            for unnormalized patterns.
%          Line width of radiation patterns can be specified.
%          Legend can be placed. To move the legend, press left mouse
%            button on the legend and drag to the desired location.
%          Grid linetype can be specified.
%          Default values are inside [], press Enter to chose default.
%          0 degrees can be at North/Top or East/Right side of plot.
%
%          Example:  radpat(a1,r1,'r-',a2,r2,'y--')
%
%          Based on polarpat.m by D. Liu, 9/13/1996.
%          T.J. Watson Research center, IBM
%          P.O.Box 218
%          Yorktown Heights, NY 10598
%          Email: dliu@watson.ibm.com
%
%          Updated by Thomas P. Montoya, SDSM&T, 1/26/2006
%          * allow up to four traces
%          * added degree symbols to plot spoke labels
%          * for plots vs. theta keep spoke labels in 0 to +180 deg
%             range and indicate that negative theta angles are for 
%             phi+180 deg and orient plot so that 0 degrees at the
%             top (North)
%          * added check on unitless inputs as to whether they are
%            proportional to power or not to get proper dB scaling

if nargin==0
  help radpat;
  return;
end

% get some parameters from user
db = input('Are input values in dB (Y/N)[Y]? ','s');
if isempty(db), db = 'Y'; end
db = upper(db);
if db == 'N'
  pwr = input('Input values proportional to power (Y/N) [Y]? ','s');
  if isempty(pwr), pwr = 'Y';  end 
  pwr = upper(pwr);
  if pwr == 'N'     % make unitless inputs proportional to power
      R1 = R1.*R1;  % to get proper dB plots
      if nargin > 3, R2 = R2.*R2;  end
      if nargin > 6, R3 = R3.*R3;  end
      if nargin > 9, R4 = R4.*R4;  end
  end
end

scale = input('Normalize to the Maximum Gain Value (Y/N)[Y]? ','s');
if isempty(scale), scale = 'Y'; end
scale = upper(scale);
maxG=0;
if scale == 'N'
  maxG = input('Maximum dB value for outer circle of plot [0]? ');
  if isempty(maxG), maxG=0;  end 
end

minG = input('Minimum dB value at plot center [-40]? ');
if isempty(minG), minG=-40; end
aminG=abs(minG); 

thetaflag = input('Are the angles theta values? (Y/N)[Y]? ','s');
if isempty(thetaflag), thetaflag = 'Y'; end
thetaflag = upper(thetaflag);

if thetaflag == 'Y'
    DIRECT='N';
else
    DIRECT=input('0 deg at North/Top or East/Right (N/E)[N]? ','s');
    if isempty(DIRECT), DIRECT='N'; end
    DIRECT=upper(DIRECT);
end

HV = input('Labels on Vertical or Horizontal axis (V/H)[V]? ','s');
if isempty(HV), HV = 'V'; end
HV = upper(HV);

LW = input('Pattern line width [1.25]: ');
if isempty(LW), LW=1.25; end

if nargin > 3
  LG = input('Legend for traces on graph (Y/N)[N]? ','s');
  if isempty(LG), LG='N'; end
  LG=upper(LG);
  if LG == 'Y'
    str1=input('Enter label for trace 1: ','s');
    str2=input('Enter label for trace 2: ','s');
    if nargin>6
      str3=input('Enter label for trace 3: ','s');
    end
    if nargin>9
      str4=input('Enter label for trace 4: ','s');
    end
    BOX = input('Put a box around the legend (Y/N)[Y]? ','s');
    if isempty(BOX), BOX='Y'; end
    BOX=upper(BOX); 
  end
end

LP = input('Line type of grid(-, --, -., :)[:]? ','s');
if isempty(LP), LP=':'; end

dr=pi/180; cc=10/log(10); er=exp(minG/cc);

% find maximum values
maxA=360; maxR=max(R1);
if nargin > 3, maxR=max([maxR, max(R2)]); end
if nargin > 6, maxR=max([maxR, max(R3)]); end
if nargin > 9, maxR=max([maxR, max(R4)]); end

if scale=='N' & maxR>maxG
  maxRdb=maxR;
  if db=='N', maxRdb=10*log10(maxR); end
  maxG=10*ceil((maxRdb-minG)/10)+minG;
end 
minmax=maxG-minG;

% normalize Ri vectors
if db == 'N'
  if scale=='Y', R1=R1/maxR;  end  
  for i=1:length(R1)
    if R1(i) < er,  R1(i)=er;  end
    R1(i)=aminG+cc*log(R1(i));
  end
  if nargin > 3
    if scale=='Y', R2=R2/maxR;  end  
    for i=1:length(R2)
      if R2(i) < er,  R2(i)=er;  end
      R2(i)=aminG+cc*log(R2(i));
    end
  end
  if nargin > 6
    if scale=='Y', R3=R3/maxR;  end  
    for i=1:length(R3)
      if R3(i) < er,  R3(i)=er;  end
      R3(i)=aminG+cc*log(R3(i));
    end
  end
  if nargin > 9
    if scale=='Y', R4=R4/maxR;  end  
    for i=1:length(R4)
      if R4(i) < er,  R4(i)=er;  end
      R4(i)=aminG+cc*log(R4(i));
    end
  end
else
  if scale=='Y', R1=R1-maxR; end
  R1=R1-minG;
  for i=1:length(R1)
    if R1(i) < 0, R1(i)=0; end
  end
  if nargin > 3 
    if scale=='Y', R2=R2-maxR; end
    R2=R2-minG;
    for i=1:length(R2)
      if R2(i) < 0, R2(i)=0; end
    end
  end
  if nargin > 6, 
    if scale=='Y', R3=R3-maxR; end
    R3=R3-minG;
    for i=1:length(R3)
      if R3(i) < 0, R3(i)=0; end
    end
  end
  if nargin > 9, 
    if scale=='Y', R4=R4-maxR; end
    R4=R4-minG;
    for i=1:length(R4)
      if R4(i) < 0, R4(i)=0; end
    end
  end
end

% get hold state
cax = newplot;
next = lower(get(cax,'NextPlot'));
hold_state = ishold;

% get x-axis text color so grid is in same color
tc = get(cax,'xcolor');

% only do grids if hold is off
if ~hold_state
% make a radial grid
   hold on;
   hhh=plot([0 maxA],[0 minmax]);
   v = [get(cax,'xlim') get(cax,'ylim')];
   ticks = length(get(cax,'ytick'));
   delete(hhh);

% check radial limits and ticks
   rmin = 0; rmax = v(4); rticks = ticks-1;
   if rticks > 9	% see if we can reduce the number
      if rem(rticks,2) == 0
	rticks = rticks/2;
      elseif rem(rticks,3) == 0
	rticks = rticks/3;
      end
   end

% define a circle
   th = 0:pi/50:2*pi;
   xunit = cos(th);
   yunit = sin(th);

   rinc = (rmax-rmin)/rticks;
   cc=minmax/(rmax-rmin);
   for i=(rmin+rinc):rinc:rmax
       if i == rmax    %force solid outer ring
           plot(xunit*rmax,yunit*rmax,'-','color',tc,'LineWidth',1.3);
       else
           plot(xunit*i,yunit*i,LP,'color',tc);
       end
     ii=i-rmax;
     if scale == 'N',  ii=minG+cc*(i-rmin);  end
     if HV == 'V'
       if ii < 0
         strr=num2str(ii);
       elseif ii == 0
         strr=[' ' num2str(ii)];
       else
         strr=[' ' num2str(ii)];
       end
       text(0.8,i+rinc/20+1,strr,'verticalalignment','middle');
     else
       text(i+rinc/20,0,num2str(ii),'verticalalignment','top',...
                                    'horizontalalignment','center');
     end
   end

   if HV == 'V'
       text(6.3,i+rinc/20+1,'dB','verticalalignment','middle',...
           'horizontalalignment','left');
   else
       text(i+rinc/20+2.2,-5,'dB','verticalalignment','top',...
                                   'horizontalalignment','center');
   end

% plot spokes
   th = (1:6)*2*pi/12; 
   if DIRECT=='N', th=th+pi/2; end
   cst = cos(th); snt = sin(th);
   cs = [-cst; cst];
   sn = [-snt; snt];
   if DIRECT=='N', cst=-cst; end
   plot(rmax*cs,rmax*sn,LP,'color',tc);

% annotate spokes in degrees
   rt = 1.1*rmax;
   for i = 1:max(size(th))
      loc = int2str(i*30);
      text(rt*cst(i),rt*snt(i),[loc,'\circ'],'horizontalalignment',...
          'center');
      if i == max(size(th))
        loc = int2str(0);
      else
        loc = int2str(180+i*30);
      end
      if thetaflag == 'Y';
        if loc == '330'
           loc = '30';
        elseif loc == '300'
          loc = '60';
        elseif loc == '270'
          loc = '90';
        elseif loc == '240'
          loc = '120';
        elseif loc =='210'
          loc = '150';
        end
      end
      text(-rt*cst(i),-rt*snt(i),[loc,'\circ'],...
          'horizontalalignment','center');
   end
   if thetaflag == 'Y';
       text(-rt*cst(max(size(th)))-0.75*rinc,...
       -rt*snt(max(size(th)))+0.75*rinc,...
       '\phi+180\circ \leftarrow \theta \rightarrow \phi',...
           'horizontalalignment','center');
   end

% set viewto 2-D
   view(0,90);
% set axis limits
   axis(rmax*[-1 1 -1.1 1.1]);
end

offset = 0;
% one curve
if nargin < 4
  if DIRECT=='N', 
      offset = pi;
      ang1 = ang1+90;
  end
  ang1 = ang1*dr;
  xx1 = R1.*cos(ang1+offset);
  yy1 = R1.*sin(ang1);
  hhh = plot(xx1,yy1,st1);
  set(hhh,'linewidth',LW)
end

% two curves
if nargin > 3 & nargin < 7
  if DIRECT=='N', 
      offset = pi;
      ang1 = ang1+90;
      ang2 = ang2+90;
  end
  ang1 = ang1*dr;
  ang2 = ang2*dr;
  xx1 = R1.*cos(ang1+offset);
  yy1 = R1.*sin(ang1);
  xx2 = R2.*cos(ang2+offset);
  yy2 = R2.*sin(ang2);
  hhh = plot(xx1,yy1,st1,xx2,yy2,st2);
  set(hhh,'linewidth',LW)
  if LG=='Y',
%      legend(str1,str2);
       hl = legend(hhh,str1,str2);
    if BOX=='N', set(hl,'Visible','off'); end  
  end
end

% three curves
if nargin > 6 & nargin < 10
  if DIRECT=='N', 
      offset = pi;
      ang1 = ang1+90;
      ang2 = ang2+90;
      ang3 = ang3+90;
  end
  ang1 = ang1*dr;
  ang2 = ang2*dr;
  ang3 = ang3*dr;
  xx1 = R1.*cos(ang1+offset);
  yy1 = R1.*sin(ang1);
  xx2 = R2.*cos(ang2+offset);
  yy2 = R2.*sin(ang2);
  xx3 = R3.*cos(ang3+offset);
  yy3 = R3.*sin(ang3);
  hhh = plot(xx1,yy1,st1,xx2,yy2,st2,xx3,yy3,st3);
  set(hhh,'linewidth',LW)
  if LG=='Y'
    hl=legend(hhh,str1,str2,str3);  
    if BOX=='N', set(hl,'Visible','off'); end  
  end
end

% four curves
if nargin > 9
  if DIRECT=='N', 
      offset = pi;
      ang1 = ang1+90;
      ang2 = ang2+90;
      ang3 = ang3+90;
      ang4 = ang4+90;
  end
  ang1 = ang1*dr;
  ang2 = ang2*dr;
  ang3 = ang3*dr;
  ang4 = ang4*dr;
  xx1 = R1.*cos(ang1+offset);
  yy1 = R1.*sin(ang1);
  xx2 = R2.*cos(ang2+offset);
  yy2 = R2.*sin(ang2);
  xx3 = R3.*cos(ang3+offset);
  yy3 = R3.*sin(ang3);
  xx4 = R4.*cos(ang4+offset);
  yy4 = R4.*sin(ang4);
  hhh = plot(xx1,yy1,st1,xx2,yy2,st2,xx3,yy3,st3,xx4,yy4,st4);
  set(hhh,'linewidth',LW)
  if LG=='Y'
    hl=legend(hhh,str1,str2,str3,str4);  
    if BOX=='N', set(hl,'Visible','off'); end  
  end
end

set(findobj('type','text'),'fontname','times')
set(findobj('type','axes'),'fontname','times')

if ~hold_state
	axis('equal');axis('off');
end

% reset hold state
if ~hold_state, set(cax,'NextPlot',next); end