load hclassic.dat
load results.dat
freq2 = hclassic(:,1);
results(1,3)=results(2,3); % to avoid NaN
qn_total = interp1(results(:,1),results(:,3),freq2);
noise_total = sqrt(hclassic(:,2).^2+qn_total.^2);

c = 299792458;
Gravity = 6.673*10^(-11);
solarmass = 1.989*10^(30);
mtotal = 1.4*solarmass+1.4*solarmass;
mreduce = 1.4*solarmass*1.4*solarmass/mtotal;
Mpc = 365.2422*24*60*60*c*3.26*10^6;

integ = 0;
for nn=1:find(freq2>1570,1);
    integ = integ + 1/(freq2(nn)^(7/3)*noise_total(nn)^2)*(freq2(nn+1)-freq2(nn));
end

range = (1/8)/sqrt(2)*2*(Gravity/c^3)^(5/6)*c*(5/96*mreduce)^(1/2)*(mtotal/pi^2)^(1/3)*sqrt(4*integ)/Mpc;

disp(['IR = ' num2str(range) 'Mpc']);