function [Sim, state] = graven(varargin)

% GRAVEN - simulate a set of GW signals for addition to a data stream
%
% [Sim, state] = graven(logFileName, gps, detId, inputFileName, ...
%           sampFreq, seed);
%   Generate signals based on information in file inputFileName
% [Sim, state] = graven(nSim, logFileName, gps, detId, inputFileName, ...
%           sampFreq, seed);
%   Generate signals based on nSim random draws from information in file
%   inputFileName
% [Sim, state] = graven(nSim, logFileName, gps, detId, simId, ampl, ...
%           startSamp, internal, external, sampFreq, seed);
%   Generate nSim signals for fixed detId, simId, and ampl; random start in
%   range startSamp, and either fixed or random internal and external
% [Sim, state] = graven(nSim, logFileName, gps, detId, inputFileName, ...
%           startSamp, external, sampFreq, seed);
%   Push a pregenerated waveform contained in inputFileName name through the
%   detector projection function and the inverse claibrator nSim times
%
% logFileName       string with name of log file (must include file 
%                   extention)
% detId             detector identifier string; type 'help getifo' for list 
%                   of currently supported detectors 
% gps               GPS start time to allow selection of calibration
% inputFileName     string with name of input file (must include file 
%                   extention) or cell containing file name of
%                   pre-generated waveform (with extention) and name of 
%                   variable or field path of timeseries 
%                       e.g. {'inputFileName','variableName'}
% nSim              positive integer indicating how many simulations to 
%                   generate
% simId             signal type string; type 'help makeh' for list of
%                   currently supported signals
% ampl              signal amplitude or amplitude range
%                   *When ampl is a scalar, the amplitude is fixed
%                   *When ampl is a 1x2 vector, the amplitude is randomly 
%                   chosen between ampl(1) and ampl(2) in a logarithmic 
%                   distribution. N.B. ampl(1)=ampl(2) is equivilent to 
%                   entering a fixed amplitude
%                   *When ampl is a 1x3 vector, the amplitdue is randomly 
%                   chosen from the discrete linear range 
%                   ampl(1):ampl(2):ampl(3)
% startSamp         1x2 vector defining range of allowed start times in
%                   samples wrt the beginning of data set
% internal          struct or vector with internal orientation angles, or 
%                   string 'RANDOM' or 'OPTIMAL'; if struct, must have 
%                   fields x and phi
% external          struct or vector with sky location & polarization angle 
%                   or string 'RANDOM' or 'OPTIMAL'; if struct, must have 
%                   fields x, phi and psi
% sampFreq          sampling frequency of detId
% seed              random number generator seed (state can be fed back in
%                   as seed)
%
% Sim(nSim)    struct array of simuluations
%  .name        simId string
%  .nstr        start sample at detector
%  .sig         signal to add starting at sample nstr
%  .gps0        signal start time in gps sec, ns at detector
%  .pkgps       signal peak time in gps sec, ns at detector
% state         state of the random number generator at the end of the code
%               (can be fed back into graven as seed)'
%
% SEE ALSO: MAKEH, MAKETT, DETPROJ, CALIBSIMFD, IFODELAY, PREGENLIST

% $Id: graven.m,v 1.38 2004/05/17 21:16:49 stuver Exp $

% DETERMINE number of input arguments
% CHECK variables for error
% DETERMINE mode of operation from number of inputs
% SET preGen flag to false
% IF mode 1 (6 input arguments)
%   ASSIGN variable names to input argumnets
%   SET random number generator seed
%   MAKE list reading line-by-line from a file
%   RETURN list
% ELSEIF mode 2 (7 input argumnets)
%   ASSIGN variable names to input arguments
%   CHECK variables for error
%   SET random number generator seed
%   MAKE list by drawing on file nSim times
%   RETURN list
% ELSEIF mode 3 (11 input arguments)
%   ASSIGN variable names to input arguments
%   SET random number generator seed
%   MAKE list based directly on input arguments
%   RETURN list
% ELSEIF mode 4 (9 input arguments)
%   ASSIGN variable names to input arguments
%   SET random number generator seed
%   MAKE list for pre-generated waveforms
% ENDIF

