anchor.diff:       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...")
  +        !call disp(" >> Assembling the residuals...")
           res  = get_residual(body1)
  -        call disp(" >> Residual assembly complete...")
  -        call disp("   R   =   ", res)
  +        !call disp(" >> Residual assembly complete...")
  +        !call disp("   R   =   ", res)
   
           ! jacobian assembly
  -        call disp(" >> Assembling the Jacobian...")
  +        !call disp(" >> Assembling the Jacobian...")
           jac  = get_jacobian(body1)
  -        call disp(" >> Jacobian assembly complete...")
  -        call disp("   JAC =   ",jac)
  +        !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...") 
  +        !call disp(" >> Calling the linear solver...") 
   
           !dq = linear_solve(jac,-res,'GMRES')
   
  -        call disp(" >> Update to residual is computed...") 
  +        !call disp(" >> Update to residual is computed...") 
   
           l2_norm   = norm2(dq); linf_norm =  maxval(abs(dq));
   
  
positive.code:   call disp("                 Rigid body dynamics                    ")
    call disp("========================================================")
  
    ! init method 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) /)
  
negative.code: 
  end subroutine residual
  
  !!$
  !!$
  !!$
  !!$ 
  !!$  ! 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"
  