function [F hrss] = cal_fluence(file, alpha, doPlot, doParseval)
% function [F hrss] = cal_fluence(file, alpha)
%   file is of the form: t hp hx
%   alpha (optional) is the scale factor
%   returns F: isotropic equivalent fluence in units of ergs/cm^2
%   and hrss = ( sum(hp^2 + hx^2) * dt )^0.5
% To convert fluence into isotropic equivalent energy, we need to know the 
% distance to the source D.  If we know that, then
%   E = 4*pi*D^2 * F
% Note that 1 Mpc = 3.1e24 cm.
% Eric Thrane
%
% Note: there are two closely related scripts cal_fluence and cluster_fluence.
% The first calculates the fluence associated with a GW waveform described by
% a time series.  The second calculates the fluence associated with a ft-map
% cluster.

% debugging options
try
  doParseval;
catch
  doParseval = false;
end
try
  doPlot;
catch
  doPlot = false;
end

% set alpha=1 if not specified
try
  alpha;
catch
  alpha=1;
  fprintf('alpha not specified; setting alpha=1\n');
end

% load waveform file and parse data
g = load(file);
t = g(:,1);
hp = g(:,2);
hx = g(:,3);

% define sampling frequency
dt = t(2)-t(1);
Fs = 1/dt;

% calculate duration in case we want to compare with 1/df
T = t(end)-t(1);

% calculate FFTs
[hp_f f] = fft_eht(hp, Fs);
[hx_f f] = fft_eht(hx, Fs);
df = f(2)-f(1);

% diagnostic plot
if doPlot
  figure;
  loglog(f, abs(hp_f), 'b', 'linewidth', 2);
  hold on;
  loglog(f, abs(hx_f), 'r--', 'linewidth', 2);
  print('-dpng', 'cal_fluence');
end

% demonstrate Parseval's theorem to show correctness of frequency domain 
% calculation
if doParseval
  X = norm(hp)^2;
  Y = norm(hp_f)^2/numel(hp_f);
  fprintf(['Parseval check: norm(hp)^2 = %1.3e, norm(hp_f)^2/numel(hp_f)' ...
    ' = %1.3e\n'], X, Y);
end

% define constants in SI units
c = 3e8;
G = 6.67e-11;

% calculate isotropic equivalent fluence
% See Eq. 7.4 in the stochastic units document:
%   sgwb/trunk/doc/TechNotes/T070045/multipole.pdf
% See also Eq. 3.12-3.13 in Thrane et al. PRD 83, 083003 (2011)
% See also Eq. 4 in Abbott et al. ApJ Lett 701, L68 (2009)
% Note: in the frequency domain: sum(hdot^2) -> sum(f^2 htilde^2) * (2*pi)^2
F = (c^3/(16*pi*G)) * (2*pi)^2 * ...
  (sum(dt * f.^2 .* (abs(hp_f).^2 + abs(hx_f).^2))) / numel(hp_f);

% SI units of fluence are J/m^2; we convert to ergs/cm^2
% conversion factor = 1e7 / 100^2 = 1e3
F = 1000 * F;

% scale F by scale factor alpha
F = alpha * F;

% calculate hrss
hrss = sum(dt*(hp.^2 + hx.^2))^0.5;

return