% Read length of varargin
inLength = length(varargin);

% Check for errors
input = varargin;
errorcheck(input);

preGen = ~true;
% Assign variable names to components of varargin & generate list
if inLength == 6
    [logFileName, gps, detId, inputFileName, sampFreq, seed] = ...
        deal(varargin{:});
    % Seed rand here
    rand('state', seed);
    % Generate simulation list
    list = makelist(inputFileName, detId);
    % list is a cell array
elseif inLength == 7
    [nSim, logFileName, gps, detId, inputFileName, sampFreq, seed] = ...
        deal(varargin{:});
    % Seed rand here
    rand('state', seed);
    % Generate simulation list
    list = makelist(nSim, inputFileName, detId);
    % list is a cell array
elseif inLength == 11
    [nSim, logFileName, gps, detId, simId, ampl, startSamp, internal, ...
            external, sampFreq, seed] = deal(varargin{:});
    % Seed rand here
    rand('state', seed);
    % Generate simulation list
    list = makelist(nSim, simId, ampl, startSamp, internal, external, detId);
    % list is a cell array
elseif inLength == 9
    [nSim, logFileName, gps, detId, inputFileName, startSamp, external, ...
            sampFreq, seed] = deal(varargin{:});    
    % Seed rand here
    rand('state', seed);
    list = pregenlist(nSim, inputFileName, startSamp, external, detId);
    % list is a cell array
else % Error
    error('graven: Number of input arguments can only be 6, 7, 9 or 11')
end

% DETERMINE number of lines in list
% INITIALIZE structures
% INITIALIZE cache for filter coefficents
% OPEN log file
% WRITE headerlines to log file

% Read length of List (each row is a simulation)
listLength = size(list, 1);

if listLength == 0
    error('graven: listing error, list has zero length')
end

% Initialize Sim, Internal and External structures
[Sim(1:listLength)] = deal(struct('name', [], 'nstr', [], 'sig', [], ...
    'gps0', [], 'pkgps' ,[]));
Internal = deal(struct('x', [], 'phi', []));
External = deal(struct('x', [], 'phi', [], 'psi', []));

% Initialize alpha and beta cache
abcache = []; 

% Open log file w/ append option to prevent over-writing
fid = fopen(logFileName, 'a');
if fid == -1
    error('graven: Failed to open log file')
end

date = clock;
% Write headerlines to log file
fprintf(fid, ['# This simulation was run on %2.0f/%2.0f/%4.0f ' ...
        'at %2.0f:%2.0f:%2.2f EDT\n'], date(2), date(3), date(1), ...
    date(4), date(5), date(6));
fprintf(fid, ['# (simId) (ampl) (startSamp1) (startSamp2) (Internal.x) '...
        '(Internal.phi) (External.x) (External.phi) (External.psi)' ...
        ' (gps) (gps0 s) (gps0 ns) (pkampl)  (pkgps s) (pkgps ns) '...
        '(detId)\n']);

% INITIALIZE FOR loop to run through the list
%   ASSIGN variable names for this simulation
%   DETERMINE number of samples to add to startSamp to account for the GW's
%             time of flight from the center of the Earth then add to 
%             startSamp
%   IF preGen flag is false
%       GENERATE metric perturbation
%       PROJECT metric perturbation into TT gauge
%   ENDIF
%   PROJECT TT gauge perturbation onto antenna pattern
%   DETERMINE the maximum strain and time of max strain
%   CALIBRATE strain into AS_Q counts
%   SAVE simulation into Sims struct
%   WRITE simulation information to log file
% ENDFOR
% CLOSE log file
% SAVE random number generator state 

RAD2DEG = 180/pi;

