Polar plot for Radiation Patern in dB format script is here

How to use openEMS. Discussion on examples, tutorials etc

Moderator: thorsten

Post Reply
Hale_812
Posts: 171
Joined: Fri 13 May 2016, 02:54

Polar plot for Radiation Patern in dB format script is here

Post by Hale_812 » Thu 12 Apr 2018, 11:52

As you know, there is a big problem plotting radiation patterns in Matlab/Octave.
Lucky, there was a good attempt, made by Steve Rickman.
Unfortunately, the script he shared in Matlab community has few flaws. The main is missing check for data boundaries, which results in radially inverted plots of deep -dB values.
There are also some cosmetic flubs, like visually interfering grid and radial tick labels, etc.
Probably, This is the third time I am fixing this script in last 6 years... but I am always loosing this file when changing place of work. So, let me share it with you, if you like it.
Capture.PNG
Screenshot. Screen output and PDF printout.
Capture.PNG (76.82 KiB) Viewed 3889 times
Please, feel free to update it for better text labels allocation etc.

The usage with OpenEMS is as follows:
figure(5);
set (gcf, "papersize", [400, 300]./72, "paperposition", [0, 0, 400, 300]./72) ; %make PDF "printout" approximately equivalent to screen output
set(gcf, "position", [430, ceil((get(0, 'Screensize')(4))-820), 400,300]); %common dimensions with readable font proportions
leg=[]; %legend collection
cmap = colormap('lines');%hsv(length(f0)); %Enforce common colormap. Sometimes Octave have wrong ideas about it, when running user-scripted plotting tools.
for iiF=1:length(f0) %multiple frequencies array
h=dirplot((nf2ff.phi .*180./pi-89.9), 20*log10(nf2ff.E_norm{iiF}./max(nf2ff.E_norm{iiF}(:)))+10*log10(nf2ff.Dmax(iiF)),[PlotMaxdB,PlotMaxdB-dBspan,dBtick]);
%The formula is from Thorsten's OpenEMS tutorial. Not sure how much correct is it, because I don't know method of nf2ff.E_norm, and Dmax calculation.
drawnow; %wait for slow drawing.
hold on;
set(h,'linewidth',2); %Because data-curve handle was returned, we can use it without paganism and witchery.
set(h,'Color',cmap(iiF,:));
leg{iiF}=['f = ' num2str(f0(iiF)/1e9,4) ' GHz'];
endfor
legend(leg, "location", "southwest");
title("H-plane pattern");
drawnow; %again, let it finish drawing.
hold off; %and finish our manipulations.

Code: Select all

function hpol = dirplot(theta,rho,line_style,params)
% DIRPLOT  Polar directivity plot.
%   A modification of The Mathworks POLAR function, DIRPLOT generates
%   directivity plots in the style commonly used in acoustic and RF work.
%   Features include:
%     1. Plots -90 to +90 or -180 to +180 degrees based on range of input
%        THETA, with 0 degrees at top center.
%     2. Produces semicircular plots when plot range is -90 to +90 degrees.
%     3. RHO is assumed to be in decibels and may include negative
%        values.
%     4. Default automatic rho-axis scaling in "scope knob" factors.
%     5. Optional PARAMS argument allows manual setting of rho-axis
%        scaling.
%   
%   DIRPLOT(THETA, RHO) makes a plot using polar coordinates of the
%   angle THETA versus the radius RHO. THETA must be in degrees, and
%   must be within the range -180 to +180 degrees. If THETA is within
%   the range -90 to +90 degrees, the plot will be semicircular. RHO is
%   assumed to be in decibels and the values may be positive or negative or
%   both. By default, with no PARAMS argument, rho-axis scaling will be determined
%   automatically using scope knob factors of 1-2-5. By default, 10
%   ticks will be plotted. Note: Like POLAR, DIRPLOT does not rescale the
%   axes when a new plot is added to a held graph.
%
%   DIRPLOT(THETA, RHO, LINE_STYLE, PARAMS) makes a plot as described above
%   using the linestyle specified in string LINE_STYLE, and using the rho-axis
%   scaling specified in vector PARAMS. Either of these optional arguments may be
%   used alone. Vector PARAMS is a 3-element row vector defined as
%   [RHOMAX RHOMIN RHOTICKS]. String LINE_STYLE is the standard MATLAB linestyle
%   string. See PLOT for a description.
%
%   HPOL = DIRPLOT(...) returns a handle to the LINE object generated by the PLOT
%   function that actually generates the plot in DIRPLOT.
% 
%   See also POLAR, PLOT, LOGLOG, SEMILOGX, SEMILOGY.
% 
%   Rev 1.0, 17 January 2002
%   Tested in MATLAB v. 6.0
%
%   Adapted from The MathWorks POLAR function by
%   Steve Rickman
%   sar@surewest.net
%   Fixed by Hale in April 2018. Thanks Steve, for initial work!

