function LCGT_ASCinspiralrange()

%% Total A2DARM plot
Larm=3000;
DARMsensitivity_data=importdata('./DARMsensitivity/bLCGTDARMsensitivity_VDRSE.txt');
freq=DARMsensitivity_data(:,1);
DARMsensitivity_data(:,2)=DARMsensitivity_data(:,2)/Larm;

pitchfile='BRSE-negative-pitch-16-Feb-2015';
yawfile='BRSE-negative-yaw-16-Feb-2015';
pitchA2DARM_data=importdata(['./results/bLCGT-',pitchfile,'/totalA2DARM.txt']);
yawA2DARM_data=importdata(['./results/bLCGT-',yawfile,'/totalA2DARM.txt']);

pitchA2DARM=resampling(pitchA2DARM_data,freq,Larm);
yawA2DARM=resampling(yawA2DARM_data,freq,Larm);
lastfreqkk=length(pitchA2DARM(:,1));

pitchA2DARM(:,2)=pitchA2DARM(:,2)*10;
yawA2DARM(:,2)=yawA2DARM(:,2)*10;

totalDARMnoise=DARMsensitivity_data;
for kk=1:lastfreqkk
    totalDARMnoise(kk,2)=sqrt(DARMsensitivity_data(kk,2)^2+pitchA2DARM(kk,2)^2+yawA2DARM(kk,2)^2);
end

ir_fundamental=calcInspiralRange(DARMsensitivity_data)
ir_angular=calcInspiralRange(totalDARMnoise)

figure();

set(gcf,'DefaultLineLineWidth',3);
set(gcf,'defaultAxesFontSize',16);
set(gcf,'defaultTextFontSize',16);

loglog(freq,DARMsensitivity_data(:,2),pitchA2DARM(:,1),pitchA2DARM(:,2),yawA2DARM(:,1),yawA2DARM(:,2),freq,totalDARMnoise(:,2),'k');
legend(['official VDRSE IR=',num2str(ir_fundamental,'%.0f'),' Mpc'],'pitch angular noise','yaw angular noise',['total IR=',num2str(ir_angular,'%.0f'),' Mpc'],1);
grid on;
ylim([1e-25,1e-11]);
set(gca,'YTick',logspace(-25,-11,15));
xlim([1e-2,3e3]);
set(gca,'XTick',logspace(-2,3,6));
xlabel('Frequency [Hz]');
ylabel('Strain [1/rtHz]');
title('Total Angular Noise Coupling to DARM (pitch and yaw)')
saveas(gcf,['./results/inspiralrange/inspiralrange',pitchfile,'-',yawfile,'.png'])


function angA2DARM=resampling(angA2DARM_data,freq,Larm);
freq2=angA2DARM_data(:,1);
angA2DARM_data(:,2)=angA2DARM_data(:,2)/Larm;
for kk=1:length(freq)
    if freq2(end) < freq(kk)
        lastfreqkk=kk-1;
        break
    end
end
angA2DARM=change_sampling(angA2DARM_data,freq(1:lastfreqkk)');

function IR=calcInspiralRange(sensitivity);
% Constants and parameters for calculating IR
c=299792458;    % speed of light [m/s] 
G=6.673e-11;    % gravitational constant [m^3/kg/s^2]
solarmass=1.989e30;  % [kg]
Mpc=3.08568025e16*1e6;  % 1 Mpc [m]
mtotal=(1.4+1.4)*solarmass;  % total mass of a binary system [kg]
mreduce=1.4*solarmass*1.4*solarmass/mtotal;  % reduced mass of a binary system [kg]

fmin=10;    % low cut off frequency of integral [Hz]
fmax=c^3/(6^(3/2)*pi*G*mtotal);    % high cut off frequency of integral [Hz]

fmax

SNR=8;  % threshold of signal to noise ratio

% IR calculation (integral)
f=sensitivity(:,1);
h=sensitivity(:,2);
for kk=1:length(f)
    if fmin < f(kk);
        kkmin=kk;
        break
    end
end
for kk=1:length(f)
    if fmax < f(kk);
        kkmax=kk;
        break
    end
end

integ=0;
for kk=kkmin:kkmax
    integ=integ+1/(f(kk)^(7/3)*h(kk)^2)*(f(kk+1)-f(kk));
end

% KAGRA definition
%IR=(1/SNR)/sqrt(2)*2*(G/c^3)^(5/6)*c*(5/96*mreduce)^(1/2)*(mtotal/pi^2)^(1/3)*sqrt(4*integ)/Mpc;

% aLIGO/AdVirgo definition (sky average)
IR=sqrt(2)*0.442478*(1/SNR)/sqrt(2)*2*(G/c^3)^(5/6)*c*(5/96*mreduce)^(1/2)*(mtotal/pi^2)^(1/3)*sqrt(4*integ)/Mpc;