for iSim = 1:listLength
    % Read list:
    simId = list{iSim, 1};
    ampl = list{iSim, 2};
    startSamp = list{iSim, 3};
    Internal.x = list{iSim, 4};
    Internal.phi = list{iSim, 5};
    External.x = list{iSim, 6};
    External.phi = list{iSim, 7};
    External.psi = list{iSim, 8};
    
    if sampFreq ~= 16384 & (upper(detId) == 'H1' | upper(detId) == 'H2' ...
            | upper(detId) == 'L1')
        warning(['graven: sampFreq does not match the sampling frequency ' ...
                'of the detId entered'])
    end
    
    % If ampl = -1, this simulation is from a pre-generated waveform
    if ampl == -1
        % Reset preGen 
        preGen = true;
    end    
    
    % startSamp is the time in samples the GW is incident on the center of 
    % the Earth.  Calculate the number of samples to ADD to this startSamp to 
    % be incident on detId
    samp = ifodelay(detId, External.x, External.phi, sampFreq);
    if isempty(samp)
        error('graven: ifodelay returned empty samp')
    end
    stimeIfo = startSamp+samp;
    
    % If this is a pre-generated waveform, the start time of the simulation
    % at the detector is rounded to the nearest integer sample since
    % inter-sample start times are not currently supported for pre-generated
    % waveforms.
    if preGen
        stimeIfo = round(stimeIfo);
        htt = simidfile(simId);
    end
    
    % Assume pre-generated waveforms are already in the TT gauge
    if ~preGen
        % Generate metric perturbation
        hij = makeh(simId, ampl, stimeIfo, sampFreq); 
        if isempty(hij)
            error('graven: makeh returned empty hij')
        end
        % Project perturbation into TT gauge
        htt = makett(hij, Internal); 
        if isempty(htt)
            error('graven: makett returned empty htt')
        end
    end
    
    % Project TT perturbation onto the detector
    h = detproj(htt, detId, External);
    if isempty(h)
        error('graven: detproj returned empty h')
    end
    
    % Find the magnitude and time of the maximum strain
    [maxStrain, maxSamp] = max(h);
    maxTimeNdx = ceil(stimeIfo)+maxSamp-1; % Time of max strain in samples
    % Convert into sec and ns
    maxTimeGps = ndx2gps(maxTimeNdx, sampFreq, gps); 
    
    
    % Convert to AS_Q counts
    [asqSeries, abcache] = calibsimfd(detId, fix(gps+(stimeIfo/sampFreq)), h, ...
        sampFreq, abcache);
    
    if isempty(asqSeries)
        error('graven: calibsimfd returned empty asqSeries')
    end
    % Save start time, time-series, GPS time and max strain time for sim    
    Sim(iSim).name = upper(simId);
    Sim(iSim).nstr = ceil(stimeIfo);
    Sim(iSim).sig = asqSeries;
    Sim(iSim).gps0 = ndx2gps(stimeIfo, sampFreq, gps);
    Sim(iSim).pkgps = maxTimeGps;
    
    % Write information for this sim to log file 
    fprintf(fid, ['%s %0.5g %9.15f %9.15f %d %d %d %d %d %d %d ' ...
            '%8.0f %d %d %8.0f %s \n'], simId, ampl, startSamp, startSamp, ...
        Internal.x, Internal.phi, External.x, External.phi, ...
        External.psi, gps, Sim(iSim).gps0.sec, Sim(iSim).gps0.ns, ...
        maxStrain, Sim(iSim).pkgps.sec, Sim(iSim).pkgps.ns, detId);
end

% Close logfile
fclose(fid);

% Save random number generator state
state = rand('state');

return

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

function errorcheck(input)

% ERRORCHECK - checks graven's input data for errors and stops the code if
%              an error is found
%
% errorcheck(input)
%
% input     the varargin input argument to graven

% DETERMINE number of input arguments
% DETERMINE mode by number of inputs
% IF mode 1  
%   ASSIGN variables 
%   CHECK variables that are not used in all modes for error
% ELSEIF mode 2
%   ASSIGN variables 
%   CHECK variables that are not used in all modes for error
% ELSEIF mode 3
%   ASSIGN variables 
%   CHECK variables that are not used in all modes for error
% ELSEIF mode 4
%   ASSIGN variables 
%   CHECK variables that are not used in all modes for error
% ENDIF
% CHECK variables that are used by all modes for error

inputLength=length(input);

