sphere3d.m

a usefull source code from internet

function [Xout, Yout, Zout, Cmap] = sphere3d(Zin,theta_min,theta_max,...
                                    phi_min,phi_max,Rho,meshscale,varargin)
%
% SPHERE3D  Plots 3D data on a spherical surface.
% 
%           SPHERE3D(Zin,theta_min,theta_max,phi_min,phi_max,Rho,meshscale) 
%           plots the 3D profile Zin as a mesh plot on a spherical surface 
%           of radius Rho, between horizontal sweep angles theta_min and 
%           theta_max and vertical phi_min and phi_max, with mesh size 
%           determined by meshscale. 
% 
%           SPHERE3D(Zin,...,meshscale,plotspec) plots the 3D profile Zin 
%           with a plot type specification. If plotspec = 'surf' a standard
%           surface is plotted, whereas 'mesh', or 'meshc' will plot mesh,
%           or mesh with contour, respectively. A special contour is
%           plotted when plotspec = 'contour'. The default is mesh if not
%           specified.   
%
%           SPHERE3D(Zin,...,meshscale,interpspec) plots the 3D profile Zin
%           with the interpolation specification, interpspec, which can be 
%           one of 'spline', 'linear', 'nearest' or 'cubic'. The default  
%           interpolation is linear if not specified.     
%
%           SPHERE3D(Zin,...,meshscale,Zscale) plots the 3D profile Zin with 
%           the data scaling factor, Zscale. This allows you to scale the 
%           peaks and troughs of the data on the surface if the radius is
%           relatively large. Zscale larger than 1 magnifies the data
%           range correspondingly. Zscale defaults to 1 if not specified.
% 
%           [Xout,Yout,Zout,Cmap] = SPHERE3D(Zin,...) returns output values
%           corresponding to the Cartesian positions (Xout,Yout,Zout) with 
%           colour map, Cmap.         
% 
% SYNTAX    sphere3D(Zin,theta_min,theta_max,phi_min,phi_max,Rho,meshscale)
%           sphere3D(Zin,theta_min,theta_max,phi_min,phi_max,Rho,meshscale,plotspec)
%           sphere3D(Zin,theta_min,theta_max,phi_min,phi_max,Rho,meshscale,interpspec)
%           sphere3D(Zin,theta_min,theta_max,phi_min,phi_max,Rho,meshscale,Zscale)
%           sphere3D(Zin,...,meshscale,plotspec,interpspec)
%           sphere3D(Zin,...,meshscale,plotspec,Zscale)
%           sphere3D(Zin,...,meshscale,interpspec,Zscale)
%           sphere3D(Zin,...,meshscale,plotspec,interpspec,Zscale)
%           [Xout,Yout,Zout,Cmap] = sphere3D(Zin,...)
%           
% INPUT     Zin         input magnitude profiles where each row in Zin is 
%                       assumed to represent the horizontal sweep information 
%                       between theta_min and theta_max at a given vertical 
%                       sweep angle phi on the sphere. Alternatively, each 
%                       column represents the data gathered between top to
%                       bottom of the sphere at given angle theta.
%                       
%                       Zin is a (M x N) matrix, where M and N are not 
%                       necessarily equal. If M is not equal to N then the
%                       data are interpolated to make them equal. The final
%                       size is determined by the larger value of (M,N),
%                       and meshscale if different from 1.
% 
%                       The N columns of Zin are assumed to be equally
%                       spaced measurements starting from top (phi_max) to 
%                       bottom of the sphere (phi_min), and at angle theta_min  
%                       and so on, to the last column at theta_max. The
%                       vertical axis of the sphere is assumed to be 
%                       parallel to the columns in the data.
%
%                       Zin(1,1) corresponds to (theta_min,phi_max) and 
%                       Zin(M,1) corresponds to (theta_min,phi_min).
%                       Zin(1,N) corresponds to (theta_max,phi_max) and 
%                       Zin(M,N) corresponds to (theta_max,phi_min).  
%                       Theta increases in the anticlockwise direction 
%                       looking from the top of the sphere, while phi 
%                       increases from lower hemisphere to the upper.
% 
%           theta_min   the lower value in radians of the angular range 
%                       (horizontal sweep) over which the data is defined. 
%                       Theta_min is a scalar quantity.
% 
%           theta_max   the upper value in radians of the angular range 
%                       (horizontal sweep) over which the data is defined. 
%                       Theta_max is a scalar quantity.
%
%                       The difference between theta_max and theta_min
%                       should be no more than 2pi.
% 
%           phi_min     the lower value in radians of the angular range 
%                       (vertical sweep) over which the data is defined. 
%                       Phi_min is a scalar quantity.
% 
%           phi_max     the upper value in radians of the angular range 
%                       (vertical sweep) over which the data is defined. 
%                       Phi_max is a scalar quantity.
%
%                       The difference between phi_max and phi_min should
%                       be no more than pi.
% 
%           Rho         the radius of the cylinder. Rho is a scalar
%                       quantity. See comments on Zscale for connection  
%                       between Rho and Zscale.
% 
%           meshscale   a scalar that determines the size of the squares 
%                       on the mesh or surf plots, and takes on integer or
%                       non-integer values greater than 0.
%
%                       If meshscale = 1, the mesh remains unchanged relative
%                       to the input grid. If meshscale = 2.15, say, the size  
%                       of the squares is increased by this factor with
%                       a consequential decrease in the dimensions of Xout,  
%                       Yout and Zout by the same factor.
%
%                       If meshscale is less than 1, a decrease in the mesh
%                       squares will follow with a consequential increase 
%                       in the dimensions of Xout, Yout and Zout by 1/meshscale.                                          
% 
%           plotspec    = 'surf'  produces a surface plot.
%                       = 'mesh'  produces a mesh plot.
%                       = 'meshc'  produces a mesh plot with countour in 
%                          the X-Y plane.
%                       = 'contour'  produces a transparent 2D contour plot
%                          on the curved surface of the cylinder.
%                       = 'off' disengages plot function.
% 
%           interpspec  = 'linear'  bilinear interpolation on Zin.
%                       = 'spline'  spline interpolation on Zin.
%                       = 'nearest'  nearest neighbour interpolation on Zin.
%                       = 'cubic'  bicubic interpolation on Zin.
%
%                       If Zin is a square matrix and meshscale = 1, no
%                       interpolation is carried out.
% 
%           Zscale      a scalar that allows the user to scale the peaks and
%                       troughs of the data so that it is more visible, 
%                       especially when the radius of the sphere is large 
%                       in comparison to the swings in the data.
%
%                       When Zscale = 1, no data scaling occurs and the object 
%                       appears in its true shape. 
%
%                       When Zscale is not 1, the maximum and minimum
%                       radial values are rescaled and may become negative.
%                       When this occurs Zscale defaults to the maximum 
%                       scaling factor possible and a warning is given.  
%                       No action is required by the user.
%
%                       When Zscale is small in relation to the radius, a 
%                       smooth spherical surface is produced with radius Rho.
%                       When used in conjunction with surf plot, the result
%                       gives the impression of a filled contour plot.
% 
% OUTPUT    Zout        output magnitude profiles defined by Zin at
%                       positions (Xout,Yout). 
%                       
%                       Zout is square with dimensions determined by the 
%                       maximum dimension of the input matrix Zin. The 
%                       dimensions of Zout are reduced or enlarged by meshscale.                      
% 
%           Xout        output X-positions corresponding to polar positions
%                       (rho,theta). Xout is square with dimensions  
%                       determined by the maximum dimension of the input 
%                       matrix Zin. The dimensions of Xout are reduced or 
%                       enlarged by meshscale.
% 
%           Yout        output Y-positions corresponding to polar positions
%                       (rho,theta). Yout is square with dimensions  
%                       determined by the maximum dimension of the input 
%                       matrix Zin. The dimensions of Yout are reduced or 
%                       enlarged by meshscale.
%
%           Cmap        colour mapping associated with (Xout,Yout,Zout). The
%                       dimension of Cmap is square and similar in size to
%                       Xout, Yout and Zout.
%
% See also POLAR3D, CYL3D, POLAR, POL2CART, SPH2CART and INTERP2

