anchor.diff:  !!$  q_dot(7:9)   = (/ 0.0_dp, 0.0_dp, 0.0_dp /)
   !!$  q_dot(10:12) = (/ 0.0_dp, 0.0_dp, 0.0_dp /)
   
  -!  CALL init_random_seed()         ! see example of RANDOM_SEED
  +  !  CALL init_random_seed()         ! see example of RANDOM_SEED
     CALL RANDOM_NUMBER(q)
     CALL RANDOM_NUMBER(q_dot)  
     q_dot = q_dot*q_dot+1.0_dp
   
     ! define the attributes of the body
     mass      = 2.0_dp
  -  ! point of interest on the body  is the point where the body axis is located
  +
  +  ! point of interest on the body
     CALL RANDOM_NUMBER(re)
     !  re        = (/ 0.0_dp, 0.0_dp, 0.0_dp /)   
   
     !-------------------------------------------------------------------!
     ! (2) create a pendulum body using the state and attr above
     !-------------------------------------------------------------------!
  -  aa = 1.0_dp
  -
  -  call disp(" >> Creating a body...")
  +
  +  !  call disp(" >> Creating a body...")
   
     body1 = create_body(mass, vector(re), q, q_dot)
   
  -  call disp(" >> Body created successfully...")
  -    
  -  ! residual assembly
  -  call disp(" >> Assembling the residuals...")
  -  res  = get_residual(body1)
  -  call disp(" >> Residual assembly complete...")
  -  call disp("   R   =   ", res)
  -
  -  ! jacobian assembly
  -  call disp(" >> Assembling the Jacobian...")
  -  jac  = get_jacobian(body1)
  -  call disp(" >> Jacobian assembly complete...")
  -  call disp("   JAC =   ",jac)
  -
  -
  -  ! implement FD to check if the jacobian is the same
  -  qtmp = q
  -  call set_state(q, q_dot, body1)
  -  dh = 1.0d-8
  -
  -  res = get_residual(body1)
  -     
  -  var: do  k = 1 , TOT_NDOF
  -
  -     ! perturb k-th variable
  -     qtmp( k ) = qtmp (k)  + dh
  -     ! updating the body state
  -     call set_state(qtmp, q_dot, body1)
  -     ! find the solution
  -     res2 = get_residual(body1)
  -     !    call disp("   R   =   ", res)
  -     
  -     dq =  (res2-res)/dh
  -     
  -     jac2 (:, k) = dq
  -     
  -     qtmp = q
  -
  -  end do var
  -  call disp("   dJAC =   ",jac2)
  - 
  -  call disp("   Error =   ",jac2-Jac)
  -
  -
  -
  -  stop"stopped"
  +  !call disp(" >> Body created successfully...")
   
     !-------------------------------------------------------------------!
     ! Settings for time-integration
  
positive.code:      newton: do k = 1, MAX_NEWTON_ITER
  
          ! (3) Assemble the residual and jacobian using the body
          
          ! call this method to update the state for every other iter
          if (k .gt. 1) then
             call set_state(q, q_dot, body1)
          end if
  
          ! residual assembly
          call disp(" >> Assembling the residuals...")
          res  = get_residual(body1)
          call disp(" >> Residual assembly complete...")
          call disp("   R   =   ", res)
  
          ! jacobian assembly
          call disp(" >> Assembling the Jacobian...")
          jac  = get_jacobian(body1)
          call disp(" >> Jacobian assembly complete...")
          call disp("   JAC =   ",jac)
  
          ! ******************************************************
          ! (3) Solve the linear system and compute update delta_q
          ! ******************************************************
  
          call disp(" >> Calling the linear solver...") 
  
          !dq = linear_solve(jac,-res,'GMRES')
  
          call disp(" >> Update to residual is computed...") 
  
          l2_norm   = norm2(dq); linf_norm =  maxval(abs(dq));
  
  
negative.code:   call disp("                 Rigid body dynamics                    ")
    call disp("========================================================")
  
    ! init routine that has sanity check on the input settings
    ! has nothing as of now, we can use to load input files, python etc
    call init() 
  
    !-------------------------------------------------------------------!
    ! (1) set the initial states and attributes of the body
    !-------------------------------------------------------------------!
  
    call disp(" >> Setting up the test problem...")
  
  !!$  ! define the initial states
  !!$  q(1:3)    = (/ 1.0_dp, 2.0_dp, 3.0_dp /)
  !!$  q(4:6)    = (/ deg2rad(20.0_dp), deg2rad(0.0_dp), deg2rad(20.0_dp) /)
  