% visualize direction cosines and pattern cuts
% modified to generate chart for radar notes (ver5.1)
clear
rad=pi/180;

% pick a point (theta,phi)
theta=input('theta (deg): ');
thetar=theta*rad;
phi=input('phi (deg): ');
phir=phi*rad;
u0=sin(thetar)*cos(phir); v0=sin(thetar)*sin(phir); w0=cos(thetar);
figure(1)
clf
plot3(u0,v0,w0,'ko')
hold on

% now plot spherical triangle solution
sa=u0; sb=v0;
cx=cos(thetar)/sqrt(1-sa^2); ss=cos(thetar)/sqrt(1-sb^2);
chi1=acos(cx); psi1=asin(ss);
%plot3(0,sin(chi1),cos(chi1),'bs')
%plot3(cos(psi1),0,sin(psi1),'bd')
% transformation matrix and great circle cut for chi1
Tchi=[1 0 0; 0 cos(chi1) -sin(chi1); 0 sin(chi1) cos(chi1)];
% cut coordinates in primed system
thr1=[0:1:90]*rad;
phr1=0;
U1=sin(thr1)*cos(phr1); V1=sin(thr1)*sin(phr1); W1=cos(thr1);
A=[U1;V1;W1];
T=transpose(Tchi);
B=T*A; U=B(1,:); V=B(2,:); W=B(3,:);
% plot3(U,V,W,'m-')
% transformation matrix and great circle cut for psi1
Tpsi=[sin(psi1) 0 -cos(psi1); 0 1 0; cos(psi1) 0 sin(psi1)];
thr1=[0:1:90]*rad;
phr1=pi/2;
U1=sin(thr1)*cos(phr1); V1=sin(thr1)*sin(phr1); W1=cos(thr1);
A=[U1;V1;W1];
T=transpose(Tpsi);
B=T*A; U=B(1,:); V=B(2,:); W=B(3,:);
% plot3(U,V,W,'g-')

%legend(['\theta = ',num2str(theta),'^o, \phi = ',num2str(phi),'^o'],...
%    '\chi sph tri','\psi sph tri','\chi rot about x','\psi rot about y')

% plot the octant
thr=[0:1:90]*rad;
st=sin(thr); ct=cos(thr);
phr=thr;
sp=st; cp=ct;
N=length(thr);
% phi=0, x-z axis
x=ct; y=zeros(size(x)); z=st;
    plot3(x,y,z),hold on
% phi=90, y-z axis
y=ct; z=st; x=zeros(size(y));
    plot3(x,y,z)
% theta=90
x=cp; y=sp; z=zeros(size(x));
    plot3(x,y,z)
% draw the axes
plot3([0 1], [0 0], [0 0], 'k-','LineWidth',2)
plot3([0 0], [0 1], [0 0], 'k-','LineWidth',2)
plot3([0 0], [0 0], [0 1], 'k-','LineWidth',2)

plot3([0 u0], [0 v0], [0 w0],'r-', 'LineWidth',2)
view(120,20)
% polar arc
U=st*cos(phir); V=st*sin(phir); W=ct;
plot3(U,V,W,'LineWidth',2)
% azimuth arc
x=sin(thetar)*cp; y=sin(thetar)*sp; z=w0*ones(size(cp));
plot3(x,y,z,'LineWidth',2)
% draw the projection box
plot3([0 0], [v0 v0], [0 w0], '--')
plot3([u0 u0], [v0 v0], [0 w0], '--')
plot3([u0 u0], [0 0], [0 w0], '--')
plot3([u0 u0], [0 v0], [0 0], '--')
plot3([0 u0], [v0 v0], [0 0], '--')
plot3([u0 u0], [0 v0], [w0 w0], '--')
plot3([0 u0], [v0 v0], [w0 w0], '--')
plot3([u0 0], [0 0], [w0 w0], '--')
plot3([0 0], [v0 0], [w0 w0], '--')
plot3([0 cos(phir)], [0 sin(phir)], [0 0], '--')
% plot3([0 0], [0 v0], [0 w0], 'g--')
% plot3([0 u0], [0 0], [0 w0], 'g--')
text(1.05,0,0,'x')
text(0,1.05,0,'y')
text(0,0,1.05,'z')


axis equal
axis square
axis off
hold off