if inputLength == 6 % mode 1
    [logFileName, gps, detId, inputFileName, sampFreq, seed] = ...
        deal(input{:});
    
    % Check inputFileName for errors since not all of the modes require
    % inputFileName
    
    % inputFileName must be a string
    if ~strcmp(class(inputFileName), 'char')
        inputFileName
        error('graven: inputFileName must be a string')
    end
    
elseif inputLength == 7 % mode 2
    [nSim, logFileName, gps, detId, inputFileName, sampFreq, seed] = ...
        deal(input{:});
    
    % Check inputFileName and nSim for errors since not all modes 
    % require them
    
    % inputFileName must be a string
    if ~strcmp(class(inputFileName), 'char')
        inputFileName
        error('graven: inputFileName must be a string')
    end
    
    % nSim must be a real positive interger
    if nSim < 1 | mod(nSim, 1) > 0 | isinf(nSim) | isnan(nSim)
        % nSim must be real for mod
        nSim
        error('graven: nSim must be a real finite integer >= 1')
    end
    
elseif inputLength == 11
    [nSim, logFileName, gps, detId, simId, ampl, startSamp, internal, ...
            external, sampFreq, seed] = deal(input{:});
    
    % Check nSim and variables that only mode 3 uses (simId, etc) for errors
    
    % nSim must be a real postive integer
    if nSim < 1 | mod(nSim, 1) > 0 | isinf(nSim) | isnan(nSim)
        % nSim must be real for mod
        nSim
        error('graven: nSim must be a real finite integer >= 1')
    end
    
    % simId must be a string
    if ~strcmp(class(simId), 'char')
        simId
        error('graven: simId must be a string')
    end
    
    % ampl must be a scalar, 1x2 or 1x3 vector with real positive numbers
    if size(ampl, 1) == 1 & size(ampl, 2) == 1 % scalar
        if ampl < 0 | isinf(ampl) | isnan(ampl) | ~isreal(ampl)
            % ampl must be a real postive number
            ampl
            error('graven: ampl must be a real positive fitite number')
        end
    elseif size(ampl, 1) == 1 & size(ampl, 2) == 2 % 1x2 vector
        if ampl < 0 | isinf(ampl(1)) | isinf(ampl(2)) | isnan(ampl(1)) ...
                | isnan(ampl(2)) | ~isreal(ampl(1)) | ~isreal(ampl(2))
            % All components of the ampl vector must be real positive numbers
            ampl
            error(['graven: ampl components must be real positive finite '...
                    'numbers'])
        end
    elseif size(ampl, 1) == 1 & size(ampl, 2) == 3 % 1x3 vector
        if ampl < 0 | isinf(ampl(1)) | isinf(ampl(2)) | isinf(ampl(3)) | ...
                isnan(ampl(1)) | isnan(ampl(2)) | isnan(ampl(3)) | ...
                ~isreal(ampl(1)) | ~isreal(ampl(2)) | ~isreal(ampl(3))
            % All components of the ampl vector must be real positive numbers
            ampl
            error(['graven: ampl components must be real positive finite '...
                    'numbers'])
        end
    else
        error('graven: ampl must be a scalar or a 1x2 or 1x3 vector')
    end
    
    % startSamp must be a 1x2 vector and contain real postive integers
    if size(startSamp, 1) ~= 1 | size(startSamp, 2) ~= 2
        startSamp
        error('graven: startSamp must be a 1x2 vector')
    elseif startSamp(1) < 1 | mod(startSamp(1), 1) > 0 | ...
            mod(startSamp(2), 1) > 0 | isinf(startSamp(1)) | ...
            isinf(startSamp(2)) | isnan(startSamp(1)) | isnan(startSamp(2))
        % startSamp must be real for mod
        % startSamp(1) and startSamp(2) must be real postive integers
        startSamp
        error(['graven: startSamp(1) and startSamp(2) must be real '...
                'finite integers >= 1'])
    end
    
    % internal must be a vector, string or struct
    if strcmp(class(internal), 'double')
        % if internal is a double, it must be a 1x2 vector
        if size(internal, 1) ~= 1 | size(internal, 2) ~= 2
            internal
            error(['graven: internal must be a 1x2 vector when not a ' ...
                    'struct or string'])
        elseif abs(internal(1)) > 1 | ~isreal(internal(1))
            % internal(1) must be real and bewteen -1 and 1
            internal
            error('graven: internal(1) must be real and bewteen -1 and 1')
        elseif internal(2) < -pi | internal(2) > pi | ~isreal(internal(2))
            % internal(2) must be real and between -pi and pi
            internal
            error('graven: internal(2) must be real and between -pi and pi')
        end
    elseif strcmp(class(internal), 'char')
        % if internal is a string, it must be either 'RANDOM' or 'OPTIMAL'
        if ~strcmp(upper(internal), 'RANDOM') & ...
                ~strcmp(upper(internal), 'OPTIMAL') 
            internal
            error(['graven: internal must be ''RANDOM'' or ''OPTIMAL'' ' ...
                    'when a string'])
        end
    elseif strcmp(class(internal), 'struct')
        % if internal is a struct is must contain fields x and phi
        if ~isfield(internal, 'x')
            internal
            error('graven: internal must contain field x when a struct')
        elseif abs(internal.x) > 1 | ~isreal(internal.x)
            % if internal.x exists, it must have a value between -1 and 1
            internal.x
            error('graven: internal.x must be real and between -1 and 1')
        end
        if ~isfield(internal, 'phi')
            internal
            error('graven: internal must contain field phi when a struct')
        elseif internal.phi < -pi | internal.phi > pi | ~isreal(internal.phi)
            % if internal.phi exits, it must be real and have a value between
            % -pi and pi
            internal.phi
            error(['graven: internal.phi must be real and betweeen -pi '...
                    'and pi'])
        end
    end
    
    % Check external angles for error
    checkexternal(external);
    
