function [t,x,f]=bhchirp(m,r,tin,fs)
%
% E. Chassande-Mottin, 22/06/99.
%
% This function evaluates the waveform of a gravitational wave
% radiated by a coalescing binary of black holes (coalescence time =
% 0). It follows the toy model given in the articles:
% W. Anderson & R. Balasubramanian, "Time-frequency detection of
% gravitational waves," gr-qc ??.
% E. Flanagan & S. Hughes, Phys. Rev D 54, 4535 (1998).
%
% Input parameters:
% m    = mass of the objects (assumed to be of equal masses), in solar masses
% r    = distance binary-earth, in Mpc
% tin  = time to ISCO time, in seconds
% fs   = sampling frequency, in Hertz
%
% Output parameters:
% t    = time vector, in seconds
% x    = waveform
% f    = instantaneous frequency of x, in Hertz
%
% Typical example: Anderson/Balasubramanian,
%
% m=30; r=100; tin=1; fs=1e3; 
% [t,x,f] = bhchirp(m,r,tin,fs);
% plot(t,x); axis([t(1) t(length(t)) -1 1]);
% xlabel('time (seconds)'); ylabel('amplitude');

% some constants
to=4.96e-6; % (s)
fo=2.01e05; % (Hz)
ro=1.49;    % (km)
Mpc=3e19;   % (km)
alpha=1/4;
beta=5/8;

% total mass

M = 2*m;

% reduced mass
mu = M/2;

% chirp mass
cm=M^(2/5)*mu^(3/5);

% normalized distance
d=r*Mpc/ro;

% last frequency reached when the system is at its I(nnermost)
% S(table) C(ircular) O(rbit)

fisco = 205*20/M;
tisco = 3*cm^(-5/3)*(fisco/100)^(-8/3);

% duration of merging

tbaam = 50*M*to;

% duration of ringdown

fqnr = 1320*(20/M);
tqnr = 3.70/fqnr;
threshold=1e-3;
tdqnr=-tqnr*log(threshold);

% sampling error
if fqnr>(fs/2),
  error('The waveform is undersampled');
end;

% time settings
ts = 1/fs;
tins = (-tin:ts:0)';
tins = tins-tins(length(tins));   % t=0 at ISCO time
xins=zeros(length(tins),1);
fins=zeros(length(tins),1);
tmer = (0:ts:tbaam)';
xmer=zeros(length(tmer),1);
fmer=zeros(length(tmer),1);
trin = (0:ts:tdqnr)';
xrin=zeros(length(trin),1);
frin=zeros(length(trin),1);

% inspiral

dchirp = 160*3^(3/8)*cm^(-5/8);
A=3.3691e-21*cm^(5/4)/r;
%A=abs(beta-1)^(1/2)*(d*beta)^((-alpha+1/2)/(beta-1))
Ains = A*(abs(tins-tisco)).^(-alpha);
phins= 2*pi*dchirp*(abs(tins-tisco)).^(beta);
fins= dchirp*beta*(abs(tins-tisco)).^(beta-1);
xins=Ains.*cos(phins);

% merger

phimer0 = 0;
dt=tbaam;
a0 = fins(length(fins));
%a1=fo*(3/8*teta_isco^(-11/8)+0.363*5/8*teta_isco^(-13/8))/(160*pi*M^2*to); % df/dt(t_isco)
a1=dchirp*beta*(beta-1)*tisco^(beta-2);
a2=(3*fqnr-3*a0-2*a1*dt)/dt^2;
a3=(-2*fqnr+2*a0+a1*dt)/dt^3;

% phimer2 = 3/(640*M^2*to^2)*(tetains(length(tetains)).^(-11/8)+1399.46*tetains(length(tetains)).^(-13/8));
% amer = fqnr/(3*tmerge.^2)-(phimer1+2*phimer2*tmerge)/(3*tmerge.^2);

% bmer = -phimer2/(3*tmerge);
phimer = phimer0+2*pi*(a0*tmer+a1/2*tmer.^2+a2/3*tmer.^3+a3/4*tmer.^4);
fmer=a0+a1*tmer+a2*tmer.^2+a3*tmer.^3;
Amer=linspace(Ains(length(Ains)),0.0505*M/d,length(tmer))';
xmer=Amer.*cos(phimer);

tmer=tmer+ts; % tmer 1st sample corresponds to time t=ts;

% ringdown

phirin=2*pi*fqnr*trin+phimer(length(phimer));
frin=fqnr*ones(length(trin),1);
Arin=0.0505*M/d*exp(-trin/tqnr);
xrin=Arin.*cos(phirin);

trin=trin+tmer(length(tmer))+ts; % trin 1st sample corresponds to time=tmer(end)+ts

% full waveform

t=[tins;tmer;trin]; % gluing all time vectors
x=[xins;xmer;xrin];
f=[fins;fmer;frin];

% keyboard
