anchor.diff:        ! ..................
         ! Translational velocity:
         ! ..................
  -      mm = maxval( abs(Mesh1_O_op%TranslationVel ) )
         call meshcopy( Mesh1_O_op, Mesh1_O, MESH_UPDATECOPY, ErrStat, ErrMsg) 
         do n1=1,15
            
            ! perturb x^S:
  -         call getRotationPerturb(LinVec_1)
  -         call getRandomVector(LinVec_1_a)
  -         call getRandomVector(LinVec_1_b)
  -         
  -         LinVec_2_a = matmul( map_mod1_mod2%dm%tv, LinVec_1 ) + matmul( map_mod1_mod2%dm%mi, LinVec_1_a ) &
  -                    + matmul( map_mod1_mod2%dm%fx_p, LinVec_1_b )! approximate delta theta^D         
  -                      
  +         
  +        !call getRandomVector(LinVec_2_b) ! u^D  ! for this, we will actually perturb u^S (LinVec_1_a) and theta^S (LinVec_1)
  +         call getRotationPerturb(LinVec_1) ! theta^S
  +         call getRandomVector(LinVec_1_a)  ! u^S
  +         call getRandomVector(LinVec_1_b)  ! v^S 
  +         call getRandomVector(LinVec_1_c)  ! omega^S 
  +         
  +         
  +            
  +         ! get the displaced position         
            do i=1,Mesh1_O%nnodes
               j = (i-1)*3 + 1
               
                  ! Mesh1_O%Orientation = theta^S|_op + delta theta^S
  -            !call SmllRotTrans( 'orientation', LinVec_1(j  ) &
  -            !                                , LinVec_1(j+1) &
  -            !                                , LinVec_1(j+2) & 
  -            !                                , Orientation, ErrStat=ErrStat, ErrMsg=ErrMsg)            
               call getLinearOrient(LinVec_1(j:j+2), Orientation)
               Mesh1_O%Orientation(:,:,i) = matmul( Mesh1_O_op%Orientation(:,:,i), Orientation )
  -
  +                        
  +               ! Mesh1_O%TranslationDisp = u^S|_op + delta u^S
  +            Mesh1_O%TranslationDisp(:,i) = Mesh1_O_op%TranslationDisp(:,i) + LinVec_1_a(j:j+2)
                           
                  ! Mesh1_O%TranslationVel = v^S|_op + delta v^S
  -            Mesh1_O%TranslationVel(:,i) = Mesh1_O_op%TranslationVel(:,i) + LinVec_1_a(j:j+2)
  -            
  +            Mesh1_O%TranslationVel(:,i) = Mesh1_O_op%TranslationVel(:,i) + LinVec_1_b(j:j+2)
  +                   
                  ! Mesh1_O%RotationVel = omega^S|_op + delta omega^S
  -            Mesh1_O%RotationVel(:,i) = Mesh1_O_op%RotationVel(:,i) + LinVec_1_b(j:j+2)
  +            Mesh1_O%RotationVel(:,i) = Mesh1_O_op%RotationVel(:,i) + LinVec_1_c(j:j+2)
               
               !IF (ErrStat /= ErrID_None) CALL WrScr("*******"//TRIM(ErrMsg))
  -                        
  -         end do
  -         
  -
  +            
  +         end do
  +         
  +         
            !M( x^S|_op + delta x^S )
  -         call TransferMotionData()
  +         call TransferMotionData() ! this sets Mesh2_I, the displaced source mesh
  +                           
   
            ! delta x^D = M( x^S|_op + delta x^S ) - x^D|_op
            do i=1,Mesh2_I%nnodes
               j = (i-1)*3 + 1
               LinVec_2(j:j+2) = Mesh2_I%TranslationVel(:,i) - Mesh2_I_op%TranslationVel(:,i)          
            end do
  -                  
  -         e=TwoNorm( LinVec_2 - LinVec_2_a )               
  -         !if (.not. equalrealNos(e,0.0_ReKi)) &
  -            write(*,ErrFmt) 'Translational Velocity: ', n1, e, e/max(mmin, TwoNorm(LinVec_2)), maxval(abs(LinVec_2)), &
  -               mm, maxval(abs(LinVec_1)), maxval(abs(LinVec_1_a)), maxval(abs(LinVec_1_b)) 
  -                              
  +         
  +         ! delta u^D
  +         do i=1,Mesh2_I%nnodes
  +            j = (i-1)*3 + 1
  +            LinVec_2_b(j:j+2) = Mesh2_I%TranslationDisp(:,i) - Mesh2_I_op%TranslationDisp(:,i)          
  +         end do
  +         
  +         LinVec_2_a = matmul( map_mod1_mod2%dm%tv_uD, LinVec_2_b ) + matmul( map_mod1_mod2%dm%tv_uS, LinVec_1_a ) &
  +                    + matmul( map_mod1_mod2%dm%mi, LinVec_1_b ) + matmul( map_mod1_mod2%dm%fx_p, LinVec_1_c )! approximate delta theta^D         
  +                               
  +         
  +         call WrErrorLine('Translational Velocity', LinVec_2, LinVec_2_a, Mesh1_O_op%TranslationVel, LinVec_1, LinVec_1_a, LinVec_1_b)
  +            
         end do      
         write(*,*) ! blank line
         
         ! ..................
         ! Rotational acceleration:
         ! ..................
  -      mm = maxval( abs(Mesh1_O_op%RotationAcc ) )
               
         call meshcopy( Mesh1_O_op, Mesh1_O, MESH_UPDATECOPY, ErrStat, ErrMsg) 
         call meshcopy( Mesh2_O_op, Mesh2_O, MESH_UPDATECOPY, ErrStat, ErrMsg) 
  
positive.code:             j = (i-1)*3 + 1
              LinVec_2(j:j+2) = Mesh2_I%RotationAcc(:,i) - Mesh2_I_op%RotationAcc(:,i)          
           end do
                    
           e=TwoNorm( LinVec_2 - LinVec_2_a )               
           !if (.not. equalrealNos(e,0.0_ReKi)) &
              write(*,ErrFmt) 'Rotational Acceleration: ', n1, e, e/max(mmin, TwoNorm(LinVec_2)), maxval(abs(LinVec_2)), &
                mm , maxval(abs(LinVec_1)) 
        end do            
        write(*,*) ! blank line
        
        ! ..................
        ! Translational acceleration:
        ! ..................
        mm = maxval( abs(Mesh1_O_op%TranslationAcc ) )
  
        call meshcopy( Mesh1_O_op, Mesh1_O, MESH_UPDATECOPY, ErrStat, ErrMsg) 
        call meshcopy( Mesh2_O_op, Mesh2_O, MESH_UPDATECOPY, ErrStat, ErrMsg) 
  
negative.code:       
        if (TestNumber == 13) then
  
           !write( 21, '(A)' )
           DO i=1,mesh2_o%nnodes
              n1 = mesh1_o%ElemTable(ELEMENT_LINE2)%Elements(Map_Mod1_Mod2%MapMotions(i)%OtherMesh_Element)%ElemNodes(1)
              n2 = mesh1_o%ElemTable(ELEMENT_LINE2)%Elements(Map_Mod1_Mod2%MapMotions(i)%OtherMesh_Element)%ElemNodes(2)
        
              WRITE (21, '(6(F15.5))' ) atan2( mesh2_i%refOrientation(1,2,i ) , mesh2_i%refOrientation(1,1,i ) ),&
                                        atan2( mesh2_i%Orientation(   1,2,i ) , mesh2_i%Orientation(   1,1,i ) ),&
                                        atan2( mesh1_o%refOrientation(1,2,n1) , mesh1_o%refOrientation(1,1,n1) ),&
                                        atan2( mesh1_o%Orientation(   1,2,n1) , mesh1_o%Orientation(   1,1,n1) ),&      
  