% crossproduct.m

in1 = input('Enter p/t of vector 1 as [p t]: ');

in2 = input('Enter p/t of vector 2 as [p t]: ');

V1 = [ cos(in1(1)*pi/180)*cos(in1(2)*pi/180)
       cos(in1(1)*pi/180)*sin(in1(2)*pi/180)
       sin(in1(1)*pi/180) ];

V2 = [ cos(in2(1)*pi/180)*cos(in2(2)*pi/180)
       cos(in2(1)*pi/180)*sin(in2(2)*pi/180)
       sin(in2(1)*pi/180) ];

C = cross(V1,V2)

mag = sqrt(C(1)^2 + C(2)^2 + C(3)^2);
C(1) = C(1)/mag;
C(2) = C(2)/mag;
C(3) = C(3)/mag;

plunge = asin(C(3))*180/pi;
trend  = atan2(C(2),C(1))*180/pi;
if (trend < 0) 
   trend = trend + 360;
end;
if (plunge < 0)
   plunge = -plunge;
   trend = trend + 180;
   if (trend > 360)
      trend = trend - 360;
   end
end   
plunge
trend