anchor.diff:     TYPE(BD_ContinuousStateType)                       :: x_tmp      ! Holds temporary modification to x
      TYPE(BD_OtherStateType     )                       :: OS_tmp     ! Holds temporary modification to x
      TYPE(BD_InputType)                                 :: u_interp   ! interpolated value of inputs
  -   TYPE(BD_InputType)                                 :: u_tmp2
      INTEGER(IntKi)                                     :: ErrStat2   ! Temporary Error status
      CHARACTER(ErrMsgLen)                               :: ErrMsg2    ! Temporary Error message
      CHARACTER(*), PARAMETER                            :: RoutineName = 'BD_GA2'
  
positive.code: 
  END SUBROUTINE BD_InitshpDerJaco
  
  
  END MODULE BeamDyn
  
negative.code:    INTEGER(IntKi),                    INTENT(  OUT)  :: ErrStat     ! Error status of the operation
     CHARACTER(*),                      INTENT(  OUT)  :: ErrMsg      ! Error message if ErrStat /= ErrID_None
  
     ! local variables
     REAL(ReKi)                                        :: temp_R(3,3) 
     REAL(ReKi)                                        :: I_minus_hA_inv(2,2) 
     REAL(ReKi)                                        :: Hub_theta_Root(3) 
     REAL(ReKi)                                        :: u_theta_pitch 
     REAL(ReKi)                                        :: denom 
     
        ! Initialize ErrStat
  
        ErrStat = ErrID_None
        ErrMsg  = ""
  
        ! Update discrete states here:
        
  ! Actuator
     IF( p%torq ) THEN
         temp_R = MATMUL(u%RootMotion%Orientation(:,:,1),TRANSPOSE(u%HubMotion%Orientation(:,:,1)))
         Hub_theta_Root = EulerExtract(temp_R)
         u_theta_pitch = -Hub_theta_Root(3)
                
  !bjj: this matrix should be stored as a parameter and calculated only at initialization:>>>>>>
         I_minus_hA_inv(1,1) =  p%pitchJ + p%pitchC*p%dt
         I_minus_hA_inv(2,1) = -p%pitchK * p%dt
         I_minus_hA_inv(1,2) =  p%pitchJ * p%dt
         I_minus_hA_inv(2,2) =  p%pitchJ
         denom               =  p%pitchJ + p%pitchC*p%dt + p%pitchK*p%dt**2
         I_minus_hA_inv      = I_minus_hA_inv / denom
  !<<<<<<<<<<<<<<
         
         xd%thetaP  = I_minus_hA_inv(1,1)*xd%thetaP + I_minus_hA_inv(1,2)*xd%thetaPD + &
                                                      I_minus_hA_inv(1,2)*(p%pitchK*p%dt/p%pitchJ)*(-Hub_theta_Root(3))
         xd%thetaPD = I_minus_hA_inv(2,1)*xd%thetaP + I_minus_hA_inv(2,2)*xd%thetaPD + &
                                                      I_minus_hA_inv(2,2)*(p%pitchK*p%dt/p%pitchJ)*(-Hub_theta_Root(3))
         
  
     ENDIF
  ! END Actuator      
        
  END SUBROUTINE BD_UpdateDiscState
  !-----------------------------------------------------------------------------------------------------------------------------------
  SUBROUTINE BD_CalcConstrStateResidual( t, u, p, x, xd, z, OtherState, Z_residual, ErrStat, ErrMsg )
  