if nargin <= 1
    error('Requires 2, 3, or 4 input arguments.')
elseif nargin == 2
    line_style = 'auto';
elseif nargin == 3 
    if isnumeric(line_style)
        params = line_style;
        line_style = 'auto';
    end
end
if exist('params')
    if length(params) ~= 3
        error('Argument PARAMS must be a 3-element vector: [RHOMAX RHOMIN RHOTICKS].')
    end
    if params(1) <= params(2)
        error('Error in PARAMS argument. RHOMAX must be greater than RHOMIN.')
    end
    if params(3) <= 0
        params(3) = 1;
        warning('Error in PARAMS argument. RTICKS set to 1.')
    end
end
if isstr(theta) | isstr(rho)
    error('THETA and RHO must be numeric.');
end
if ~isequal(size(theta),size(rho))
    error('THETA and RHO must be the same size.');
end
if (max(theta) - min(theta)) < 6.3
    warning('THETA must be in degrees');
end
if min(theta) >= 0
    warning('Plot is -90 to +90 or -180 to +180 degrees');
end
if max(abs(theta)) > 180
    error('Plot is -90 to +90 or -180 to +180 degrees');
end

% Get range of theta and set flag for full or half plot.
if (max(theta)-min(theta)) > 180 | max(theta) > 90
    fullplot = 1;
else
    fullplot = 0;
end

% Translate theta degrees to radians
theta = theta*pi/180;

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

if hold_state & exist('params')
    warning('Plot is held. New plot parameters ignored')
end

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

% Hold on to current Text defaults, reset them to the
% Axes' font attributes so tick marks use them.

% ARNAUT: made a little change to fontsize
fAngle  = get(cax, 'DefaultTextFontAngle');
fName   = get(cax, 'DefaultTextFontName');
fSize   = get(cax, 'DefaultTextFontSize');
fWeight = get(cax, 'DefaultTextFontWeight');
fUnits  = get(cax, 'DefaultTextUnits');
set(cax, 'DefaultTextFontAngle',  get(cax, 'FontAngle'), ...
    'DefaultTextFontName',   get(cax, 'FontName'), ...
    'DefaultTextFontSize',   get(cax, 'FontSize') + 5, ...
    'DefaultTextFontWeight', get(cax, 'FontWeight'), ...
    'DefaultTextUnits','data')

