Draw a circle for arbitrary orientation on spherical surface

I would like to draw a circular loop on spherical surface for a fixed orientation of theta(i.e. polar angle on sphere) and variying the azimuthal angle. The following code only generates the circle which are parallel to the equiatorial plane, but I need arbitrary orientation of ploar angle on the sphere. Pl somebody help me.

clear; clc;
N=10;
[X,Y,Z]=sphere(N);
C=zeros(N+1,N+1);
x=7;
y=1;
r=10;
for i=1:N+1
    for j=1:N+1
        d=sqrt(((i-x)^2)+((j-y)^2));
        if (d<=r)
            C(i,j)=1;
        end
    end
end
figure
surf(X,Y,Z,C)
axis equal

Answer 1

[Jim Riggs]

My approach is to define the circle around the X-axis, then rotate it into the desired position through a Yaw-Pitch rotation.

% Specify the user inputs
Rmag= 10;   % Sphere radius
radius = 1; % circle radius
psi = 30*pi/180;    % yaw rotation angle
theta = -45*pi/180; % pitch rotation angle (negative rotation is up)

% Define vectors and Calculate the YAW-PITCH transformation matrix
alpha = asin(radius/Rmag);  

YAW = [cos(psi), -sin(psi), 0; sin(psi), cos(psi), 0; 0, 0, 1]; % Planar YAW rotation
PITCH = [cos(theta), 0, sin(theta); 0,1,0; -sin(theta), 0, cos(theta)]; % Planar PITCH rotation
YP = YAW*PITCH; % YAW-PITCH rotation matrix

% Rc = Column Vector pointing to circle on X-Axis (start point)
Rc = [Rmag*cos(alpha); 0; radius];  
R = YP*[Rmag; 0; 0]; % Vector pointing to Circle center 

% Now sweep the Rc vector around the X-axis to generate the circle 
% This is done by adding a planar ROLL rotation to YP
clear C
C = [];  % C is the vector containing the circle X-Y-Z cordinates
for phi = 0:pi/50:2*pi
    ROLL = [1, 0, 0;  0, cos(phi), -sin(phi); 0, sin(phi), cos(phi) ]; 
    YPR = YP*ROLL;  % 3-dimentional transform
    Rnew = YPR*Rc;
    C = [C, Rnew];
end

% plot the result
[Sx,Sy,Sz]=sphere(50);

figure;
mesh(Rmag.*Sx, Rmag.*Sy, Rmag.*Sz);  % draw the sphere
hold on;
plot3([0, R(1)], [0, R(2)], [0,R(3)], 'r'); % Vector to Circle center
plot3(C(1,:), C(2,:),C(3,:), 'b'); % Circle 
axis equal

[Adam Danz]

+1; nice answer, Jim Riggs. I started on this yesterday and got sidetracked, then forgot about it.

Here are some minor fixes & suggestions

1. There's a variable name Radius that should be lower case; otherwise it throws an "unrecognized variable" error.
2. Replace * with .* in this line: YP*[Rmag, 0, 0], otherwise there should be an "incorrect dimensions" error.
3. clear C: no need for this.
4. preallocate the C matrix (see section of code below)
5. add terminal semicolon to supress output: Rnew = YPR*Rc
6. Instead of using phi as your loop variable and building the C matrix using C = [C, Rnew]; create a vector named phi and loop through indices of phi. That will allow you to build the C matrix using indexing.

phi = 0:pi/50:2*pi; 
C = nan(3, numel(phi));
for i = 1:numel(phi)
    ROLL = [1, 0, 0;  0, cos(phi(i)), -sin(phi(i)); 0, sin(phi(i)), cos(phi(i)) ]; 
    YPR = YP*ROLL;  % 3-dimentional transform
    Rnew = YPR*Rc;
    C(:,i) = Rnew;
end

[Jim Riggs]

Thanks Adam. I appreciate the comments.

(I do not have Matlab on this machine, so I can't actually test run my code when I make a post!)

[AVM]

@Jim & @ Adam: Lot of thanks to both of you. This is really amazing. The thing that I also need is to plot just one circle on the spherical surface. By using mesh command, the sphere itself get full of meshline so that the circular loop is not visible. Pl advise me, what command has to be modified so that I can get desired thing. I need a just one circular loop on a clear spherical surface.

[Jim Riggs]

1) Make sure that your circle is calculated corectly.

2) If you still can't see the circle, that means that it is inside the mesh.

Try making the size of the sphere slightly smaller, (or slightly increase Rmag when calculating the circle)

Then the circle should be visible.

[EDIT] Also, you can make the mesh transparent with the following;

h = mesh(Rmag.*Sx, Rmag.*Sy, Rmag.*Sz);  % draw the sphere
h.FaceColor = 'none';

or, you can make the face transparent by changing the "Alpha" value (1 = opaque (default) , <1 = transparent)

h.FaceAlpha = '0.85';

[AVM]

@Jim: Thanks a lot for your reply. It is appearing. It's really great.

[Adam Danz]

Also, make sure the circle isn't appearing on the other side of the sphere that is not visible from the view point.

There are two solutions to that

1) use view(az,el) to rotate the plot.

2) make the sphere partially transparent.