anchor.diff:                             RETURN
                           ENDIF
   
  +                        !--------------------------
  +                        ! Phase shift due to offset
  +                        !--------------------------
  +                        if (p%NBodyMod == 2) then
  +                           !> The phase shift due to an (x,y) offset for second order difference frequencies is of the form
  +                           !! \f$  exp[-\imath ( k(\omega_1) ( X cos(\Beta(w_1)) + Y sin(\beta(w_1)) )
  +                           !!                  - k(\omega_2) ( X cos(\Beta(w_2)) + Y sin(\beta(w_2)) ) ) ]\f$
  +                           !  NOTE: the phase shift applies to the aWaveElevC of the incoming wave.  Including it here instead
  +                           !        of above is mathematically equivalent, but only because each frequency has only one wave
  +                           !        direction associated with it through the equal energy approach used in multidirectional waves.
  +
  +                           WaveNmbr1   = WaveNumber ( REAL(Omega1,SiKi), InitInp%Gravity, InitInp%WtrDpth )    ! SiKi returned
  +                           WaveNmbr2   = WaveNumber ( REAL(Omega2,SiKi), InitInp%Gravity, InitInp%WtrDpth )    ! SiKi returned
  +                           TmpReal1    = WaveNmbr1 * ( InitInp%PtfmRefxt(1)*cos(InitInp%WaveDirArr(J+K)*D2R) + InitInp%PtfmRefyt(1)*sin(InitInp%WaveDirArr(J+K)*D2R) )
  +                           TmpReal2    = WaveNmbr2 * ( InitInp%PtfmRefxt(1)*cos(InitInp%WaveDirArr(K)*D2R)   + InitInp%PtfmRefyt(1)*sin(InitInp%WaveDirArr(K)*D2R)   )
  +
  +                           ! Set the phase shift for the set of difference frequencies
  +                           PhaseShiftXY = CMPLX( cos(TmpReal1 - TmpReal2), -sin(TmpReal1 - TmpReal2) )
  +
  +                           ! For similicity, apply to the QTF_Value (mathematically equivalent to applying to the wave elevations)
  +                           QTF_Value = QTF_Value*PhaseShiftXY
  +
  +                        endif ! Phaseshift for NBodyMod==2
  +
                              ! Calculate this value and add it to what we have so far.
                           TmpHMinusC = TmpHMinusC + aWaveElevC1 * CONJG(aWaveElevC2) * QTF_Value
   
  
positive.code: 
           ! Now loop through all the dimensions and perform the calculation
        DO IBody=1,p%NBody
           DO ThisDim=1,6
              Idx = (IBody-1)*6+ThisDim
  
  
negative.code:       REAL(ReKi)                                         :: Omega1               !< First  wave frequency
        REAL(ReKi)                                         :: Omega2               !< Second wave frequency
        REAL(SiKi),    ALLOCATABLE                         :: MnDriftForce(:)      !< Mean drift force (first term).  MnDrift_InitCalc routine will return this.
        REAL(SiKi)                                         :: RotateZdegOffset     !< Offset to wave heading (NBodyMod==2 only)
        REAL(SiKi)                                         :: RotateZMatrixT(2,2)  !< The transpose of rotation in matrix form for rotation about z (from global to local)
        COMPLEX(SiKi)                                      :: PhaseShiftXY         !< The phase shift offset to apply to the body
        REAL(SiKi)                                         :: WaveNmbr1            !< Wavenumber for this frequency
        REAL(SiKi)                                         :: WaveNmbr2            !< Wavenumber for this frequency
  
           ! Interpolation routine indices and value to search for, and smaller array to pass
        INTEGER(IntKi)                                     :: LastIndex4(4)        !< Last used index for searching in the interpolation algorithms.  First  wave freq
  