procedure circle3d { xc yc zc nx ny nz radius }
{
expression mag {pow((nx*nx+ny*ny+nz*nz),0.5)}
expression nx {nx/mag}
expression ny {ny/mag}
expression nz {nz/mag}
if_then_else {abs(nx) > abs(ny)}
{
expression ux {-nz}
expression uy {0}
expression uz {nx}
expression mag {pow((ux*ux+uy*uy+uz*uz),0.5)}
expression ux {ux/mag}
expression uy {uy/mag}
expression uz {uz/mag}
expression vx {ny*uz - uy*nz}
expression vy {-(nx*uz - ux*nz)}
expression vz {nx*uy - ux*ny}
}
{
expression ux {0}
expression uy {nz}
expression uz {-ny}
expression mag {pow((ux*ux+uy*uy+uz*uz),0.5)}
expression ux {ux/mag}
expression uy {uy/mag}
expression uz {uz/mag}
expression vx {ny*uz - uy*nz}
expression vy {-(nx*uz - ux*nz)}
expression vz {nx*uy - ux*ny}
}
expression mag {pow((vx*vx+vy*vy+vz*vz),0.5)}
expression vx {vx/mag}
expression vy {vy/mag}
expression vz {vz/mag}
ang3d_draw_parametric_curve t
   {0; t<=6.28+0.005; t+0.05}
   {xc+radius*cos(t)*ux+radius*sin(t)*vx; yc+radius*cos(t)*uy+radius*sin(t)*vy; zc+radius*cos(t)*uz+radius*sin(t)*vz}
}


procedure arc3d { xc yc zc ux uy uz vx vy vz radius}
{
expression mag {pow((ux*ux+uy*uy+uz*uz),0.5)}
expression ux {ux/mag}
expression uy {uy/mag}
expression uz {uz/mag}
expression mag {pow((vx*vx+vy*vy+vz*vz),0.5)}
expression vx {vx/mag}
expression vy {vy/mag}
expression vz {vz/mag}
expression theta_stop {acos(ux*vx+uy*vy+uz*vz)}
expression v_dot_u {vx*ux+vy*uy+vz*uz}
expression v_perp_x {vx-(v_dot_u*ux)}
expression v_perp_y {vy-(v_dot_u*uy)}
expression v_perp_z {vz-(v_dot_u*uz)}
expression vx {v_perp_x}
expression vy {v_perp_y}
expression vz {v_perp_z}
expression mag {pow((vx*vx+vy*vy+vz*vz),0.5)}
expression vx {vx/mag}
expression vy {vy/mag}
expression vz {vz/mag}
ang3d_draw_parametric_curve t
   {0; t<=theta_stop; t+0.005}
   {xc+radius*cos(t)*ux+radius*sin(t)*vx; yc+radius*cos(t)*uy+radius*sin(t)*vy; zc+radius*cos(t)*uz+radius*sin(t)*vz}
}

% Setup the 3D Cartesian System & View
ang3d_picture 0 0 100 100
ang3d_origin 60 50 0.75 0.75
ang3d_unit 5
ang3d_scale 1 1 1
ang3d_axes_drawing_range -12 12 -12 12 -12 12
%ang3d_drawsystem_p 2 10 10 10 1 2

% Draw the 1st Circle FIXED
number x0 0		% center of circle
number y0 0
number z0 0
number n1 0		% normal to circle
number n2 0
number n3 1
color 255 0 0
call circle3d { x0 y0 z0 n1 n2 n3 5 }

% Draw the 2nd Circle INPUT
number n1 0		% normal to circle
number n2 1
number n3 0
call circle3d { x0 y0 z0 n1 n2 n3 5 }

% Draw the 3rd Circle OUTPUT
number n1 1		% normal to circle
number n2 0
number n3 0
call circle3d { x0 y0 z0 n1 n2 n3 5 }
color 0 0 0

% Draw the 1st Arc FIXED LINK
% Specify a vector u from the center towards the start of the arc
number u1 0		
number u2 1
number u3 0
% Specify a vector v from the center towards the end of the arc
% The magnitudes of u & v don't matter since the radius of the 
% arc is explicitly prescribed in the call to the arc3d routine
number v1 1
number v2 0
number v3 0
double
call arc3d { x0 y0 z0 u1 u2 u3 v1 v2 v3 5 }
normal

% Draw the 2nd Arc INPUT LINK
number u1 1		% start of arc
number u2 0
number u3 0
number v1 1		% end of arc
number v2 0
number v3 1
double
call arc3d { x0 y0 z0 u1 u2 u3 v1 v2 v3 5 }
normal

% Draw the 3rd Arc COUPLER LINK
number u1 1	% start of arc
number u2 0
number u3 1
number v1 0		% end of arc
number v2 1
number v3 2
% Draw the 4th Circle COUPLER
expression n1 {u2*v3-u3*v2}	% normal to circle
expression n2 {-1*(u1*v3-v1*u3)}
expression n3 {u1*v2-v1*u2}
color 255 0 0
call circle3d { x0 y0 z0 n1 n2 n3 5 }
color 0 0 0
double
call arc3d { x0 y0 z0 u1 u2 u3 v1 v2 v3 5 }
normal

% Draw the 3rd Arc OUTPUT LINK
number u1 0	% start of arc
number u2 1
number u3 2
number v1 0		% end of arc
number v2 1
number v3 0
double
call arc3d { x0 y0 z0 u1 u2 u3 v1 v2 v3 5 }
normal

ang3d_point center x0 y0 z0
ang3d_point O 7.5 0 0
ang3d_point C 0 7.5 0
ang3d_point in_mov_axis 1 0 1
ang3d_point out_mov_axis 0 1 2
expression mag {pow((1*1+0*0+1*1),0.5)}
expression in_x {7.5*1/mag}
expression in_y {7.5*0/mag}
expression in_z {7.5*1/mag}
ang3d_point A in_x in_y in_z
expression mag {pow((0*0+1*1+2*2),0.5)}
expression out_x {7.5*0/mag}
expression out_y {7.5*1/mag}
expression out_z {7.5*2/mag}
ang3d_point B out_x out_y out_z
drawvector center O
drawvector center C
drawvector center A 
drawvector center B

mark_l O
mark_r C 
mark_l A
mark_r B
ang3d_point fixed_label 2.9 4.8 0.0
printat fixed_label { \gamma }
ang3d_point input_label 5.8 0.9 4.0
printat input_label { \alpha }
ang3d_point coupler_label 6 6 10.9
printat coupler_label { \eta }
ang3d_point output_label 1.3 5.8 6.0
printat output_label { \beta }