
function pflag = plotBox(origin,vec1,vec2,vec3,a,b,c)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Function to plot a box of dimensions a, b, c
% centered at a given origin with edges defined
% by frame vectors vec1, vec2, vec3.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%set colors and line width for plot

color1 = 'g' ; color2 = 'r' ; color3 = 'b' ; color4 = 'k' ;

lwidth = 2 ;


%draw frame vectors as colored line segments

x(1) = origin(1) ; y(1) = origin(2) ; z(1) = origin(3) ;
x(2) = origin(1) + (a/2)*vec1(1) ; 
y(2) = origin(2) + (a/2)*vec1(2) ; 
z(2) = origin(3) + (a/2)*vec1(3) ;
line(x,y,z,'Color',color1,'LineWidth',lwidth)

x(1) = origin(1) ; y(1) = origin(2) ; z(1) = origin(3) ;
x(2) = origin(1) + (b/2)*vec2(1) ; 
y(2) = origin(2) + (b/2)*vec2(2) ; 
z(2) = origin(3) + (b/2)*vec2(3) ;
line(x,y,z,'Color',color2,'LineWidth',lwidth)

x(1) = origin(1) ; y(1) = origin(2) ; z(1) = origin(3) ;
x(2) = origin(1) + (c/2)*vec3(1) ; 
y(2) = origin(2) + (c/2)*vec3(2) ; 
z(2) = origin(3) + (c/2)*vec3(3) ;
line(x,y,z,'Color',color3,'LineWidth',lwidth)


%draw box face: x = a/2.

