function calcTranslateStep(direction) {
    // calculate the translateStep along each direction
    // direction parameter is the direction choosen by the user
    // to translate the box
    
    //  3------1  7------5
    //  |      |  |	 |
    //  |      |  |	 |
    //  4------2  8------6
    //  => x <=> 7--3
    //  => y <=> 3--1
    //  => z <=> 4--3
  
    if (direction=='x') {
    
     // first get the vector along x direction of user box (7-->3)
        vX = V3.x - V7.x
	vY = V3.y - V7.y 
        vZ = V3.z - V7.z } // end of if direction=='x'
    
    

    if (direction=='y') {
    
     // first get the vector along y direction of user box (3-->1)
        vX = V1.x - V3.x
	vY = V1.y - V3.y 
        vZ = V1.z - V3.z } // end of if direction=='y'
    
                        
    
    if (direction=='z') {
    
     // first get the vector along z direction of user box (4-->3)
        vX = V3.x - V4.x
	vY = V3.y - V4.y 
        vZ = V3.z - V4.z  } // end of if direction=='z'


    normv = sqrt(vX*vX + vY*vY + vZ*vZ);
    vX = vX/normv
    vY = vY/normv
    vZ = vZ/normv
    	 
    translateStepX = translateStep * vX   
    translateStepY = translateStep * vY   
    translateStepZ = translateStep * vZ   
    
     
    }// end of phiAxis

    function drawPhiAxis(ON) {
      // draw or hide the phi axis; V1 must be defined before
      if ON {
          vect1 = V1;
	  vect2 = point(0,0,1);
	  u = cross(vect1,vect2);
	  normU = sqrt(u.x*u.x + u.y*u.y + u.z*u.z);
          u = u/normU;
          u = @{u*rho};
  	  pt9  =point(bbic.x-u.x,bbic.y-u.y, bbic.z-u.z);
	  pt10 =point(bbic.x+u.x,bbic.y+u.y, bbic.z+u.z); 	  
	  draw phiAxis  @pt9 @pt10;
	  phiAxisDisplay=true;
	    }
      else  {draw phiAxis  OFF; phiAxisDisplay=false;}	    

                             }// end of phiAxis
// -----------------------------------------------------------

    function calculateRhoTetaPhi(xx, yy, zz) {
      // calculate rho teta phi from new reference vector (xx,yy,zz)
      // angles in degrees here
	rho = sqrt(xx*xx + yy*yy + zz*zz);
	
	
	if (yy>=0) {teta= acos(xx/sqrt(xx*xx+yy*yy))} else {teta= 360.0 - acos(xx/sqrt(xx*xx+yy*yy))};
         
        
        // 0 < phi < 180 in our system
	phi = acos(zz/rho); 
	                                 
					 } // end of calculateRhoTetaPhi
					 


// -----------------------------------------------------------

    function setTranslationStep(boxVect) {

      diagonal=boxVect.distance({0 0 0});
      translateStep=translateRatio*diagonal;;
     }

// -----------------------------------------------------------

// -----------------------------------------------------------

    function setScaleStep() {

scaleStep=scaleRatio*0.1;
     }

// -----------------------------------------------------------




    function updateBoxVerticesForAScale(sign) {

//  update the 8 vertices by shifting them along (pt5, pt1) vector

Vx = (V1.x - V5.x)*sign*scaleStep;
Vy = (V1.y - V5.y)*sign*scaleStep;
Vz = (V1.z - V5.z)*sign*scaleStep;

// ********** V1,V5 **********
V1x = V1.x + Vx;
V1y = V1.y + Vy;
V1z = V1.z + Vz;

V1 = point(V1x, V1y, V1z);

V5x = V5.x - Vx;
V5y = V5.y - Vy;
V5z = V5.z - Vz;

V5 = point(V5x, V5y, V5z);

// ********** V3, V7 **********
V3x = V3.x + Vx;
V3y = V3.y + Vy;
V3z = V3.z + Vz;

V3 = point(V3x, V3y, V3z);

V7x = V7.x - Vx;
V7y = V7.y - Vy;
V7z = V7.z - Vz;

V7 = point(V7x, V7y, V7z);

// ********** V2,6 **********
V2x = V2.x + Vx;
V2y = V2.y + Vy;
V2z = V2.z + Vz;

V2 = point(V2x, V2y, V2z);

V6x = V6.x - Vx;
V6y = V6.y - Vy;
V6z = V6.z - Vz;

V6 = point(V6x, V6y, V6z);

// ********** V4,8 **********
V4x = V4.x + Vx;
V4y = V4.y + Vy;
V4z = V4.z + Vz;

V4 = point(V4x, V4y, V4z);

V8x = V8.x - Vx;
V8y = V8.y - Vy;
V8z = V8.z - Vz;

V8 = point(V8x, V8y, V8z);


// update rho, phi and teta
calculateRhoTetaPhi(V1.x,V1.y,V1.z)


// update translation step because depends on dimension "a" that changed
setTranslationStep(V1)

}
// -----------------------------------------------------------

    function updateBoxVerticesForBScale(sign) {

//  update the 8 vertices by shifting them along (pt3, pt1) vector

Vx = (V1.x - V3.x)*sign*scaleStep;
Vy = (V1.y - V3.y)*sign*scaleStep;
Vz = (V1.z - V3.z)*sign*scaleStep;

// ********** V1,V3 **********
V1x = V1.x + Vx;
V1y = V1.y + Vy;
V1z = V1.z + Vz;

V1 = point(V1x, V1y, V1z);

V3x = V3.x - Vx;
V3y = V3.y - Vy;
V3z = V3.z - Vz;

V3 = point(V3x, V3y, V3z);

// ********** V5, V7 **********
V5x = V5.x + Vx;
V5y = V5.y + Vy;
V5z = V5.z + Vz;

V5 = point(V5x, V5y, V5z);

V7x = V7.x - Vx;
V7y = V7.y - Vy;
V7z = V7.z - Vz;

V7 = point(V7x, V7y, V7z);

// ********** V6,8 **********
V6x = V6.x + Vx;
V6y = V6.y + Vy;
V6z = V6.z + Vz;

V6 = point(V6x, V6y, V6z);

V8x = V8.x - Vx;
V8y = V8.y - Vy;
V8z = V8.z - Vz;

V8 = point(V8x, V8y, V8z);

// ********** V2,4 **********
V2x = V2.x + Vx;
V2y = V2.y + Vy;
V2z = V2.z + Vz;

V2 = point(V2x, V2y, V2z);

V4x = V4.x - Vx;
V4y = V4.y - Vy;
V4z = V4.z - Vz;

V4 = point(V4x, V4y, V4z);

// update rho, phi and teta
calculateRhoTetaPhi(V1.x,V1.y,V1.z)


// update translation step because depends on dimension "b" that changed
setTranslationStep(V1)

} // end of  updateBoxVerticesForBScale