%           Written by JM DeFreitas, QinetiQ Ltd, Winfrith Technology
%           Centre, Dorchester DT2 8XJ, UK. [e-mail].
%
%           Released 27 September 2005. (Beta Release).
%
%           Terms and Conditions of Use
% 
%           1.	This function is made available to Matlab� users under the 
%               terms and conditions set out in the Matlab Exchange by 
%               The Mathworks, Inc.
%           2.	Where the use of SPHERE3D is to be cited, the following is recommended:
%               J M De Freitas. SPHERE3D: A Matlab Function to Plot 3-Dimensional 
%               Data on a Spherical Surface. 
%               QinetiQ Ltd, Winfrith Technology Centre, Winfrith, 
%               Dorchester DT2 8XJ. UK. 15 September 2005. 
%           3.	No offer of warranty is made or implied by the author and 
%               use of this work means that the user has agreed to take full 
%               responsibility for its use.
%
%           Edited by K. Van Caekenberghe on 17/03/2009
if (nargin < 7)
    disp('SPHERE3D Error: Too few input arguments.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
elseif (nargin > 10)
    disp('SPHERE3D Error: Too many input arguments.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
end

[p,q] = size(theta_min);
if (((p ~= 1)|(q ~= 1))|~isreal(theta_min))|ischar(theta_min)
    disp('SPHERE3D Error: theta_min must be scalar and real.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
end
[p,q] = size(theta_max);
if (((p ~= 1)|(q ~= 1))|~isreal(theta_max))|ischar(theta_max)
    disp('SPHERE3D Error: theta_max must be scalar and real.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
end
if theta_max <= theta_min
    disp('SPHERE3D Error: theta_max less than or equal theta_min.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
end
if abs(theta_max - theta_min) > 2*pi
    disp('SPHERE3D Error: range of theta greater than 2pi.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
end

[p,q] = size(phi_min);
if (((p ~= 1)|(q ~= 1))|~isreal(phi_min))|ischar(phi_min)
    disp('SPHERE3D Error: phi_min must be scalar and real.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
end
[p,q] = size(phi_max);
if (((p ~= 1)|(q ~= 1))|~isreal(phi_max))|ischar(phi_max)
    disp('SPHERE3D Error: phi_max must be scalar and real.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
end
if phi_max <= phi_min
    disp('SPHERE3D Error: phi_max less than or equal phi_min.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
end
if abs(phi_max - phi_min) > pi
    disp('SPHERE3D Error: range of phi greater than pi.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
end

[p,q] = size(Rho);
if (((p ~= 1)|(q ~= 1))|~isreal(Rho))|(ischar(Rho)|Rho < 0)
    disp('SPHERE3D Error: Rho must be scalar, positive and real.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
end

[p,q] = size(meshscale);
if (((p ~= 1)|(q ~= 1))|~isreal(meshscale))|ischar(meshscale)
    %disp('SPHERE3D Warning: mesh scale must be scalar and real.');
    meshscale = 1;
end
if (meshscale <= 0)
    %disp('SPHERE3D Warning: mesh scale must be scalar and positive.');
    meshscale = 1;
end

% Set up default plot and interpolation specifications.
str1 = 'mesh';
str2 = 'linear';
str3 = 1.0;
if length(varargin) == 3
% Sort out plot,interpolation and data scaling specification if three variables given.   
    str1 = [varargin{1}(:)]';
    str2 = [varargin{2}(:)]';
    str3 = [varargin{3}(:)]';
    g1 = (~isequal(str1,'mesh')&~isequal(str1,'surf'))&~isequal(str1,'off');
    g2 = (~isequal(str1,'meshc'));
    g5 = (~isequal(str1,'contour'));
    g3 = (~isequal(str2,'cubic')&~isequal(str2,'linear'));
    g4 = (~isequal(str2,'spline')&~isequal(str2,'nearest'));
    g6 = ~isnumeric(str3);
    if (g1&g2&g5)
        %disp('SPHERE3D Warning: Incorrect plot specification. Default to mesh plot.');
        str1 = 'mesh';
    end
    if (g3&g4)
        %disp('SPHERE3D Warning: Incorrect interpolation specification.'); 
        %disp('Default to linear interpolation.');
        str2 = 'linear';
    end  
    if (g6)
        %disp('SPHERE3D Warning: Incorrect data scaling factor.'); 
        disp('Default to Zscale = 1.');
        str3 = 1.0;
    end  
elseif length(varargin) == 2
% Sort out plot,interpolation or data scaling specification if two variables given.   
    str1 = [varargin{1}(:)]';
    str2 = [varargin{2}(:)]';
    g1 = (~isequal(str1,'mesh')&~isequal(str1,'surf'))&~isequal(str1,'off');
    g2 = (~isequal(str1,'meshc'));
    g3 = (~isequal(str1,'contour'));
    g4 = (~isequal(str1,'cubic')&~isequal(str1,'linear'));
    g5 = (~isequal(str1,'spline')&~isequal(str1,'nearest'));  
    if (g1&g2&g3&g4&g5)
        %disp('SPHERE3D Warning: Incorrect plot or interpolation specification.'); 
        %disp('Default to mesh plot and linear interpolation.');
        str1 = 'mesh';
        str2 = 'linear';
    end   
    g6 = (~isequal(str2,'cubic')&~isequal(str2,'linear'));
    g7 = (~isequal(str2,'spline')&~isequal(str2,'nearest'));  
    g8 = ~isnumeric(str2);
    if (g6&g7&g8)
        %disp('SPHERE3D Warning: Incorrect interpolation specification or data scaling factor.'); 
        %disp('Default to linear interpolation and Zscale = 1.');
        str2 = 'linear';
        str3 = 1.0;
    end
    temp2 = str2;
    if isequal(str1,'cubic')
        str2 = str1;
        str1 = 'mesh';
    elseif isequal(str1,'linear')
        str2 = str1;
        str1 = 'mesh';
    elseif isequal(str1,'spline')
        str2 = str1;
        str1 = 'mesh';
    elseif isequal(str1,'nearest')
        str2 = str1;
        str1 = 'mesh';   
    end 
    if isnumeric(temp2)
        str3 = temp2;
        if ~(g1&g2&g3)
            str2 = 'linear';
        else
            str1 = 'mesh';
        end
    end

elseif length(varargin) == 1
% Sort out plot, interpolation specification or data scaling factor from single variable input.    
    str1 = [varargin{1}(:)]';
    g1 = (~isequal(str1,'mesh')&~isequal(str1,'surf'))&~isequal(str1,'off');
    g2 = (~isequal(str1,'meshc'));
    g5 = (~isequal(str1,'contour'));
    g3 = (~isequal(str1,'cubic')&~isequal(str1,'linear'));
    g4 = (~isequal(str1,'spline')&~isequal(str1,'nearest'));  
    g6 = ~isnumeric(str1);
    if (g1&g2)&(g3&g4&g5&g6)
        disp('SPHERE3D Error: Incorrect plot,interpolation specification or data scaling factor.');
        Xout = []; Yout = []; Zout = []; Cmap = [];
        return
    elseif isequal(str1,'cubic')
        str2 = str1;
        str1 = 'mesh';
    elseif isequal(str1,'linear')
        str2 = str1;
        str1 = 'mesh';
    elseif isequal(str1,'spline')
        str2 = str1;
        str1 = 'mesh';
    elseif isequal(str1,'nearest')
        str2 = str1;
        str1 = 'mesh';
    elseif isequal(str1,'off')
        str2 = 'linear'; 
    elseif isnumeric(str1)
        str3 = str1;
        str1 = 'mesh';