/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
/* [ Created with wxMaxima version 14.12.1 ] */

/* [wxMaxima: input   start ] */
reset()$
kill(all)$
load(vect) $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Calculation of the torque applied by the TM
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
/* Young modulus of piano wire */
E_piano : 19e10 $

/* Acceleration due to gravity in m/s^2 */
g : 9.81 $

/* Diameter of the piano wire in m */
diameter_piano_TM : 200e-6 $

/* Mass of TM in kg */
mass_TM : 10.7 $

/* Just an unitary vector point downwards */
j : [ 0 , -1 , 0 ] $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
/* 2nd moment of area */
moment2ndarea_piano : float(%pi) * diameter_piano_TM ^ 4 / 64 $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
/* Position of bending poiint from the clamping point */
y_bend : sqrt( E_piano * moment2ndarea_piano / ( g * mass_TM /4 ) ) $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
/* In such a case the TM moves sideways by the following amount */
/* theta : float( 2 * %pi/180 ) ; */
x_TM : 2 * y_bend * sin( theta );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
/* Nominal position of the T center of mass with */
/* respect to the bending point of the rod holdeing the IM */
Rcm_TM : [ x_TM ,  -587e-3 , 0 ] + [ 0 , -1e-3 , 0 ] $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
( Rcm_TM ~ j ) * mass_TM * g  $
Torque_TM : express( % ) $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Torque_TM_Z : Torque_TM[3]  ;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Torque applied by the RM.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
/* Young modulus of tunsgten wire */
E_tungsten : 411e9 $

/* Diameter of the tungsten wire in m */
diameter_tungsten_RM : 600e-6 $

/* Mass of RM in kg */
mass_RM :  10.7 + 0.745  $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
/* 2nd moment of area */
moment2ndarea_RM : float(%pi) * diameter_tungsten_RM ^ 4 / 64 $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
/* Position of bending point from the clamping point */
y_bend_tungsten : sqrt( E_tungsten * moment2ndarea_RM / ( g * mass_RM /4 ) );
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
/* In such a case the RM moves sideways by the following amount */
/* theta : float( 2 * %pi/180 ) ; */
x_RM : 2 * y_bend_tungsten * sin( theta ) $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
/* Nominal position of the the center of mass with */
/* respect to the bending point of the rod holding the IM */
Rcm_RM : [ x_RM ,  -587e-3 , 0 ] + [ 0 , -1e-3 , 0 ] $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
( Rcm_RM ~ j ) * mass_RM * g  $
Torque_RM : express( % ) $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
Torque_RM_Z : Torque_RM[3]  ;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Calculation of the torque applied by the IM
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
/**/
Mass_IM_total : 26.464 $
/**/
Mass_X : 0.879 $
/**/
Mass_Z : 0.891 $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
/* Position of the center of mass of the payload WITHOUT the */
/* movable bodies with respect to the */
/* suspension point above the IM */
Rcm_IM_reduced:  ( Mass_IM_total * [ 0 , -1e-3 , 0 ] -
    Mass_X * ( [ 0 , 29e-3 , 0 ] + [ 0 , -1e-3 , 0 ] ) -
    Mass_Z * ( [ 0 , -36e-3 , 0 ] + [ 0 , -1e-3 , 0 ] ) )
    / Mass_IM_reduced $

/**/
Rcm_X : [ x , 29e-3 , 0 ] + [ 0 , -1e-3 , 0 ] $

/**/
Rcm_Z :  [ 0 , -36e-3 , 0e-3 ] + [ 0 , -1e-3 , 0 ]  $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
/* The center of mass is calculated */
/* It should be [ 0 -1e-3 0 ] when x and z above are zero */
Rcm : ( Mass_IM_reduced * Rcm_IM_reduced +
        Mass_X * Rcm_X + Mass_Z * Rcm_Z ) / Mass_IM_total $
norm_Rcm : sqrt( Rcm. Rcm ) $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
alpha : acos( Rcm . j / norm_Rcm ) $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
/* Torque applied by the IM */
Torque_IM_Z: norm_Rcm * Mass_IM_total * g * sin( alpha - theta ) ;
/* [wxMaxima: input   end   ] */

/* [wxMaxima: comment start ]
Total torque and final value of theta.
   [wxMaxima: comment end   ] */

/* [wxMaxima: input   start ] */
/* tilt : find_root( Torque_TM_Z + Torque_IM_Z + Torque_RM_Z = 0 ,
                            theta , 0 , %pi/4 ) ;  */

/* alpha - tilt must be greater then zero */
alpha - tilt $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
tilt * 180 / float( %pi ) $
/* [wxMaxima: input   end   ] */

/* [wxMaxima: input   start ] */
plot2d( find_root( Torque_TM_Z + Torque_IM_Z + Torque_RM_Z = 0 ,
  theta , 0 , 5* %pi/180 )  ,[x,0, 10e-3 ],
 [xlabel,"Displacement of mass (m)"], [ylabel,"Roll (rad)"]);
/* [wxMaxima: input   end   ] */

/* Maxima can't load/batch files which end with a comment! */
"Created with wxMaxima"$