anchor.diff:  
                           ! Copy this value difference frequency information over to the array we will take the IFFT.  Divide
                           ! by two for the single sided FFT given in the documentation.
  -                     TmpComplexArr(J) = TmpHMinusC / 2.0_SiKi
  +                     TmpComplexArr(J,ThisDim) = TmpHMinusC / 2.0_SiKi
   
                     ELSE     ! outside the frequency range, so
   
  -                     TmpComplexArr(J) = CMPLX(0.0_SiKi,0.0_SiKi,SiKi)
  +                     TmpComplexArr(J,ThisDim) = CMPLX(0.0_SiKi,0.0_SiKi,SiKi)
   
                     ENDIF    ! frequency check
   
   
                  ENDDO
  -
  -                  ! Now we apply the FFT to the result of the sum
  -               CALL ApplyFFT_cx(  TmpDiffQTFForce(:),  TmpComplexArr(:), FFT_Data, ErrStatTmp )
  -               CALL SetErrStat(ErrStatTmp,'Error occured while applying the FFT to the second term of the difference QTF.', &
  -                              ErrStat,ErrMsg,RoutineName)
  -               IF ( ErrStat >= AbortErrLev ) THEN
  -                  IF (ALLOCATED(TmpData4D))        DEALLOCATE(TmpData4D,STAT=ErrStatTmp)
  -                  IF (ALLOCATED(DiffQTFForce))     DEALLOCATE(DiffQTFForce,STAT=ErrStatTmp)
  -                  IF (ALLOCATED(TmpDiffQTFForce))  DEALLOCATE(TmpDiffQTFForce,STAT=ErrStatTmp)
  -                  IF (ALLOCATED(TmpComplexArr))    DEALLOCATE(TmpComplexArr,STAT=ErrStatTmp)
  -                  RETURN
  -               END IF
  -
  -
  -                  ! Now we multiply the result by 2 and save it to the DiffQTFForce array and add the MnDrift term
  -               DO K=0,InitInp%NStepWave-1  ! bjj: added the "-1" here because TmpDiffQTFForce(InitInp%NStepWave) is not set and DiffQTFForce(InitInp%NStepWave,Idx) gets overwritten next, anyway
  -                  DiffQTFForce(K,Idx) = 2.0_SiKi * TmpDiffQTFForce(K) + MnDriftForce(Idx)
  -               ENDDO
  -
  -                  ! Copy the last first term to the first so that it is cyclic
  -               DiffQTFForce(InitInp%NStepWave,Idx) = DiffQTFForce(0,Idx)
  -
               ENDIF    ! Load component to calculate
  +         ENDDO ! ThisDim -- The 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
  +            TmpComplexArr(J,1:2) = MATMUL(RotateZMatrixT, TmpComplexArr(J,1:2))
  +            TmpComplexArr(J,4:5) = MATMUL(RotateZMatrixT, TmpComplexArr(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(  TmpDiffQTFForce(:),  TmpComplexArr(:,ThisDim), FFT_Data, ErrStatTmp )
  +            CALL SetErrStat(ErrStatTmp,'Error occured while applying the FFT to the second term of the difference QTF.', &
  +                           ErrStat,ErrMsg,RoutineName)
  +            IF ( ErrStat >= AbortErrLev ) THEN
  +               IF (ALLOCATED(TmpData4D))        DEALLOCATE(TmpData4D,STAT=ErrStatTmp)
  +               IF (ALLOCATED(DiffQTFForce))     DEALLOCATE(DiffQTFForce,STAT=ErrStatTmp)
  +               IF (ALLOCATED(TmpDiffQTFForce))  DEALLOCATE(TmpDiffQTFForce,STAT=ErrStatTmp)
  +               IF (ALLOCATED(TmpComplexArr))    DEALLOCATE(TmpComplexArr,STAT=ErrStatTmp)
  +               RETURN
  +            END IF
  +
  +
  +               ! Now we multiply the result by 2 and save it to the DiffQTFForce array and add the MnDrift term
  +               ! NOTE: phase shift and orientations on the MnDriftForce term have already been applied
  +               ! NOTE: the "-1" since TmpDiffQTFForce(InitInp%NStepWave) is not set and DiffQTFForce(InitInp%NStepWave,Idx) gets overwritten
  +            DO K=0,InitInp%NStepWave-1
  +               DiffQTFForce(K,Idx) = 2.0_SiKi * TmpDiffQTFForce(K) + MnDriftForce(Idx)
  +            ENDDO
  +
  +               ! Copy the last first term to the first so that it is cyclic
  +            DiffQTFForce(InitInp%NStepWave,Idx) = DiffQTFForce(0,Idx)
   
            ENDDO ! ThisDim -- The current dimension
         ENDDO    ! IBody -- This WAMIT body
  
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:                            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
  
  