% only do grids if hold is off
if ~hold_state    
    % make a radial grid
    hold on;
    if ~exist('params')
        rticks = 10; % default ticks
        lims = findscale(rho,rticks); % get click, rmax, rmin
        click = lims(1); rmax = lims(2); rmin = lims(3);
        rngdisp = rmax - rmin;
    else
        rmax = params(1); rmin = params(2); rticks = params(3);
        rngdisp = rmax - rmin;
        click = rngdisp/rticks;        
    end
   
    set(cax,'userdata',[rngdisp rmax rmin]); % save variables for added plots
 
    % define a circle
    th = 0:pi/50:2*pi;
    xunit = cos(th);
    yunit = sin(th);
    % now really force points on x/y axes to lie on them exactly
    inds = 1:(length(th)-1)/4:length(th);
    xunit(inds(2:2:4)) = zeros(2,1);
    yunit(inds(1:2:5)) = zeros(3,1);
    % plot background if necessary
    if ~isstr(get(cax,'color')),
        patch('xdata',xunit*rngdisp,'ydata',yunit*rngdisp, ...
            'edgecolor',tc,'facecolor',get(gca,'color'),...
            'handlevisibility','off');
    end
    
    % draw radial circles
    % angles for text labels
    c88 = cos(80*pi/180);
    s88 = sin(80*pi/180);
    c92 = -cos(149*pi/180);
    s92 = -sin(149*pi/180);
    
    for i=click:click:rngdisp
        tickt = i+rmin;
        if abs(tickt) < .001
            tickt = 0;
        end
        ticktext = ['' num2str(tickt)];
        hhh = plot(xunit*i,yunit*i,ls,'color',tc+.5,'linestyle',':','linewidth',0.5,...
            'handlevisibility','off');
        if not(fullplot)
	        if i < rngdisp
	            text((i-.4)*c88,(i-.4)*s88, ...
	                ticktext,'rotation', 10,'verticalalignment','middle',...
	                'horizontalalignment','right','handlevisibility','off','fontsize',12)
	        else
	            text((i-.4)*c88,(i-.4)*s88, ...
	                [ticktext,' dB'],'rotation', 10,'verticalalignment','middle',...
	                'horizontalalignment','center','handlevisibility','off','fontsize',12)
	        end
	end
        if fullplot
            if i < rngdisp
                text((i-1.5)*c92,(i-1.5)*s92, ...
                    ticktext,'rotation', 59 , 'verticalalignment','middle',...
                    'horizontalalignment','right','handlevisibility','off','fontsize',12)
            else
                text((i-1.5)*c92,(i-1.5)*s92, ...
                    [ticktext,' dB'],'verticalalignment','middle','rotation', 59,...
                    'horizontalalignment','center','handlevisibility','off','fontsize',12)
            end            
        end
    end
    set(hhh,'linestyle','-','color','black') % Make outer circle solid
    
    % plot spokes at 10 degree intervals
    th = (0:18)*2*pi/36;
    lss=repmat(":",length(th),1);
    lss(find(mod((0:18).*2,3)==0))=ls;
    cst = cos(th); snt = sin(th);
    cs = [-cst; cst];
    sn = [-snt; snt];
    h=plot(rngdisp*cs,rngdisp*sn,ls,'color',tc+.5,'linestyle',ls,'linewidth',.5,...
        'handlevisibility','off');
    hold off;
    for iih=1:length(h)
    	set(h(iih),'linestyle',lss(iih));
    endfor;
    hold on;

    % label spokes in 30 degree intervals
    rt = 1.15*rngdisp;
    for i = 1:3:19
        text(rt*cst(i)*1.03,rt*snt(i),[int2str(90-(i-1)*10),'\deg '],...
            'horizontalalignment','center',...
            'handlevisibility','off'); 
    end
    if fullplot
        for i = 3:3:6
            text(-rt*cst(i+1)*1.03,-rt*snt(i+1),[int2str(-90-i*10),'\deg '],...
                'horizontalalignment','center',...
                'handlevisibility','off'); 
        end
        for i = 9:3:15
            text(-rt*cst(i+1)*1.03,-rt*snt(i+1),[int2str(270-i*10),'\deg '],...
                'horizontalalignment','center',...
                'handlevisibility','off'); 
        end        
    end
    
    % set view to 2-D
    view(2);
    % set axis limits
    if fullplot
        axis(rngdisp*[-1 1 -1.15 1.15]);
    else
        axis(rngdisp*[-1 1 0 1.15]);        
    end
end

if hold_state
    v = get(cax,'userdata');
    rngdisp = v(1);
    rmax = v(2);
    rmin = v(3);
end
        
% Reset defaults.
set(cax, 'DefaultTextFontAngle', fAngle , ...
    'DefaultTextFontName',   fName , ...
    'DefaultTextFontSize',   fSize, ...
    'DefaultTextFontWeight', fWeight, ...
    'DefaultTextUnits',fUnits );

%limit data at min and max for avoiding plotting to negative, and out of bounds.
rho_plot=rho;
rho_plot(find(rho<rmin))=rmin;
rho_plot(find(rho>rmax))=rmax;


% transform data to Cartesian coordinates.
% Rotate by pi/2 to get 0 degrees at top. Use negative
% theta to have negative degrees on left.
xx = (rho_plot+rngdisp-rmax).*cos(-theta+pi/2);
yy = (rho_plot+rngdisp-rmax).*sin(-theta+pi/2);

% plot data on top of grid
if strcmp(line_style,'auto')
    q = plot(xx,yy);
else
    q = plot(xx,yy,line_style,'LineWidth',1.5);
