anchor.diff:  
     ! output filename
     character                          :: filename   
  +
  +  ! logical for success of newton iterations
  +  logical                            :: newton_success = .false.
   
     call disp("========================================================")
     call disp("'                 Rigid body dynamics                  '")
  
positive.code:   call disp(" >> Setting initial state of the body...")
  
    call random_seed(); call random_number(q); call random_number(q_dot);
    q_dot = q_dot
  
    ! mass of the body
    mass = 2.0_dp
    
    ! used to calculate the inertial properties J and C
    re   = (/ 0.1_dp, 0.2_dp, 0.3_dp /) 
  
  
negative.code: 
    ! import the necessary modules
    use global_constants, only: sp, dp, NUM_SPAT_DIM, TOT_NDOF
    use global_variables, only: ABS_TOL, REL_TOL, start_time, dT, &
         & end_time, aa, bb, MAX_NEWTON_ITER, fcnt, init, time
    use utils, only:deg2rad
    use types, only: body, vector
    use solver_utils, only: get_updated_q, get_updated_q_dot,&
         & get_updated_q_double_dot, get_approx_q, get_approx_q_dot
    use bodies, only: create_body, set_state
    use residual, only: get_residual
    use jacobian, only: get_jacobian
    use finite_diff, only: finite_difference2
    use linear_system, only: direct_solve
    use bodies, only: print_body, create_body
    use dispmodule, only: disp
    use filehandler
  
    implicit none
  