clf 
theta45=asin(sqrt(2)/sqrt(3))*180/pi;
sim_thetas=[0,45,theta45,45,90,90,90];
sim_phis=[0,0,45,90,0,45,90];

[xx,yy,zz]=sphere(120);
zz=zz*.7;
s=surf(xx,yy,zz);
s.LineStyle='none';
s.FaceColor='white';
hold on;

R = 1;
xo = 0;
yo = 0;
zo = 0;
num_angles = 3600;
latitudes = -90:15:90; % latitudes in degrees
longitudes = -180:30:180;
theta = linspace(0,2*pi,num_angles)';
% lat_angle = latitides*pi/180; 
x = cosd(latitudes);
z = .7 * sind(latitudes);
X = cos(theta) * x + xo;
Y = sin(theta) * x + yo;
Z = repmat(z, size(X,1),1) + zo;
%clf
plot3(X,Y,Z,'-k');
x=cos(theta);
y=0*theta;
z=.7*sin(theta);
%keyboard
for phi=longitudes
    X=rotz(phi)*([x,y,z].');
    plot3(X(1,:),X(2,:),X(3,:),'-k');
end
 %keyboard  

camlight

windmag=.2; %Length of wind vector
bmag = .2; %Length of B field vector
origindistance=2; %Distance from ellipse origin to wind/bmag schematic

for i=1:numel(sim_thetas)
    sim_theta=sim_thetas(i);
    sim_phi=sim_phis(i);

    RSim=rotz(sim_phi)*roty(sim_theta);
    sSim=RSim*[0;0;1]; %Transformation used to get spin axes
    windSim=[1;0;0];
    bSim=[0;1;0];
    
    wind=RSim.'*windSim;
    b=RSim.'*bSim;
    
    % Now calculate angle to rotate wind vector around so it is in xz plane
    phi=atan2(wind(2),wind(1))*180/pi;
    
    % Rotate wind and magnetic field vectors around z axis so that wind vector is in xz plane
    wind=rotz(-phi)*wind;
    b=rotz(-phi)*b;
    
    % Origin of wind/bfield schematic is in opposite direction of wind
    origin=-origindistance*wind;
    
    hold on;
    plot3(origin(1)+[0,wind(1)*windmag], origin(2)+[0,wind(2)*windmag], origin(3)+[0,wind(3)*windmag],'-k','LineWidth',2);
    
    % Now calculate elevation of wind vector (or azimuthal angle)
    windtheta=atan2(-wind(1),-wind(3))*180/pi;
    
    % Generate circle in xy plane at origindistance above ellipse (in z direction)
    theta=linspace(0,2*pi,500);
    z=origindistance+0*theta;%cos(theta)*bmag;
    y=sin(theta)*bmag;
    x=cos(theta)*bmag;
    
    % Rotate circle around y axis so it's origin coincides with wind/bfield schematic
    X=roty(windtheta)*[x;y;z];
    p=fill3(X(1,:),X(2,:),X(3,:),'g');
    p.FaceAlpha=.2;
    plot3(origin(1)+[0,b(1)*bmag],origin(2)+[0,b(2)*bmag],origin(3)+[0,b(3)*bmag],'r','LineWidth',2);
    %plot3(origin(1)+[0,0],origin(2)+[0,0],origin(3)+[0,bmag],'--k','LineWidth',2)
    scw=cross([0,0,1],wind);
    scw=scw./sqrt(sum(scw.^2))*bmag;
    plot3(origin(1)+[-scw(1),scw(1)],origin(2)+[-scw(2),scw(2)], origin(3)+[-scw(3),scw(3)],'--k');
    scw=cross(scw,wind);
    plot3(origin(1)+[-scw(1),scw(1)],origin(2)+[-scw(2),scw(2)], origin(3)+[-scw(3),scw(3)],'--k');
    plot3(origin(1)+[-scw(1),scw(1)],origin(2)+[-scw(2),scw(2)], origin(3)+[-scw(3),scw(3)],'--k');
end
axis equal;