end
if nargout > 0
    hpol = q;
end
set(gca,'dataaspectratio',[1 1 1]), axis off; set(cax,'NextPlot',next);
set(get(gca,'xlabel'),'visible','on')
set(get(gca,'ylabel'),'visible','on')

% Subfunction finds optimal scaling using "scope knob"
% factors of 1, 2, 5. Range is limited to practical
% decibel values.
function lims = findscale(rho, rticks)
	clicks = [.001 .002 .005 .01 .02 .05 .1 ...
              .2 .5 1 2 5 10 20 50 100 200 500 1000];
	lenclicks = length(clicks);
	rhi = max(rho);
	rlo = min(rho);
	rrng = rhi - rlo;
	rawclick = rrng/rticks;
	n = 1;
	while clicks(n) < rawclick
        n = n + 1;
        if n > lenclicks
            close;
            error('Cannot autoscale; unrealistic decibel range.');
        end
	end
	click = clicks(n);
	
	m = floor(rhi/click);
	rmax = click * m;
	if rhi - rmax ~= 0
        rmax = rmax + click;
	end	
	rmin = rmax - click * rticks;
	
	% Check that minimum rho value is at least one tick
	% above rmin. If not, increase click value and
	% rescale.
    if rlo < rmin + click
        if n < lenclicks
            click = clicks(n+1);
        else
            error('Cannot autoscale; unrealistic decibel range.');
        end
        
        m = floor(rhi/click);
        rmax = click * m;
        if rhi - rmax ~= 0
            rmax = rmax + click;
        end
        rmin = rmax - click * rticks;
    end
    lims = [click rmax rmin];
Don't panic, no endfunction is used, because it is the end of the file. Matlab guys often do such weird things.
Last edited by Hale_812 on Fri 13 Apr 2018, 02:43, edited 2 times in total.

thorsten
Posts: 1393
Joined: Mon 27 Jun 2011, 12:26

Re: Polar plot for Radiation Patern in dB format script is h

Post by thorsten » Thu 12 Apr 2018, 20:24

Hi, sounds very nice.
What are the main advantages in comparison to the plotFFdB function?

Hale_812
Posts: 171
Joined: Fri 13 May 2016, 02:54

Re: Polar plot for Radiation Patern in dB format script is h

Post by Hale_812 » Fri 13 Apr 2018, 04:42

Andvantages are:
1)Sector lines
2)Angle and dB labels
3)Correct handling of handles results in correct colors in the legend. Sorry Thorsten, there is a flaw in your code making all the multiplot legend 'gray'.
4)Feeding angular range is possible. It is useful when you have a non standard model direction, but interested in readable patterns with ZERO in zenith.
5)Grid handles are hidden, only plot handle is returned.
6)Plotting over does not repeat drawing of grid lines. Scaling is remembered.
7)Plotting half circle pattern is possible, when feeding angles between (-90;90). (needs rewriting....maybe later)

Disadvantages:
1)Unclear copyright. Seems like wild coding, so use it at your risk. I do not recommend inserting it into OpenEMS as is. It should be rewritten in the same style, probably.
2)dB label still can not be switched off, need more properties scripting (which is boring as hell, when you make it right, and failproof)
3)Octave does not rotate screen fonts, so they become overlapped at 60 degree. Printed fonts are rotated as you see.
4*)Octave (gnuplot?) inserts weird space before \deg, °, or ^o symbols, when printing to PDF. I don't know why, but it causes headache when trying to meet both, onscreen and printout formats.

*this is a bug of Octave 4.2.2, and should be fixed in 4.4.4, according to developers.

PaulUK
Posts: 76
Joined: Wed 28 Oct 2015, 14:23

Re: Polar plot for Radiation Patern in dB format script is here

Post by PaulUK » Sat 24 Nov 2018, 10:04

Thank you for doing this, I'll give it a try.

Hale_812
Posts: 171
Joined: Fri 13 May 2016, 02:54

Re: Polar plot for Radiation Patern in dB format script is here

Post by Hale_812 » Fri 04 Jan 2019, 06:09

>Noob question:
you should refer to any linux forums, or man. in short, the latter is a better serviced envelope to the former. And the former itself is a satellite tool to apt. so the usage depends on how tricky you are in linux environment.

Post Reply