graven

 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
                           (gps < 0 use |gps| but no 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']
                - + polarization only
            ['inputFileName,variableName']
                - + polarization only 
            ['inputFileName;inputFileName']
                - + and x polarization (in that order), ; separator
            ['nputFileName,variableName;inputFileName,variableName']
                - + and x polarization (in that order), ; separator
            {'inputFileName'}
                - + polarization only
            {'inputFileName','variableName'}
                - + polarization only
            { 'inputFileName','variableName','inputFileName','variableName'}
                - + and x polarization (in that order)
 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 assuming a sample
                   rate of 16384 sample/sec.  If detID = 'VIRGO' | 'TAMA', 
                   sampFreq will automatically be adjusted for a sampling
                   frequency of 20 kHz.
 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
  .startSamp   start sample at detector
  .h           strain signal to add starting at sample startSamp
  .AS_Q        AS_Q signal to add starting at sample startSamp
  .DARM_ERR    DARM_ERR signal to add starting at sample startSamp
  .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)

 For more complete documentation, please see LIGO-T040020.

 SEE ALSO: MAKEH, MAKETT, DETPROJ, CALIBSIMFD, IFODELAY, PREGENLIST

0001 function [Sim, state] = graven(varargin)
0002 
0003 % GRAVEN - simulate a set of GW signals for addition to a data stream
0004 %
0005 % [Sim, state] = graven(logFileName, gps, detId, inputFileName, ...
0006 %           sampFreq, seed);
0007 %   Generate signals based on information in file inputFileName
0008 % [Sim, state] = graven(nSim, logFileName, gps, detId, inputFileName, ...
0009 %           sampFreq, seed);
0010 %   Generate signals based on nSim random draws from information in file
0011 %   inputFileName
0012 % [Sim, state] = graven(nSim, logFileName, gps, detId, simId, ampl, ...
0013 %           startSamp, internal, external, sampFreq, seed);
0014 %   Generate nSim signals for fixed detId, simId, and ampl; random start in
0015 %   range startSamp, and either fixed or random internal and external
0016 % [Sim, state] = graven(nSim, logFileName, gps, detId, inputFileName, ...
0017 %           startSamp, external, sampFreq, seed);
0018 %   Push a pregenerated waveform contained in inputFileName name through the
0019 %   detector projection function and the inverse claibrator nSim times
0020 %
0021 % logFileName       string with name of log file (must include file
0022 %                   extention)
0023 % detId             detector identifier string; type 'help getifo' for list
0024 %                   of currently supported detectors
0025 % gps               GPS start time to allow selection of calibration
0026 %                           (gps < 0 use |gps| but no calibration)
0027 % inputFileName     string with name of input file (must include file
0028 %                   extention) or cell containing file name of
0029 %                   pre-generated waveform (with extention) and name of
0030 %                   variable or field path of timeseries
0031 %       e.g. ['inputFileName']
0032 %                - + polarization only
0033 %            ['inputFileName,variableName']
0034 %                - + polarization only
0035 %            ['inputFileName;inputFileName']
0036 %                - + and x polarization (in that order), ; separator
0037 %            ['nputFileName,variableName;inputFileName,variableName']
0038 %                - + and x polarization (in that order), ; separator
0039 %            {'inputFileName'}
0040 %                - + polarization only
0041 %            {'inputFileName','variableName'}
0042 %                - + polarization only
0043 %            { 'inputFileName','variableName','inputFileName','variableName'}
0044 %                - + and x polarization (in that order)
0045 % nSim              positive integer indicating how many simulations to
0046 %                   generate
0047 % simId             signal type string; type 'help makeh' for list of
0048 %                   currently supported signals
0049 % ampl              signal amplitude or amplitude range
0050 %                   *When ampl is a scalar, the amplitude is fixed
0051 %                   *When ampl is a 1x2 vector, the amplitude is randomly
0052 %                   chosen between ampl(1) and ampl(2) in a logarithmic
0053 %                   distribution. N.B. ampl(1)=ampl(2) is equivilent to
0054 %                   entering a fixed amplitude
0055 %                   *When ampl is a 1x3 vector, the amplitdue is randomly
0056 %                   chosen from the discrete linear range
0057 %                   ampl(1):ampl(2):ampl(3)
0058 % startSamp         1x2 vector defining range of allowed start times in
0059 %                   samples wrt the beginning of data set assuming a sample
0060 %                   rate of 16384 sample/sec.  If detID = 'VIRGO' | 'TAMA',
0061 %                   sampFreq will automatically be adjusted for a sampling
0062 %                   frequency of 20 kHz.
0063 % internal          struct or vector with internal orientation angles, or
0064 %                   string 'RANDOM' or 'OPTIMAL'; if struct, must have
0065 %                   fields x and phi
0066 % external          struct or vector with sky location & polarization angle
0067 %                   or string 'RANDOM' or 'OPTIMAL'; if struct, must have
0068 %                   fields x, phi and psi
0069 % sampFreq          sampling frequency of detId
0070 % seed              random number generator seed (state can be fed back in
0071 %                   as seed)
0072 %
0073 % Sim(nSim)    struct array of simuluations
0074 %  .name        simId string
0075 %  .startSamp   start sample at detector
0076 %  .h           strain signal to add starting at sample startSamp
0077 %  .AS_Q        AS_Q signal to add starting at sample startSamp
0078 %  .DARM_ERR    DARM_ERR signal to add starting at sample startSamp
0079 %  .gps0        signal start time in gps sec, ns at detector
0080 %  .pkgps       signal peak time in gps sec, ns at detector
0081 % state         state of the random number generator at the end of the code
0082 %               (can be fed back into graven as seed)
0083 %
0084 % For more complete documentation, please see LIGO-T040020.
0085 %
0086 % SEE ALSO: MAKEH, MAKETT, DETPROJ, CALIBSIMFD, IFODELAY, PREGENLIST
0087 
0088 % $Id: graven.m,v 1.52 2005/04/28 19:55:02 stuver Exp $
0089 
0090 % DETERMINE number of input arguments
0091 % CHECK variables for error
0092 % DETERMINE mode of operation from number of inputs
0093 % SET preGen flag to false
0094 % IF mode 1 (6 input arguments)
0095 %   ASSIGN variable names to input argumnets
0096 %   SET random number generator seed
0097 %   MAKE list reading line-by-line from a file
0098 %   RETURN list
0099 % ELSEIF mode 2 (7 input argumnets)
0100 %   ASSIGN variable names to input arguments
0101 %   CHECK variables for error
0102 %   SET random number generator seed
0103 %   MAKE list by drawing on file nSim times
0104 %   RETURN list
0105 % ELSEIF mode 3 (11 input arguments)
0106 %   ASSIGN variable names to input arguments
0107 %   SET random number generator seed
0108 %   MAKE list based directly on input arguments
0109 %   RETURN list
0110 % ELSEIF mode 4 (9 input arguments)
0111 %   ASSIGN variable names to input arguments
0112 %   SET random number generator seed
0113 %   MAKE list for pre-generated waveforms
0114 % ENDIF
0115 
0116 LIGOGEO = 16384; % Sample rate for LIGO and GEO600
0117 VIRGOTAMA = 2e+4; % Sample rate for VIRGO and TAMA
0118 
0119 % Read length of varargin
0120 inLength = length(varargin);
0121 
0122 % Check for errors
0123 input = varargin;
0124 [gps,skipCalibFlag] = errorcheck(input);
0125 
0126 preGen = ~true;
0127 
0128 % Assign variable names to components of varargin & generate list
0129 if inLength == 6
0130     [logFileName, gps, detId, inputFileName, sampFreq, seed] = ...
0131         deal(varargin{:});
0132     % Seed rand here
0133     rand('state', seed);
0134     % Generate simulation list
0135     list = makelist(inputFileName, detId);
0136     % list is a cell array
0137 elseif inLength == 7
0138     [nSim, logFileName, gps, detId, inputFileName, sampFreq, seed] = ...
0139         deal(varargin{:});
0140     % Seed rand here
0141     rand('state', seed);
0142     % Generate simulation list
0143     list = makelist(nSim, inputFileName, detId);
0144     % list is a cell array
0145 elseif inLength == 11
0146     [nSim, logFileName, gps, detId, simId, ampl, startSamp, internal, ...
0147             external, sampFreq, seed] = deal(varargin{:});
0148     % Seed rand here
0149     rand('state', seed);
0150     % Generate simulation list
0151     list = makelist(nSim, simId, ampl, startSamp, internal, external, detId);
0152     % list is a cell array
0153 elseif inLength == 9
0154     [nSim, logFileName, gps, detId, inputFileName, startSamp, external, ...
0155             sampFreq, seed] = deal(varargin{:});    
0156     % Seed rand here
0157     rand('state', seed);
0158     list = pregenlist(nSim, inputFileName, startSamp, external, detId);
0159     % list is a cell array
0160 else % Error
0161     error('graven: Number of input arguments can only be 6, 7, 9 or 11')
0162 end
0163 
0164 % IF detector is not part of LIGO
0165 %   SET flag to produce output in strain (for now)
0166 % ENDIF
0167 oldSkipFlag = skipCalibFlag;
0168 testDet = upper(detId);
0169 switch testDet
0170     case 'H1'
0171         skipCalibFlag = oldSkipFlag;
0172     case 'H2'
0173         skipCalibFlag = oldSkipFlag;
0174     case 'L1'
0175         skipCalibFlag = oldSkipFlag;
0176     otherwise
0177         skipCalibFlag = true;
0178 end
0179 
0180 % DETERMINE number of lines in list
0181 % INITIALIZE structures
0182 % OPEN log file
0183 % WRITE headerlines to log file
0184 
0185 % Read length of List (each row is a simulation)
0186 listLength = size(list, 1);
0187 
0188 if listLength == 0
0189     error('graven: listing error, list has zero length')
0190 end
0191 
0192 % Initialize Sim, Internal and External structures
0193 [Sim(1:listLength)] = deal(struct('name', [], 'startSamp', [], 'h', [], ...
0194     'AS_Q', [], 'DARM_ERR', [], 'gps0', [], 'pkgps' ,[]));
0195 Internal = deal(struct('x', [], 'phi', []));
0196 External = deal(struct('x', [], 'phi', [], 'psi', []));
0197 
0198 % Open log file w/ append option to prevent over-writing
0199 fid = fopen(logFileName, 'a');
0200 if fid == -1
0201     error('graven: Failed to open log file')
0202 end
0203 
0204 date = clock;
0205 % Write headerlines to log file
0206 fprintf(fid, ['# This simulation was run on %2.0f/%2.0f/%4.0f ' ...
0207         'at %2.0f:%2.0f:%2.2f EDT\n'], date(2), date(3), date(1), ...
0208     date(4), date(5), date(6));
0209 fprintf(fid, ['# (simId) (ampl) (startSamp1) (startSamp2) (Internal.x) '...
0210         '(Internal.phi) (External.x) (External.phi) (External.psi)' ...
0211         ' (gps) (gps0 s) (gps0 ns) (pkampl)  (pkgps s) (pkgps ns) '...
0212         '(detId)\n']);
0213 
0214 % INITIALIZE FOR loop to run through the list
0215 %   ASSIGN variable names for this simulation
0216 %   DETERMINE number of samples to add to startSamp to account for the GW's
0217 %             time of flight from the center of the Earth then add to
0218 %             startSamp
0219 %   IF preGen flag is false
0220 %       GENERATE metric perturbation
0221 %       PROJECT metric perturbation into TT gauge
0222 %   ENDIF
0223 %   PROJECT TT gauge perturbation onto antenna pattern
0224 %   DETERMINE the maximum strain and time of max strain
0225 %   CALIBRATE strain into AS_Q counts
0226 %   SAVE simulation into Sims struct
0227 %   WRITE simulation information to log file
0228 % ENDFOR
0229 % CLOSE log file
0230 % SAVE random number generator state
0231 
0232 RAD2DEG = 180/pi;
0233 
0234 for iSim = 1:listLength
0235     % Read list:
0236     simId = list{iSim, 1};
0237     % Extra white space must be removed from simId
0238     space = find(simId == ' ');
0239     simId(space) = [];
0240     ampl = list{iSim, 2};
0241     startSamp = list{iSim, 3};
0242     % Adjust startSamp to appropriate sampling frequency for TAMA & VIRGO
0243     if sum(strcmp(upper(detId),{'VIRGO', 'TAMA'})) > 0
0244        startSamp = ((startSamp-1)*(VIRGOTAMA/LIGOGEO))+1;
0245     end
0246     Internal.x = list{iSim, 4};
0247     Internal.phi = list{iSim, 5};
0248     External.x = list{iSim, 6};
0249     External.phi = list{iSim, 7};
0250     External.psi = list{iSim, 8};
0251     
0252     % If ampl = -1, this simulation is from a pre-generated waveform
0253     if ampl == -1
0254         % Reset preGen
0255         preGen = true;
0256     end    
0257     
0258     % startSamp is the time in samples the GW is incident on the center of
0259     % the Earth.  Calculate the number of samples to ADD to this startSamp to
0260     % be incident on detId
0261     samp = ifodelay(detId, External.x, External.phi, sampFreq);
0262     if isempty(samp)
0263         error('graven: ifodelay returned empty samp')
0264     end
0265     stimeIfo = startSamp+samp;
0266     
0267     % If this is a pre-generated waveform, the start time of the simulation
0268     % at the detector is rounded to the nearest integer sample since
0269     % inter-sample start times are not currently supported for pre-generated
0270     % waveforms.
0271     if preGen
0272         stimeIfo = round(stimeIfo);
0273         htt = simidfile(simId);
0274     end
0275     
0276     % Assume pre-generated waveforms are already in the TT gauge
0277     if ~preGen
0278         id = simparse(simId);
0279         if ~strcmp(upper(id), 'CHIRPLET')
0280             % Generate metric perturbation
0281             hij = makeh(simId, ampl, stimeIfo, sampFreq);
0282             if isempty(hij)
0283                 error('graven: makeh returned empty hij')
0284             end
0285             % Project perturbation into TT gauge
0286             htt = makett(hij, Internal);
0287             if isempty(htt)
0288                 error('graven: makett returned empty htt')
0289             end
0290         elseif strcmp(upper(id), 'CHIRPLET')
0291             % Make the htt for the chirplet waveform
0292             htt = makechirplet(simId, ampl, startSamp, sampFreq)
0293         end
0294     end
0295 
0296     % Project TT perturbation onto the detector
0297     h = detproj(htt, detId, External);
0298     if isempty(h)
0299         error('graven: detproj returned empty h')
0300     end
0301     
0302     % Find the magnitude and time of the maximum strain
0303     [maxStrain, maxSamp] = max(h);
0304     maxTimeNdx = ceil(stimeIfo)+maxSamp-1; % Time of max strain in samples
0305     % Convert into sec and ns
0306     maxTimeGps = ndx2gps(maxTimeNdx, sampFreq, gps); 
0307     
0308     
0309     % IF not calibrating
0310     %    SET AS_Q to raw data
0311     % ELSE
0312     %    Convert to AS_Q counts
0313     % ENDIF
0314     if(skipCalibFlag == true)
0315         asqSeries = zeros(size(h'));
0316         darmSeries = asqSeries;
0317     else
0318         [asqSeries] = calibsimfd(detId, fix(gps+(stimeIfo/sampFreq)), h, ...
0319         sampFreq, 'AS_Q');
0320         if isempty(asqSeries)
0321             error('graven: calibsimfd returned empty asqSeries')
0322         end
0323         runName = getcalibrun(gps);
0324         if strcmp(upper(runName), 'S4')
0325             [darmSeries] = calibsimfd(detId, fix(gps+(stimeIfo/sampFreq)), h, ...
0326             sampFreq, 'DARM_ERR');
0327             if isempty(darmSeries)
0328                 error('graven: calibsimfd returned empty darmSeries')
0329             end
0330         else  % Return zeroes for the DARM_ERR time series if not S4
0331             darmSeries = zeros(size(asqSeries));
0332         end
0333     end
0334     % Save start time, time-series, GPS time and max strain time for sim
0335     Sim(iSim).name = upper(simId);
0336     Sim(iSim).startSamp = ceil(stimeIfo);
0337     Sim(iSim).h = h';
0338     Sim(iSim).AS_Q = asqSeries;
0339     Sim(iSim).DARM_ERR = darmSeries;
0340     Sim(iSim).gps0 = ndx2gps(stimeIfo, sampFreq, gps);
0341     Sim(iSim).pkgps = maxTimeGps;
0342     
0343     % Write information for this sim to log file
0344     fprintf(fid, ['%s %0.5g %9.15f %9.15f %d %d %d %d %d %d %d ' ...
0345             '%8.0f %d %d %8.0f %s \n'], simId, ampl, startSamp, startSamp, ...
0346         Internal.x, Internal.phi, External.x, External.phi, ...
0347         External.psi, gps, Sim(iSim).gps0.sec, Sim(iSim).gps0.ns, ...
0348         maxStrain, Sim(iSim).pkgps.sec, Sim(iSim).pkgps.ns, detId);
0349 end
0350 
0351 % Close logfile
0352 fclose(fid);
0353 
0354 % Save random number generator state
0355 state = rand('state');
0356 
0357 return
0358 
0359 % ~~~ END MAIN ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0360 
0361 function [gpsVal,skipCalibFlag] = errorcheck(input)
0362 
0363 % ERRORCHECK - checks graven's input data for errors and stops the code if
0364 %              an error is found
0365 %
0366 % [gpsVal,skipCalibFlag] = errorcheck(input)
0367 %
0368 % input     the varargin input argument to graven
0369 %
0370 %  gpsVal GPS value (positive integer)
0371 %   skipCalibFlag (true = do not calibrate sim)
0372 %    - set if input GPS < 0
0373 
0374 % DETERMINE number of input arguments
0375 % DETERMINE mode by number of inputs
0376 % IF mode 1
0377 %   ASSIGN variables
0378 %   CHECK variables that are not used in all modes for error
0379 % ELSEIF mode 2
0380 %   ASSIGN variables
0381 %   CHECK variables that are not used in all modes for error
0382 % ELSEIF mode 3
0383 %   ASSIGN variables
0384 %   CHECK variables that are not used in all modes for error
0385 % ELSEIF mode 4
0386 %   ASSIGN variables
0387 %   CHECK variables that are not used in all modes for error
0388 % ENDIF
0389 % CHECK variables that are used by all modes for error
0390 gpsVal = 0;
0391 skipCalibFlag = false;
0392 
0393 inputLength=length(input);
0394 
0395 if inputLength == 6 % mode 1
0396     [logFileName, gps, detId, inputFileName, sampFreq, seed] = ...
0397         deal(input{:});
0398     
0399     % Check inputFileName for errors since not all of the modes require
0400     % inputFileName
0401     
0402     % inputFileName must be a string
0403     if ~strcmp(class(inputFileName), 'char')
0404         inputFileName
0405         error('graven: inputFileName must be a string')
0406     end
0407     
0408 elseif inputLength == 7 % mode 2
0409     [nSim, logFileName, gps, detId, inputFileName, sampFreq, seed] = ...
0410         deal(input{:});
0411     
0412     % Check inputFileName and nSim for errors since not all modes
0413     % require them
0414     
0415     % inputFileName must be a string
0416     if ~strcmp(class(inputFileName), 'char')
0417         inputFileName
0418         error('graven: inputFileName must be a string')
0419     end
0420     
0421     % nSim must be a real positive interger
0422     if nSim < 1 | mod(nSim, 1) > 0 | isinf(nSim) | isnan(nSim)
0423         % nSim must be real for mod
0424         nSim
0425         error('graven: nSim must be a real finite integer >= 1')
0426     end
0427     
0428 elseif inputLength == 11
0429     [nSim, logFileName, gps, detId, simId, ampl, startSamp, internal, ...
0430             external, sampFreq, seed] = deal(input{:});
0431     
0432     % Check nSim and variables that only mode 3 uses (simId, etc) for errors
0433     
0434     % nSim must be a real postive integer
0435     if nSim < 1 | mod(nSim, 1) > 0 | isinf(nSim) | isnan(nSim)
0436         % nSim must be real for mod
0437         nSim
0438         error('graven: nSim must be a real finite integer >= 1')
0439     end
0440     
0441     % simId must be a string
0442     if ~strcmp(class(simId), 'char')
0443         simId
0444         error('graven: simId must be a string')
0445     end
0446     
0447     % ampl must be a scalar, 1x2 or 1x3 vector with real positive numbers
0448     if size(ampl, 1) == 1 & size(ampl, 2) == 1 % scalar
0449         if ampl < 0 | isinf(ampl) | isnan(ampl) | ~isreal(ampl)
0450             % ampl must be a real postive number
0451             ampl
0452             error('graven: ampl must be a real positive fitite number')
0453         end
0454     elseif size(ampl, 1) == 1 & size(ampl, 2) == 2 % 1x2 vector
0455         if ampl < 0 | isinf(ampl(1)) | isinf(ampl(2)) | isnan(ampl(1)) ...
0456                 | isnan(ampl(2)) | ~isreal(ampl(1)) | ~isreal(ampl(2))
0457             % All components of the ampl vector must be real positive numbers
0458             ampl
0459             error(['graven: ampl components must be real positive finite '...
0460                     'numbers'])
0461         end
0462     elseif size(ampl, 1) == 1 & size(ampl, 2) == 3 % 1x3 vector
0463         if ampl < 0 | isinf(ampl(1)) | isinf(ampl(2)) | isinf(ampl(3)) | ...
0464                 isnan(ampl(1)) | isnan(ampl(2)) | isnan(ampl(3)) | ...
0465                 ~isreal(ampl(1)) | ~isreal(ampl(2)) | ~isreal(ampl(3))
0466             % All components of the ampl vector must be real positive numbers
0467             ampl
0468             error(['graven: ampl components must be real positive finite '...
0469                     'numbers'])
0470         end
0471     else
0472         error('graven: ampl must be a scalar or a 1x2 or 1x3 vector')
0473     end
0474     
0475     % startSamp must be a 1x2 vector and contain real postive integers
0476     if size(startSamp, 1) ~= 1 | size(startSamp, 2) ~= 2
0477         startSamp
0478         error('graven: startSamp must be a 1x2 vector')
0479     elseif startSamp(1) < 1 | mod(startSamp(1), 1) > 0 | ...
0480             mod(startSamp(2), 1) > 0 | isinf(startSamp(1)) | ...
0481             isinf(startSamp(2)) | isnan(startSamp(1)) | isnan(startSamp(2))
0482         % startSamp must be real for mod
0483         % startSamp(1) and startSamp(2) must be real postive integers
0484         startSamp
0485         error(['graven: startSamp(1) and startSamp(2) must be real '...
0486                 'finite integers >= 1'])
0487     end
0488     
0489     % internal must be a vector, string or struct
0490     if strcmp(class(internal), 'double')
0491         % if internal is a double, it must be a 1x2 vector
0492         if size(internal, 1) ~= 1 | size(internal, 2) ~= 2
0493             internal
0494             error(['graven: internal must be a 1x2 vector when not a ' ...
0495                     'struct or string'])
0496         elseif abs(internal(1)) > 1 | ~isreal(internal(1))
0497             % internal(1) must be real and bewteen -1 and 1
0498             internal
0499             error('graven: internal(1) must be real and bewteen -1 and 1')
0500         elseif internal(2) < -pi | internal(2) > pi | ~isreal(internal(2))
0501             % internal(2) must be real and between -pi and pi
0502             internal
0503             error('graven: internal(2) must be real and between -pi and pi')
0504         end
0505     elseif strcmp(class(internal), 'char')
0506         % if internal is a string, it must be either 'RANDOM' or 'OPTIMAL'
0507         if ~strcmp(upper(internal), 'RANDOM') & ...
0508                 ~strcmp(upper(internal), 'OPTIMAL') 
0509             internal
0510             error(['graven: internal must be ''RANDOM'' or ''OPTIMAL'' ' ...
0511                     'when a string'])
0512         end
0513     elseif strcmp(class(internal), 'struct')
0514         % if internal is a struct is must contain fields x and phi
0515         if ~isfield(internal, 'x')
0516             internal
0517             error('graven: internal must contain field x when a struct')
0518         elseif abs(internal.x) > 1 | ~isreal(internal.x)
0519             % if internal.x exists, it must have a value between -1 and 1
0520             internal.x
0521             error('graven: internal.x must be real and between -1 and 1')
0522         end
0523         if ~isfield(internal, 'phi')
0524             internal
0525             error('graven: internal must contain field phi when a struct')
0526         elseif internal.phi < -pi | internal.phi > pi | ~isreal(internal.phi)
0527             % if internal.phi exits, it must be real and have a value between
0528             % -pi and pi
0529             internal.phi
0530             error(['graven: internal.phi must be real and betweeen -pi '...
0531                     'and pi'])
0532         end
0533     end
0534     
0535     % Check external angles for error
0536     checkexternal(external);
0537     
0538 elseif inputLength == 9
0539     [nSim, logFileName, gps, detId, inputFileName, startSamp, external, ...
0540             sampFreq, seed] = deal(input{:});
0541     
0542     % inputFileName must be a string
0543     if ~strcmp(class(inputFileName), 'char') & ~iscell(inputFileName)
0544         inputFileName
0545         error('graven: inputFileName must be a string')
0546     end
0547     
0548     % startSamp must be a 1x2 vector and contain real postive integers
0549     if size(startSamp, 1) ~= 1 | size(startSamp, 2) ~= 2
0550         startSamp
0551         error('graven: startSamp must be a 1x2 vector')
0552     elseif startSamp(1) < 1 | mod(startSamp(1), 1) > 0 | ...
0553             mod(startSamp(2), 1) > 0 | isinf(startSamp(1)) | ...
0554             isinf(startSamp(2)) | isnan(startSamp(1)) | isnan(startSamp(2))
0555         % startSamp must be real for mod
0556         % startSamp(1) and startSamp(2) must be real postive integers
0557         startSamp
0558         error(['graven: startSamp(1) and startSamp(2) must be real '...
0559                 'finite integers >= 1'])
0560     end
0561     
0562     % Check external angles for error
0563     checkexternal(external);
0564     
0565     % nSim must be a real positive interger
0566     if nSim < 1 | mod(nSim, 1) > 0 | isinf(nSim) | isnan(nSim)
0567         % nSim must be real for mod
0568         nSim
0569         error('graven: nSim must be a real finite integer >= 1')
0570     end
0571 end
0572 
0573 % Check variables that all modes use for errors
0574 
0575 % logFileName must be a string
0576 if ~strcmp(class(logFileName), 'char')
0577     logFileName
0578     error('graven: logFileName must be a string')
0579 end
0580 
0581 % gps must be a real non-zero integer
0582 if gps == 0 | mod(gps, 1) > 0 | isinf(gps) | isnan(gps)
0583     % gps must be real for mod
0584     error('graven: gps must be a real finite integer')
0585 end
0586 % IF GPS < 0
0587 %   SET skip calibration flag
0588 %   SET output gps to absolute value
0589 % ELSE
0590 %   CLEAR skip calibration flag
0591 %   SET output gps to current value
0592 % ENDIF
0593 if (gps < 0)
0594     skipCalibFlag = true;
0595     gpsVal = abs(gps);
0596 else
0597     skipCalibFlag = false;
0598     gpsVal = gps;
0599 end
0600 
0601 % detId must be a string
0602 if ~strcmp(class(detId), 'char')
0603     detId
0604     error('graven: detId must be a string')
0605 end
0606 
0607 % sampFreq must be a real positive number
0608 if sampFreq <= 0 | isinf(sampFreq) | isnan(sampFreq) | ~isreal(sampFreq)
0609     sampFreq
0610     error('graven: sampFreq must be a real finite number > 0')
0611 end
0612 
0613 if sampFreq ~= 16384 & (strcmp(upper(detId),'H1') || strcmp(upper(detId),'H2') ...
0614         || strcmp(upper(detId),'L1') || strcmp(upper(detId),'GEO'))
0615     warning(['graven: sampFreq does not match the sampling frequency ' ...
0616         'of the detId entered'])
0617 %     fprintf(fid, ['sampFreq does not match the sampling frequency ' ...
0618 %         'of the detId entered'])
0619 elseif sampFreq ~= 20000 & (strcmp(upper(detId),'VIRGO') ...
0620         || strcmp(upper(detId),'TAMA'))
0621     warning(['graven: sampFreq does not match the sampling frequency ' ...
0622         'of the detId entered'])
0623 %     fprintf(fid, ['sampFreq does not match the sampling frequency ' ...
0624 %         'of the detId entered'])
0625 end
0626 
0627 % seed must be a number
0628 if ~strcmp(class(seed), 'double')
0629     seed
0630     error('graven: seed must be a number')
0631 end
0632 
0633 return
0634 
0635 % ~~~ END ERRORCHECK ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0636 
0637 function checkexternal(external)
0638 
0639 % external must be a vector, string or struct
0640 if strcmp(class(external), 'double')
0641     % if external is a double, it must be a 1x3 vector
0642     if size(external, 1) ~= 1 | size(external, 2) ~= 3
0643         external
0644         error(['graven: external must be a 1x3 vector when not a ' ...
0645                 'struct or string'])
0646     elseif abs(external(1)) > 1 | ~isreal(external(1))
0647         % external(1) must be real and bewteen -1 and 1
0648         external
0649         error('graven: external(1) must be real and bewteen -1 and 1')
0650     elseif external(2) < -pi | external(2) > pi | ~isreal(external(2))
0651         % external(2) must be real and between -pi and pi
0652         external
0653         error('graven: external(2) must be real and between -pi and pi')
0654     elseif external(3) < 0 | external(3) > 2*pi | ~isreal(external(3))
0655         % external(3) must be real and between 0 and 2*pi
0656         external
0657         error('graven: external(3) must be real and between 0 and 2*pi')
0658     end
0659 elseif strcmp(class(external), 'char')
0660     % if external is a string, it must be either 'RANDOM' or 'OPTIMAL'
0661     if ~strcmp(upper(external), 'RANDOM') & ...
0662             ~strcmp(upper(external), 'OPTIMAL') 
0663         external
0664         error(['graven: external must be ''RANDOM'' or ''OPTIMAL'' ' ...
0665                 'when a string'])
0666     end
0667 elseif strcmp(class(external), 'struct')
0668     % if external is a struct is must contain fields x, phi and psi
0669     if ~isfield(external, 'x')
0670         external
0671         error('graven: external must contain field x when a struct')
0672     elseif abs(external.x) > 1 | ~isreal(external.x)
0673         % if external.x exists, it must be real and have a value between
0674         % -1 and 1
0675         external.x
0676         error('graven: external.x must be real and between -1 and 1')
0677     end
0678     if ~isfield(external, 'phi')
0679         external
0680         error('graven: external must contain field phi when a struct')
0681     elseif external.phi < -pi | external.phi > pi | ~isreal(external.phi)
0682         % if external.phi exits, it must be real and have a value between
0683         % -pi and pi
0684         external.phi
0685         error(['graven: external.phi must be real and betweeen -pi '...
0686                 'and pi'])
0687     end
0688     if ~isfield(external, 'psi')
0689         external
0690         error('graven: external must contain field psi when a struct')
0691     elseif external.psi < 0 | external.psi > (2*pi) | ...
0692             ~isreal(external.psi)
0693         % if external.psi exists, it must be real and have a value
0694         % between 0 and 2*pi
0695         external.psi
0696         error('graven: external.psi must be real and between 0 and 2*pi')
0697     end
0698 end
0699 return
0700 
0701 % ~~~ END CHECKEXTERNAL ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0702 
0703 function htt = simidfile(simId)
0704 
0705 % SIMIDFILE - retrieve pre-generated waveform when simId is a file name
0706 %
0707 % htt = simidfile(simId);
0708 %
0709 % simId     file name of the pre-generated waveform, (needs to be massaged
0710 %           into useable form)
0711 %
0712 % htt       2xtime array containing + polarization waveform in the first row
0713 %           and x polarization waveform in the second row (pre-generated
0714 %           waveforms are presummed to only contain + polarization)
0715 
0716 % simId is the file containing the pre-generated waveform but its
0717 % form may need to be massaged to be used properly:
0718 
0719 % DETERMINE class of simId
0720 % IF simId is a cell,
0721 %   SEND simId to loadpregendata
0722 % ELSEIF simId is a string
0723 %   PROCESS the simId string into a cell
0724 %   SEND simId to loadpregendata
0725 % ELSE
0726 %   ERROR
0727 % ENDIF
0728 % RETURN htt
0729 
0730 if strcmp(class(simId), 'cell')
0731     % simId does not need to be massaged
0732     if length(simId) == 1 | length(simId) == 2
0733         htt = loadpregendata(simId);
0734     elseif length(simId) == 4
0735         plusPolFile = {simId{1}, simId{2}};
0736         hPlus = loadpregendata(plusPolFile);
0737         crossPolFile = {simId{3}, simId{4}};
0738         hcross = loadpregendata(crossPolFile);
0739         htt = [hPlus(1,:), hCross(1,:)];
0740     else
0741         simId
0742         error(['graven -> simidfile: if simId is a cell, it must have ' ...
0743             'length of 1, 2 or 4'])
0744     end
0745 elseif strcmp(class(simId), 'char')
0746     % simId needs to be massaged
0747     % ; is the separator between polarization files
0748     % + polarization is assumed to be listed first, then x
0749     G = find(simId == ';');
0750     if isempty(G)
0751         plusPolFile = simId(1:end);
0752     elseif length(G) == 1
0753         plusPolFile = simId(1:G-1);
0754         crossPolFile = simId(G+1:end);
0755     else
0756         simId
0757         error('graven -> simidfile: more than 2 pregenerated polarizations specified')
0758     end
0759     Iplus = find(plusPolFile == ',');
0760     if length(Iplus) == 0
0761         plusFile = plusPolFile(1:end);
0762     elseif length(Iplus) == 1
0763         plusFile{1,1} = plusPolFile(1:Iplus(1)-1);
0764         plusFile{1,2} = plusPolFile(Iplus(1)+1:end);
0765     else
0766         plusPolFile
0767         error('graven -> simidfile: Invalid plusPolFile')
0768     end
0769     if exist('crossPolFile') == 1
0770         Icross = find(crossPolFile == ',');
0771         if length(Icross) == 0
0772             crossFile = crossPolFile(1:end);
0773         elseif length(Icross) == 1
0774             crossFile{1,1} = crossPolFile(1:Icross(1)-1);
0775             crossFile{1,2} = crossPolFile(Icross(1)+1:end);
0776         else
0777             crossPolFile
0778             error('graven -> simidfile: Invalid crossPolFile')
0779         end
0780         hCross = loadpregendata(crossFile);
0781         % Load the pre-generated waveform
0782         hPlus = loadpregendata(plusFile);
0783         htt = [hPlus(1,:); hCross(1,:)];
0784     else
0785         htt = loadpregendata(plusFile);
0786     end
0787 else
0788     simId
0789     error('graven -> simidfile: simId must be cell or string')
0790 end
0791 
0792 return
0793 
0794 % ~~~ END SIMIDFILE ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
0795 
0796 function htt = loadpregendata(inputFileName)
0797 
0798 % LOADPREGENDATA - load pre-generated waveform data from file
0799 %
0800 % htt = loadpregendata(inputFileName);
0801 %
0802 % inputFileName     1x1 or 1x2 cell containing a string with name of input
0803 %                   file (must include file extention) and an optional
0804 %                   field path for *.ilwd files
0805 %
0806 % htt               2xtime array containing + polarization waveform in the
0807 %                   first row and x polarization waveform in the second row
0808 %                   (pre-generated waveforms are presummed to only contain +
0809 %                   polarization)
0810 
0811 % PARSE inputFileName
0812 % IF ext is 'ilwd'
0813 %   LOAD ilwd-file
0814 %   GET field path
0815 %   GET appropriate field
0816 %   PUT waveform into htt
0817 % ELSEIF ext is 'txt'
0818 %   LOAD txt-file
0819 %   PUT waveform into htt
0820 % ELSE
0821 %   ERROR
0822 % ENDIF
0823 % RETURN htt
0824 
0825 if ischar(inputFileName)
0826     inputFileName = {inputFileName};
0827 end
0828 [dirPath, fileName, ext] = fileparts(inputFileName{1});%fileext(inputFileName);
0829 dirPath = [dirPath filesep];
0830 if isempty(fileName) | isempty(ext)
0831     error('graven: fileext returned empty file or ext')
0832 end
0833 % Assumed that the pre-generated waveform contains only + polarization
0834 if strcmp(ext, '.ilwd')
0835     if exist([dirPath fileName ext]) ~= 2
0836         file = [dirPath fileName ext];
0837         error('graven: File does not exist on search path')
0838     end
0839     % load waveform
0840     waveform = ilwdread([dirPath fileName ext]);
0841     % get field path
0842     if length(inputFileName) == 2
0843         ilwdString = inputFileName{2};
0844         if strcmpi(ilwdString, 'default');
0845             % Set default ilwdstring to match format on Burst Sim page:
0846             % http://www.ligo.caltech.edu/~ajw/bursts/burstsim.html
0847             ilwdString = 'm1.real_8';
0848         end
0849     else
0850         % Set default ilwdstring to match format on Burst Sim page:
0851         % http://www.ligo.caltech.edu/~ajw/bursts/burstsim.html
0852         ilwdString = 'm1.real_8';
0853     end
0854     % Since getfield only will handle one field at a time, the ilwdString
0855     % must be used one field at a time.  The following lines of code handle
0856     % this:
0857     % Parse the ilwd-file at the '.'
0858     I = find(ilwdString == '.');
0859     if ~isempty(I) & I(end) == length(ilwdString)
0860         ilwdString
0861         error('graven: Invalid ilwdString')
0862     end
0863     % Loop through each field using getfield
0864     nFields = length(I)+1;
0865     for iFields = 1:nFields
0866         if iFields == 1 & nFields ~= 1
0867             temp = getfield(waveform, ilwdString(1:I(iFields)-1));
0868         elseif iFields == 1 & nFields == 1
0869             temp = getfield(waveform, ilwdString);
0870         elseif iFields == nFields & nFields ~= 1
0871             temp = getfield(temp, ilwdString(I(iFields-1)+1:end));
0872         else
0873             temp = getfield(temp, ilwdString(I(iFields-1)+1:I(iFields)-1));
0874         end
0875     end
0876     if ~strcmp(class(temp), 'double')
0877         ilwdString
0878         error('graven: Retrieved data not class double')
0879     end
0880     htt(1,:) = temp(:)';
0881     % The following line will not compile but will work in the MATLAB
0882     % environment
0883     %   eval(['htt(1,:) = waveform.' ilwdString '(:)'';'])
0884     htt(2,:) = zeros(1,length(htt(1,:)));
0885 elseif strcmp(ext, '.txt')
0886     if exist([dirPath fileName ext]) ~= 2;
0887         file = [dirPath fileName ext];
0888         error('graven: File does not exist on search path')
0889     end
0890     % Load waveform
0891     waveform = load([dirPath fileName ext]);
0892     htt(1,:) = waveform(:)';
0893     htt(2,:) = zeros(1, length(waveform));
0894     % The following code will not compile but can be used in the MATLAB
0895     % environment:
0896     % elseif strcmp(ext, 'mat') % Cannot contain struct, only double arrays
0897     %     if exist([dirPath fileName '.' ext]) ~= 2
0898     %         file = [dirPath fileName '.' ext];
0899     %         error('graven: File does not exist on search path')
0900     %     end
0901     %     eval(['load ' dirPath fileName '.' ext]) % inputFileName{1}
0902     %     if length(inputFileName) == 2
0903     %         eval(['waveform = ' inputFileName{2} ';'])
0904     %     else
0905     %         % Set default to assume that the variable name matched the file
0906     %         % name
0907     %         eval(['waveform = ' fileName ';'])
0908     %     end
0909     %     htt(1,:) = waveform(:)';
0910     %     htt(2,:) = zeros(1, length(waveform));
0911 else
0912     ext
0913     error('graven: inputFileName must have extention *.ilwd or *.txt')
0914 end
0915 
0916 return