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(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) = TmpHPlusC
  -
  +                     Term2ArrayC(J,ThisDim) = TmpHPlusC
   
                     ENDIF    ! Check on the limits
   
                  ENDDO       ! Second term calculation -- frequency step on the sum frequency
  -
  -
  -
  -                  ! Now we apply the FFT to the result of the sum.
  -               CALL ApplyFFT_cx(  Term1Array(:),  Term1ArrayC(:), 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(:), 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)
  -
               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
  
positive.code: 
  
           ! Setup the arrays holding the SumQTF terms, both the complex frequency domain and real time domain pieces
        ALLOCATE( Term1ArrayC( 0:InitInp%NStepWave2), STAT=ErrStatTmp )
        IF (ErrStatTmp /= 0) CALL SetErrStat(ErrID_Fatal,' Cannot allocate array for the first term of one load component of the full sum '// &
                                               'QTF 2nd order force in the frequency domain.',ErrStat, ErrMsg, RoutineName)
        ALLOCATE( Term2ArrayC( 0:InitInp%NStepWave2), STAT=ErrStatTmp )
        IF (ErrStatTmp /= 0) CALL SetErrStat(ErrID_Fatal,' Cannot allocate array for the second term of one load component of the full sum '// &
                                               'QTF 2nd order force in the frequency domain.',ErrStat, ErrMsg, RoutineName)
        ALLOCATE( Term1Array( 0:InitInp%NStepWave), STAT=ErrStatTmp )
  
negative.code:                   ! Set an initial search index for the 4D array interpolation
                 LastIndex4 = (/0,0,0,0/)
  
                    ! 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
  