function hij = makeh(simId, ampl, startSamp, sampFreq)

% MAKEH - generates metric perturbation for GW source
%
% hij = makeh(simID, ampl, startSamp, sampFreq);
%
% simId         string containing signal information
% ampl          fixed signal amplitude
% startSamp     start time of (in samples) of signal; mostly interested in 
%               fractional part to properly line up signal samples with data 
%               samples
% sampFreq      sampling frequency
%
% hij           metric perturbation
%
% Currently supported signals are:
%       Gaussians                         'G'
%       sine-Gaussians                    'SG'
%       cosine-Gaussian                   'CG' 
%       BH ringdowns                      'BH'
%       sine                              'SIN'
%       cosine                            'COS'
%       sinc                              'SINC'
%       linear chirp                      'CHIRP'
%       square wave                       'SQUARE'
%       Gaussian modulate square wave     'SQUAREG'
%       sawtooth wave                     'SAWTOOTH'
%       Gaussian modulated sawtooth wave  'SAWTOOTHG'
%
% SEE ALSO: GRAVEN, SIMPARSE

% $Id: makeh.m,v 1.23 2004/05/17 21:21:01 stuver Exp $

hij = [];

% GET signal information from simId
% CHECK for errors
% GENERATE time vector

% Make simId all upper case for ease in using strcmp
simId = upper(simId);
% Get signal information from simId
[id, f0, tau, durationSec] = simparse(simId);

% Check that simparse did not return an empty value
if isempty(id) | isempty(f0) | isempty(tau) | isempty(durationSec)
    [id, f0, tau, durationSec]
    error('makeh: simparse returned an empty value')
end

% Check that sampFreq is > 0
if sampFreq <= 0 | ~isreal(sampFreq) | isnan(sampFreq) | isinf(sampFreq)
    sampFreq
    error('makeh: sampFreq must the > 0, real and finite')
end

% Check that durationSec is > 0
if durationSec <= 0 | ~isreal(durationSec) | isnan(durationSec) ...
        | isinf(durationSec)
    durationSec
    error('makeh: durationSec must be > 0, real and finite')
end

% Deterimine length of time vector in samples (floor trucated)
durationSamp = floor(durationSec*sampFreq);
% Check that durationSamp is not < 1
if durationSamp < 1
    durationSamp
    error('makeh: durationSamp is < 1')
end
% Time vector = times at which IFO samples waveform
ti = ((1:durationSamp)./sampFreq)+((1-mod(startSamp,1))/sampFreq);
% N.B. This accounts for intersample start times

% Check to make sure ti is not empty
if isempty(ti)
    ti
    error('makeh: ti is empty')
end

% GENERATE needed base functions

% The logical(sum(strcmp(id, {...}))) makes it so that if any one the
% strings match id, then the appropriate base is generated
if sum(strcmp(id, {'G', 'SG', 'CG', 'SQUAREG', 'SAWTOOTHG'})) > 0
    % A gaussian base needs to be generated
    
    % Check that tau is > 0
    if tau <= 0 | ~isreal(tau) | isnan(tau) | isinf(tau)
        tau
        error('makeh: tau must be > 0, real, and finite for Gaussians')
    end
    
    % Generate the Gaussian waveform that is also the base for SG and CG
    % Offset times for Gaussians
    tg = ti-(ti(end)/2);
    % Generate Gaussian
    gauss = exp(-(tg/tau).^2); 
end

if sum(strcmp(id, {'SIN', 'SG'})) > 0
    % A sine wave base needs to be generated
    
    % Check that f >= 0 for periodic function
    checkperiodic(f0);
    
    sine = ampl.*sin(2*pi*f0*ti);
end

if sum(strcmp(id, {'COS', 'CG'})) > 0
    % Check that f >= 0 for periodic function
    checkperiodic(f0);

    % A cosine wave base needs to be generated
    cosine = ampl.*cos(2*pi*f0*ti);
end

if sum(strcmp(id, {'SQUARE', 'SQUAREG'})) > 0
    % Check that f >= 0 for periodic function
    checkperiodic(f0);

    % A square wave base needs to be generated
    squ = ampl*square(2*pi*f0*ti);
