anchor.diff:                    ! 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
  
positive.code:                         RETURN
                       ENDIF
  
                          ! Set the value of the first term in the frequency domain
                       Term1ArrayC(2*J) = aWaveElevC1 * aWaveElevC1 * QTF_Value
  
  
                    ENDIF    ! Check on the limits
                 ENDDO       ! First term calculation
  
  
  
                 !---------------------------------------------------------------------------------
  
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)) )
                             !!                  - 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(K)*D2R)   + InitInp%PtfmRefyt(1)*sin(InitInp%WaveDirArr(K)*D2R)   )
                             TmpReal2    = WaveNmbr2 * ( InitInp%PtfmRefxt(1)*cos(InitInp%WaveDirArr(J-K)*D2R) + InitInp%PtfmRefyt(1)*sin(InitInp%WaveDirArr(J-K)*D2R) )
  
                             ! Set the phase shift for the set of sum frequencies
                             PhaseShiftXY = CMPLX( cos(TmpReal1 + TmpReal2), -sin(TmpReal1 + TmpReal2) )
  
                             QTF_Value = QTF_Value*PhaseShiftXY
  
                          endif ! Phaseshift for NBodyMod==2
  
                             ! Set the value of the first term in the frequency domain.
                          TmpHPlusC = TmpHPlusC + aWaveElevC1 * aWaveElevC2 * QTF_Value
  
                       ENDDO
  
                          ! Save the value from the summation.
                       Term2ArrayC(J,ThisDim) = TmpHPlusC
  
                    ENDIF    ! Check on the limits
  
                 ENDDO       ! Second term calculation -- frequency step on the sum frequency
              ENDIF    ! Load component to calculate
           ENDDO ! ThisDim -- current dimension
  
  
           !----------------------------------------------------
           ! Rotate back to global frame
           !----------------------------------------------------
  
              ! Set rotation
              ! NOTE: RotateZMatrixT is the rotation from local to global.
           RotateZMatrixT(:,1) = (/  cos(InitInp%PtfmRefztRot(IBody)), -sin(InitInp%PtfmRefztRot(IBody)) /)
           RotateZMatrixT(:,2) = (/  sin(InitInp%PtfmRefztRot(IBody)),  cos(InitInp%PtfmRefztRot(IBody)) /)
  
              ! Loop through all the frequencies
           DO J=1,InitInp%NStepWave2
  
                 ! Apply the rotation to get back to global frame -- term 1
              Term1ArrayC(J,1:2) = MATMUL(RotateZMatrixT, Term1ArrayC(J,1:2))
              Term1ArrayC(J,4:5) = MATMUL(RotateZMatrixT, Term1ArrayC(J,4:5))
  
                 ! Apply the rotation to get back to global frame -- term 2
              Term2ArrayC(J,1:2) = MATMUL(RotateZMatrixT, Term2ArrayC(J,1:2))
              Term2ArrayC(J,4:5) = MATMUL(RotateZMatrixT, Term2ArrayC(J,4:5))
  
           ENDDO    ! J=1,InitInp%NStepWave2
  
  
  
           !----------------------------------------------------
           ! Apply the FFT to get time domain results
           !----------------------------------------------------
  
           DO ThisDim=1,6
              Idx = (IBody-1)*6+ThisDim
  
                 ! Now we apply the FFT to the result of the sum.
              CALL ApplyFFT_cx(  Term1Array(:),  Term1ArrayC(:,ThisDim), FFT_Data, ErrStatTmp )
              CALL SetErrStat(ErrStatTmp,'Error occured while applying the FFT to the first term of the Sum QTF.', &
                             ErrStat,ErrMsg,RoutineName)
              IF ( ErrStat >= AbortErrLev ) THEN
                 IF (ALLOCATED(TmpData4D))        DEALLOCATE(TmpData4D,STAT=ErrStatTmp)
                 RETURN
              END IF
  
                 ! Now we apply the FFT to the result of the sum.
              CALL ApplyFFT_cx(  Term2Array(:), Term2ArrayC(:,ThisDim), FFT_Data, ErrStatTmp )
              CALL SetErrStat(ErrStatTmp,'Error occured while applying the FFT to the second term of the Sum QTF.', &
                             ErrStat,ErrMsg,RoutineName)
              IF ( ErrStat >= AbortErrLev ) THEN
                 IF (ALLOCATED(TmpData4D))        DEALLOCATE(TmpData4D,STAT=ErrStatTmp)
                 RETURN
              ENDIF
  
                 ! Now we add the two terms together.  The 0.5 multiplier on is because the double sided FFT was used.
              DO J=0,InitInp%NStepWave-1  !bjj: Term1Array and Term2Array don't set the last element, so we can get over-flow errors here. SumQTFForce(InitInp%NStepWave,Idx) gets overwritten later, so Idx'm setting the array bounds to be -1.
                 SumQTFForce(J,Idx) = 0.5_SiKi*(REAL(Term1Array(J) + 2*Term2Array(J), SiKi))
              ENDDO
  
                 ! Copy the last first term to the first so that it is cyclic
              SumQTFForce(InitInp%NStepWave,Idx) = SumQTFForce(0,Idx)
  
  
           ENDDO ! ThisDim -- current dimension
        ENDDO    ! IBody -- current WAMIT body
  