elseif inputLength == 9
    [nSim, logFileName, gps, detId, inputFileName, startSamp, external, ...
            sampFreq, seed] = deal(input{:});
    
    % inputFileName must be a string
    if ~strcmp(class(inputFileName), 'char') & ~iscell(inputFileName)
        inputFileName
        error('graven: inputFileName must be a string')
    end
    
    % startSamp must be a 1x2 vector and contain real postive integers
    if size(startSamp, 1) ~= 1 | size(startSamp, 2) ~= 2
        startSamp
        error('graven: startSamp must be a 1x2 vector')
    elseif startSamp(1) < 1 | mod(startSamp(1), 1) > 0 | ...
            mod(startSamp(2), 1) > 0 | isinf(startSamp(1)) | ...
            isinf(startSamp(2)) | isnan(startSamp(1)) | isnan(startSamp(2))
        % startSamp must be real for mod
        % startSamp(1) and startSamp(2) must be real postive integers
        startSamp
        error(['graven: startSamp(1) and startSamp(2) must be real '...
                'finite integers >= 1'])
    end
    
    % Check external angles for error
    checkexternal(external);
    
    % nSim must be a real positive interger
    if nSim < 1 | mod(nSim, 1) > 0 | isinf(nSim) | isnan(nSim)
        % nSim must be real for mod
        nSim
        error('graven: nSim must be a real finite integer >= 1')
    end
end

% Check variables that all modes use for errors

% logFileName must be a string
if ~strcmp(class(logFileName), 'char')
    logFileName
    error('graven: logFileName must be a string')
end

% gps must be a real positive integer
if gps < 0 | mod(gps, 1) > 0 | isinf(gps) | isnan(gps)
    % gps must be real for mod
    gps
    error('graven: gps must be a real positive finite integer')
end

% detId must be a string
if ~strcmp(class(detId), 'char')
    detId
    error('graven: detId must be a string')
end

% sampFreq must be a real positive number
if sampFreq <= 0 | isinf(sampFreq) | isnan(sampFreq) | ~isreal(sampFreq)
    sampFreq
    error('graven: sampFreq must be a real finite number > 0')
end

% seed must be a number
if ~strcmp(class(seed), 'double')
    seed
    error('graven: seed must be a number')
end

return

% ~~~ END ERRORCHECK ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

function checkexternal(external)

