Dear salhus and dav-og
Many thanks for your work. I have some questions about infinite frequency added mass, especially for the calculation results by AQWA.

• I found both WAMIT and AQWA calcuate IRF and infinite frequency added mass with $R_{ij}(t)=\frac{2}{\pi}\int_0^\infty\frac{B_{ij}(\omega)}{\omega}\sin(\omega t)d\omega$ and $A_{ij}(\infty)=A_{ij}(\omega)-\int_{0}^{\infty}R_{ij}(t)\cos(\omega t)dt$ (1), but radiationIRF.m uses the equations $R_{ij}(t)=\frac{2}{\pi}\int_{0}^{\infty}B(\omega)\cos(\omega t)d\omega$ and $A_{ij}(\infty)=A_{ij}(\omega)+\frac{1}{\omega}\int_{0}^{\infty}R_{ij}(t)\sin(\omega t)dt$ proposed in State-Space Realization of the Wave-Radiation Force within FAST. Are they identical? Moreover, I found additional $\omega$ term in radiationIRF.m: $R_{ij}(t)=\frac{2}{\pi}\int_{0}^{\infty}B(\omega)\cos(\omega t)\omega d\omega$. Is it a mistake, or an equation from other reference?

% WAMIT and AQWA
hydro.ra_K(i,j,:) = (2/pi)*trapz(w,ra_B.*(sin(w.*t(:))./w), 2);
ra_Ainf_temp(k,1)  = ra_A(k) - trapz(t,ra_K.*cos(w(k).*t.'));

% radiationIRF.m
hydro.ra_K(i,j,:) = (2/pi)*trapz(w,ra_B.*(cos(w.*t(:)).*w), 2);
ra_Ainf_temp(k,1)  = ra_A(k) + (1./w(k))*trapz(t,ra_K.*sin(w(k).*t.'));

• For AQWA, wave frequency range is highly depend on the mesh quality, and generally the maximum frequency in AQWA is much smaller than that in WAMIT. Therefore, IRF and infinite frequency added mass are calcualted based on truncated frequency range $[{\Omega_{1}},{\Omega_{2}}]$, where $\Omega_{1}\simeq0$. Truncated IRF is denoted by $R_{ij}^{r}(t)=\frac{2}{\pi}\int_{\Omega_{1}}^{\Omega_{2}}\frac{B_{ij}(\omega)}{\omega}\sin(\omega t)d\omega$. I attempt to correct the IRF and infinite frequency added mass within $[{\Omega_{2}},\infty]$ using the method proprosed in Computation of wave effects using the panel method: $R_{ij}^{c}(t)\simeq-\frac{4}{\pi^{2}\Omega_{2}}\left[\cos\Omega_{2} t+\Omega_{2} t\sin(\Omega_{2} t)\right]\int_{0}^{\Omega_{2}}B_{i j}(\omega)d\omega$, $A_{ij}(\infty)=A_{ij}(\omega)-\int_{0}^{\infty}[R_{ij}^{r}(t)+R_{ij}^{c}(t)]\cos(\omega t)dt$ (2), but results are unsatisfactory. I therefore compared the infinite frequency added mass calculated by WAMIT (pink), AQWA (blue), old radiationIRF(AQWA) (green), new radiationIRF(AQWA) with interpolated added mass (grey), AQWA results calculated by equation (1) (orange), and AQWA results calculated by equation (2) (yellow)

% Computation of wave effects using the panel method
hydro.ra_K(i,j,:) = (2/pi)*trapz(w,ra_B.*(sin(w.*t(:))./w), 2);
ra_K_correct(i,j,:) = (-4/(pi^2*wMax))*(cos(wMax*t)+wMax*t.*sin(wMax*t))*trapz(w,ra_B);
ra_Ainf_temp(k,1)  = ra_A(k) - trapz(t,ra_K.*cos(w(k).*t.'));

• The correction results calculated based on equation (2) (yellow) is obviously wrong. I presume it is caused by formula itself, as the formula assumes that ${\Omega_{2}}$ is large enough, or is caused by difference in dimension.
• I think AQWA has calculated infinite frequency added mass as ADDEDMASS_HF (blue), with similar accuracy compared with new radiationIRF.m (grey). Therefore, there is no need to calculate infinite frequency added mass for AQWA in radiationIRF.m

if isempty(hydro.Ainf) == 1 || isfield(hydro,'Ainf') == 0 || strcmp(hydro(F).code,'WAMIT') || strcmp(hydro(F).code,'AQWA')==0

@hachikoi1,

Thanks for the discussion!

re: Kernel differences,

I am assuming you are using a sufficiently small $\Delta\omega$ and by WAMIT generated IRFs you mean using the f2t utility.

Here are some thoughts on the kernel differences you pointed out. I will respond to the comparisons in a separate comment if the comments that follow do not make any difference to your calculations/comparisons.
In the WAMIT manual, notice that the second part of the radiation IRF Eq. $13.1$ is written for accelerations,

Usually the IRF is represented using the velocity [see 1, 2, 3, and 4], and that should explain the differences between the kernels.

So it really matters, if the IRF kernel was written so as to be convoluted with device velocity or accelerations.
For your reference, I have some personal scripts (that are not part of the WEC-Sim repositories and come with no warranty or support from WEC-Sim) that can run and extract the f2t results from WAMIT (all operated from within MATLAB but you can use the cmd for the same). They work for the $.1, .2$ results in the JR and KR formats, refer here

re: $\omega$ term in radiationIRF.m

Note, the radiation damping terms are normalized by $\omega\rho$. While $\rho$ is a scalar constant and can be multiplied later, the $omega$ term does need to be in the kernel.

[1] Falnes, J., 2002, Ocean Waves and Oscillating Systems: Linear Interactions Including Wave-Energy Extraction, Cambridge University Press, Cambridge.
[2] Korde, U. A., and Ringwood, J., 2016, Hydrodynamic Control of Wave Energy Devices, Cambridge University Press.
[3] Cruz, J., ed., 2008, Ocean Wave Energy: Current Status and Future Prepectives [i.e. Perspectives], Springer, Berlin.
[4] Folley, M., 2016, Numerical Modelling of Wave Energy Converters State-of-the-Art Techniques for Single Devices and Arrays, Academic Press is an imprint of Elsevier, London, UK.

** Run F2T **

`
fid=fopen('inputs.f2t','w');
fprintf(fid,'Frequency to Time Domain\n');
fprintf(fid,'cylR1D1\n');
fprintf(fid,'0,0,1 \n'); % (IRAD IDIFF NUMHDR)
fprintf(fid,'0,0,1 \n'); % (INUMOPT5 INUMOPT6 IPEROUT)
fprintf(fid,'0.01,2000 \n') ; % (DT NT time step and number of time steps)
fprintf(fid,'0 \n') ; % (IOUTFCFS, output both cosine and sine transforms)
fprintf(fid,'1,9.8 ULEN GRAV\n') ;
status=fclose(fid);
%%

system('C:\WAMITv7\utils\f2t\f2t.exe')

`

** Read Radiation IRF from _IR in .1 format **

`clc;clear;close all

mode = 3;
fname = '';

n_angle = 1;
n_body = 1; %number of bodies
modes_tmp = [1 2 3 4 5 6];
modes = [];
for i = 1:n_body
modes = [modes modes_tmp];
end
n_modes = length(modes);

fileID = fopen([fname '_IR.1']);
irf_rad_out = textscan(fileID,'%f %f %f %f %f','HeaderLines',1);
fclose(fileID);

irf_rad.data = cell2mat(irf_rad_out);

n_modes = find(and(irf_rad.data(2:end,2)==irf_rad.data(1,2), ...
irf_rad.data(2:end,3)==irf_rad.data(1,3)),1);

I = irf_rad.data(1:n_modes,2);
J = irf_rad.data(1:n_modes,3);

I0 = find(irf_rad.data(:,1)>0,1,'first');

nt = size(irf_rad.data(I0:end,:),1)/n_modes;

irf_tmp = zeros(6, 6, nt);
irf_time = zeros(6,6,nt);

for ii = 1:n_modes
ind = ii + I0-1;
irf_tmp(I(ii),J(ii),:) = irf_rad.data(ind:n_modes:end,5);
irf_time(I(ii),J(ii),:) = irf_rad.data(ind:n_modes:end,1);
end

irf.Radirf_dot = squeeze(irf_tmp(mode,mode,:));
irf.radt = squeeze(irf_time(mode,mode,:));
irf.dt = abs(abs(irf.radt(2)) - abs(irf.radt(1)));
irf.Radirf = (diff(irf.Radirf_dot)./irf.dt );
irf.RadirfTime = irf.radt(2:end);
%%

% RHO = 1000;
% %%
% dt=0.1;Tend=20;
% tnum=0:dt:Tend;dw=w(2)-w(1);
% Bmode=squeeze(B(mode,mode,:));iB_irf=zeros(size(tnum));
% for i = 1:length(tnum)
% for j = 1:length(w)
% iB_irf(i) = iB_irf(i) + (2/(RHOpiw(j))).*Bmode(j).*sin(w(j)*tnum(i)).*dw ;
% end
% end
% %%
% B_irf=zeros(size(tnum));
%
% for i = 1:length(tnum)
% for j = 1:length(w)
% B_irf(i) = B_irf(i) + (2/pi).*Bmode(j).*cos(w(j)*tnum(i)).*dw ;
% end
% end
clearvars -except irf
%%
figure(1)
plot(irf.radt,irf.Radirf_dot *1000,'-k','LineWidth',2)
grid on
hold on
% plot(tnum,iB_irf,'-.r','LineWidth',2)
legend('f2t','numerical integration')
figure(2)
plot(irf.RadirfTime,irf.Radirf *1000)
grid on
hold on
% plot(tnum,B_irf,'-.r','LineWidth',2)
legend('f2t','numerical integration')`

** Read Radiation IRF in _JR .1 format **

`clc;clear;close all
mm = 3;
n=mm;m=mm;
fname = '';

load(strcat(sprintf('%s',fname),'.mat'))
n_angle = 1;
n_body = 1; %number of bodies
fileID = fopen([fname '_JR.1']);
irf_rad_out = textscan(fileID,'%f','delimiter',' ','HeaderLines',1);
fclose(fileID);

JR1irf.data = rmmissing(cell2mat(irf_rad_out));
% For both cosine transform and sine transform
JR1irf.mode = 2*((n-1)6 + m + 0);
JR1irf.spacer = 2(6n_body)^2 + 1;
% For only cosine trasnform or sine trasnform
% JR1irf.mode = 1((n-1)6 + m + 1);
% JR1irf.spacer = 1(6*n_body)^2 + 1;
% Write Data
JR1irf.time_for_dot = JR1irf.data(1:JR1irf.spacer:end);
JR1irf.nmdata_dot = JR1irf.data(JR1irf.mode :JR1irf.spacer:end);
JR1irf.dt = abs(abs(JR1irf.time_for_dot(2)) ...
- abs(JR1irf.time_for_dot(1)));
JR1irf.radIRF = diff(real(JR1irf.nmdata_dot))./JR1irf.dt;
JR1irf.radIRF = JR1irf.radIRF(2:end);
JR1irf.time = JR1irf.time_for_dot(3:end);
figure()
plot(JR1irf.time_for_dot,JR1irf.nmdata_dot)

figure()
plot(JR1irf.time,JR1irf.radIRF *1000)

clearvars -except JR1irf`

** Read Radiation IRF in _KR .1 format **

`clc;clear;close all
mm = 3;
n=mm;m=mm;
fname = '';

load(strcat(sprintf('%s',fname),'.mat'))
n_angle = 1;
n_body = 1; %number of bodies
fileID = fopen([fname '_KR.1']);
irf_rad_out = textscan(fileID,'%f','delimiter',' ','HeaderLines',1);
fclose(fileID);

KR1irf.data = rmmissing(cell2mat(irf_rad_out));
% For both cosine transform and sine transform
KR1irf.mode = 2*((n-1)6 + m + 0);
KR1irf.spacer = 2(6n_body)^2 + 1;
% For only cosine trasnform or sine trasnform
% KR1irf.mode = 1((n-1)6 + m + 1);
% KR1irf.spacer = 1(6*n_body)^2 + 1;
% Write Data
KR1irf.time = KR1irf.data(1:KR1irf.spacer:end);
KR1irf.nmdata = KR1irf.data(KR1irf.mode :KR1irf.spacer:end);

figure()
plot(KR1irf.time(3:end),real(KR1irf.nmdata(3:end)))

clearvars -except KR1irf

`

** Read Excitation Force IRF in _IR .2 format **

`clc;clear;close all

mode = 3;
fname = '';

%%

n_angle = 1;
n_body = 1; %number of bodies
modes_tmp = [1 2 3 4 5 6];
modes = [];
for i = 1:n_body
modes = [modes modes_tmp];
end

fileID = fopen([fname '_IR.2']);
irf_rad_out = textscan(fileID,'%f %f %f %f %f','HeaderLines',1);
fclose(fileID);
irf_IR2.data = cell2mat(irf_rad_out);
n_modes = find(and(irf_IR2.data(2:end,2)==irf_IR2.data(1,2), ...
irf_IR2.data(2:end,3)==irf_IR2.data(1,3)),1);

I0 = find(irf_IR2.data(:,1)==0,1,'first');
nt = size(irf_IR2.data(I0:end,:),1)/n_modes;

%% Package it all up

IR2.time = irf_IR2.data(mode:length(modes):end,1);
IR2.Eirf = irf_IR2.data(mode:length(modes):end,4);
IR2.Ephase = irf_IR2.data(mode:length(modes):end,5);

close all
clearvars -except IR2

%%
figure()
plot(IR2.time,IR2.Eirf10009.81,'-k','LineWidth',2)
grid on
figure()
plot(IR2.time,IR2.Ephase,'--b','LineWidth',2)
grid on
hold on
legend('Eirf f2t IR2','Ephase f2t IR2')`

** Read Excitation Force IRF in JR .2 format **

`clc;clear;close all

mode = 3;
fname = '';

%%

n_angle = 1;
n_body = 1; %number of bodies
fileID = fopen([fname '_JR.2']);
irf_rad_out = textscan(fileID,'%f %f %f %f %f %f %f','HeaderLines',1);
fclose(fileID);
irf_JR2.data = cell2mat(irf_rad_out);
n_modes = find(and(irf_JR2.data(2:end,2)==irf_JR2.data(1,2), ...
irf_JR2.data(2:end,3)==irf_JR2.data(1,3)),1);
I0 = find(irf_JR2.data(:,1)==0,1,'first');
nt = size(irf_JR2.data(I0:end,:),1)/n_modes;

%% Package it all up

JR2.time = irf_JR2.data(:,1);
JR2.Firf = irf_JR2.data(:,mode+1);

close all
clearvars -except JR2

%%
figure(1)
plot(JR2.time,JR2.Firf10009.81,'-k','LineWidth',2)
grid on
legend('Firf f2t JR2')
`

Dear salhus
Thank you for your sharing. I understand the differences between WAMIT/AQWA with WEC-Sim.
Have a nice day!
best regards,

Thanks @hachikoi1 !
Let me know if you would like to explore this more or happen to find any other issues.
Cheers,
Sal

Dear salhus
Your reply was clear and I no longer have any doubts.
Best regards,