sphere3d.m

a usefull source code from internet

    end 
end;

% Check if data scaling factor is acceptable.
Zscale = str3;
[p,q] = size(Zscale);
if (((p ~= 1)|(q ~= 1))|~isreal(Zscale))|(ischar(Zscale)|Zscale < 0)
    disp('SPHERE3D Error: Zscale must be scalar, positive and real.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return
end

% Check if dimensions of input data are acceptable.
[r,c] = size(Zin');
if (r < 5)&(c < 5)
    disp('SPHERE3D Error: Input matrix dimensions must be greater than (4 x 4).');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return    
end

% Check if input data has two rows or columns or less.
if (r < 3)|(c < 3)
    disp('SPHERE3D Error: One or more input matrix dimensions too small.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return    
end

% Setup input magnitude matrix. 
Zin = flipud(Zin);
[r,c] = size(Zin);

% Check if meshscale is compatible with dimensions of input data.
scalefactor = round(max(r,c)/meshscale);
if scalefactor < 3
    disp('SPHERE3D Error: mesh scale incompatible with dimensions of input data.');
    Xout = []; Yout = []; Zout = []; Cmap = [];
    return    
end

% Set up meshgrid corresponding to larger matrix dimension of Zin
% for interpolation if required.
n = meshscale;
if r > c
    L = r;
    L2 = fix(L/n)*n;
    step = r/(c-1);
    [X1,Y1] = meshgrid(0:step:r,1:r);
    if n < 1
        [X,Y] = meshgrid(0:n:(L2-n),0:n:(L2-n));
    else
        [X,Y] = meshgrid(1:n:L2,1:n:L2); 
    end
    T = interp2(X1,Y1,Zin,X,Y,str2);
elseif c > r
    L = c;
    L2 = fix(L/n)*n;
    step = c/(r-1);
    [X1,Y1] = meshgrid(1:c,0:step:c);
    if n < 1
        [X,Y] = meshgrid(0:n:(L2-n),0:n:(L2-n));
    else
        [X,Y] = meshgrid(1:n:L2,1:n:L2); 
    end
    T = interp2(X1,Y1,Zin,X,Y,str2);
else
    L = r;
    L2 = fix(L/n)*n;
    [X1,Y1] = meshgrid(1:r,1:r);
    if n < 1
        [X,Y] = meshgrid(0:n:(L2-n),0:n:(L2-n));
    else
        [X,Y] = meshgrid(1:n:L2,1:n:L2); 
    end
    T = interp2(X1,Y1,Zin,X,Y,str2);
end

[p,q] = size(T);
L2 = max(p,q);

% Set up angles theta
angl = theta_min:abs(theta_max-theta_min)/(L2-1):theta_max;
for j = 1:L2  
    theta(j,:) = ones([1 L2])*angl(j);
end
theta = theta';

% Set up angles phi
angl2 = phi_min:abs(phi_max-phi_min)/(L2-1):phi_max;
for j = 1:L2  
    phi(j,:) = ones([1 L2])*angl2(j);
end


% Set up autoscaling if required.
% In some interpolation methods NaNs are returned; these are tracked and
% replaced by 0's, but may not always be sufficient. Try a different
% interpolation method.
Zscale = str3;
[ind] = find(isnan(T));
T(ind) = 0;
Zavg = mean(mean(T));
MaxZin = max(max(T-Zavg));
MinZin = min(min(T-Zavg));
Rho_max = Rho + Zscale*MaxZin;
Rho_min = Rho + Zscale*MinZin;
if Rho_min < 0
    Zscale = -Rho/MinZin;
    Rho_min = Rho + Zscale*MinZin;
    %disp(['SPHERE3D Warning: Data scaling changed to ',num2str(Zscale),...
    %    ' to avoid negative radial values.']);
end

% Convert to Cartesian coordinates
for j = 1:L2  
    [Xout(j,:) Yout(j,:) Zout(j,:)] = ...
    sph2cart(theta(j,:),phi(j,:),Rho*ones([1 L2])+Zscale*T(j,:));
end

% Set up color mapping for output
nColor = 50;
[Cmap,Rmax,Rmin] = RadialColorMap(Xout,Yout,Zout,nColor); 

% Plot the data
switch str1;
case 'mesh'   
    colormap([0 0 0]);
    mesh(Xout,Yout,Zout);
    axis vis3d
    axis off;
    grid off;
    set(gca,'DataAspectRatio',[1 1 1]);
    set(gca,'XDir','rev','YDir','rev');
    
case 'meshc'
    Zoutc = Zout + 0.3*max(max(Zout));
    colormap([0 0 0]);
    meshc(Xout,Yout,Zoutc);
    axis off;
    grid off;
    set(gca,'DataAspectRatio',[1 1 1])    
    set(gca,'XDir','rev','YDir','rev'); 
    
case 'surf'   
    surf(Xout,Yout,Zout,Cmap);
    set(gca,'DataAspectRatio',[1 1 1]);
    axis vis3d
    axis off;
    grid off;
    set(gca,'XDir','rev','YDir','rev');
    cbar = colorbar;
    numYTicks = length(get(cbar,'YTickLabel'));
    gap = (MaxZin - MinZin)/(numYTicks-1);
    k = 0:numYTicks-1;
    if gap < 5
        minLabel = round(MinZin*100)/100;
        stepLabel = k.*round(gap*100)/100 + minLabel;
    else
        minLabel = round(MinZin);
        stepLabel = k.*round(gap) + minLabel;
    end    
    colorbar('YTickLabel',stepLabel); 
    
case 'contour' 
    C = contourc(T',25);
    [row,col] = size(C);
    mag(1) = C(1,1);    % magnitude at nLevel = 1
    num = C(2,1);       % number of (x,y)contour pairs at nLevel = 1
    finished = 0;       %
    sumPairs = 0;       % total number of (x,y) contour pairs at all levels
    nLevel = 1;         % nth contour level

    % retrieve contour coordinates (x,y) from contourc function output
    while ~finished
        x(1:num,nLevel) = C(1,2:num+1);
        y(1:num,nLevel) = C(2,2:num+1);
        C(1,:) = shift(C(1,:),-(num+1));
        C(2,:) = shift(C(2,:),-(num+1));
        sumPairs = sumPairs + num + 1;
        nLevel = nLevel + 1;
        num = C(2,1);
        mag(nLevel) = C(1,1);
        if sumPairs >= col
            finished = 1;
        end
    end
    % Replace zeros with NaNs in x and y
    [ind] = find(x == 0);
    x(ind) = NaN;
    [ind] = find(y == 0);
    y(ind) = NaN;

    % set up equivalent phi angles from contour x and radius Rho
    alpha = (x/Rho);
    minalpha = min(min(alpha));
    alpha = alpha - minalpha;
    maxalpha = max(max(alpha));
    alpha = alpha*(phi_max - phi_min)/maxalpha + phi_min;
    minalpha = min(min(alpha));
    maxalpha = max(max(alpha));
    
    % set up equivalent theta angles from contour y and radius Rho
    gamma = (y/Rho);
    mingamma = min(min(gamma));
    gamma = gamma - mingamma;
    maxgamma = max(max(gamma));
    gamma = gamma*(theta_max - theta_min)/maxgamma + theta_min;
    mingamma = min(min(gamma));
    maxgamma = max(max(gamma));
    
    % convert to cartesian
    Xsph = Rho*cos(alpha).*cos(gamma); 
    Ysph = Rho*cos(alpha).*sin(gamma);
    Zsph = Rho*sin(alpha);

    % draw contour lines one at a time with preset colour in rgbCol
%     figure
    nCol = (nLevel-2:-1:0)'/nLevel;
    rgbCol = hsv2rgb(0.8*[nCol ones(nLevel-1,2)]);
    maxZ = max(max(Zin));
    minZ = min(min(Zin));
    caxis([minZ, maxZ]);
    caxis manual
    colorbar
    [r,c] = size(Xsph);
    for i = 1:nLevel-1
        xx(1:r) = Xsph(:,i);
        yy(1:r) = Ysph(:,i);
        zz(1:r) = Zsph(:,i);
        line(xx,yy,zz,'Color',[rgbCol(i,1) rgbCol(i,2) rgbCol(i,3)]);
    end
    xmin = -Rho; xmax = Rho;
    ymin = -Rho; ymax = Rho;
    zmin = -Rho; zmax = Rho;
    axis([xmin xmax ymin ymax zmin zmax])
    axis vis3d
    set(gca,'DataAspectRatio',[1 1 1])
    set(gca,'XDir','rev','YDir','rev')
    axis off;
    grid off;

% set up sphere plot borders
    nPoints = 80;
    j = 1:nPoints;
    theta1 = theta_min + (j-1)*(theta_max - theta_min)/(nPoints - 1);
    phi1 = phi_min + (j-1)*(phi_max - phi_min)/(nPoints - 1);
    XborderRight = Rho*cos(phi1)*cos(theta_max); 
    YborderRight = Rho*cos(phi1)*sin(theta_max); 
    ZborderRight = Rho*sin(phi1); 
    XborderLeft = Rho*cos(phi1)*cos(theta_min); 
    YborderLeft = Rho*cos(phi1)*sin(theta_min); 
    ZborderLeft = Rho*sin(phi1);
    XborderTop = Rho*cos(phi_max)*cos(theta1); 
    YborderTop = Rho*cos(phi_max)*sin(theta1); 
    ZborderTop = Rho*ones(size(phi1))*sin(phi_max);
    XborderBot = Rho*cos(phi_min)*cos(theta1); 
    YborderBot = Rho*cos(phi_min)*sin(theta1); 
    ZborderBot = Rho*ones(size(phi1))*sin(phi_min);

% plot borders on sphere
    BorderLine(XborderRight,YborderRight,ZborderRight,Rho);
    BorderLine(XborderTop,YborderTop,ZborderTop,Rho);
    BorderLine(XborderBot,YborderBot,ZborderBot,Rho); 
    BorderLine(XborderLeft,YborderLeft,ZborderLeft,Rho);
end %switch

%
% Local Functions
%
    function BorderLine(Xborder,Yborder,Zborder,Rho)
        % Draws a border line defined by (Xborder,Yborder,Zborder).
        % Used to plot borders of contour on spherical surface.
        %
        hold on
        line(Xborder,Yborder,Zborder,'Color',[.5 .5 .5]);
        xlo = -Rho; xhi = Rho;
        ylo = -Rho; yhi = Rho;
        zlo = -Rho; zhi = Rho;
        axis([xlo xhi ylo yhi zlo zhi]);
        axis vis3d;
        set(gca,'XDir','rev','YDir','rev');
        axis off;
        grid off;
        hold off;
        return

    function [A] = shift(A,n)
        %
        % SHIFT     moves each element down by n steps along vector A,
        %           treating the positions j as though they are on a circle,
        %           so that A(j) moves to position(j+n), and A(L-n+1:L)occupy
        %           positions(1:n), where L is the number of elements in A.
        %           If n is positive the shift is in the forward direction
        %           whereas if negative, the shift is backwards.
        %
        % Usage:    [A] = shift(A,n)
        %
        % Input:    A,  array whose elements will be shifted
        %           n,  number of shifts to perform
        %               n can be negative as well.
        %
        % Output:   A,  array with shifted elements.
        %
        change = 0;
        [r,c] = size(A);
        if (r > 1)&(c == 1)
            A = A'; change = 1;
        end
        [r,c] = size(A);
        if (r == 1)&(c >= 1)
            [r,c] = size(n);
            if (r == 1)&(c == 1)
                n = round(n);
                L = length(A);
                n = rem(n,L);
                if n < 0 n = L + n; end
                T = A;
                A(n+1:L) = T(1:L-n);
                A(1:n) = T(L-n+1:L);
            else
                %disp('Shift warning: number of shifts cannot be a vector or matrix');
            end
        else
            %disp('Shift warning: input must be row or column vector');
        end
        if change == 1 A = A'; end;
        return      
        
   function [Cmap,Rmax,Rmin] = RadialColorMap(X,Y,Z,nColor)
       % RadialColorMap     maps a colour scale to the radial information
       %                    formed from coordinates X, and Y. The number of
       %                    colors formed is nColor.
       %
       % Usage              [Cmap,Rmax,Rmin] = RadialColorMap(X,Y,nColor)
       %
       % Input              X       a square matrix of X-coordinates
       %                    Y       a square matrix of Y-coordinates
       %                    Z       a square matrix of Z-coordinates
       %                    nColor  number of colours to use
       %
       % Output             Cmap    the colour map
       %                    Rmax    maximum radius of input data
       %                    Rmin    minimum radius of input data
       %
       R = sqrt(X.^2 + Y.^2 + Z.^2);
       [r,c] = size(R);
       Rmax = max(max(R));
       Rmin = min(min(R));
       Col = (Rmax - Rmin)/(nColor-1);
       for i = 1:r
           for j = 1:r
               for k = 1:nColor
                   if (R(i,j)>= Rmin + (k-1)*Col)&(R(i,j) < Rmin + k*Col)
                       Cmap(i,j) = Rmin + (k-1)*Col;
                   end
               end
           end
       end 
       return
            
return