% external must be a vector, string or struct
if strcmp(class(external), 'double')
    % if external is a double, it must be a 1x3 vector
    if size(external, 1) ~= 1 | size(external, 2) ~= 3
        external
        error(['graven: external must be a 1x3 vector when not a ' ...
                'struct or string'])
    elseif abs(external(1)) > 1 | ~isreal(external(1))
        % external(1) must be real and bewteen -1 and 1
        external
        error('graven: external(1) must be real and bewteen -1 and 1')
    elseif external(2) < -pi | external(2) > pi | ~isreal(external(2))
        % external(2) must be real and between -pi and pi
        external
        error('graven: external(2) must be real and between -pi and pi')
    elseif external(3) < 0 | external(3) > 2*pi | ~isreal(external(3))
        % external(3) must be real and between 0 and 2*pi
        external
        error('graven: external(3) must be real and between 0 and 2*pi')
    end
elseif strcmp(class(external), 'char')
    % if external is a string, it must be either 'RANDOM' or 'OPTIMAL'
    if ~strcmp(upper(external), 'RANDOM') & ...
            ~strcmp(upper(external), 'OPTIMAL') 
        external
        error(['graven: external must be ''RANDOM'' or ''OPTIMAL'' ' ...
                'when a string'])
    end
elseif strcmp(class(external), 'struct')
    % if external is a struct is must contain fields x, phi and psi
    if ~isfield(external, 'x')
        external
        error('graven: external must contain field x when a struct')
    elseif abs(external.x) > 1 | ~isreal(external.x)
        % if external.x exists, it must be real and have a value between
        % -1 and 1
        external.x
        error('graven: external.x must be real and between -1 and 1')
    end
    if ~isfield(external, 'phi')
        external
        error('graven: external must contain field phi when a struct')
    elseif external.phi < -pi | external.phi > pi | ~isreal(external.phi)
        % if external.phi exits, it must be real and have a value between
        % -pi and pi
        external.phi
        error(['graven: external.phi must be real and betweeen -pi '...
                'and pi'])
    end
    if ~isfield(external, 'psi')
        external
        error('graven: external must contain field psi when a struct')
    elseif external.psi < 0 | external.psi > (2*pi) | ...
            ~isreal(external.psi)
        % if external.psi exists, it must be real and have a value 
        % between 0 and 2*pi
        external.psi
        error('graven: external.psi must be real and between 0 and 2*pi')
    end
end
return

% ~~~ END CHECKEXTERNAL ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

function htt = simidfile(simId)

% SIMIDFILE - retrieve pre-generated waveform when simId is a file name
%
% htt = simidfile(simId);
%
% simId     file name of the pre-generated waveform, (needs to be massaged
%           into useable form)
%
% htt       2xtime array containing + polarization waveform in the first row
%           and x polarization waveform in the second row (pre-generated 
%           waveforms are presummed to only contain + polarization)

% simId is the file containing the pre-generated waveform but its
% form may need to be massaged to be used properly:

% DETERMINE class of simId
% IF simId is a cell, 
%   SEND simId to loadpregendata
% ELSEIF simId is a string
%   PROCESS the simId string into a cell
%   SEND simId to loadpregendata
% ELSE
%   ERROR
% ENDIF
% RETURN htt

if strcmp(class(simId), 'cell')
    % simId does not need to be massaged
    htt = loadpregendata(simId);
elseif strcmp(class(simId), 'char')
    % simId needs to be massaged
    I = find(simId == ',');
    if length(I) == 0
        file = simId(1:end);
    elseif length(I) == 1
        file{1,1} = simId(1:I(1)-1);
        file{1,2} = simId(I(1)+1:end);
    else
        simId
        error('graven -> simidfile: Invalid simId')
    end
    % Load the pre-generated waveform
    htt = loadpregendata(file);
else
    simId
    error('graven -> simidfile: simId must be cell or string')
end

return

% ~~~ END SIMIDFILE ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

function htt = loadpregendata(inputFileName)

% LOADPREGENDATA - load pre-generated waveform data from file
%
% htt = loadpregendata(inputFileName);
%
% inputFileName     1x1 or 1x2 cell containing a string with name of input 
%                   file (must include file extention) and an optional
%                   field path for *.ilwd files 
%
% htt               2xtime array containing + polarization waveform in the
%                   first row and x polarization waveform in the second row 
%                   (pre-generated waveforms are presummed to only contain + 
%                   polarization)

