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