Reference for data functions
GTIoverlap
Synopsis
computes the overlap of an interval [t1, t2] with a set of good time invervals
Usage
Double_Type GTIopverlap(Double_Type t1, Double_Type t2, Struct_Type GTI
See also: getCommonGTIs
atime_calc_beta
Synopsis
converts a pulse phaseshift into arrival times
Usage
(Double_Type[] arrtimes, errors) = atime_calc(Struct_Type atime, eph[, orb]);
or (Double_Type[] arrtimes, errors) = atime_calc(Double_Type[] t0, phi, error, Struct_Type eph[, orb]);
Qualifiers
none
Description
Takes the structure determined by 'atime_det' or an array of times, phase shifts and their errors and returns the pulse arrival times t = t0 + phi * p(t0) [+ z(t0)/c] To calculate the pulse period at each time a pulse ephemeris structure is needed (see 'check_pulseperiod_orbit_struct'). If an optional structure containing orbital parameters is given, the Doppler shift delay is included as well.
See also: atime_det, check_pulseperiod_orbit_struct, pulseperiod
atime_dataind
Synopsis
returns the datasets, which contain arrival times
Usage
Integer_Type[] atime_dataind();
Qualifiers
none
Description
Loops over all defined datasets and checks, which are containing arrival times. Therefore the function 'atime_metavalid' is used. If no dataset is found, the function returns NULL.
See also: atime_metavalid, define_atime
atime_det
Synopsis
Determines the arrival times from a lightcurve using a pulse pattern
Usage
Struct_Type atime_det(Struct_Type lc, Struct_Type[] pattern, Double_Type t0, Struct_Type ephemeris[, Struct_Type orbit]);
or Struct_Type atime_det(String_Type lc, Struct_Type[] pattern, Double_Type t0, Double_Type[] period[, Struct_Type orbit]);
Qualifiers
time - name of the time field in the FITS-file, see fits_read_lc for details indiv - determine individual pulses. If no value is assigned, all arrival times are re- turned. Otherwise a given integer sets the number of pulses to average movem - if 'indiv' is greater one, so it is aver- aged over a number of arrival times, a moving mean is used to get all individual pulses. Note that the pulses are then not statistically independent! ccfint - interpolates the cross correlation of the pulse pattern and the actual analyzed pulse to increase the accuracy. Only works well for clear signals! The number of interpolated bins between each original bin is set to the assigned number (default = 4) varerr - estimates the error by a local standard deviation. Therefore the variation of the given number of arrival times is used (default = 10) mcerr - estimates the error of each arrival time by performing Monte Carlo simulations. If an integer is assigned it sets the number of runs (default = 10000). Has a higher priority than 'varerr' mcgaus - fit a gaussian to the Monte Carlo distri- bution, if used for error estimation getpat - if multiple patterns are given, the index of the used one is stored in this variable, which has to be given as a reference (@var) match - reference to a variable (&var) where the matching pulses are saved as structures similar to the reference profile ccflim - matches below the given cross-correlation values are skipped (default = 0), allowed range is -1 to 1. chatty - boolean value for output messages (default = 1). If set to 2, also echo result of Monte Carlo error estimation debug - plot the cross correlation and the match- ing pulse, echo the found phase shift and sleep the given seconds (default = 1) or, if set to 'user' wait until a key is pressed
Description
Using one or more pulse pattern the arrival times of pulses in a lightcurve are determined using phase connection. The pattern must be given as a structure with the fields Double_Type bin_lo Double Type bin_hi Double_Type value Double_Type error, as, for example, returned by the 'epfold' function. The lightcurve can be passed by a structure with the fields 'time', 'rate', 'error' and 'fracexp', or by the filename to a FITS-file. To transform the phase shift to time, the pulse period at the actual position in the lightcurve is needed. For a first guess this may be a constant value (in seconds) or an array of two elements defining a range of periods (min/max values, in seconds). In the latter case the pulse period is searched by using the 'getVarPeriod' function and additional qualifiers are passed. However, the correct way is to pass the pulse ephemeris and additional orbital parameters to calculate the pulse poriod more precisely. The corresponding structures are described in the 'check_pulseperiod_orbit_struct' function. The values may be found by a fit to the arrival times using a constant pulse period. This leads to slightly different results, hence the determination of the arrival times and their analysis is an iterative procedure. If the 'indiv' qualifier is omitted the given lightcurve is folded and the returned arrival time corresponds to the phase shift of the given pattern to the resulting profile. If the method for estimating the uncertainty of the arrival times is not specified by qualifiers, the default uncertainty is set to one phase bin. The number of bins are derived from the given pulse pattern. The returned structure has the following fields: arrtimes - array of determined arrival times (MJD) error - their uncertainties (days) arrnum - relative pulse number of each arrival time to a specific pulse. Dramatically increases the speed of a fit. numref - index of the reference pulse, should be zero. Is updated if arrival times are merged. reft0 - the given reference time t0
See also: atime_merge, save_atime, define_atime, arrtimes, pfold
atime_det_beta
Synopsis
Determines the arrival times from a lightcurve using a pulse pattern
Usage
Struct_Type atime_det_beta(Struct_Type lc, Struct_Type[] pattern, Double_Type t0[, Double_Type period]);
or Struct_Type atime_det_beta(String_Type lc, Struct_Type[] pattern, Double_Type t0[, Double_Type period]);
Qualifiers
time - name of the time field in the FITS-file, see fits_read_lc for details indiv - determine individual pulses. If no value is assigned, all arrival times are re- turned. Otherwise a given integer sets the number of pulses to average movem - if 'indiv' is greater one, so it is aver- aged over a number of arrival times, a moving mean is used to get all individual pulses. Note that the pulses are then not statistically independent! ccfint - interpolates the cross correlation of the pulse pattern and the actual analyzed pulse to increase the accuracy. Only works well for clear signals! The number of interpolated bins between each original bin is set to the assigned number (default = 4) varerr - estimates the error by a local standard deviation. Therefore the variation of the given number of arrival times is used (default = 10) mcerr - estimates the error of each arrival time by performing Monte Carlo simulations. If an integer is assigned it sets the number of runs (default = 10000). Has a higher priority than 'varerr' mcgaus - fit a gaussian to the Monte Carlo distri- bution, if used for error estimation ccflim - matches below the given cross-correlation values are skipped (default = 0), allowed range is -1 to 1. numdif - maximum allowed deviation from expected pulse number (default = 0.5 phases) match - reference to a variable (&var) where the matching pulses are saved as structures similar to the reference profile slope - reference to a variable (&var) where the slope of the lightcurve is saved, which is calculated by an interpolation of the mean count rate in the single pulses skipsl - do not take the slope of the lightcurve into account mmnorm - by default the pulse pattern and the pulses are renormalized by f = (f - mean(f)) / sdev(f) If this qualifier is set, the min/max- normalizaton is used instead f = (f - max(f)) / (max(f) - min(f)) chatty - boolean value for output messages (default = 1). If set to 2, also echo result of Monte Carlo error estimation debug - plot the cross correlation and the match- ing pulse, echo the found phase shift and sleep the given seconds (default = 1) or, if set to 'user' wait until a key is pressed
Description
Using one or more pulse pattern the arrival times of pulses in a lightcurve are determined using phase connection. The patterns must be given as a structure with the fields Double_Type bin_lo Double Type bin_hi Double_Type value Double_Type error, as, for example, returned by the 'epfold' function. If multiple patterns are given, the best matching one is used to determine the arrival times. The lightcurve can be passed by a structure with the fields 'time', 'rate', 'error' and 'fracexp', or by the filename to a FITS-file. To determine individual pulses as enabled by the 'indiv' qualifier an approximate pulse period must be given (in seconds) to cut the lightcurve into segments were the pulses are searches. Instead of the pulse period an array of two elements defining a range of periods (min/max values, in seconds) may also be given. In that case the 'getVarPeriod' function is used to determine the pulse period and additional qualifiers are passed. If the 'indiv' qualifier is omitted the given lightcurve is folded and the returned phase shift corresponds to the phase shift of the used pattern to the resulting profile. If the method for estimating the uncertainty of the phase shifts is not specified by qualifiers, the default uncertainty is set to one phase bin. The number of bins are derived from the given pulse pattern. The returned structure has the following fields: t0 - array of first time bins of the pulses determined from the lightcurve (MJD) phi - array of phase shifts of the pulses with respect to the used pattern error - uncertainties of 'phi' arrnum - relative pulse number of each arrival time to a specific pulse. Dramatically increases the speed of a fit. numref - index of the reference pulse, should be zero. Is updated if arrival times are merged. reft0 - the given reference time t0 p0 - the given pulse period usedpat - index of the used pulse pattern in case of multiple ones given The arrival times can be calculated by t = t0 + phi * p0 as a first guess, since p0 is a function of time. This fact is handled properly by the fit function 'arrtimes'. The function 'atime_calc' calculates the arrival times according to the above equation with respect to a given pulse ephemeris and orbital parameters.
See also: atime_calc, atime_merge, save_atime, define_atime, arrtimes, pfold
atime_get_ephemeris
Synopsis
retrieves the actual model parameters of a given dataset of pulse arrival times
Usage
(Struct_Type, Struct_Type) atime_get_ephemeris(Integer_Type idx);
Qualifiers
none
Description
Returnes two structures containing the model parameters of a dataset, which contains pulse arrival times. The first structure describes the pulse ephemeris: ppuls - pulse period (s) pdot - first derivative (s/s) p2dot - second derivative (s/s^2) p3dot - third derivative (s/s^3) tpuls0 - reference time (MJD) The orbital parameters are stored in the second structure: porb - orbital period (d) torb0 - time of periastron passage (MJD) asini - projected semi major axis (lts) ecc - eccentricity omega - angle of periastron (degrees) If porb or asini is zero, the returned orbital structure will be set to NULL.
See also: arrtimes, check_pulseperiod_orbit_struct
atime_get_t0
Synopsis
returns the reference time of the arrival times
Usage
Double_Type atime_get_t0(Integer_Type[] dataid);
Qualifiers
none
Description
This function returns the emitting time of the pulse of number 0 as used in the polynomial to calculate pulse arrival times.
See also: arrtimes, pulse_time
atime_load_ephemeris
Synopsis
returns previously ssaved pulse ephemeris and orbital parameters
Usage
atime_load_ephemeris(Integer_Type idx, String_Type filename);
or (Struct_Type, Struct_Type) atime_load_ephemeris(String_Type filename);
Qualifiers
none
Description
The pulse ephemeris and orbital parameters, which have been previously saved by using the function 'atime_save_ephemeris', are either assigned to the given dataset 'idx' or returned as structures. The first one defines the pulse ephemeris and the second one the orbit. The following commands are equal: (eph, orb) = atime_load_ephemeris(filename); (eph, orb, ) = evalfile(filename);
See also: atime_save_ephemeris, atime_get_ephemeris, evalfile
atime_merge
Synopsis
merges two or more structures of arrival times into one
Usage
Struct_Type atime_merge(Struct_Type atime1, Struct_Type atime2);
or Struct_Type atime_merge(Struct_Type[] atime);
Qualifiers
nosort - do not sort the merged arrival times
Description
Either both given structures or an array of structures containing arrival times are merged. The fields of the structures must be equal to the ones described in 'atime_det'. The relative pulse numbers and their references are updated to the indices of the new array of arrival times. In addition the arrival times are sorted by time unless the 'nosort' qualifier is given.
See also: atime_det
atime_merge_beta
Synopsis
merges two or more structures of arrival times into one
Usage
Struct_Type atime_merge(Struct_Type atime1, Struct_Type atime2);
or Struct_Type atime_merge(Struct_Type[] atime);
Qualifiers
nosort - do not sort the merged arrival times
Description
Either both given structures or an array of structures containing arrival times are merged. The fields of the structures must be equal to the ones described in 'atime_det'. The relative pulse numbers and their references are updated to the indices of the new array of arrival times. In addition the arrival times are sorted by time unless the 'nosort' qualifier is given.
See also: atime_det
atime_metavalid
Synopsis
checks if a given dataset of arrival times contains valid metadata
Usage
Integer_Type[] atime_metavalid(Integer_Type index);
Qualifiers
none
Description
If a dataset of arrival times is defined using the 'define_atime' function, there are also a lot of addiotional needed informations defined in the metadata. These are used by fitting the data by the 'arrtimes' fit function. This function can be used to find all defined data containing arrival times, which is implemented in 'atime_dataind'.
See also: define_atime, arrtimes, get_dataset_metadata, atime_dataind
atime_save_ephemeris
Synopsis
saves pulse ephemeris and orbital parameters into as an S-lang script
Usage
atime_save_ephemeris(Integer_Type idx, String_Type filename);
or atime_save_ephemeris(Struct_Type eph[, Struct_Type orb], String_Type filename)
Qualifiers
none
Description
The fitted pulse ephemeris and orbital paramters of the given dataset or the structures themself are saved into a file. This file is an S-lang script and can be called directly or loaded by using 'atime_load_ephemeris'.
See also: atime_load_ephemeris, atime_get_ephemeris
atime_set_ephemeris
Synopsis
sets the actual model parameters of a given dataset of pulse arrival times
Usage
atime_set_ephemeris(Integer_Type[] idx, Struct_Type eph[, Struct_Type orb]);
or atime_set_ephemeris(Integer_Type[] idx, Double_Type eph[, Struct_Type orb]);
Qualifiers
sett0 - also set the reference time of the pulse ephemeris, which is used internaly
Description
The given structure(s) containing the pulse ephemeris and the orbital parameters are used to set the model parameters of dataset 'idx'. Instead of the pulse ephemeris the pulse period may be given only. Other- wise the structure has to contain: ppuls - pulse period (s) pdot - first derivative (s/s) p2dot - second derivative (s/s^2) p3dot - third derivative (s/s^3) tpuls0 - reference time (MJD) The orbital parameters are stored in the second structure: porb - orbital period (d) torb0 - time of periastron passage (MJD) asini - projected semi major axis (lts) ecc - eccentricity omega - angle of periastron (degrees) By default the reference time of the pulse ephemeris is NOT set, because this should be fixed and set to the same time as used for the pulse pattern. This time therefore defines the arrival time of the pulse of number n=0. Only change this time if you know what effects will result (see arrtimes function)!
See also: arrtimes, atime_get_ephemeris
atime_shiftref
Synopsis
switches the shift of the reference pulse number
Usage
atime_shiftref(Integer_Type[] dataid[, Integer_Type boolean]);
Qualifiers
none
Description
While fitting pulse arrival times a reference pulse number may be used. The reference pulse is the first one of the dataset, but while fitting it is by default shifted to the pulse nearest to the fixed reference time 'tpuls0' of the pulse ephe- meris. With this function this shift can be turned off (boolean=0) or on (boolean=1, default).
See also: arrtimes, atime_useref, atime_det
atime_sim
Synopsis
simulates pulse arrival times
Usage
Struct_Type atime_sim([Integer_Type n,] Double_Type tmin, tmax, Struct_Type eph[, orb]);
Qualifiers
scramble - 1 sigma noise of the simulated arrival times in seconds (default = 0) dn - the increment of the arrival times if they are not randomly generated (default = 1)
Description
Simulates the orbit and the pulsation of a pulsar and returns the pulse arrival times in the time range given by 'tmin' and 'tmax'. If the number of arrival times 'n' is given, then n randomly selected arrival times are returned. The returned structure is equal to the atime_det function. The uncertainties of the simulated arrival times are calculated using the 'scramble' value. If no value is given, the uncertainties are set to 0.1% of the pulse period to avoid ISIS from treating the arrival times as counts, resulting in poisson errors. The structures describing the pulse ephemeris and the orbit must fullfil the conditions described in the function check_pulseperiod_orbit_struct.
NOTE Due to speed reasons the returned arrival times might be beyond the given time interval or not all pulses in this interval are returned.
See also: atime_det, check_pulseperiod_orbit_struct, arrtimes
atime_sim_orbitimpact
Synopsis
determines the impact of uncertainties of orbital parameters on the pulse ephemeris
Usage
atime_sim_orbitimpact(Double_Type tmin, tmax, Struct_Type eph, orb, dorb);
Qualifiers
tpd - number of simulated arrival times per day (default: 10) pars - parameters of the pulse ephemeris which are used to fit the simulated data (default: ["ppuls", "pdot", "p2dot", "p3dot"]) range - allowed fitting ranges of the parameters given by the 'pars' qualifier in accordant units. The ranges have to be given as an array of alternate min/max values for each parameter in 'pars'
Description
An uncertainty of an orbital parameter may lead to a wrong pulse ephemeris after a successful fit of pulse arrival times. This function deter- mines the impact of uncertainties on the pulse ephemeris by simulating arrival times in the range 'tmin' and 'tmax' and fitting them within the given orbital uncertainties. The returned structure itself contains structures for each orbital parameter, which contain the fitted pulse ephemeris and therefore the impact of this parameter. The highest phase shift left in the residuals is also returned for each orbital parameter. The given orbital and pulse ephemeris structure must fullfil the conditions described in the check_pulseperiod_orbit_struct function. The uncertainties of the parameters are itself given by a structure, which may either contain double types for the errors or arrays with two elements for the lower and upper errors.
See also: atime_sim, check_pulseperiod_orbit_struct, arrtimes
atime_sim_residuals
Synopsis
simulates arrival times to calculate how the residuals depend on parameter variations
Usage
atime_sim_residuals(Double_Type tmin, tmax, Struct_Type eph, orb, String_Type psfile);
Qualifiers
n - number of variations for each parameter and sign (-> 2*n variatons, default: 3) tpd - number of simulated arrival times per day (default: 10) mres - maximum allowed residuals (=yrange) in units of pulse phase (default: 0.45), must be in the range of 0.0-0.45 minr - smallest allowed relative variation of a parameter, used to determine the overall variation of this parameter to produce nice plots (default: 1e-10) pars - only simulate this given parameters (default: all except tpuls0 and pporb)
Description
This function determines the dependency of the residuals on the given orbit and pulse ephemeris. Therefore pulse arrival times are simulated and the fit parameters are variied in the way, that the shape of the resulting residuals can be seen well. The resulting plots are stored in a Post- Script file specified by the filename 'psfile'. The labels show the absolute variation of the the parameter. The simulated times will be in the range of 'tmin' to 'tmax' (in MJD). The orbital parameter and the pulse ephemeris structure must be in the same form described in the function check_pulseperiod_orbit_struct.
See also: atime_sim, check_pulseperiod_orbit_struct, arrtimes
atime_useref
Synopsis
switches the use of a reference pulse number
Usage
atime_useref(Integer_Type[] dataid[, Integer_Type boolean]);
Qualifiers
none
Description
While fitting pulse arrival times the pulse number of each pulse must be determined. By using the atime_det function to determine the arrival times each pulse gets a pulse number relative to the first pulse found in the input lightcurve. Once the pulse number of this reference pulse is deter- mined during a fit, the numbers of all other pulses are also set. Using this function the use of rela- tive pulse numbers during a fit can be turned off (boolean=0) or on (boolean=1, default).
See also: arrtimes, atime_det, atime_shiftref
atime_xinclude
Synopsis
Specifies the datasets, which should be taken into account in the fit
Usage
atime_xinclude(Integer_Type[] index);
Qualifiers
none
Description
During the fit only the given datasets containing arrival times are taken into account. The remain- ing ones are noticed such that no bins are used. Hence the fit function and parameters are not affected by changed dataset indices.
See also: xnotice_atime, define_atime
check_pulseperiod_orbit_struct
Synopsis
checks if the pulse period- and/or orbit-structure are valid
Usage
Integer_Type check_pulseperiod_orbit_struct(Struct_Type structure);
or Struct_Type check_pulseperiod_orbit_struct(Double_Type pulseperiod);
Description
This function performs a checks whether the given structure matches one of the following definitions of a pulse period- or orbital parameter structure:
-
pulseperiod struct { % given as period over time Double_Type[] time, % in MJD Double_Type[] period % in seconds }
-
pulseperiod struct { % given as taylor coefficients Double_Type t0, % reference time in MJD Double_Type p0[, % pulse period at t0 in seconds Double_Type pdot[, % first derivative in s/s Double_Type p2dot[, % second derivative in s/s^2 ... Double_Type pNdot]]] % higher orders }
-
orbit struct { Double_Type tau or t90, % time of periastron passage % or mean longitude of 90 degrees in MJD Double_Type porb, % orbital period in days Double_Type asini, % projected semi-major axis in lt-s Double_Type ecc[, % eccentricity Double_Type omega] % longitude of periastron in degrees % required if ecc > 0 }
On a successful match the number of the definition (1-3) is returned, 0 otherwise. Note that additional field names are allowed.
If an input pulse period (as a number) is given instead a structure, a new pulse period structure (type 2) is returned with t0 set to 0 and the pulse period as given.
Note that units cannot checked by this function, but are a suggestion. Read the help of any function using the same definition on their required units.
See also: pulseperiod
colacal
Synopsis
Timing Tools: Coherence and Lag Calculation
Usage
([6 outputs]) = colacal ([10 required inputs]);
Description
Input:
freq - Fourier Frequency Array
cpd - CPD Array
noicpd - Poisson Noise Contribution to CPD
lopsd - Non-noise-corrected PSD (low channel)
hipsd - Non-noise-corrected PSD (high channel)
noilopsd - Noise-Contribution to low PSD
noihipsd - Noise-Contribution to high PSD
siglopsd - Noise-corrected PSD (low channel)
sighipsd - Noise-corrected PSD (high channel)
alln - number of averaged segments in each frequency bin
Output:
rawcof - non-noise-corrected coherence function [Vaughan & Nowak, 1997, ApJ, 474, L43 (Eqn. 2)]
cof - noise-corrected coherence function [Vaughan & Nowak, 1997, ApJ, 474, L43 (Eqn. 8, Part 1)]
errcof - one-sigma uncertainty of cof [Vaughan & Nowak, 1997, ApJ, 474, L43 (Eqn. 8, Part 2)]
lag - time lag [Nowak et al., 1999, ApJ, 510, 874 (Sect. 4)]
errlag - one-sigma uncertainty of lag [Nowak et al., 1999, ApJ, 510, 874 (Eqn. 16)]
sigcpd - noise-corrected cross-power-density [sigcpd = cpd - noicpd]
color_color_data
Synopsis
extracts light curves in 3 consecutive bands and the corresponding colors from an event list
Usage
Struct_Type dat = color_color_data(Struct_Type evts, Double_Type E0, E1, E2, E3);
Qualifiers
- dt: time resolution [default: 100]
- field: field used for selection of the bands [default: "pi"]
Description
dat.a
= lightcurve in "E0 <= field < E1" band
dat.b
= lightcurve in "E1 <= field < E2" band
dat.c
= lightcurve in "E2 <= field < E3" band
dat.
x_
y = x/y color for x, y = a | b | c
See also: histogram
data_map_function
Synopsis
executes a function on specific datasets
Usage
data_map_function(String_Type filter, Ref_Type function, arguments...; qualifiers);
Description
Loops over all defined datasets and matches each data information against the given filter. In the simplest way, this filter is a string respresenting the wanted instrument as specified in the field of 'get_data_info'. If the 'fi' qualifier is provided, the given filter is interpreted as if-statement, which has to be fullfiled. Thereby, the associated data information can be accessed via the 'info' variable. The given function is applied to each matching dataset in the form function(dataset, arguments...; qualifiers);
Example
% set the energy ranges for all RXTE-PCA sepctra data_map_function("PCA", &xnotice_en, 3.5, 50);
% apply grouping to all Suzaku-PIN spectra
data_map_function(is_substr(info.file, "hxd_pin")
,
&group; min_sn = 40, fi);
See also: get_data_info, eval, array_map
dcf
Synopsis
compute the Edelson & Krolik discrete correlation function and returns a structure with the correlation
Usage
dcf(time_a, val_a, tima_be, val_b);
Qualifiers
- err_a [=0]: array containing the measurement errors of val_a. If given, this is taken into account by using eq. (3) of Edelson & Krolik. NOTE: There are lots of problems with interpreting the DCF computed this way instead of using Edelson & Krolik, eq. (2), and it is NOT recommended to give erra and errb (see, e.g., White & Peterson, 1994, PASP, 106, 879)
- err_b [=0]: as err_a except for val_b
- minlag [=mininum time difference between time_a and time_b]: minumum lag to consider
- maxlag [=maximum time difference between time_a and time_b]: maximum lag to consider
- numf [=int(min([length(val_a) ,length(val_b) ])/10]: number of lag bins for which to compute the DCF
- minpt [=10]: minimum number of data points for a DCF value to be
Description
This is the imported version of the IDL subroutine dcf.pro
It will compute the Edelson & Krolik discrete correlation function.
define_atime
Synopsis
defines a dataset of arrival times
Usage
Integer_Type define_atime(Struct_Type atime[, Integer_Type divide]);
Qualifiers
modnum - a reference to a variable, where to store the pulse numbers found by the model noff - do not add 'arrtimes' to the actual fit- function correspondig to the dataset
Description
Takes the given arrival times structure to define an ISIS dataset. This structure may be created by 'atime_det' or loaded by 'load_atime'. Accordingly all usual ISIS routines, which operate on datasets, e.g. fitting, can be used. The fit function 'arrtimes' handle arrival times including orbital motion and pulse ephemeris. If not disabled by the 'noff' qualifier the actual fit function is auto- matically extended by the resulting dataset using the 'arrtimes' model. It also takes the qualifier 'modnum' into account to return the pulse numbers determined by the model. The reference to the pulse numbers is stored into the metadata of the dataset. If the second parameter 'divide' is given the dataset is defined several times accordingly to the given integer. The datasets are then noticed automatically such that the whole data is divided into the given number of parts of equal length. Each part may be then fitted individually. The returned integer corresponds to ISIS dataset index.
See also: arrtimes, atimes_det, load_atime, atime_metavalid, xnotice_atime
dof
Synopsis
Number of degrees of freedom
Usage
Double_Type = dof (hist_index);
##### Description
Use this function to retrieve number of
degrees of freedom.
EXAMPLE
isis>xray = load_data("data.pha"); isis>variable num = dof(1);
See also: num_bin
e_bv
Synopsis
Calculates the E(B_V) color excess after Predehl & Schmitt (1995)
Usage
Double_Type = e_bv(Double_Type N_H);
Qualifiers
R_V - Scalar specifying the ratio of total to selective extinction R(V) = A(V) / E(B - V). If not specified, then R = 3.1 Extreme values of R(V) range from 2.3 to 5.3
Description
From a given hydrogen absorption column density N_H in the diredtory of the object, the color excess E(B-V) is calculated after Predehl & Schmitt (1995): E(B-V) = N_H/(1.79e21*R_V).
EXAMPLE Calculate the color excess for Cen A (RA 13h25m27.6s DEC -43d01m09s) for N_H = 8.09e20 and an "average" reddening of for the diffuse interstellar medium (R(V) = 3.1).
isis> N_H = 8.09e20; isis> ebv = e_bv(N_H); isis> print(ebv);
See also: fm_unred;
energyflux
Synopsis
evaluates the energy flux of the current fit model in a given energy range
Usage
Double_Type energyflux(Double_Type Emin, Emax)
Qualifiers
- factor [=
1.001
]: step size of the logarithmic energy grid - Emin: minimum energy of (extended) grid
- Emax: maximum energy of (extended) grid
- cgs: return flux in erg/cm^2/s
Description
This function calculates the energy flux of the current best fit model
in keV/cm^2/s by integrating over the model on a logarithmic energy grid.
The logarithmic step size of the model is given by the factor
qualifier.
In other words, energyflux returns
int_{Emin}^{Emax} E * S_E(E) dE
(in keV/s/cm^2)
where S_E(E)
(in 1/s/cm^2/keV) is defined by fit_fun
.
If the fit-function is a convolution model, e.g., Compton reflection,
it may be necessary to use an extended grid defined by the Emin
and Emax
qualifiers.
Note: While physicists generally use the term flux density'' for the quantity returned by this function, with the exception of radio astronomy astronomers generally call the returned value the
flux''.
The change in function name is based on the experience that most
users of isisscripts did not recognize that this function is what
they required.
See also: eval_fun, eval_fun_keV, luminosity
energyfluxdensity
Synopsis
deprecated function
Usage
do not use
Description
This function had the wrong name. Please use the energyflux function instead.
See also: energyflux
epferror
Synopsis
estimate epoch folding error with monte-carlo simulation approach.
Usage
(Double_Type err) = epferror(Double_Type time, rate, period[, rate_error]);
Qualifiers
- pstart: start period for epoch folding period search (default 0.5*period)
- pstop: stop period for epoch folding period search (default 1.5*period)
- sampling: periods per peak for epoch folding
- ntrial: number of MC loops (default 20)
- nbins: number of bins for the pulse profile (default 20)
- dt: exposure of every time bin.
- chatty: chattiness of the function.
- fchatty: chattiness piped to epfold.
Description
This routine tries to estimate the error of a previous received period using the epoche folding approach. It is adopted from the IDL script with the same name, but does not yet allow for GTI or Poisson statistics. It follows the following basic points: 1.) calculate a mean profile with given period. 2.) compute the intensity for all times applying the period multiplied profile. 3.) simulate an error for all times (assuming normal distribution with sigma = error, or, if not given sqrt(rate)). 4.) Perform epoch folding for that simulated lightcurve. 5.) Go to step 2.) Ntrial times, sum up the maximum of epoche folding found. 6.) Compute the standard deviation of the Ntrial maxima obtained and take these as the error.
NOTE: The processing may last a long time (be prepared to wait hours to weeks)!
References: Davies, S.R., 1990, MNRAS 244, 93 Larsson, S., 1996, A&AS 117, 197 Leahy, D.A., Darbro, W., Elsner, R.F., et al., 1983, ApJ 266, 160-170 Schwarzenberg-Czerny, A., 1989, MNRAS 241, 153
epfold
Synopsis
peforms epoch folding on a lightcurve in a given period interval
Usage
(Struct_Type res) = epfold(Double_Type t, r, pstart, pstop);
or (Struct_Type res) = epfold(Double_Type t, pstart, pstop) ; (event data)
Qualifiers
- nbins: number of bins for the pulse profile
- exact: calculate the pulse profile in a more exact way, see description of pfold (not recommed as it takes a very long time!).
- dt:exposure of every lightcurve time bin, should be given to ensure correct results.
- sampling: how many periods per peak to use (default=10)
- nsrch: how many periods to search in a linear grid (default not set)
- dp: delta period of linear period grid (default not set)
- lstat: use L-statistics instead of chi^2 statistics
- chatty: set the verbosity, chatty-1 is piped to pfold (default=0)
- gti:GTIs for event data, given as struct{start=Double_Type, stop=Double_Type}
Description
Performs epoch folding on a given lightcurve between the periods pstart and pstop. The function was adopted from the IDL routine of the same name. GTI correction only implemented for event-data yet.
By default, periods are sampled according to the triangular rule for estimating the period error, using "sampling" periods per peak. If a linear grid is to be used, either "dp" for a given distance between two consecutive periods or "nsrch" for a given number of periods can be given. These qualifiers are mutually exclusive.
The returned structure "res" contains four fields: "p" for the evaluated period and "stat" for the value of the statistic used. Additionally the field "nbins" contains the number of phase bins used and "badp" contains indices for res.p marking periods where the pulse profile showd empty phase bins. Values of res.stat[res.badp] should be taken with great care!
Compared to the similar function sitar_epfold_rate, epfold.sl uses the chi^2 statistic as a default and is based on pfold.sl, which can take errors of the lightcurve into account. If the qualifier "lstat" is given, the statistic is switched to the l-stat statistic as in sitar_epfold_rate, but errors of the lightcurve are no longer taken into account. lstat is not available for event-data.
If the "exact" qualifier is given, the function takes the exposure time of every time bin into account in the sense that, that a given time bin may overlap over two phase bins. The corresponding exposure time in every phase bin is reduced accordingly.
NOTE: the "exact" qualifiers slows the script considerably down, depending on the length of the lightcurve and the number of bins for up to a factor of >100!
If "exact" is not given, the script is ~10% slower than sitar_epfold_rate.
The script is still in the develepment phase, please report any bugs or missing features to Felix (felix.fuerst@sternwarte.uni-erlangen.de).
NOTE: Please read Davies, S.R., 1990, MNRAS 244, 93 (L-statistics!) Larsson, S., 1996, A&AS 117, 197 Leahy, D.A., Darbro, W., Elsner, R.F., et al.,1983, ApJ 266, 160-170 Schwarzenberg-Czerny, A., 1989, MNRAS 241, 153
epfoldpdot
Synopsis
peforms epoch folding on a lightcurve in a given period interval
Usage
(Struct_Type res) = epfoldpdot(Double_Type t, r, pstart, pstop);
or (Struct_Type res) = epfoldpdot(Double_Type t, pstart, pstop) ; (event data)
Qualifiers
- nbins: number of bins for the pulse profile
- exact: calculate the pulse profile in a more exact way, see description of pfold (not recommed as it takes a very long time!).
- dt:exposure of every lightcurve time bin, should be given to ensure correct results.
- sampling: how many periods per peak to use (default=10)
- nsrch: how many periods to search in a linear grid (default not set)
- dp: delta period of linear period grid (default not set)
- lstat: use L-statistics instead of chi^2 statistics
- chatty: set the verbosity, chatty-1 is piped to pfold (default=0)
- gti:GTIs for event data, given as struct{start=Double_Type, stop=Double_Type}
- pdstart: start p-dot for grid search (default=0)
- pdstop: stop p-dot for grid search (default=0)
- pdnsrch: search point for p-dot grid, (default = nsrch (if defined, otherwise 10))
Description
Performs epoch folding on a given lightcurve between the periods pstart and pstop and period derivatives (p-dot) pdstart and pdstop. The function was adopted from the IDL routine of the same name. GTI correction only implemented for event-data yet.
By default, periods are sampled according to the triangular rule for estimating the period error, using "sampling" periods per peak. If a linear grid is to be used, either "dp" for a given distance between two consecutive periods or "nsrch" for a given number of periods can be given. These qualifiers are mutually exclusive.
P-dot is always sampled on a linear grid with pdnsrch points.
The routine uses the S-lang function 'parallel_map' to use Isis_Slaves.num_slaves cores on your machine for parallel computing. On a MacBook Pro with quadcorse i7 this gives a speed improvement of a factor ~>2 over single core calculations.
The returned structure "res" contains five fields: "p" for the evaluated period, pd for the evaluated period derviates, and "stat" for the value of the statistic used. "stat" is a 2D array of dimnesion (np, npd). Additionally the field "nbins" contains the number of phase bins used and "badp" contains indices for res.p marking periods where the pulse profile showd empty phase bins. Values of res.stat[res.badp] should be taken with great care!
Compared to the similar function sitar_epfold_rate, epfold.sl uses the chi^2 statistic as a default and is based on pfold.sl, which can take errors of the lightcurve into account. If the qualifier "lstat" is given, the statistic is switched to the l-stat statistic as in sitar_epfold_rate, but errors of the lightcurve are no longer taken into account. lstat is not available for event-data.
If the "exact" qualifier is given, the function takes the exposure time of every time bin into account in the sense that, that a given time bin may overlap over two phase bins. The corresponding exposure time in every phase bin is reduced accordingly.
NOTE: the "exact" qualifiers slows the script considerably down, depending on the length of the lightcurve and the number of bins for up to a factor of >100!
NOTE: "exact" qualifier is currenlty untested and may lead to erronoeus results!
If "exact" is not given, the script is ~10% slower than sitar_epfold_rate.
The script is still in the develepment phase, please report any bugs or missing features to Felix (felix.fuerst@fau.de).
NOTE: Please read Davies, S.R., 1990, MNRAS 244, 93 (L-statistics!) Larsson, S., 1996, A&AS 117, 197 Leahy, D.A., Darbro, W., Elsner, R.F., et al.,1983, ApJ 266, 160-170 Schwarzenberg-Czerny, A., 1989, MNRAS 241, 153
eval_fun2_keV
Synopsis
evaluate a fit-function on a user-defined energy-grid
Usage
flux = eval_fun2_keV(handle, E_bin_lo, E_bin_hi[, params[, args]]);
Description
The fit-function given by handle
(S_E(E)
) is evaluated
on an arbitrary grid defined by E_bin_lo
and E_bin_hi
(in keV):
flux = int_{E_bin_lo}^{E_bin_hi} S_E(E) dE
The unit of flux
is ph/s/cm^2/bin.
See also: eval_fun2, eval_fun_keV
eval_fun_keV
Synopsis
evaluate the fit-function on a user-defined energy-grid
Usage
Double_Type flux[] = eval_fun_keV(Double_Type E_bin_lo[], E_bin_hi[]);
Description
The currently defined fit-function S_E(E)
is evaluated
on an arbitrary grid defined by E_bin_lo
and E_bin_hi
(in keV):
flux = int_{E_bin_lo}^{E_bin_hi} S_E(E) dE
Note that flux is bin integrated, i.e., the unit of flux
is
ph/s/cm^2. To plot the flux density (ph/s/cm^2/keV), you need to
divide this quantity by the bin width de=E_bin_hi - E_bin_lo.
See also: eval_fun
factorized_arf_rmf
Synopsis
changes an ARF/RMF pair into factorized ARF/normalized RMF
Usage
(newARF, newRMF) = factorized_arf_rmf(ARF, RMF);
See also: factor_rsp
fake_pulsar_lightcurve
Synopsis
creates a synthetic lightcurve of a pulsar
Usage
Struct_Type synthetic_pulsar_lightcurve(
Struct_Type lightcurve or Double_Type[] time,
Double_Type or Struct_Type period[, Struct_Type orbit
[, Struct_Type profile[, Struct_Type or Double_type fluxevolution
[, Ref_Type noise_fun]]]]
);
Qualifiers
lcdt - time resolution of the input lightcurve in days (default: difference of first two time bins) interpol - function reference used to to time-grid interpolations (default: &interpol_points) pfold - structure of qualifiers to be passed to 'pfold' (default: struct { nbins = 32, dt = ..., t0 = ..., pdot = ... }) fluxlc - structure of qualifiers to be passed to 'pulse2pulse_flux_lc' (default: NULL) tophase - structure of qualifiers to be passed to 'pulseperiod2phase' (default: struct { t0 = ... }) chatty - show or hide output messages (default: 1)
Description
This function fakes a pulsar's lightcurve including the following aspects:
- longterm flux evolution (on timescales larger than the pulse period)
- lightcurve modulation by pulse profile
- pulse period change including orbital motion
- gaussian or user-defined observation noise
Two modi are possible: a) providing an observed initial lightcurve from which all needed aspects are derived. A pulse period or its evolution is mandatory. Certain aspects can be overwritten by user input. The resulting faked and the input lightcurve have the same time-grid. b) providing the output time-grid. The to be included aspects have to be given explicitely.
See also: check_pulseperiod_orbit_struct,pulse2pulse_flux_lc,pulseperiod2phase,pfold
filter_gti
Usage
Integer_Type ind[] = filter_gti(Double_Type time[], Struct_Type gti);
or filter_gti(Double_Type time_lo[], time_hi[], Struct_Type gti);
Qualifiers
- minfracexp: minimum fractional exposure a time bin has to have, otherwise it is considered bad (lightcurves only, default: 1e-4)
- fracexp: if set to a reference, returns the fractional exposure of each time bin (lightcurves only)
- exposure: if set to a reference, returns the livetime in each time bin (lightcurves only)
- indarray: return an array of index arrays instead (in case of events only)
Description
For a lightcurve given by time_lo
and time_hi
, this function
returns all indices ind
to the time bins, for which the fractional
exposure time as defined by the good time intervals defined by gti
is at least fracexp
.
For a list of events measured at times time
, return a list of indices
to all events that were measured during the given set of good time intervals.
The good time intervals are defined by a struct {start,stop}
where
start
and stop
define the start and stop times.
The code assumes that
- time and time_hi are in ascending order
- the gti does not contain any overlapping intervals
- gti, time, and time_hi have the same unit
Warning: the function applies a time sorting to the gti if necessary. This will MODIFY THE INPUT since structures are passed as references in S-Lang!
find_peak
Synopsis
Written for epoch folding purposes, it finds peaks of certain characteristics
Usage
Struct_Type find_peak(Struct_Type input)
Description
The input should be a Struct_Type of the form: { p, stat, [ lc ], [ nbins_epfold ] , [ expectation ] , [ sigma ] } The individual fields are: * p, stat: arrays of doubles, the peak is looked for in stat * lc: struct { time, rate }, the original lightcurve * nbins_epfold: necessary in order to obtain an error estimate * expectation: approximate location of the peak in p * sigma: approximate error the expectation value has lc and nbins_epfold reflect the original application of this function, but it can be used without these arguments to find peaks in other curves. The output is a Struct_Type of the form: { period, error, [ profile ], err_area, err_theory, width_lo, width_hi, statvsp, flags, [ bayes ], [ fit ], badness } flags is a structure of indicators regarding the quality of the output, badness is a coarse estimate how well the process worked. A value of 0 is desirable, 1-2 could be acceptable, 3-4 is usually an indication of bad input data and for larger values something has seriously gone wrong. Fields marked [] will not be output if the qualifiers chosen require it. In the default configuration, this function works the following way:
(1) A bayesian block analysis finds possible maxima * if the block analysis fails, a smoothing algorithm is used to remove noise in the curve, since the block analysis is quite susceptible to noise (2) The maxima are ranked considering: * the width (broad is an advantage, but very broad maxima in comparison to dp = p^2 / T resp sigma are discarded) * the distance to the expected value, where anything within sigma is ranked approximately equal (3) Maxima with a quality exceeding the tolerance level are taken and the accurate peak is found: * a first estimate is a weighted mean of points within the maximum block * the accurate value is determined by fitting a given xy-function to the peak. The proportion of dp used for choosing the range for the fit is controlled by the qualifier fit_part. If you do not expect considerable secondary maxima, this variable can be set as big as 1. However, if secondary peaks are expected to occur, a value this large can make the fit useless. (4) For each possible peak the lightcurve is folded and the resulting pulse profile is examined. The peak corresponding to the best overall result is taken.
If the fit_fct qualifier is used, the function has to have the following parameters (in this specific order): (1) peak (2) measure for the height (3) measure for the broadness of the peak (4) absolute offset (5) measure for the asymmetry of the peak All qualifiers can also be passed using a structure named 'find_peak_qualifiers'. If this qualifier is given, its content will overwrite all other qualifiers. If using this structure, all qualifiers are expected to have values (i.e. only_max = 1). This is helpful if you want to make things easier to read in scripts.
Qualifiers
- only_max: neither block analysis nor fit are performed
- blocks_but_only_max: block analysis is performed but only the highest bin in the best block taken
- weighting: weights for properties of a peak [area, sharpness, distance to expectation, profile quality] (default: [1,1,1,1])
- no_fit: nlock analysis is performed but only the first estimate taken
- fit_fct: an xyfit function for the fit, passed as a string (default: "sqrsinc")
- fit_part: measure for the part of the peak used for the fit (default: .3)
- tolerance: determines how many peaks are passed to step (4), in [0,1) (default: .99)
- smoothing: how strong the smoothing algorithm works (default: 10.)
- ncp_prior: argument for block analysis. If a string is given, the default block analysis routine will be used (default: 100)
- fp_rate: argument for block analysis, see its help
- no_profile: No pfold will be executed. Corresponds to tolerance -> 1 and less computation
- nbins_pfold: argument for pfold, see its help
- pfold_not_exact: changes how pfold is executed. Recommended as long as the exact-bug is not fixed
- exact: passes the information that the statistic was calculated using the exact qualifier
- flag_info: detailed information regarding the content of the flags structure is given, NOTHING else is done
- error_info: detailed information about the calculation of errors is given, NOTHING else is done
- chatty: boolean value (default: 1);
See also: pulseperiod_search, epfold, pulseperiod_epfold, split_and_epfold_lc, bayesian_blocks, pfold
fits_save_data_model_flux
Synopsis
saves ISIS spectral data into a FITS file
Usage
fits_save_data_model_flux([filename[, ids]]);
flux2lum
Synopsis
Calculates the source luminosity
Usage
Double_Type lum = flux2lum (Double_Type flux, Double_Type z);
Qualifiers
- silent: If set, the program will not display adopted cosmological parameters at the terminal.
- h0 [=70]: Hubble parameter in km/s/Mpc
- omega_m [=0.3]: Matter density, normalized to the closure density, default is 0.3. Must be non-negative.
- omega_lambda [=0.7]: Cosmological constant, normalized to the critical density.
- omega_k [=0]: curvature constant, normalized to the critical density. Default is 0, indicating a flat universe.
- q0 [=-0.55]: Deceleration parameter, numeric scalar = -R*(R'')/(R')^2
Description
This function calculates the luminosity of a source
given a flux
in erg/s/cm^2 as well as a redshift z
using
the cosmological parameters specified by the qualifiers
(which are passed to the function cosmo_param
). The
distance is calculated with the function lumdist
.
Flux
and z
can be scalars or vectors. Use
energyflux to calculate the flux.
Example
variable flux=3e-13; variable z=0.0705; variable l = flux2lum(flux,z);
See also: energyflux, cosmo_param, lumdist
fm_unred
Synopsis
Deredden a flux vector using the Fitzpatrick (1999) parameterization
Usage
Double_Type = fm_unred(Double_Type[] wave, Double_Type[] flux, Double_Type ebv);
Qualifiers
N_H - Hydrogen absorption column R_V - Scalar specifying the ratio of total to selective extinction R(V) = A(V) / E(B - V). If not specified, then R = 3.1 Extreme values of R(V) range from 2.3 to 5.3 LMC2 - If set, then the fit parameters are set to the values determined for the LMC2 field (including 30 Dor) by Misselt et al. Note that neither 'AVGLMC' or 'LMC2' will alter the default value of 'R_V' which is poorly known for the LMC. AVGLMC - If set, then the default fit parameters c1,c2,c3,c4,gamma,x0 are set to the average values determined for reddening in the general Large Magellanic Cloud (LMC) field by Misselt et al. (1999, ApJ, 515, 128) extcurve - If set to a variable as Ref_Type, the E(wave-V)/E(B-V) extinction curve is returned, interpolated onto the input wavelength vector gamma - Width of 2200 A bump in microns (default = 0.99) x0 - Centroid of 2200 A bump in microns (default = 4.596) c1 - Intercept of the linear UV extinction component (default = 2.030 - 3.007 * c2) c2 - Slope of the linear UV extinction component (default = -0.824 + 4.717 / R_V) c3 - Strength of the 2200 A bump (default = 3.23) c4 - FUV curvature (default = 0.41) wilm - If N_H value has been determined using wilm abundances and Verner cross-sections use this updated correlation (Nowak et al., 2012) ngc3227 - Use the reddening curve for NGC3227 after Crenshaw, D.~M., Kraemer, S.~B., Bruhweiler, F.~C., & Ruiz, J.~R. 2001, ApJ, 555, 633 provided locally in /home/beuchert/work/ngc3227/scripts/reddening/
Description
The unreddened flux vector of the input 'flux' is calculated on the wavelength vector 'wave' using the Fitzpatrick (1999) parameterization. The scalar 'ebv' is the color excess E(B-V). If a negative EBV is supplied, then fluxes will be reddened rather than dereddened. If the scalar is not known the hydrogen absorption column in the directory of the object must be given as a qualifier and E(B-V) is calculated by E(B-V) = N_H/(1.79e21*R_V).
The R-dependent Galactic extinction curve is that of Fitzpatrick & Massa (Fitzpatrick, 1999, PASP, 111, 63; astro-ph/9809387 ). Parameterization is valid from the IR to the far-UV (3.5 microns to 0.1 microns). UV extinction curve is extrapolated down to 912 Angstroms. This Function is adopted from the IDL-Function fm_unred.
The five input qualifiers 'gamma', 'x0', 'c1', 'c2', 'c3' and 'c4' allow the user to customize the adopted extinction curve. For example, see Clayton et al. (2003, ApJ, 588, 871) for examples of these parameters in different interstellar environments.
EXAMPLE Determine how a flat spectrum (in wavelength) between 1200 A and 6500 A is altered by a reddening of E(B-V) = 0.5. Assume an "average" reddening for the diffuse interstellar medium (R(V) = 3.1)
isis> wave = [1200:6500:#500]; %Create a wavelength vector from 1200 to 6500 Angstrom with 500 steps between isis> flux = 1.*ones(500); %Create a "flat" flux vector isis> ebv = 0.5; %Value for E(B-V) isis> variable var; %Referenz for extcurve isis> funred = fm_unred(wave, flux, ebv; extcurve=&var); %Redden flux vector isis> plot(wave,var); %Plots the extinctioncurve versus the wavelength
See also: e_bv;
foucalc
Synopsis
calculates power-, cross power-, coherence- and timelag-spectra for timing analysis
Usage
Struct_Type foucalc(Struct_Type lc, Integer_Type dimseg)
Qualifiers
- verbose
- normtype: normalization type of PSD data, can be
"Miyamoto"
[default],"Leahy"
, or"Schlittgen"
- normindiv
- avgbkg: array of average background rates for each energy band
- numinst [
=1
]: number of activated PCUs on XTE, required for noise correction - deadtime [
=1e-5
]: detector deadtime in seconds - fmin: minimum frequency used for RMS calculation
- fmax: maximum frequency used for RMS calculation
- RMS: reference to a variable to store the RMS of each light curve in the [fmin, fmax] band. Can only be used with normtype="Miyamoto".
- avgrate: reference to a variable to store the average rate of each light curve
- compact: compact the output structure, only keep the most important quantities
- noCPD: do not calculate cross power densities and derived quantities
Description
lc
contains (properly segmented) light curves in several energy bands.
The best performance is achieved with a common structure of arrays,
lc = struct { time=[t_1, t_2, ...], rate1=[r1_1, r1_2, ...], rate2=... };
.
[However, one can also use an array of structures (with an enormous overhead),
lc = [ struct { time=t_1, rate1=r1_1, rate2=...}, struct { time=t_2, rate1=r1_2, rate2=... }, ... ];
.]
dimseg
is the segmentsize used for the FFTs, whihch should therefore be a power of 2.
The returned structure contains the following fields:
-
freq
: the frequency grid -
numavgall
: The bini
in all power spectra has been averaged overnumavgall[i]
original bins. Here,numavgall[i] = numseg
, as no frequency rebinning has been performed. -
Power spectra for each energy band:
* rawpsd
, errpsd = rawpsd/sqrt(numseg)
, noipsd
:
the raw power spectrum (from makepsd
),
its error from the average over segments,
and the noise level (from psdcorr_zhang
)
* sigpsd = rawpsd - noipsd
:
the signal power spectrum
* rawnormpsd
, errnormpsd
, noinormpsd
, effnoinormpsd = noinormpsd / sqrt(numseg)
:
the normalized power spectrum,
its error, the noise level,
and the effective noise level in the normalized power spectrum
* signormpsd = rawnormpsd - noinormpsd
:
the signal in the normalized power spectrum
- Cross power spectra, coherence and time lag functions for each pair of energy bands:
* realcpd
, imagcpd
:
real and imaginary part of the cross power density
* noicpd = ( sigpsd_lo * noipsd_hi + sigpsd_hi * noipsd_lo + noipsd_lo * noipsd_hi ) / numseg
* rawcof
, cof
, errcof
:
non-noise-corrected (raw) and noise-corrected coherence function and its one-sigma uncertainty
* lag
, errlag
:
time lag and its one-sigma uncertainty
foufreq
Synopsis
Timing Tools: Calculation of the Fourier Frequency Array
Usage
Double_Type freq[] = foufreq(Double_Type time[]);
Description
Calculates the Fourier frequency array corresponding to a given equally-binned time array.
getCommonGTIs
Synopsis
returns a new set of GTIs which are both contained in gti1 and gti2
Usage
Struct_Type gti = getCommonGTIs(Struct_Type gti1, Struct_Type gti2);
Description
Good time intervals (GTIs) are stored in { start, stop }
-structs.
getVarPeriod
Synopsis
Determines the period of a variation in a lightcurve
Usage
Double_Type getVarPeriod(Struct_Type lc, Double_Type minperiod, Double_Type maxperiod)
Qualifiers
plot - plot chi-square distribution and modelled gaussians. Pause the assigned number in seconds between each step (default 0.3) nogauss - skip the gaussian fit
piped to epfold dt - exposure of every lightcurve time bin, should be given to ensure correct results. sampling - how many periods per peak to use (default=10) nsrch - how many periods to search in a linear grid (default not set) dp - delta period of linear period grid (default not set) lstat - use L-statistics instead of chi^2 statistics nbins - number of bins for the pulse profile exact - calculate the pulse profile in a more exact way, see description of pfold (not recommed as it takes a very long time!).
Description
**** This function is deprecated!! It will be removed in 2017 Mai. If you heavily depend on this function, please write an email to matthias.kuehnel@sternwarte.uni-erlangen.de **** Performs an Epoch Folding of the lightcurve in the given period range (same time unit as used in the lightcurve) and returns the period corresponding to the maximum. If not skipped by the 'nogauss' qualifier, a gaussian it fitted to the maximum to get its position more accurate. The gaussian parameters and used data- points are applied automatically and improved step by step. The steps can be monitored using the 'plot' qualifier.
See also: epfold
get_all_data
Synopsis
Get a list of all data-set indices
Usage
List_Type = get_all_data;
Description
This function returns a list of data indices, which have been obtained via all_data. The information is saved in a List, not an Array, which allows plotting of all datasets.
See Example
EXAMPLE
isis>variable d = get_all_data; isis>plot_data(d);
See also: all_data
get_bin_corr_factor
Synopsis
reads the bin correction factor set by set_bin_corr_factor
Usage
Double_Type[nbins] corr_factor = set_bin_corr_factor(Integer_Type data_id)
See also: load_fermi,set_bin_corr_factor
get_combined_data_model_residuals
Synopsis
returns these 3 structures for a combination of data sets
Usage
(d, m, r) = get_combined_data_model_residuals(Integer_Type id[]);
Description
All d
, m
and r
are bin_lo, bin_hi, value, err
structures.
d
and m
are the sum of data and model counts of all data sets,
rebinned to the grid of the first data set id[0].
d.err
is calculated from quadratic error propagation.
r.value = (d.value-m.value)/d.err
contains the residuals
unless the ratio
qualifier is set (see below).
Qualifiers
- ratio: use ratio-residuals
r.value = d.value/m.value
See also: get_data_counts, get_model_counts, rebin
get_count_rate
Synopsis
Get background-subtracted count rates
Usage
Double_Type = get_count_rate (hist_index);
##### Description
Use this function to retrive background subtracted
countrates (default), countrates or background rates
for spectra that have been loaded.
WARNING: This implementation for faked
datasets is experimental!!
##### Qualifiers
tex: additional TeXoutput of CRs, choose this if you need
errors
bkg: 0 (default: background subtracted CR), 1 (count rates),
2 (background count rate)
fake: Set this qualifier if you have a faked spectrum and a
faked background (which has not been deleted!);
Assumption: hist_index(bkg) = hist_index(dset) +1 !!
EXAMPLE
isis>xray = load_data("data.pha"); isis>variable cr = get_count_rate (1 ;tex, bkg=1);
See also: get_data_counts; TeX_value_pm_error;
get_data_counts_with_tot_err
Synopsis
returns spectral data, taking systematic errors into account
Usage
Struct_Type data = get_data_counts_with_tot_err(Integer_Type id);
Description
data.err = sqrt( stat_err^2 + [sys_err_frac * data.value]^2 );
See also: get_data_counts, get_sys_err_frac
get_flux_corrected_convolved_model_flux
Synopsis
computes the model flux which can be compared with a flux-corrected spectrum
Usage
Struct_Type get_flux_corrected_convolved_model_flux(Integer_Type hist_index)
Description
The flux-corrected convolved model flux is computed as
int K(R(h,E), A(E), S(E)) dE
F(h) = --------------------------------
,
int K(R(h,E), A(E), 1 ) dE
where S(E)
is the flux model, A(E)
is the effective area (the ARF),
R(h,E)
is the redistribution function (the RMF),
and the kernel K
defaults to K(R,A,S) = R*A*S
.
The flux-corrected convolved model fluxes is defined in such a way
that a folded model is "unfolded" in the same way as the data are
by virtue of flux_corr
. It can thus be compared with flux-corrected data.
See also: flux_corr, get_model_counts, get_model_flux, get_convolved_model_flux
get_modelcomponents
Synopsis
load binned model component values
Usage
Struct_Type[] mc = get_modelcomponent( Integer_Type hist_index);
or
```c
Struct_Type[] mc = get_modelcomponent( );
##### Qualifiers
* fct [=NULL]: NULL: Using 'get_model' (default)
"counts": Using 'get_model_counts'
"flux": Using 'get_model_flux'
##### Description
For each model component (see 'get_fun_components') within the given
dataset 'hist_index' this funtion returns a stucture with three array
fields, bin_lo [Angstrom], bin_hi [Angstrom] and value, containging
only the contribution of this component to the overall model.
If 'hist_index' isn ot given all included datasets are used!
The model contribution of each component is obtained by setting all
the norm parameters (get_par_info: is_a_norm=1) of all other components
to zero and evaluating the model. NOTE THAT for each component an
'eval_counts' is called!
The 'fct' qualifier determines which function should be used to
obtain the model components.
__See also__: get_fun_components, get_model, get_model_counts, get_model_flux
#### get_ratio
##### Synopsis
calculates data/model ratios
##### Usage
```c
Struct_Type rat[] = get_ratio(Integer_Type id[]);
Description
Each rat[i]
is a { bin_lo, bin_hi, value, err }
structure
where rat[i].value
contains the (data counts)/(model counts)
ratios obtained for the data set id[i]
.
If only a scalar id
is given, only a single structure is returned.
See also: get_data_counts, get_model_counts, get_residuals
get_residuals
Synopsis
calculates (data-model)/error residuals
Usage
Struct_Type res = get_residuals(Integer_Type id[]);
Description
Each res[i]
is a { bin_lo, bin_hi, value, err }
structure
where res[i].value
contains the (data counts - model counts)/error
residuals obtained for the data set id[i]
.
If only a scalar id
is given, only a single structure is returned.
Qualifiers
- noticed: restrict to noticed bins
- keV: convert Angstrom-bins to keV-bins
See also: get_data_counts, get_model_counts, get_ratio
get_selected_data_flux_en
Synopsis
returns the flux of a data set in energy units
Usage
flux = get_selected_data_flux_en(id[, Emin, Emax[, alpha]]);
Description
flux
is a { bin_lo, bin_hi, value, err}
structure
containing energy bins of data set id
and the spectral photon flux density (in 1/s/cm^2/keV, unless alpha!=0
).
If Emin
and Emax
are used, only values in this energy range are considered.
If alpha
!=0, the flux values are multiplied with E^alpha
.
Example
hplot( get_selected_data_flux_en(1, 4, 20, 1) );
% plots the spectral energy flux density of data set 1 in the 4--20 keV range
get_warmabs_model
Synopsis
computes the contribution of separate elements in a warmabs-model
Usage
Struct_Type data = get_warmabs_model();
group_noticed_data
Synopsis
groups previously noticed spectral bins by an integer factor
Usage
group_noticed_data(Integer_Type id, Integer_Type factor);
See also: group_data
ignore_large_bins
Synopsis
ignores those bins of a spectral dataset exceeding a maximal size
Usage
ignore_large_bis(Integer_Type id[], Double_Type maxsize);
Qualifiers
- verbose
- unit[=
"A"
]: maxsize may be in A or keV
lc_from_events
Synopsis
extracts light curves from an event list
Usage
Struct_Type dat = lc_from_events(Struct_Type evts [, Array_Type bands]);
Qualifiers
- dt: time resolution in sec [default: 100]
- back = Struct_Type b_evts : substract background events)
- gti = struct { Double_Type start, stop } : GTIs for events
- minfracexp: minimum allowed fractional exposure of a time bin, requires GTIs to be set
Description
This function extracts light curves from an event list given as the
following structure:
Double_Type[] time - the event arrival times (in seconds)
Double_Type[] pi - the associated event energies (in eV)
If bands
is given, an array of structures containing light curves of each band
is returned. The energy bands have to be given in eV. To receive
two lightcurves with energies between 1-2keV and 2-3keV, e.g.,
If GTIs are provided, the fractional exposure 'fracexp' will be added to the output structure.
dat = lc_from_events(evts, [1000,2000,3000]);
See also: color_color_data, histogram
loadDataset
Synopsis
loads and assigns one spectrum and its detector response
Usage
Integer_Type id = loadDataset(phaFile[, rmfFile[, arfFile[, backFile[, row]]]]);
Description
loadHETGlc_sum
Synopsis
ADD Chandra HETG lightcurves
Usage
Struct_Type lcsum = loadHETGlcs( [lcs] );
Qualifiers
- dir [="."] : path to the extracted spectra
- heg [=1] : load HEG spectra (=0: don't load)
- meg [=1] : load MEG spectra (=0: don't load)
- order [=[-1,1]]: diffraction orders to load. Sign matters.
- energyband [=[500,10000]]: extracted energy bands. See description below for details.
- unit [="eV"] : Unit for the energyband: eV, keV, A
Description
Sum multiple Chandra HETG lightcurves extracted with the (new) Remeis extraction scripts and return a structure with the result.
The new (summer 2014) Remeis Chandra gratings extraction scripts do not longer specify observation / extraction specific information in the file name, but only the arm, diffraction order, and energyband of the lightcurve, e.g., heg_m1_500-1000.lc or meg_p1_3500-10000.lc, and only extract lightcurves on a per arm/order basis.
Often, we are more interested in the combined lightcurve of multiple arms / orders. loadHETGlc_sum can be used in 2 ways:
- call the function with no arguments and use the qualifiers to specify which lightcurves you are interested in. The qualifiers are passed to loadHETGlcs to load the lightcurves first.
- call the function with 1 argument: the lightcurves produced by a previous call of loadHETGlcs. In this case, all lightcurves contained in the structure will be summed, and the qualifiers will have no effect [in the future, they could be used to choose a subset of the given lightcurves, but this is not implemented yet].
The returned structure has the same fields as simply reading a single lightcurve file would produce. If reading of the lightcurves fails, the return value is -1.
The function assumes that the lightcurves have not been tempered with, i.e., that they are correctly (Chandra) formatted lightcurves and have the same length, time grid, etc. No checks are performed. Fields like time and exposure are taken from the first lightcurve in the structure, since they are expected to be the same for various extractions of the same ObsID.
The energyband can be specified in multiple ways: * a list or array of limits: [0.5,1.5,3,10] keV. In this case, the energybands for the lightcurves are assumed to be 0.5-1.5, 1.5-3, and 3-10 keV, i.e., there are no gaps. * a list or an array of energy pairs: {[0.5,1.5],[3,10]} or [{0.5,1.5},{3,10}] keV. Then the bands are taken as they are, i.e., it is possible to have gaps between energy bands.
See also: loadHETGlcs, loadHETGspec
loadHETGlcs
Synopsis
Load Chandra HETG lightcurves
Usage
lcs = loadHETGlcs();
Qualifiers
- dir [="."] : path to the extracted spectra
- heg [=1] : load HEG spectra (=0: don't load)
- meg [=1] : load MEG spectra (=0: don't load)
- order [=[-1,1]]: diffraction orders to load. Sign matters.
- energyband [=[500,10000]]: extracted energy bands. See description below for details.
- unit [="eV"] : Unit for the energyband: eV, keV, A
Description
Load Chandra HETG lightcurves extracted with the (new) Remeis extraction scripts and return a structure with the results.
The new (summer 2014) Remeis Chandra gratings extraction scripts do not longer specify observation / extraction specific information in the file name, but only the arm, diffraction order, and energyband of the lightcurve, e.g., heg_m1_500-1000.lc or meg_p1_3500-10000.lc.
loadHETGlcs reads the corresponding lightcurve for each specified arm, diffraction order, and enegeryband and stores it as a field in a structure. The returned structure is then in the form (in case of the defaults): struct{ dir = ".", orders = [-1,1], arms = ["heg","meg"], energyband = [500, 10000], unit = "eV", heg_m1_500_10000 = Struct_Type, % content of the fits file heg_p1_500_10000 = Struct_Type, meg_m1_500_10000 = Struct_Type, ... } The function returns -1 if both arms are ignored, i.e., no lightcurves chosen.
The energyband can be specified in multiple ways: * a list or array of limits: [0.5,1.5,3,10] keV. In this case, the energybands for the lightcurves are assumed to be 0.5-1.5, 1.5-3, and 3-10 keV, i.e., there are no gaps. * a list or an array of energy pairs: {[0.5,1.5],[3,10]} or [{0.5,1.5},{3,10}] keV. Then the bands are taken as they are, i.e., it is possible to have gaps between energy bands.
See also: loadHETGlc_sum, loadHETGspec
loadHETGspec
Synopsis
Load locally extracted Chandra HETG spectra and responses
Usage
Struct_Type ids = loadHETGspec();
Qualifiers
- dir [="."] : path to the extracted spectra
- heg [=1] : load HEG spectra (=0: don't load)
- meg [=1] : load MEG spectra (=0: don't load)
- order [=[-1,1]]: diffraction orders to load. Sign matters.
Description
Load Chandra HETG spectra extracted with the (new) Remeis extraction scripts and return a structure with the dataset indices.
The new (summer 2014) Remeis Chandra gratings extraction scripts do not longer specify observation / extraction specific information in the file name, but only the arm and diffraction order of the spectrum, e.g., heg_m1.pha or meg_p1.pha. The old extraction scripts used to generate a loaddata0.sl file for convenience. With the more generic file names, this function replaces the load script. To ensure backwards compatibility with dataset-index sensitive par files, the function loads the arms and orders in the exact same order as the old load scripts used to: heg_m1 heg_p1 meg_m1 meg_p1 heg_m2 heg_p2 meg_m2 meg_p2 ...
The returned structure is then in the form (in case of the defaults): struct{ dir = ".", orders = [-1,1], arms = ["heg","meg"], heg_m1 = 1, heg_p1 = 2, meg_m1 = 3, ... } The function returns -1 if both arms are ignored, i.e., no spectra chosen.
Negative and positive diffraction orders can be chosen independently of each other. However, all selected orders are applied to all selected grating arms: You can load HEG without MEG and vice versa. But if you load them both in a single call of loadHETGspec, the same orders will be loaded for both of them.
If you would like to load configurations like [heg_m1, meg_p1], where different orders are loaded for each grating arms, you have to call the function multiple times or load the spectra by hand.
See also: loadHETGlcs, loadHETGlc_sum
load_1storderHETGS_datasets
Synopsis
loads MEG+-1 & HEG+-1 Chandra HETGS spectra
Usage
Integer_Type ids[] = load_1storderHETGS_datasets(String_Type specpath, RMFpath);
load_HETGS_datasets
Synopsis
loads Chandra HETGS spectra from a type II pha file
Usage
Integer_Type ids[] = load_HETGS_datasets(String_Type specpath, RMFpath[, Integer_Type ms[]]);
Description
ids
is an array containing negative/positive MEG/HEG spectra
load_atime
Synopsis
loads a structure containing arrival times from a FITS-file
Usage
Struct_Type load_atime(String_Type filename, String_Type extname);
Qualifiers
none
Description
The arrival times, which were stored previously by 'save_atime', are loaded and returned as a structure. This structure has the fields described in the 'atime_det' function. The extension table containing the arrival times must be specified as well as the FITS-file itself.
See also: save_atime, atime_det
load_data_combined
Synopsis
reads and combines the given rows from a Fits Type II pha file
Usage
Integer_Type load_data_combined(String_Type, pha_filename, Integer_Type[] rows);
Qualifiers
are passed to load_data
Description
A combination of load_data and combine_datasets, with the exception that the rows appear as a single dataset.
See also: load_data, combine_datasets
load_data_integral
Synopsis
reads an INTEGRAL spectrum from a PHA file and performs the necessary tweaks
Usage
Integer_Type id = load_data_integral(String_Type pha_filename);
Description
ISIS's load_data often does not work with INTEGRAL spectra for the following reasons:
-
The extension of the PHA file is not called 'SPECTRUM'.
-
The extensions of the RMF file are not called 'MATRIX' and 'EBOUNDS'.
-
The extension of the ARF file is not called 'SPECRESP'.
-
The columns ENERG_LO and ENERG_HI of the ARF file contain only zeros.
The function load_data_integral tries to make the appropriate changes (modify the extension names, write a new ARF file) and finally uses load_data.
Qualifiers
- verbose [=1]: show changes
See also: load_data
load_fermi
Synopsis
loads fermi data produced by the "Fermi-SED Script"@remeis
Usage
Integer_Type data_id = load_fermi(String_Type filename);
Qualifiers
- rmf_factor[=10]: define a fine RMF to properly evaluate the model
Description
Loads fermi data produced by the "Fermi-SED Script"@remeis. Note that special care is taken here, as the energy bins are relatively large. Therefore the function set_bin_corr_factor is used to calculate a correction factor. This factor is used automatically for loading the data. In order to plot the data properly, plot_data/unfold has to be used. WARNING: Use clear_all to delete the correction factors (as well as the data). This is necessary as otherwise the correction factors are applied to a newly loaded dataset with the same ID.
See also: load_fermi_catalog
load_fermi_catalog
Synopsis
loads Fermi data from a catalog for a given source
Usage
Integer_Type data_id = load_fermi_catalog(String_Type sourcename);
WARNING: Use clear_all to delete the correction factors (as well
as the data). This is necessary as otherwise the correction factors
are applied to a newly loaded dataset with the same ID.
Qualifiers
- list: return list of available sources
- ul: Load includes upper limits
- catalog: specify the catalog (default: 2FGL_point_source@Remeis)
- fgl3: Use this when reading in the 3FGL catalog, either with the catalog qualifier or without
See also: load_fermi
load_piled1storderHETGS_datasets
Synopsis
loads 1st order Chandra-HETGS data and sets up the simple_gpile model
Usage
Integer_Type ids[] = load_1storderHETGS_dataset(specpath, RMFpath);
load_piledHETGS_datasets
Synopsis
loads Chandra HETGS spectral and sets up the simple_gpile model
Usage
Integer_Type ids[] = load_1storderHETGS_dataset(specpath, RMFpath);
See also: load_HETGS_datasets, use_simple_gpile
load_radio2
Synopsis
reads and loads radio data
Usage
Integer_Type data_id = load_radio (String_Type filename);
Description
Takes data an ASCII file , comprised of three columns: Frequency[Hz], Flux density [mJy], Error [mJy], and loads the data as a dataset. WARNING: The standard bin width is set to 10%
It can be changed with the qualifier binwidth (given in fraction of the frequency). Since a log frequency grid is used, it will not be exactly the size given. WARNING: Bins cannot overlap, as grid needs to be continuous!
Qualifiers
- binwidth: specify bin width
- use_struct: use an ISIS structure with the following fields: freq[Hz], bandwidth[Hz], flux[mJy], err[mJy]; with full bandwidth given
Example
isis> variable data_file = "/some/path/radio_file.txt"; isis> variable radio_data = load_radio2(data_file);
example using a textfile that includes bin widths:
isis> variable abc = ascii_read_table ("radio_dat.txt",[{"%F","freq"},{"%F", "bandwidth"},{"%F","flux"},{"%F","err"}]); isis> variable radio_data = load_radio2(abc; use_struct);
See also: read_radio, load_radio
logging
Synopsis
write a log file
Usage
logging ("filename.log", "Log message");
Description
This function writes output to a log file (which will be created if it does not exist. It can also handle input in the sprintf format (see example). By default all input to logging is also printed to stdout. Repeated calls to the function add the new message at the end of the file. It is therefore recommended to specify the date and time in the name of the log file.
Qualifiers
- v: verbosity: [default: 1] print output also to stdout, set this qualifier to 0 for no output to stdout
Example
isis> variable cutime = strftime("%Y-%m-%d_%H:%M:%S"); isis> variable l = cutime+".log"; isis> logging(l, "* Found %.u observations for %s", length(src_info), src); isis> logging(l, "***Error in function %s", _function_name() );
See also: strftime, _function_name()
luminosity
Synopsis
computes a source luminosity assuming the current fit-function and a distance
Usage
Double_Type luminosity(Double_Type Emin, Emax, d);
Qualifiers
- factor [=
1.001
]: step of logarithmic energy grid - Emin: minimum energy of (extended) grid
- Emax: maximum energy of (extended) grid
Description
returns 4pi (d kpc)^2 * int_Emin^Emax E*S_E(E) dE
(in erg/s)
Use flux2lum to calculate luminosities for extragalactic sources.
See also: energyflux, flux2lum
make_fine_rmf
Synopsis
creates a fine RMF
Usage
make_fine_rmf (Data_Id [Integer], Factor [Integer]);
Description
Creates a fine RMF. Therefore each data bin will be split into a certain number of bins, specified by the variable "Factor".
See also: load_slang_rmf
makepsd
Synopsis
Calculate the power density spectrum of a single energy band.
Usage
(psd, nipsd) = makepsd(rate, timeseg, dimseg);
Description
Input: rate - input rate array timeseg - real time of one segment of rate dimseg - number of bins in each segment of rate
Output: psd - unnormalized, averaged PSD nipsd - averaged PSD, individual segments normalized
Qualifiers
- avgbkg: Average background count rate for Miyamoto normalization
- normtype: "Miyamoto", "Leahy", "Schlittgen" normalization type
normpsd
Usage
psd_norm = normpsd(psd, normtype, avrate, avbackrate, timeseg, dimseg);
Description
psd
is the unnormalized power spectrum.
avrate
is the average countrate of corresponding lightcurve.
avbackrate
is the average background countrate of corresponding lightcurve.
timeseg
is the real time of one corresponding lightcurve segment.
dimseg
is the number of bins in one lightcurve segment.
String_Type normtype
can be:
- "Miaymoto" (default) [Miyamoto et al. (1991), ApJ 383, 784]
psd *= 2 * timeseg / (dimseg * [avrate-avbackrate]^2)
- "Leahy" [Leahy et al. (1983), ApJ 266, 160]
psd *= 2 * timeseg / (dimseg^2 * avrate)
- "Schlittgen" [Schlittgen, H.J., Streitberg, B. (1995), Zeitreihenanalyse, R. Oldenbourg]
psd /= dimseg
notice_human2isis
Synopsis
create isis notice_list from str = notice_isis2human
Usage
Array_Type notice_list = notice_human2isis(String_Type str);
Description
This routine converts a notice_list string created with notice_isis2human back to the ISIS conform binning array (see e.g. get_data_info(idx).notice_list). It can be used to restore the otcing of data directly, by using it as argument in the function notice_list(...).
See also: fits_save_fit, get_data_info, notice_isis2human
notice_isis2human
Synopsis
converts isis notice_list in a human readable string
Usage
String_Type str = notice_isis2human(Array_Type notice_list);
Description
This routine converts the notice_list as used by ISIS (see e.g. get_data_info(idx).notice_list) to a string which is easy to read.
The routine fits_save_fit() uses this string to save the noticed bins of the data. It can be converted back by using notice_human2isis().
See also: fits_save_fit, get_data_info, notice_human2isis
notice_listFromString
Synopsis
creates a isis notice_list from a string containing the noticed bins separated with ",".
Usage
Array_Type notice = notice_listFromString(String_Type str);
See also: notice_list, notice_listToString, get_data_info, fits_save_fit
notice_listToString
Synopsis
converts an isis notice_list to a string
Usage
String_Type str = notice_listToString(Array_Type notice_list);
See also: notice_list, notice_listFromString, get_data_info, fits_save_fit
nuFnu
Synopsis
changes an Energy/Flux Spectrum to a Frequency/Flux*Frequency Spectrum
Usage
Struct_Type = nuFnu (hist_index, E_min(keV), E_max(keV));
Description
Use this function on a dataset that has been modeled. This function uses get_data_flux/eval_fun_keV to load spectral data into isis and convert energy bins in a given range to frequency (s^-1) bins as well as Flux to flux*frequency. If no range specified: E_min=0.5 keV, E_max=10 keV. Values and Errors are given in erg s^-1 cm^-2.
WARNING: Currently only frequency has been implemented as x-unit.
For the evaluation of a deabsorbed component use the qualifier "deabs". deabs will try to calculate the factor between absorbed and deabsorbed model components and correct the data.
WARNING: This does not work when a convolution normalization (e.g., cflux, enflux, phflux) is part of the model!
The "numbin" qualifier allows to specify the number of output bins. Currently only a logarithmic grid has been implemented. If numbin is not specified the flux for the current grid will be returned.
The "group" qualifier allows to specify the S/N that will be used for the flux calculation. Flux calculations are only correct for large number of bins. The final grid can be smaller (qualifier: numbin).
WARNING: Currently tested only on low counts spectra. S/N might have to be increased for large number of counts.
Qualifiers
- deabs: evaluate deabsorbed flux
- ff: define deabsorbed fit function
- numbin: define number of bins
- group: define S/N for grouping
Example
isis>xray = load_data("data.pha"); isis>load_par("x.par"); isis>() = fit_counts; isis>new=nuFnu(xray,2,10); %change spectrum (2-10 keV) into frequency spectrum
%Plot data: isis>hplot(new.lo, new.hi, new.val);
% Evaluate Function for deabsorbed model % (WARNING: DOES NOT WORK IF ENFLUX IS USED IN FIT FUNCTION)
isis>xray = load_data("data.pha"); isis>fit_fun("tbnew_simple_z(1)*powerlaw(1)"); isis>() = fit_counts; isis>new = nuFnu(xray,2,10; deabs, ff="powerlaw(1)"); isis>%Calculate a deabsorbed powerlaw
See also: get_data_flux, rebin_mean, eval_fun_keV, group
num_bin
Synopsis
Number of noticed bins
Usage
Double_Type = num_bin (hist_index);
##### Description
Use this function to retrieve number of
noticed bins.
##### Qualifiers
dof: gives degrees of freedom instead of number of
bins
EXAMPLE
isis>xray = load_data("data.pha"); isis>variable num = num_bin(1);
See also: dof
oplotGTIs
Synopsis
visualizes Good Time Intervals
Usage
oplotGTIs(Struct_Type gti[, Double_Type offset]);
partial_correlation
Synopsis
Tests two luminosities for partial correlation due to redshift
Usage
partial_correlations(String_Type filename);
Description
This function tests if a correlation of two parameters is due to the redshift. Description of method and code adapted from Akritas & Siebert, 1996, MNRAS, 278, 919. Works with FITS and ASCII files.
The input file needs 6 columns with luminosity1, UL, luminosity2, UL, redshift, UL. The upper limit columns should have a 1 for a detection, and a 0 for an upper limit.
Examples
partial_correlation("data.txt")
partial_correlation("data.fits")
See also:
pfold
Synopsis
folds a lightcurve on a given period
Usage
(Struct_Type pp) = pfold(Double_Type t, r, p);
or (Struct_Type pp) = pfold(Double_Type t, r, p, e);
or (Struct_Type pp) = pfold(Double_Type t, p) ; (event data)
Qualifiers
- nbins: number of bins for the pulse profile
- exact: calculate the pulse profile in a more exact way, see description.
- returnlc: returns the LC and the corresponding phases and phase bins for every LC bin. Used mainly for debugging.
- dt:exposure of every lightcurve time bin, should be given to ensure correct results.
- t0: set reference time.
- pdot:first derivative of pulse period p (default=0).
- pddot:second derivative of pulse period p (default=0).
- gti:GTIs for event data, given as struct{start=Double_Type, stop=Double_Type}
Description
Calculates the pulse profile of a given lightcurve, given as time t and rate r for the period p. The function was adopted from the IDL routine of the same name. The GTI filtering is only implemented for eventdata.
The data arrays time (t), rate (r), and error (r), can also be given as qualifier (t, r, e, respectively), which allows pfold to be run with "array_map_qualifiers" on an array of periods (p).
Compared to the similar function sitar_pfold_rate, pfold.sl calculates the errors pp.error (if given as a forth argument) with correct error propagtion, if e is given. When no errors are given, the variance in each bin is used as in sitar_pfold_rate. If the qualifier "dt" is given, the mean count rate in every phase bin is calculated by weighing the countrates by the exposure time, unlike sitar_pfold_rate.
If the "exact" qualifier is given, the function takes the exposure time of every time bin into account in the sense that, that a given time bin may overlap over two phase bins. The corresponding exposure time in every phase bin is reduced accordingly.
NOTE: the "exact" qualifiers slows the script down, depending on the length of the lightcurve and the number of bins for up to a factor of >10.
If "exact" is not given, the script is ~10% slower than sitar_pfold_rate.
If p is a function of time described by a taylor expansion (p(t) = p0 + pdot/2*t + pddot/6*t^2 + ...), the first 2 derivatives "pdot" and "pddot" can be set as qualifieres.
%
phasebin_GTIs
Synopsis
computes the GTI time in phase bins
Usage
Double_Type dt[] = phasebin_GTIs(tstart, tstop, T0, P, n);
photon_flux
Synopsis
Usage
String_Type = photon_flux (hist_index, E_min, E_max);
Description
This function returns the integrated photon flux of a defined Energy range. If no energy range specified E_min=0.5, E_max=10.
See also: get_data_flux;
plotGTI
Synopsis
visualizes Good Time Intervals
Usage
plotGTI(Struct_Type gti[, Double_Type offset]);
Qualifiers
- overplot
plot_atime
Synopsis
plots the residuals of a pulse arrival times fit
Usage
plot_atime([Integer_Type[] index]);
Qualifiers
t0 - reference time of xticks col - array of colors for datasets noerrbars - do not plot errorbars connect_points - see 'connect_points' style - see 'point_style' xunit - used time unit (d or s) yunit - used residual unit (phi or s) xlunit - label of xunit ylunit - label of yunit xrng - time range yrng - residual range (default -0.5 to 0.5) ploteph - overplot the pulse ephemeris plotorb - overplot the orbital motion plotpnum - overplot the modelled pulse numbers ('modnum' reference must be set, see 'define_atime') plotdif - overplot the meassured pulse period ('modnum' reference must be set, see 'define_atime') y2rng - yrange of overplot y2rel - if pulse period tics are displayed wrong, you should switch to relative values and milliseconds extcol - color of overplot extcolY - yaxis color of overplot oplot - do not erase plot window lshift - if 'peph' and 'porb' set, shift the labels printed at the ephemeris by the given value in y-direction. May be a two value array [eph,orb]. Values are interpreted as relative coordinates.
Description
By default plots the residuals of all datasets containing arrival times. The residuals are calculated by residuals = data - model Using the optional first parameter the indices of the datasets to be plotted can be specified. Also works if some datasets are excluded from the fit (see 'atime_xinclude'). The residuals in units of pulse phase or seconds are plotted against the meassured pulse arrival times in days or seconds. The units can be set by using the accordant qualifiers. Further informations can be overplotted, like the pulse period in the binary barycentre.
See also: arrtimes, atime_dataind, atime_xinclude, plot_data
plot_combined_data_model_residuals
Synopsis
plots data and model counts, and residuals for a combination of data sets
Usage
plot_combined_data_model_residuals([Integer_Type id[]]);
Qualifiers
- dcol: color of data [default=1]
- mcol: color of model [default=2]
Description
If no indices id
of data sets are specified, all_data
are used.
The combined data, model and residuals are plotted in a two panel multiplot.
See also: get_combined_data_model_residuals
psdcorr_zhang
Synopsis
Timing Tools: Poisson Noise and Deadtime Correction (Zhang)
Usage
(freq, psd) = psdcorr_zhang(totrate, tseg, dimseg);
Description
Inputs: totrate - total countrate of all instruments tseg - realtime length of the lc segments used for psd calculation dimseg - bincount of the lc segments used for psd calculation
Outputs: freq - fourier frequency array noipsd - array of observational noise of the psd
See also: Zhang et al. (1995), ApJ, 449, 930
pulsar_GTI
Synopsis
Calculates GTIs for pulse phase resolved spectroscopy
Usage
pulsar_GTI(Double_Type tstart, tstop, String_Type satellite, basefilename, Double_Type t0, p, phase_lo, phase_hi);
Qualifiers
- pdot: first derivative of the pulse period in s/s (default: 0)
- pddot: second derivative of the pulse period in s/s^2 (default: 0)
- MJD: If set, tstart, tstop and t0 are assumed to in MJD
- local: create GTIs in satellite's local time system. Requires qualifiers 'nobarevt' and 'barevt' to be set.
- barevt: event file in barycentered time frame
- nobarevt: event file in local time frame
- all other qualifiers are passed to BinaryPos
Description
This function calculates the GTIs for pulse phase resolved spectroscopy. Input arguments are
tstart - start time of the time interval to be covered (typically the start of the observation) tstart - end time of the time interval to be covered (typically the end of the observation) satellite - the name of the satellite (needed for setting the reference MJD correctly) basefilename - the file created will be named 'basefilename.gti' t0 - the reference time (time of phase zero) p - the pulse period phase_lo - the lower phase boundary phase_hi - the upper phase boundary
Upon qualifier request, tstart, tstop, and t0 can be given in MJD instead of seconds. The pulse period is always in seconds.
If the pulse period and the reference time are corrected for the binary
motion, the orbital parameters should be provided by qualifiers such that
the GTIs can be transformed into the observers time system. The
qualifiers are passed and equal to the BinaryPos
function.
If GTIs in the satellite's local reference time system are needed (e.g., for Suzaku-XIS), an event file in both local and barycenterd time can be passed via qualifiers and the GTIs are interpolated to the local time system. Header keywords are set accordingly.
Example
pulsar_GTI (2.7082e8, 2.7097e8, "nustar", "phase_0.25-0.35", 2.708242e8, 443.07, 0.25, 0.35);
creates a GTI file named 'phase_0.25-0.35.gti', ranging from 2.7082e8--2.7097e8 seconds in NuSTAR's mission specific time system for the pulsar 4U 1907+09 with pulse period 443.07 seconds for the pulse phase interval 0.25--0.35 where phase 0.0 is set to be at 2.708242e8 seconds.
See also: MJDref_satellite, BinaryPos, pulse_time
pulse2pulse_flux_lc
Synopsis
averages the given lightcurve over the pulsar's period
Usage
Struct_Type = pulse2pulse_flux_lc(
Struct_Type lightcurve,
Double_Type or Struct_Type period,
[, Struct_Type orbit]
);
Qualifiers
remap - interpolate the resulting flux lightcurve on the input time grid (default: no) interpol - function reference to perform the time interpolation (default: &interpol_points) dphitol - minimum phase coverage (dphi) a pulse in the lc has to have at least (default: .95) gaptol - minimum time difference between bins which defines a gap (default: 2*min(diff(lc.time))) t0 - time of pulse phase zero (can be provided using the ephemeris structure as well; default: first time bin of the lightcurve) phase - array of the pulse phases corresponding to the lightcurve time array (default: its calculated using 'pulseperiod2phase'). If the qualifier is set to a reference the calculated phases are assigned to the given variable.
Description
The underlying slope in a lightcurve of a pulsar might affect any timing analysis of its pulsations due to pulse to pulse variations of the luminosity. This functions averages the count rates over each pulse to get this luminosity dependance, which is afterwards interpolated to the binning of the input lightcurve.
The input lightcurve must be of struct { Double_Type[] time, rate, error }
The input pulse period may be either a single number implying a constant pulse period or a structure containing the pulse ephemeris as defined in 'check_pulseperiod_orbit_struct'. In addition, the orbital parameters might be also given as a structure to take also the Doppler shift of the pulse period into account. This is necessary only in case of a non binary corrected lightcurve.
The output lightcurve is a struct { Double_Type[] time (center of pulse), rate, error, dphi (phase coverage of each pulse, <1 for a gap) }
By default, bad sampled pulses are ignored in the final lightcurve, which would lead to wrong mean count rates otherwise.
Note, that all parameters must have the same time unit, if applicable.
Note further, that the lightcurve's time grid has to be evenly spaced!
See also: pulseperiod2phase, interpol_points
pulse_time
Synopsis
returns the pulse arrival time of the given pulse number with respect to the pulse ephemeris and orbit
Usage
Double_Type pulse_time(Integer_Type[] number, Struct_Type ephemeris[, Struct_Type orbit]);
or Double_Type pulse_time(Integer_Type[] number, Double_Type pulseperiod[, Struct_Type orbit]);
Qualifiers
MJD - the values of the pulse ephemeris are given in days (default: seconds) eph - may be set to the pulse ephemeris structure orb - may be set to the orbital structure dphi - an additive constant phase shift
Description
Calculates the expected pulse arrival time of the given pulse number, which may be a single value or an array. The structures containing the pulse ephemeris and the orbital parameters must follow the conditions decribed in 'check_pulseperiod_orbit_struct'. Instead of the pulse ephemeris the pulse period may be given only. These structures may be also passed to the function by qualifiers, hence the structure parameters can be omitted. The equation to calculate the arrival time is similar to Hilditch eq. 3.53, including the terms of the fourth order (p3dot). The numerical calcu- lation is optimized using the Horner schema. For details see 'arrtimes'.
See also: check_pulseperiod_orbit_struct, pulseperiod, arrtimes
pulseperiod
Synopsis
calculates the pulse period at the given time depending on the given pulse period evolution
Usage
Double_Type[] pulse_period(Double_Type[] time, Struct_Type pulseperiod[, Struct_Type orbit]);
Qualifiers
interpol : reference to a function to interpolate the pulse period evolution on the requested time grid (default: &interpol_points) sameunit : set if the unit of 'time' is in the same unit as the pulse period or set to the conversion factor from days to 'time'
Description
Calculates the expected pulse period at the given time (in days), which may be a single value or an array. The structures containing the pulse period and the optional orbital parameters must follow the conditions described in 'check_pulseperiod_orbit_struct'.
If the pulse period is given as evolution (time vs. period) the period at the requested time is calculated by interpolation. Otherwise the pulse period is calculated by a 'taylor' series using a given period and its derivatives at a reference time in. In any case the requested time has to be in MJD and the period (and its Nth derivative) in seconds (/seconds^N)
Finally, if any orbital parameters are provided, the returned period gets modified by the binary motion.
See also: check_pulseperiod_orbit_struct, taylor, radial_velocity_binary
pulseperiod2phase
Synopsis
calculates the pulse phase from a given pulse period
Usage
Double_Type[] pulseperiod2phase(
Double_Type[] time,
Struct_Type pulseperiod[, Struct_Type orbit]
);
Qualifiers
interpol : interpolation method to re-map the period- onto the input time- grid (default: &interpol_points) getphi : variable reference to return the phase on the given period over time grid (method a), see below). sameunit : set if the unit of 'time' is in the same unit as the pulse period or set to the conversion factor from days to 'time'
Description
The pulse period p(t) of a pulsar and its pulse phase phi(t) are connected by dphi(t) / dt = 1. / p(t) Thus, from a given pulse period the phase can be calculated by integration, which this function provides.
The pulse period may be given as a) the pulse period over time struct { time, period } b) the taylor coefficients of the pulse ephemeris via struct { p0[, t0, pdot[, p2dot]] } as defined in 'check_pulseperiod_orbit_struct'.
It is assumed that the unit of the pulse period is a factor of 86400 larger than the time unit of the light curve, i.e., light curve in days and pulse period in seconds. If both have the same unit (regardless whether its seconds or days) use the sameunit-qualifier.
In case a), the phase is numerically integrated using the trapez method on the full provided period time grid. Finally, the calculated phase is re-mapped onto the input time grid by interpolation.
In case b), the phase is calculated analytically by a taylor-series: phi(t) = (t-t0)/p0 - pdot/p0^2*(t-t0)^2 + ...
Note, that the analytical calculation only supports a pulse ephemeris up the the second order (p2dot). If you need higher orders, you can first calculate the pulse period over time via a taylor-series and finally use method a).
In case orbital elements are provided, the additional phase shift is calculated after Hilditch Eq 3.43: \delta \phi = z(t)/c (f0 \dot{z}(t)/c - f(t-t0)) with the projected position z(t) of the neutron star and its spin frequency evolution f(t) = 1 / p(t), where p(t) is the given pulse period evolution.
See also: pulse_period, taylor, BinaryPos, radial_velocity_binary
pulseperiod_epfold
Synopsis
UNDER DEVELOPMENT; performs an automatic pulse period search on the given light curve(s)
Usage
Struct_Type pulseperiod_epfold(Struct_Type[] lc, Double_Type p0);
Qualifiers
nbins - profile bins (default: 32) fracexp - minimum fracexp all bins should have (default: 1.0) dpscale - period search range relative to formal resolution of epfold (default: 3) dpmin - minimum period search range (default: 0.01 s) gapscale - factor for maximum allowed gap length relative to formal resolution o f epfold (default: .5) goodness - threshold for the goodness of any signal (default: 3) pbins - mininum number of consecutive bins in epfold with 'goodness' to define a signal (default: 3) chatty - chattiness (default: 1) plotlc - reference to function(Struct_Type[] lc) plotepf - reference to function(Struct_Type epf, Double_Type median, norm)
Description
UNDER DEVELOPMENT, USE WITH CAUTION! Send questions or bugs to matthias.kuehnel@sternwarte.uni-erlangen.de
lc[] = struct { time, rate, error, fracexp } with time in MJD
p0 = pulse period in seconds
returns struct { time (mean in MJD) period (in s) error (in s) epfold (structure returned by epfold) lc (input light curves, but maybe splitted) }
See also: epfold
pulseperiod_magic
Synopsis
UNDER HEAVY DEVELOPMENT; USAGE ON YOUR OWN RISK
Usage
Struct_Type pulseperiod_magic(Struct_Type[] lc, Double_Type p0, dpmax);
Qualifiers
see pulseperiod_epfold ccflim - cross-correlation threshold (default: .7)
Description
UNDER DEVELOPMENT, USE WITH CAUTION! Send questions or bugs to matthias.kuehnel@sternwarte.uni-erlangen.de
lc[] = struct { time, rate, error, fracexp } with time in MJD
p0 = pulse period in seconds dpmax = maximum allowed period derivative
returns struct { time (mean in MJD) period (in s) error (in s) }
See also: pulseperiod_epfold, pulseperiod_phase_connect, bayesian_blocks
pulseperiod_search
Synopsis
Looks for periodic signal in a lightcurve.
Usage
Struct_Type pulseperiod_search( lc [, p0 [, sigma ] ] )
;
Description
The input lightcurve should be a structure of the form lc = struct { time, rate, [ error ], [ fracexp ] } If error resp. fracexp are not given, they will be filled with sqrt(rate) resp. ones. All fields are arrays of doubles of the same length. p0 and sigma determine an initial guess for the period and an approximate error. The output is a structure of the form { ep, pp, lc } Each field is a list, in which each element corresponds to one segment of the lightcurve (multiple elements if the lightcurve contains considerable gaps). ep[i] contains epoch folding information, pp[i] the folded pulse profile and lc[i] the lightcurve under consideration. The function works the following way: (the functions used are given)
(1) If no p0 is given, a Fourier method is used to find an initial guess for the period. A splitting is performed first if necessary * split_and_epfold_lc (; split_only), foucalc (2) The input lightcurve is split if it contains gaps. * split_and_epfold_lc (; split_only) (3) A correction for any underlying variation in flux of the lightcurve is performed. * pulse2pulse_flux_lc (4) For each segment an epochfolding is performed. * split_and_epfold_lc (5) The peak in the statistics of epfold is found. * find_peak
Qualifiers for the individual functions can be passed as structures with the names given below. See their helps for further information.
Qualifiers
- compact: output contains only the essential information
- epfold_fourier_qu: contains qualifiers for step (1)
- fourier_qu: contains qualifiers for step (1)
- epfold_split_qu: contains qualifiers for step (2)
- pulse2pulse_slope_qu: contains qualifiers for step (3)
- epfold_qu: contains qualifiers for step (4)
- find_peak_qu: contains qualifiers for step (5)
- no_slope_correction: step (3) is not done
- no_splitting: step (2) is not done
- exact: forces epfold to use exact
- not_exact: forces epfold not to use exact
- fourier: fourier is used, even if p0 is given -- can be used to pass sigma
- chatty: boolean value (default: 1)
See also: foucalc, split_and_epfold_lc, pulse2pulse_flux_lc, find_peak
pulseperiod_transform
Synopsis
transforms a period evolution given as taylor coefficients to a new t0
Usage
Struct_Type pulse_transform(Double_Type new_t0, Struct_Type pulseperiod);
Qualifiers
sameunit : set if the unit of 'time' is in the same unit as the pulse period or set to the conversion factor from days to 'time'
Description
The structure containing the pulse period and its derivatives at a certain time t0 must fullfil the conditions given in 'check_pulseperiod_orbit_struct'. This description of the pulse period evolution is transformed to a new given t0. In case this time is given as an array an array of transformed structures is returned.
Note, that the uncertainty of the transformed pulse ephemeris scales with (t - new_t0)^N * pNdot with the highest order N of the taylor series.
See also: check_pulseperiod_orbit_struct
pulseprofile_compose
Synopsis
composes a pulse profile from sine and cosine functions
Usage
Struct_Type pulseprofile_compose(Struct_Type decomposition);
Description
Inverse of pulseprofile_decompose, see its help for details.
Qualifiers
none
See also: pulseprofile_decompose, pfold
pulseprofile_decompose
Synopsis
decomposes a pulse profile into sine and cosine functions
Usage
Struct_Type pulseprofile_decompose(
Struct_Type profile[, Integer_Type[] a_orders, b_orders]
);
Qualifiers
amin/amax - array of min/max values for the allowed range of coefficient a_n during the fit (same order as a_orders). Needs a_orders and b_orders to be specified (default: NULL) bmin/bmax - same as amin/amax for b_n phi0rng - range of allowed phase offets [min,max] (default: [-1,1]) amplrng - range of allowed amplitudes [min,max] (default: [-_Inf,+_Inf]) initpars - array of initial parameters in the form [phi0[, ampl[, a_n..., b_n...]]] The values for a_n and b_n only apply if a_orders and b_orders are specified (default: [0, mean(profile.value), 0..., 0...])
Description
The given pulse profile, F, of the form struct { Double_Type[] bin_lo, bin_hi, value } is fitted by a series of sine and cosine functions,
F(phi) = ampl + \sum_n^N a_n*sin(2PI*(phi+phi_0)*n)
- b_n*cos(2PI*(phi+phi_0)*n)
where phi is the phase bin, ampl is the mean flux, N is the highest order to use, phi_0 is a phase offset, and a_n and a_b are the coefficients of the sine and cosine, respectively. By default the series is calculated up to the highest possible order, which is related to the number of phase bins (Nyquist frequency). It is also possible to specify the orders, which has to be taken into account during the fitting, for the sine and cosine function (a_orders and b_orders, respectively). The returned structure contains the resulting fit parameters,
struct { Integer_Type nbins, % number of phase bins Double_Type phi0, % phase offset Double_Type ampl, % mean amplitude (zero order) Double_Type[] a, % sine coefficients Double_Type[] b % cosine coefficients }
where the indices of the arrays a and b specify the order, n, of the coefficient (starting with n=1). This structure can be passed to 'pulseprofile_compose' to calculate the modelled pulse profile from the coefficients.
Example
% synthetic example profile variable nbins = 32; variable prof = struct { bin_lo, bin_hi, value }; (prof.bin_lo, prof.bin_hi) = linear_grid(0, 1, nbins); prof.value = 2 + 8*sin((prof.bin_lo)*2*PI)^2; prof.value += 4*sin((prof.bin_lo)*PI); hplot(prof);
% decomposition using the first order sine and second % order cosine coefficients only variable decomp = pulseprofile_decompose(prof, [1], [2]); ohplot(pulseprofile_compose(decomp));
Qualifiers
maxord: highest order to use (if no specific orders for a or b are given; default: nbins/2-1)
See also: pfold, pulseprofile_compose
pulseprofile_energy_interpolate
Synopsis
interpolates a pulsephase-energy-histogram to a finer grid
Usage
Double_Type[egrid\*e_interp,phigrid\*p_interp] pulseprofile_energy_interpolate(
Struct_Type map, Integer_Type p_interp, Integer_Type e_interp
);
Description
Example
See also: pulseprofile_energy_map
pulseprofile_energy_lag
Synopsis
calculates a lag (=shift) in a pulsephase-energy-histogram
Usage
Double_Type[] pulseprofile_energy_lag(
Struct_Type map, 'p' or 'e'[, Double_Type[] reference]
);
Description
Example
See also: pulseprofile_energy_map, CCF_1d
pulseprofile_energy_map
Synopsis
sorts the given events into a pulsephase-energy-histogram
Usage
Double_Type[egrid,phigrid] pfold_event_energy_map(
Double_Type[] events, energies, Double_Type pulse_period;
gti = Struct_Type, egrid = Struct_Type, pgrid = Struct_Type
);
Description
Example
See also: pfold
pulseprofile_energy_normalize
Synopsis
normalize a given pulse profile, spectrum, or pulseprofile-energy-map
Usage
Double_Type[] pulseprofile_energy_normalize(
Double_Type[] profile_or_spectrum,
String_Type method
);
or
```c
Double_Type[egrid,pgrid] pulseprofile_energy_normalize(
Double_Type[egrid,pgrid] map,
String_Type method,
Char_Type dimension
);
##### Description
A function to normalize either a 1D array representing a pulse profile or
spectrum, or a 2D pulse profile energy map. In the latter case the dimension
the normalization should apply to has to be specified, as the normalization
is performed in 1D, i.e., either for each phase normalize the energies
(dim = 'e') or for each energy normalize the pulse profile (dmi = 'p').
Following normalization methods are available:
"sdev": Substracts the mean value and devides by the standard deviation:
isis> mom = moment( value );
isis> normvalue = ( value - mom.ave ) / mom.sdev;
"minmax": Substracts the min. value and deviedes by the min-max range.
isis> normvalue = ( value - min(value) ) / (max(value)-min(value));
##### Example
__See also__: pulseprofile_energy_map
#### pulseprofile_phase_connect
##### Synopsis
UNDER DEVELOPMENT; returns the phase shift between two pulse profiles
##### Usage
```c
(Double Type phi, ccf) = pulseprofile_phase_connect(Struct_Type prof, ref);
Qualifiers
shift - automatically shifts the given profile in order to match the phase of the reference profile (caution: overwrites the input!)
Description
UNDER DEVELOPMENT, USE WITH CAUTION! Send questions or bugs to matthias.kuehnel@sternwarte.uni-erlangen.de
Determines the phase shift of the pulse profile 'prof' with respect to the reference profile 'ref' by a cross-correlation. The cross-correlation is interpolated to enhance the precision beyond the pulse profile binning. Both pulse profiles need to have the same number of bins.
The input structures have to have the same fields as returned by 'epfold'.
The resulting phase shift 'phi' and the maximum value of the cross-correlation 'ccf' are returned.
See also: pfold, CCF_1d
rebinGroup
Synopsis
rebins a {bin_lo, bin_hi, value, err} struct
ure by a grouping factor
Usage
Struct_Type s_ = rebinGroup(Struct_Type s, Integer_Type n);
See also: rebin
rebin_fouquan
Synopsis
rebins a Fourier quantity from timing analysis with foucalc
Usage
Struct_Type rebin_fouquan(freq, fouquan, numseg)
Description
freq
is an array of the original Fourier frequencies.
fouquan
is the array of original Fourier quantities.
numseg
is the number of seqments used to calculate fouquant
.
The output structure has the following fields, which are arrays:
-
freq_lo
andfreq_hi
define the new frequency bin. -
freq
is the average original frequency in each bin. -
n
is the number of original frequencies, and -
n_tot
is the total number of values (including segmentation) contributing to each bin. -
value
is the average Fourier quantity. -
error = value/sqrt(n_tot)
is the 1 sigma error for then_tot
values. -
sigma
is the standard deviation of the original quanities in each bin
Qualifiers
- logfreq
- linfreq
- newfreq: array of new frequencies
- nofreq: no rebinning
- verbose
rebin_human2isis
Synopsis
create isis binning from str = rebin_human2isis
Usage
Array_Type binning = rebin_human2isis(String_Type str);
Description
This routine converts a binning string created with rebin_isis2human back to the ISIS conform binning array (see e.g. get_data_info(idx).rebin). It can be used to restore the binning of data directly.
See also: fits_save_fit, get_data_info, rebin_isis2human
rebin_isis2human
Synopsis
converts isis binning in a human readable string
Usage
String_Type str = rebin_isis2human(Array_Type binning);
Description
This routine converts the binning as used by ISIS (see e.g. get_data_info(idx).rebin) to a string which is easy to read.
The routine fits_save_fit() uses this string to save the binning of the data. It can be converted back by using rebin_human2isis().
See also: fits_save_fit, get_data_info, rebin_human2isis
rebin_lc
Synopsis
rebins a lightcurve structure to a new time resolution
Usage
Struct_Type new_lc = rebin_lc(Struct_Type lc, Double_Type new_dt);
Qualifiers
- time[=
"time"
]: field name in the structure containing the time - rate[=
"rate"
]: field name(s) in the structure containing the rate(s) - error[=
"error"
]: field name(s) containig the corresponding error(s) - float: type cast rate(s) and error(s) to Float_Type
- verbose: shows the assumed time resolution of the initial light curve
Description
The original light curve lc
may contain discontinuities.
new_lc.rate[i]
contains the average lc.rate
where new_lc.time[i] <= lc.time < new_lc.time[i]+new_dt
.
error^2
is rebinned accordingly.
The structure fields "time, rate, error
" may have arbitrary names,
but these must then be specified by the according qualifiers.
new_lc
will also have these field names. In addition, a "time_hi
"
field is added, which contains the upper boundary of the time bins.
Example
lc1 = struct { time=[1:10], rate=[1:10]^2, , error=[1:10] };
LC1 = rebin_lc(lc1, 2);
lc2 = struct { t=[1:10], r1=[1:10]^2, r2=[1:10]^3, e1=[1:10], e2=[1:10]^1.5 };
LC2 = rebin_lc(lc2, 2; time="t", rate=["r1", "r2"], error=["e1", "e2"]);
See also: rebin
rebin_to_energy_grid
Synopsis
rebin data to a given energy grid
Usage
rebin_to_energy_grid(Dataset ID, bin_lo, bin_hi);
Description
This routine bins the data to a given energy grid. Bins outside the given energy grid are left unbinned.
See also: rebin_data,rebin_satellite, rebin_human2isis
rebin_to_instrument
Synopsis
rebin data from one instrument to another instrument
Usage
rebin_to_instrument(Dataset id_reference, Dataset id_to_be_rebinned);
Description
This routine tries its best to rebin the given dataset with the ID "id_to_be_rebinned" to the dataset with id "id_reference".
See also: rebin_to_energy_grid,rebin_data
resolution_rebin
Synopsis
Rebin spectra based on approximated detector resolution.
Usage
Array_Type resolution_rebin(dataset_id, instrument; oversampling, min_counts);
Qualifiers
- sampling: Float_Type, oversampling factor of detector resolution/bin width, default: 3
- min_sn: Float_Type, minimum signal-to-noise ratio, default: 5
Description
This function rebins a dataset to fulfill the following two criteria:
1: Each bin contains enough counts so that - assuming Poissonian noise - Signal-to-Noise is larger than the given minimum.
2: If criterion 1 is fulfilled, the bin size is chosen as a fraction of the detector resolution, an oversampling of 3 is recommended (Kaastra 2016).
Required arguments are the dataset ID and a string that specifies the used instrument. Supported are:
XMM EPICpn: "epn"
Suzaku XIS: "xis"
NuSTAR FPM: "fpm"
Examples
Load and plot EPICpn data with a histogram oversampling detector resolution by 5:
variable EPN_data = load_data("src_s.pha");
resolution_rebin(EPN_data,"epn"; sampling = 5);
plot_data(EPN_data; dsym=0);
See also: rebin_data, group
save_atime
Synopsis
saves a structure of arrival times into a FITS-file
Usage
save_atime(String_Type filename, Struct_Type atime, String_Type extname[, String_Type[] comments]);
Qualifiers
obj - observed object, written into FITS header sat - used satellite, written into extension header nfold - number of arrival times, which were merged during the determination. Corresponds to the 'indiv' qualifier for the 'atime_det' function. Written into extension header hfits - structure for additional FITS-header fields hext - structure for additional extension-header fields newfile - if the given filename exists a new file is created instead of updating the existing one newext - if the given extension already exists in the FITS-file a further one with the same name is added instead of updating the existing one
Description
A structure of arrival times, e.g. as returned by the 'atime_det' function, is save into the given FITS-file creating or updating the extension of the geiven name. Most of the qualifiers allow to store addiotional informations into the header of the FITS-file or the extension. The optional fourth parameter can be used to write additional comments into the FITS-file.
See also: atime_det, atime_merge, load_atime
save_atime_beta
Synopsis
saves a structure of arrival times into a FITS-file
Usage
save_atime(String_Type filename, Struct_Type atime, String_Type extname[, String_Type[] comments]);
Qualifiers
obj - observed object, written into FITS header sat - used satellite, written into extension header nfold - number of arrival times, which were merged during the determination. Corresponds to the 'indiv' qualifier for the 'atime_det' function. Written into extension header hfits - structure for additional FITS-header fields hext - structure for additional extension-header fields newfile - if the given filename exists a new file is created instead of updating the existing one newext - if the given extension already exists in the FITS-file a further one with the same name is added instead of updating the existing one
Description
A structure of arrival times, e.g. as returned by the 'atime_det' function, is save into the given FITS-file creating or updating the extension of the geiven name. Most of the qualifiers allow to store addiotional informations into the header of the FITS-file or the extension. The optional fourth parameter can be used to write additional comments into the FITS-file.
See also: atime_det, atime_merge, load_atime
savitzky_golay_coefficients
Synopsis
Calculate the Savitzky-Golay coefficients used for data smoothing
Usage
Double_Type[] c = savitzky_golay_coefficients(positive Integer_Type nl, nr, p; qualifiers)
Description
Calculate the Savitzky-Golay coefficients used for data smoothing. The arguments
nl'/
nr' are hereby the data points to the left/right while p' is the polynomial degree (all have to be positive integer numbers; otherwise their absolute value is rounded to the next nearest integer). The quanity
p' must not exceed nl'+
nr'!
The returned array of coefficients is ordered as [c_-nl, ..., c_0, ..., c_nr].
Notes
Requires GSL module.
Qualifiers
- derivative [=0]: Savitzky-Golay coefficients used for calculating the numerical derivative of the corresponding order (must not exceed `p', i.e., the order of the polynomial).
Example
% coefficients to smooth data: c = savitzky_golay_coefficients(15,15,4); % -> coefficients to compute first derivative: c = savitzky_golay_coefficients(15,15,4; derivative=1);
See also: savitzky_golay_smoothing
savitzky_golay_smoothing
Synopsis
Smooth noisy data by using a Savitzky-Golay filter
Usage
Double_Type[] smoothed_data = savitzky_golay_smoothing(Double_Type[] data,
positive Integer_Type nl, nr, p; qualifiers)
Description
The idea of Savitzky-Golay filtering is to approximate the given data within
a moving window (ranging from nl' data points to the left to
nr' data points
to the right) by a polynomial of order `p'. The respective polynomials are -
in principle - determined by least-squares fits to the window data points.
Luckily, Savitzky & Golay found a way to replace the fitting and evaluating
of the polynomial by taking just linear combinations of neighboring data points
tremendously speeding up the smoothing process. However, their method, which
is implemented here, is valid only for regularly spaced data points. For more
details on the method, see ``Numerical Recipes'', Third Edition, Section 14.9.
Note that the input parameters nl',
nr', and p' have to be positive integer numbers (otherwise their absolute value is automatically rounded to the next nearest integer). The quantity
p' must not exceed the minimum window semi-
length min(nl',
nr') and nl'+
nr' has to be smaller than or equal to the
number of total data points.
Practical hint: Best results are obtained when the full width of the degree 4 Savitzky-Golay filter is between 1 and 2 times the full width at half maximum of the desired features in the data.
Notes
Requires GSL module.
Qualifiers
- derivative [=0]: Apart from smoothing noisy data, the Savitzky-Golay method
is also capable to compute numerical derivatives from the fitted polynomials.
In order to calculate the
k'-th derivative of the data, set this qualifier equal to
k' and divide the returned array by the data stepsize to the power ofk'. Note that
k' has to be a positive integer (otherwise it is replaced by the next nearest integer of its absolute value). Note thatk' must not exceed
p'. - periodic: Set this qualifer to assume periodic boundary conditions for your
data avoiding special treatment of data points close to the edges and thus
speeding up the computation. In this case, the condition
p'<min(
nl',nr') is replaced by
p'<nl'+
nr'.
Example
x = [0:20:#1000]; data = sin(x); data = data + 0.1*grand(length(x)); smoothed_data = savitzky_golay_smoothing(data,100,100,4); first_derivative = savitzky_golay_smoothing(data,100,100,4; derivative=1)/(20./1000.); second_derivative = savitzky_golay_smoothing(data,100,100,4; derivative=2)/(20./1000.)^2;
See also: savitzky_golay_coefficients
scargle
Synopsis
Computes the lomb-scargle periodogram of an unevenly sampled lightcurve
Usage
Struct_Type res = scargle (t, c); % where t and c contain time and counts of the lc
Qualifiers
- fmin: minimum frequency to be used (NOT ANGULAR FREQ!), has precede over pmin, pmax
- fmax: maximum frequency to be used (NOT ANGULAR FREQ!), has precede over pmin, pmax
- pmin: minimum period to be used
- pmax: maximum period to be used
- omega: array of angular frequencies for which the PSD values are desired; if set, value for numf will be reset to length of omega during code
- noise: for the normalization of the periodogram. if not set, equal to the variance of the original lc.
- numf: number of independent frequencies
- old: if set computing the periodogram according to J.D. Scargle, 1982, ApJ 263, 835 if not set, computing the periodogram with the fast algorithm of W.H. Press and G.B. Rybicki, 1989, ApJ 338, 277.
- nu: if set, output structure also contains frequency
- om: if set, output structure also contains angluar frequency
Description
(transcribed from IDL-program scargle.pro)
The Lomb Scargle PSD is computed according to the definitions given by Scargle, 1982, ApJ 263, 835, and Horne and Baliunas, 1986, ApJ 302, 757. Beware of patterns and clustered data points as the Horne results break down in this case! Read and understand the papers and this code before using it! For the fast algorithm read W.H. Press and G.B. Rybicki, 1989, ApJ 338, 277.
The code is still stupid in the sense that it wants normal frequencies, but returns angular frequency...
The transcribed version is version 1.7, 2000.07.28
Unlike the IDL function, the isis code returns a structure containing the power spectral density (psd) and period and, if specified by qualifiers, also the frequency (nu) and angular frequency (om).
Example
variable t = [0:100:0.1]; variable c = sin(t)+((rand(1000)/(2^32-1)*0.2)+0.9); variable res = scargle(t, c; nu); plot(res.period,res.psd); plot(res.nu,res.psd);
set_bin_corr_factor
Synopsis
sets a simple bin correction factor (used in plot_data/unfold)
Usage
set_bin_corr_factor(Integer_Type data_id, Double_Type[nbins] corr_factor);
See also: load_fermi,get_bin_corr_factor,plot_unfold
simfit_namespace
Synopsis
implements all SimFit-functions into a namespace
Usage
simfit_namespace(Struct_Type SimFit);
Qualifiers
name - name of the new namespace (default: simfit) chatty - chattiness of this function (default: 1)
Description
Takes a simultaneous fit structure as input and implements all available functions defined in there into a new namespace. In this way the user no longer has to type the structure's name in front of the functions, but use the functions directly.
Note that in this namespace the ISIS intrinsic functions, such as 'fit_fun', 'set_par', etc. are being overwritten by the SimFit-versions.
Furthermore, 'eval_counts' is defined as a combination of 'eval_groups' and 'eval_global', and 'fit_counts' now points to 'fit_smart'.
To access the original ISIS functions just switch the namespace back to 'isis' using 'use_namespace' or access the functions via isis->function_name
WARNING: there is no check implemented yet that the given structure actually is a SimFit-structure, so every structure containing reference to functions may be passed, which might HARM your ISIS-session or machine!
FINALLY, if there will be any error the namespace is most likely set to 'isis' afterwards (we tried to catch errors in the SimFit-functions but this does not work for, e.g., syntax errors within the shell).
See also: simultaneous_fit, use_namespace
simulate_data
Synopsis
fakes a spectrum
Usage
Integer_Type id = simulate_data(String_Type RSPfile, Double_Type exposure);
Description
The currently defined fit-function is used to fake the spectrum.
See also: fakeit
simultaneous_fit
Synopsis
initializes a structure to handle simultaneous/combined fits of a lot of data
Usage
Struct_Type simultaneous_fit();
Description
The advantage of fitting a lot of data at the same is that fit parameters, which seem to be the same for all data, can be fitted properly. That results in a reduced number of free parameters for individual observations or data groups, and thus constrains parameters even better at low data quality.
The disadvantages are, that a) a lot of painful work has to be done on tieing and freezing parameters and b) a fit using 'fit_counts' will take hours to days. There are several functions implemented within the returned structure, which can help to solve that issues.
To access the help of a function, simply call it with the 'help' qualifier (like for the xfig functions).
References for simultaneous fits in ISIS are Kuehnel et al. 2015, Acta Polytechnica 55(2), 123127 Kuehnel et al. 2016, Acta Polytechnica 56(1), 4146
NOTE: The simultaneous fits are still in development.
spearmanrho
Synopsis
Timing Tools: Spearman's Rank Correlation Coefficient
Usage
(rho) = spearmanrho (a,b);
Description
rho = 1 - frac {6 sum d_i^2} {n(n^2 -1)}
d_i = a_i - b_i : difference between the ranks of corresponding values n : number of values
Alternative: The statistics module (require("stats");
) includes
several functions for determining (rank) correlations coefficients,
such as spearman_r
, which provides the p-value (as well as the
correlation coefficient rho
).
split_and_epfold_lc
Synopsis
Takes lightcurve and pulse period, splits if necessary and executes epfold for each segment.
Usage
Struct_Type split_and_epfold_lc(Struct_Type OR Array_Type lc, Double_Type p0), all times in the same units
Qualifiers
- exact: Forces execution of pfold with exact qualifier
- not_exact: Forces execution of pfold without exact qualifier
- exact_threshold: see description (default: 6.)
- split_only: Only the splitting is performed
- nbins: Sampling of pfold (default: min(32, p0 / dt), dt time resolution
- nsrch: Number of trial periods for epfold (default: 1e4)
- grid_steps: Period resolution for epfold (default: not set)
- search_range: Period interval considered left/right of p0 in terms of dp=p0^2/T (default: 2.)
- dpmin: For very long lightcurves, this can replace dp (default: p0e-3)
- gap_scale: the smaller, the higher is the sensitivity for gaps (default: .5)
- fracexp: Only timebins with fracexp >= this value will be considered (default: 1.)
- chatty: chattiness (default: 1)
Description
The input lightcurve(s) should be a structure or an array of structures of the form: lc = struct { time, rate, fracexp, error } where each field is an array of doubles. The output is a structure of the form: out = struct { time, [ epfold ], lc, dp, exact } Each field is an array, where each element corresponds to one of the segments the split has produced. The individual fields are: * time: average time of segment * epfold: struct { p, stat, nbins, badp } (not output if 'split_only' qualifier has been used) * lc: of the same form as input structure lc * dp: either p^2 / T or dpmin * exact: 1 if exact qualifier was used for epfold, 0 otherwise Remarks regarding qualifiers: * by default, exact is only chosen if there are less than exact_threshold pulses in the timespan considered. However, as of Sep 16, there was a bug in pfold, so that exact is ONLY used if exact_threshold the 'exact_threshold' qualifier exists. When this bug is fixed, please change the constant at the beginning of the script * nsrch and grid_steps are mutually exclusive. If you set nsrch to 0, the default epfold routine will be used * All qualifiers can also be passed using a structure named 'epfold_qualifiers'. If this qualifier is present, its content will overwrite all other qualifiers given. If passed using this structure, all qualifiers must be assigned values (e.g. exact = 0, split_only = 1). Makes scripts easier to read.
See also: pfold, epfold, pulseperiod_epfold, pulseperiod_search, find_peak
string2sys_err_array
Synopsis
create sys_err_array from str = sys_err_array2string
Usage
Array_Type sys_err = string2sys_err_array(String_Type str);
Description
This routines converts a string created by the function "sys_err_array2string" back to an array, which can then be applied to the data by the function set_sys_err_frac.
See also: fits_save_fit, set/get_sys_err_frac, sys_err_array2string
sum_lc
Synopsis
combines light curves from one structure
Usage
Struct_Type sum_lc(Struct_Type lc, Array_Type inds)
Qualifiers
- time: name of the time field (default:
"time"
) - rate: format string for the rate field (default:
"rate%S"
) - error: format string for the error field (default:
"error%S"
)
Description
The rates in the fields of lc
(named by inds
via the rate
format statement) are summed up.
The corresponding errors (named by the error
format statement)
are added in quadrature.
The returned structure has the fields time, rate, error
.
Examples
% 1. Assume you have lc = struct { time, rate1, error1, rate2, error2, rate3, error3 }; % Then slc = sum_lc(lc, [1:3]); % is equivalent to slc = struct { time = lc.time, rate = lc.rate1 + lc.rate2 + lc.rate3, error = sqrt( lc.error1^2 + lc.error2^2 + lc.error3^2 ) };
% 2. Assume you have lc = struct { time, rate_a, err_a, rate_b, err_b, rate_c, err_c }; % Then slc = sum_lc(lc, ["a", "b", "c"]; rate="rate_%s", error="err_%s"); % is equivalent to slc = struct { time = lc.time, rate = lc.rate_a + lc.rate_b + lc.rate_c, error = sqrt( lc.err_a^2 + lc.err_b^2 + lc.err_c^2 ) };
sys_err_array2string
Synopsis
creates a string from an Array returned by "get_sys_err_frac"
Usage
String_Type str = sys_err_array2string(Array_Type sys_err);
Description
This routine creates a String from the Array of systematic errors apllied to the data. Such an array is returned by the function "get_sys_err_frac".
See also: fits_save_fit, set/get_sys_err_frac, string2sys_err_array
xfigplot_pulseprofile_energy_map
Synopsis
returnes an Xfig plot object containing the pulse-profile-energy-map
Usage
Xfig_Object xfigplot_pulseprofile_energy_map(Struct_Type map);
Description
Struct_Type map = struct{ pgrid = struct{ bin_lo = Double_Type[np], bin_hi = Double_Type[np] }, egrid = struct{ bin_lo = Double_Type[ne], bin_hi = Double_Type[ne] }, map = Double_Type[ne,np] };
Example
See also: pulseprofile_energy_map
xnotice_atime
Synopsis
notices a range of defined arrival times
Usage
xnotice_atime(Integer_Type index[, Double_Type time_lo, Double_Type time_hi]);
or xnotice_atime(Integer_Type index, Double_Type time);
or xnotice_atime(Double_Type time);
Qualifiers
none
Description
There are three different ways to notice the time range of arrival times defined in a dataset:
- Like 'xnotice' set the time range of dataset 'index' explicitly from 'time_lo' to 'time_hi'.
- Find the border of the time range of dataset 'index' nearest to the given time. Then set the found border to this time. The border of a possible adjacent dataset is changed respect- ively.
- Find that dataset, which time range has the nearest border to the given time and proceed like option two. If the dataset only is given, then the full time range is used.
See also: define_atime, xnotice
z2fold
Synopsis
peforms Z^2_m search On a lightcurve or eventdata in a given period interval
Usage
(Struct_Type res) = z2fold(Double_Type t, r, pstart, pstop);
or (Struct_Type res) = z2fold(Double_Type t, pstart, pstop) ; (event data)
Qualifiers
%* sampling: how many periods per peak to use (default=10)
- nsrch: how many periods to search in a linear grid (default not set)
- dp: delta period of linear period grid (default not set)
- m: number of harmonics used (default =2)
- chatty: set the verbosity, (default=0)
- gti:GTIs for event data, given as struct{start=Double_Type, stop=Double_Type}
Description
Performs Z^2_m test on a given lightcurve or event data between the periods pstart and pstop. The function is based on epfold.sl. GTI correction only implemented for event-data yet.
By default, periods are sampled according to the triangular rule for estimating the period error, using "sampling" periods per peak. If a linear grid is to be used, either "dp" for a given distance between two consecutive periods or "nsrch" for a given number of periods can be given. These qualifiers are mutually exclusive.
The returned structure "res" contains four fields: "p" for the evaluated period and "stat" for the value of the statistic used.
See also: zsquare, epfold, pfold
zsquare
Synopsis
calculate Z^2 statistics for pulse period search
Usage
Double_Type Z2 = zsquare(Double_Type times, Double_Type testperiod
Description
Calculates the Z^2 statistics of an event list or light-curve for a given period. By itself not very useful, but is used in z2fold to search for periodicities using the Z^2 statistics.
Qualifiers
- lc: =Double_Type rate; allows the use of a lightcurve instead of event files
- m: = Integer_Type; determines the highest order of summations (default 2)