% PARSE inputFileName
% IF ext is 'ilwd'
%   LOAD ilwd-file
%   GET field path
%   GET appropriate field
%   PUT waveform into htt
% ELSEIF ext is 'txt'
%   LOAD txt-file
%   PUT waveform into htt
% ELSE
%   ERROR
% ENDIF
% RETURN htt

[fileName, ext, dirPath] = fileext(inputFileName);
if isempty(fileName) | isempty(ext)
    error('graven: fileext returned empty file or ext')
end
% Assumed that the pre-generated waveform contains only + polarization
if strcmp(ext, 'ilwd')
    if exist([dirPath fileName '.' ext]) ~= 2
        file = [dirPath fileName '.' ext];
        error('graven: File does not exist on search path')
    end
    % load waveform
    waveform = ilwdread([dirPath fileName '.' ext]);
    % get field path
    if length(inputFileName) == 2
        ilwdString = inputFileName{2};
        if strcmp(lower(ilwdString), 'default');
            % Set default ilwdstring to match format on Burst Sim page:
            % http://www.ligo.caltech.edu/~ajw/bursts/burstsim.html
            ilwdString = 'm1.real_8';
        end
    else
        % Set default ilwdstring to match format on Burst Sim page:
        % http://www.ligo.caltech.edu/~ajw/bursts/burstsim.html
        ilwdString = 'm1.real_8';
    end
    % Since getfield only will handle one field at a time, the ilwdString
    % must be used one field at a time.  The following lines of code handle
    % this:
    % Parse the ilwd-file at the '.'
    I = find(ilwdString == '.');
    if ~isempty(I) & I(end) == length(ilwdString)
        ilwdString
        error('graven: Invalid ilwdString')
    end
    % Loop through each field using getfield
    nFields = length(I)+1;
    for iFields = 1:nFields
        if iFields == 1 & nFields ~= 1
            temp = getfield(waveform, ilwdString(1:I(iFields)-1));
        elseif iFields == 1 & nFields == 1
            temp = getfield(waveform, ilwdString);
        elseif iFields == nFields & nFields ~= 1
            temp = getfield(temp, ilwdString(I(iFields-1)+1:end));
        else
            temp = getfield(temp, ilwdString(I(iFields-1)+1:I(iFields)-1));
        end
    end
    if ~strcmp(class(temp), 'double')
        ilwdString
        error('graven: Retrieved data not class double')
    end
    htt(1,:) = temp(:)';
    % The following line will not compile but will work in the MATLAB
    % environment
    %   eval(['htt(1,:) = waveform.' ilwdString '(:)'';'])
    htt(2,:) = zeros(1,length(htt(1,:)));
elseif strcmp(ext, 'txt')
    if exist([dirPath fileName '.' ext]) ~= 2
        file = [dirPath fileName '.' ext];
        error('graven: File does not exist on search path')
    end
    % Load waveform
    waveform = load([dirPath fileName '.' ext]);
    htt(1,:) = waveform(:)';
    htt(2,:) = zeros(1, length(waveform));
    % The following code will not compile but can be used in the MATLAB
    % environment:
    % elseif strcmp(ext, 'mat') % Cannot contain struct, only double arrays
    %     if exist([dirPath fileName '.' ext]) ~= 2
    %         file = [dirPath fileName '.' ext];
    %         error('graven: File does not exist on search path')
    %     end
    %     eval(['load ' dirPath fileName '.' ext]) % inputFileName{1}
    %     if length(inputFileName) == 2
    %         eval(['waveform = ' inputFileName{2} ';'])
    %     else 
    %         % Set default to assume that the variable name matched the file
    %         % name
    %         eval(['waveform = ' fileName ';'])
    %     end
    %     htt(1,:) = waveform(:)';
    %     htt(2,:) = zeros(1, length(waveform));
else
    ext
    error('graven: inputFileName must have extention *.ilwd or *.txt')
end

return