end

if sum(strcmp(id, {'SAWTOOTH', 'SAWTOOTHG'})) > 0
    % Check that f >= 0 for periodic function
    checkperiodic(f0);

    % A sawtooth wave base needs to be generated
    saw = ampl*sawtooth(2*pi*f0*ti);
end

% GENERATE desired metric perturbation based on simId
% RETURN appropriate timeseries

switch id
    case 'G'
        hij = make33time(gauss);
    case 'SIN'
        hij = make33time(sine);
    case 'SG'
        hij = make33time(sine.*gauss);
    case 'COS'
        hij = make33time(cosine);
    case 'CG'
        hij = make33time(cosine.*gauss);
    case 'SINC'
        func = ampl*sinc(2*pi*f0*ti);
        hij = make33time(func);
    case 'CHIRP'
        % Generate a chirp signal that goes from 0 to f0 frequency
        func = ampl*chirp(ti, 0, ti(end), f0);
        hij = make33time(func);
    case 'SAWTOOTH'
        hij = make33time(saw);
    case 'SAWTOOTHG'
        hij = make33time(saw.*gauss);
    case 'SQUARE'
        hij = make33time(squ);
    case 'SQUAREG'
        hij = make33time(squ.*gauss);
    otherwise
        if strcmp(id, 'BH')
            % Construct damped, periodic source for BH ringdown
            % hij = ampl*exp(-t/tau)*[sin, cos, 0; cos, sin, 0; 0, 0, 0]
            
            % GENERATE periodic part of metric perturbation in 3x3xtime 
            % array format
            % MAKE BH ringdown
            
            % Initialize metric perturbation
            hij = zeros(3, 3, length(ti));
            
            % Generate vector of zeros in time dimension for use in 
            % constructing the metric perturbation
            oh = zeros(1, 1, length(ti));
            % Generate vector of the sine function in the time dimension for 
            % use in constructing the metric perturbation
            xx = reshape(sin(2*pi*f0*ti), [1, 1, length(ti)]);
            % Generate vector of the cosine function in the time dimension 
            % for use in constructing the metric perturbation
            xy = reshape(cos(2*pi*f0*ti), [1, 1, length(ti)]);
            
            % Assemble the part of the metric perturbation that deals with 
            % periodic functions
            periodic = [xx, xy, oh; xy, xx, oh; oh, oh, oh]; 
            % This is in full 3x3xtime array format
            
            % Multiply each time step of the periodic portion of the metric 
            % of the coresponding time step of the exponential and the 
            % amplitude
            for k = 1:length(ti)
                hij(:,:,k) = (ampl*exp(-ti(k)/tau)).*periodic(:,:,k);
            end
            
        else    
            id
            error(['makeh: Code currently only supports simId=''G'''...
                    ', ''SG'', ''CG'', ''BH'', ''SIN'', ''COS'','...
                    ' ''SINC'', ''SQUARE'', ''SQUAREG'', ''SAWTOOTHG'','...
                    ' ''SAWTOOTH'' & ''CHIRP'''])
        end
end

return

% ~~~ END MAIN ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

function checkperiodic(f0)

% CHECKPERIODIC - check the frequency of a periodic function for errors
%
% checkperiodic(f0)
%
% f0    if f0 is < 0, imaginary, NaN or Inf, the program errors out

% Check that F0 is > 0, real, and finite
if f0 < 0 | ~isreal(f0) | isnan(f0) | isinf(f0)
    f0
    error(['makeh: f0 must be positive, real and finite for periodic ' ...
            'functions'])
end

return

% ~~~ END CHECKPERIODIC ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

function hij = make33time(func)

% MAKE33TIME - put a timeseries vector into a 3x3xtime array 
%
% hij = make33time(func)
%
% func  function timeseries to be copied into the 3x3xtime array
%
% hij   3x3xtime array

% Put in full 3x3xtime array format
hij = repmat(reshape(func, [1, 1, length(func)]), [3, 3, 1]);
% reshape puts the elements of the func vector in the time dimension
%   of the 3x3xtime array 
% repmat replicates the time dimension in the 3x3xtime array 

return