anchor.diff:  
                     ! Set an initial search index for the 4D array interpolation
                  LastIndex4 = (/0,0,0,0/)
  -
  -                  ! Initialize the temporary arrays for each term to zero.
  -               Term2ArrayC = CMPLX(0.0_SiKi,0.0_SiKi,SiKi)
  -
   
                     ! Check the limits for the high frequency cutoff. If WvHiCOffS is less than the
                     ! maximum frequency possible with the value of WaveDT (omega_max = pi/WaveDT = NStepWave2*WaveDOmega),
  
positive.code:                   ! Set an initial search index for the 4D array interpolation
                 LastIndex4 = (/0,0,0,0/)
  
                    ! Initialize the array to zero
                 Term1ArrayC = CMPLX(0.0_SiKi,0.0_SiKi,SiKi)
  
  
                    ! The limits look a little funny.  But remember we are placing the value in the 2*J location,
                    ! so we cannot overun the end of the array, and the highest frequency must be zero.  The
                    ! floor function is just in case (NStepWave2 - 1) is an odd number
  
negative.code:                         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)) )
                          !!                  1 k(\omega_2) ( X cos(\Beta(w_2)) + Y sin(\beta(w_2)) ) ) ]\f$.
                          !! For the first term, \f$ \omega_1 = \omega_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
                          TmpReal1    = WaveNmbr1 * ( InitInp%PtfmRefxt(1)*cos(InitInp%WaveDirArr(J)*D2R) + InitInp%PtfmRefyt(1)*sin(InitInp%WaveDirArr(J)*D2R) )
  
                          ! Set the phase shift for the set of sum frequencies
                          PhaseShiftXY = CMPLX( cos(TmpReal1 + TmpReal1), -sin(TmpReal1 + TmpReal1) )
  
                          ! 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
  
                          ! Set the value of the first term in the frequency domain
                       Term1ArrayC(2*J,ThisDim) = aWaveElevC1 * aWaveElevC1 * QTF_Value
  
  
                    ENDIF    ! Check on the limits
                 ENDDO       ! First term calculation
  
  
                 !---------------------------------------------------------------------------------
  