anchor.diff:     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:
  -
  -!      xd%DummyDiscState = 0.0
  -
  +      
  +! 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 )
  
positive.code: 
  END SUBROUTINE BD_InitshpDerJaco
  
  
  END MODULE BeamDyn
  
negative.code:       CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )   
     CALL BD_CopyInput(u, u_tmp, MESH_NEWCOPY, ErrStat2, ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        
     ! Actuator
     IF( p%torq ) THEN
         CALL PitchActuator_SetBC(p, u_tmp, xd)
     ENDIF
     ! END Actuator
     
     CALL BD_CopyInput(u_tmp, u_tmp2, MESH_NEWCOPY, ErrStat2, ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        
     CALL BD_CrvExtractCrv(p%GlbRot,temp_glb,ErrStat2,ErrMsg2)
        CALL SetErrStat( ErrStat2, ErrMsg2, ErrStat, ErrMsg, RoutineName )
        if (ErrStat >= AbortErrLev) then
  