anchor.diff:                          Omega1 = (J + K) * InitInp%WaveDOmega        ! the mth frequency -- \mu^- + n = m
                           Omega2 = K * InitInp%WaveDOmega              ! the nth frequency
   
  -
                              ! Find the Wave amplitudes 1 and 2
                           aWaveElevC1 = CMPLX( InitInp%WaveElevC0(1,J+K), InitInp%WaveElevC0(2,J+K), SiKi)  / InitInp%NStepWave2
                           aWaveElevC2 = CMPLX( InitInp%WaveElevC0(1,K),   InitInp%WaveElevC0(2,K),   SiKi)  / InitInp%NStepWave2
   
                              ! Set the (omega1,omega2,beta1,beta2) point we are looking for.
                           Coord4 = (/ REAL(Omega1,SiKi), REAL(Omega2,SiKi), InitInp%WaveDirArr(J+K), InitInp%WaveDirArr(K) /)
  +
  +                           ! Apply local Z rotation to heading angle (degrees) to put wave direction into the local (rotated) body frame
  +                        Coord4(3) = Coord4(3) - RotateZdegOffset
  +                        Coord4(4) = Coord4(4) - RotateZdegOffset
   
                              ! get the interpolated value for F(omega1,omega2,beta1,beta2)  --> QTF_Value
                           CALL WAMIT_Interp4D_Cplx( Coord4, TmpData4D, DiffQTFData%Data4D%WvFreq1, DiffQTFData%Data4D%WvFreq2, &
  
positive.code: 
           ! Wave information and QTF temporary
        COMPLEX(SiKi)                                      :: QTF_Value            !< Temporary complex number for QTF
        COMPLEX(SiKi), ALLOCATABLE                         :: TmpComplexArr(:)     !< Temporary complex array for frequency domain of one complete load component
        COMPLEX(SiKi)                                      :: TmpHMinusC           !< Temporary variable for holding the current value of \f$ H^- \f$
        COMPLEX(SiKi)                                      :: aWaveElevC1          !< Wave elevation of the first  frequency component.  NStepWave2 factor removed.
        COMPLEX(SiKi)                                      :: aWaveElevC2          !< Wave elevation of the second frequency component.  NStepWave2 factor removed.
  
negative.code: 
                          ! 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,ThisDim) = TmpHMinusC / 2.0_SiKi
  
                    ELSE     ! outside the frequency range, so
  
                       TmpComplexArr(J,ThisDim) = CMPLX(0.0_SiKi,0.0_SiKi,SiKi)
  
                    ENDIF    ! frequency check
  
  
                 ENDDO
              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
  