x(1) = origin(1) + (a/2)*vec1(1) + (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(1) = origin(2) + (a/2)*vec1(2) + (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(1) = origin(3) + (a/2)*vec1(3) + (b/2)*vec2(3) - (c/2)*vec3(3) ;
x(2) = origin(1) + (a/2)*vec1(1) - (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(2) = origin(2) + (a/2)*vec1(2) - (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(2) = origin(3) + (a/2)*vec1(3) - (b/2)*vec2(3) - (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)

x(1) = origin(1) + (a/2)*vec1(1) - (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(1) = origin(2) + (a/2)*vec1(2) - (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(1) = origin(3) + (a/2)*vec1(3) - (b/2)*vec2(3) - (c/2)*vec3(3) ;
x(2) = origin(1) + (a/2)*vec1(1) - (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(2) = origin(2) + (a/2)*vec1(2) - (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(2) = origin(3) + (a/2)*vec1(3) - (b/2)*vec2(3) + (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)

x(1) = origin(1) + (a/2)*vec1(1) - (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(1) = origin(2) + (a/2)*vec1(2) - (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(1) = origin(3) + (a/2)*vec1(3) - (b/2)*vec2(3) + (c/2)*vec3(3) ;
x(2) = origin(1) + (a/2)*vec1(1) + (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(2) = origin(2) + (a/2)*vec1(2) + (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(2) = origin(3) + (a/2)*vec1(3) + (b/2)*vec2(3) + (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)

x(1) = origin(1) + (a/2)*vec1(1) + (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(1) = origin(2) + (a/2)*vec1(2) + (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(1) = origin(3) + (a/2)*vec1(3) + (b/2)*vec2(3) + (c/2)*vec3(3) ;
x(2) = origin(1) + (a/2)*vec1(1) + (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(2) = origin(2) + (a/2)*vec1(2) + (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(2) = origin(3) + (a/2)*vec1(3) + (b/2)*vec2(3) - (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)


%draw box face: x = -a/2.

x(1) = origin(1) - (a/2)*vec1(1) + (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(1) = origin(2) - (a/2)*vec1(2) + (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(1) = origin(3) - (a/2)*vec1(3) + (b/2)*vec2(3) - (c/2)*vec3(3) ;
x(2) = origin(1) - (a/2)*vec1(1) - (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(2) = origin(2) - (a/2)*vec1(2) - (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(2) = origin(3) - (a/2)*vec1(3) - (b/2)*vec2(3) - (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)

x(1) = origin(1) - (a/2)*vec1(1) - (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(1) = origin(2) - (a/2)*vec1(2) - (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(1) = origin(3) - (a/2)*vec1(3) - (b/2)*vec2(3) - (c/2)*vec3(3) ;
x(2) = origin(1) - (a/2)*vec1(1) - (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(2) = origin(2) - (a/2)*vec1(2) - (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(2) = origin(3) - (a/2)*vec1(3) - (b/2)*vec2(3) + (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)

x(1) = origin(1) - (a/2)*vec1(1) - (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(1) = origin(2) - (a/2)*vec1(2) - (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(1) = origin(3) - (a/2)*vec1(3) - (b/2)*vec2(3) + (c/2)*vec3(3) ;
x(2) = origin(1) - (a/2)*vec1(1) + (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(2) = origin(2) - (a/2)*vec1(2) + (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(2) = origin(3) - (a/2)*vec1(3) + (b/2)*vec2(3) + (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)

x(1) = origin(1) - (a/2)*vec1(1) + (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(1) = origin(2) - (a/2)*vec1(2) + (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(1) = origin(3) - (a/2)*vec1(3) + (b/2)*vec2(3) + (c/2)*vec3(3) ;
x(2) = origin(1) - (a/2)*vec1(1) + (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(2) = origin(2) - (a/2)*vec1(2) + (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(2) = origin(3) - (a/2)*vec1(3) + (b/2)*vec2(3) - (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)

%finish box face: y = b/2.

x(1) = origin(1) + (a/2)*vec1(1) + (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(1) = origin(2) + (a/2)*vec1(2) + (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(1) = origin(3) + (a/2)*vec1(3) + (b/2)*vec2(3) + (c/2)*vec3(3) ;
x(2) = origin(1) - (a/2)*vec1(1) + (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(2) = origin(2) - (a/2)*vec1(2) + (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(2) = origin(3) - (a/2)*vec1(3) + (b/2)*vec2(3) + (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)


x(1) = origin(1) + (a/2)*vec1(1) + (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(1) = origin(2) + (a/2)*vec1(2) + (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(1) = origin(3) + (a/2)*vec1(3) + (b/2)*vec2(3) - (c/2)*vec3(3) ;
x(2) = origin(1) - (a/2)*vec1(1) + (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(2) = origin(2) - (a/2)*vec1(2) + (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(2) = origin(3) - (a/2)*vec1(3) + (b/2)*vec2(3) - (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)


%finish box face: y = -b/2.

x(1) = origin(1) + (a/2)*vec1(1) - (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(1) = origin(2) + (a/2)*vec1(2) - (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(1) = origin(3) + (a/2)*vec1(3) - (b/2)*vec2(3) + (c/2)*vec3(3) ;
x(2) = origin(1) - (a/2)*vec1(1) - (b/2)*vec2(1) + (c/2)*vec3(1) ;
y(2) = origin(2) - (a/2)*vec1(2) - (b/2)*vec2(2) + (c/2)*vec3(2) ;
z(2) = origin(3) - (a/2)*vec1(3) - (b/2)*vec2(3) + (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)


x(1) = origin(1) + (a/2)*vec1(1) - (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(1) = origin(2) + (a/2)*vec1(2) - (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(1) = origin(3) + (a/2)*vec1(3) - (b/2)*vec2(3) - (c/2)*vec3(3) ;
x(2) = origin(1) - (a/2)*vec1(1) - (b/2)*vec2(1) - (c/2)*vec3(1) ;
y(2) = origin(2) - (a/2)*vec1(2) - (b/2)*vec2(2) - (c/2)*vec3(2) ;
z(2) = origin(3) - (a/2)*vec1(3) - (b/2)*vec2(3) - (c/2)*vec3(3) ;
line(x,y,z,'Color',color4,'LineWidth',lwidth)


%box faces z = c/2 and z = -c/2 are automatically complete


pflag=1 ;   %return value to indicate successful completion