// -----------------------------------------------------------

    function updateBoxVerticesForCScale(sign) {

//  update the 8 vertices by shifting them along (pt2, pt1) vector

Vx = (V1.x - V2.x)*sign*scaleStep;
Vy = (V1.y - V2.y)*sign*scaleStep;
Vz = (V1.z - V2.z)*sign*scaleStep;

// ********** V3,V4 **********
V3x = V3.x + Vx;
V3y = V3.y + Vy;
V3z = V3.z + Vz;

V3 = point(V3x, V3y, V3z);

V4x = V4.x - Vx;
V4y = V4.y - Vy;
V4z = V4.z - Vz;

V4 = point(V4x, V4y, V4z);

// ********** V1, V2 **********
V1x = V1.x + Vx;
V1y = V1.y + Vy;
V1z = V1.z + Vz;

V1 = point(V1x, V1y, V1z);

V2x = V2.x - Vx;
V2y = V2.y - Vy;
V2z = V2.z - Vz;

V2 = point(V2x, V2y, V2z);

// ********** V7,8 **********
V7x = V7.x + Vx;
V7y = V7.y + Vy;
V7z = V7.z + Vz;

V7 = point(V7x, V7y, V7z);

V8x = V8.x - Vx;
V8y = V8.y - Vy;
V8z = V8.z - Vz;

V8 = point(V8x, V8y, V8z);

// ********** V5,6 **********
V5x = V5.x + Vx;
V5y = V5.y + Vy;
V5z = V5.z + Vz;

V5 = point(V5x, V5y, V5z);

V6x = V6.x - Vx;
V6y = V6.y - Vy;
V6z = V6.z - Vz;

V6 = point(V6x, V6y, V6z);

// update rho, phi and teta
calculateRhoTetaPhi(V1.x,V1.y,V1.z)


// update translation step because translationStep depends on dimension "c" that changed
setTranslationStep(V1)

} // end of updateBoxVerticesForCScale



// -----------------------------------------------------------
    function updateBoxVerticesForPhi() {

//
// rotates the 8 box vertices
//suppose that PHI rotation angle is now  PHI (has already been updated)

// calculate the rotation axis from cross product between reference boxvector and z axis
vect1 = V1;
if     ( (teta>0.0)  && (teta <180.0)) {vect2 = {0  0  1};} \
elseif ( (teta>180.0) &&(teta <360.0)) {vect2 = {0  0  -1};};

u = cross(vect1,vect2);
normU = sqrt(u.x*u.x + u.y*u.y + u.z*u.z);
u = u/normU;

// given an angle and a rotation axis, 
//calculate the new vertices from quaternions  
rotat = -1.0*rotAngle;
quatROT = quaternion(u, rotat);

// apply this quaternion to the 8 vertices
// remark: since V1 and vect2 are vectors expressed in
//         the original frame (centered on 0,0,0), 
//         and that V1,V2,... V8, are vectors (not points)
//         then, the userBox is first rotated in the original frame
//         next, in drawbox function, V1,V2, ... V8 will serve to define
//         pt1,pt2,...pt8 to define the rotated user box in the userBox frame.
//         (accounting for possible translations)
//         in other words: since bbic is not orginally at 0,0,0
//         the resulting rotation vector "u" will be finally translated at the userBoxcenter
//         (in drawPhiAxis function)
V1 = quatROT%V1;
V2 = quatROT%V2;
V3 = quatROT%V3;
V4 = quatROT%V4;
V5 = quatROT%V5;
V6 = quatROT%V6;
V7 = quatROT%V7;
V8 = quatROT%V8;

// update rho, phi and teta
calculateRhoTetaPhi(V1.x,V1.y,V1.z)


} // end of updateBoxVerticesForPhi

// -----------------------------------------------------------
    function updateBoxVerticesForTeta() {


//suppose that TETA rotation angle is now  TETA (has already been updated)

// the TETA case is easier than the PHI one, because
// the axis is known: it is the z axis
u = {0 0 1} 

// given an angle and a rotation axis, 
//calculate the new vertices from quaternions    

quatROT = quaternion(u, rotAngle);
// apply this quaternion to the 8 vertices
V1 = quatROT%V1;
V2 = quatROT%V2;
V3 = quatROT%V3;
V4 = quatROT%V4;
V5 = quatROT%V5;
V6 = quatROT%V6;
V7 = quatROT%V7;
V8 = quatROT%V8;

// update rho, phi and teta
calculateRhoTetaPhi(V1.x,V1.y,V1.z)

} // end of updateBoxVerticesForTeta