Reference for misc functions
AS
Synopsis
Evaluate equations of motion, total energy, or circular velocity derived from a revised Allen & Santillan potential
Usage
AS(Double_Types t, m[6,n]; qualifiers)
Description
Evaluate the equations of motion, the total energy, or the circular velocity at time 't' derived from the revised Galactic gravitational potential by Allen & Santillan (see Model I in Irrgang et al., 2013, A&A, 549, A137) using either cylindrical coordinates (r [kpc], phi [rad], z [kpc]) and their canonical momenta vr [kpc/Myr], Lz [kpc^2/Myr], vz [kpc/Myr]) or cartesian coordinates (x [kpc], y [kpc], z [kpc], vx [kpc/Myr], vy [kpc/Myr], vz [kpc/Myr]), see qualifier 'coords'. Conservation of angular momentum Lz is implemented in the equations of motion for cylindrical coordinates only. The total energy E_total [kpc^2/Myr^2] is not used to integrate the equations of motion although being a conserved quantity, too. Therefore, conservation of energy, i.e., of E_total, is a measure for the precision of the numerical methods applied.
For computing orbits with n different initial conditions, the input parameter m is a [6,n]-matrix with (qualifier("coords")=="cyl") or (qualifier("coords")=="cart") m[0,*] = r; m[0,*] = x; m[1,*] = phi; m[1,*] = y; m[2,*] = z; m[2,*] = z; m[3,*] = vr; m[3,*] = vx; m[4,*] = Lz; m[4,*] = vy; m[5,*] = vz; m[5,*] = vz; If the qualifier 'eomecd' is set to "eom", the function returns a [6,n]-matrix delta with delta[0,*] = vr; delta[0,*] = vx; delta[1,*] = Lz/r^2; % = vphi delta[1,*] = vy; delta[2,*] = vz; delta[2,*] = vz; delta[3,*] = -d/dr (Potential(r,z) + Lz^2/r^2); delta[3,*] = -d/dx Potential(x,y,z); delta[4,*] = 0; % -d/dphi Potential(r,z) delta[4,*] = -d/dy Potential(x,y,z); delta[5,*] = -d/dz Potential(r,z); delta[5,*] = -d/dz Potential(x,y,z); If the qualifier 'eomecd' is set to "energy", the function returns a [1,n]-array storing the total energy for each orbit: E_total(r,z,vr,vz,Lz) = Double_Type[n] = 0.5*(vr^2+vz^2+Lz^2/r^2) + Potential(r,z) or E_total(x,y,z,vx,vy,vz) = Double_Type[n] = 0.5*(vx^2+vy^2+vz^2) + Potential(x,y,z) If the qualifier 'eomecd' is set to "circ", the function returns a [1,n]-array storing the circular velocity for each orbit: v_circ(r,z) = Double_Type[n] = sqrt( r * d/dr Potential(r,z) ) If the qualifier 'eomecd' is set to "sgcd", the function returns the Sun-Galactic center distance found to fit best to this potential.
Qualifiers
- coords [
="cyl"
]: Use cylindrical ("cyl") or cartesian ("cart") coordinates. - eomecd [
="eom"
]: Return equations of motion ("eom"), total energy ("energy"), circular velocity ("circ"), or Sun-Galactic center distance ("sgcd"). - Mb [
=409
]: Mass of bulge in Galactic mass units, see Irrgang et al. 2013. - Md [
=2856
]: Mass of disc in Galactic mass units, see Irrgang et al. 2013. - Mh [
=1018
]: Mass scale factor of halo in Galactic mass units, see Irrgang et al. 2013. - bb [
=0.23
]: Bulge scale length, see Irrgang et al. 2013. - ad [
=4.22
]: Disc scale length, see Irrgang et al. 2013. - bd [
=0.292
]: Disc scale length, see Irrgang et al. 2013. - ah [
=2.562
]: Halo scale length, see Irrgang et al. 2013. - exponent [
=2
]: Exponent in the halo mass distribution, see Irrgang et al. 2013. - cutoff [
=200
]: Halo cutoff, see Irrgang et al. 2013.
Example
delta = AS(0, m); energy = AS(0, m; eomecd="energy"); v_circ = AS(0, m; eomecd="circ"); sgcd = AS(; eomecd="sgcd");
See also: orbit_calculator, MN_NFW, MN_TF, plummer_MW
Chandra_display_mask
Synopsis
shows the content of a Chandra mask file (msk1.fits)
Usage
Chandra_display_mask([String_Type filename]);
Description
filename
can be a globbing expression.
If it is omitted, acis*msk1.fits*
is assumed.
The mask file contains the valid part of the CCD, i.e., the portion for which events can be telemetered.
EchoMap
Synopsis
reprocesses a signal on a given surface
Usage
Struct_Type EchoMap(
Struct_Type surface, Struct_Type signal
Vector_Type observer[, Double_Type length]
);
Qualifiers
c - defines the speed of light (default: 1) noFiniteC - infinite speed of light noDoppler - disable Doppler Shift powRedist - reference to a function returning the re- distributed power of all surface elements (i.e., all parameters are arrays!) in direction to the observer. The parameters passed are * the power to be redistributed * the cosine of the angle between the element and the observer * the surface element as structure with the corresponding fields as below * wether the passed power is the surface's intrinsic one (1) or incident power (2) * a reference to a variable, which may be set to 1 to trigger re-calculating the powers (i.e. the given function is called again for the two cases *after* it was called for both cases in the first place; might be useful if, e.g., the incident power changes properties of the surface) This function is called twice: 1) for the intrinsic power of a surface element and 2) the reprocessed relative (!!) power. The last parameter specifies these two cases. If you need to pass additional qualifiers you may set the powRedistQual-qualifier to a structure holding the qualifiers (default: power*cosine_to_observer) response - response function of a surface element, see EchoMap_makelc for details reproment - instead of returning the echo signal the reproment structure is returned (see text)
all qualifiers are passed to EchoMap_makelc
Description
This function propagates the given 'signal' on an object defined by a 'surface'. There the signal is absorbed and re-emitted (reprocessed) by each surface element. The resulting response signal, as seen in direction to an 'observer', is returned by calling EchoMap_makelc after all geometrical have been calculated (called reproment here). If the corresponding qualifier is set, the latter call is skipped and the reproment structure is returned instead. This structure has the following fields: visible_obs - boolean value if the surface is visible visible_sig by the observer and the signal origin lt_surf - light travel times from the signal to the lt_obs surface and from there to the observing plane doppler - Doppler factor in direction to the observer flux_int - intrinsic flux of the surface as seen by the observer flux_rec - received signal flux of the surface flux_emi - emitted relative flux of the surface as seen by the observer
The following effects are taken into account during the calculation of the response signal:
- light travel time between each surface element to the signal source and to the observing plane
- Doppler shift by a moving surface
- effective areas and projection effects
- power redistribution function at the surface (qualifier 'powRedist')
- power response function of the surface (qualifier 'response', see EchoMap_makelc for details)
Restrictions:
- sqrt(surface.A) << vector_norm(surface.r-signal.origin) (small surface elements relative to distance to source)
- surface.v << c
- vector_norm(observer) != 1
- signal fields have to be sorted by time
The 'surface' structure contains the surface elements, defined by the fields Vector_Type[] r - position vector to each element Double_Type[] A - area of each element Vector_Type[] n - normal vector of each element Vector_Type[] v - velocity vector of each element (optional) Double_Type[] L - intrinsic power of each element, which will be added to the response signal (optional, in energy/time) The 'signal' is treated as power (energy/time) as a function of time and defined by the fields Double_Type[] time - lower time bin (in seconds) Double_Type[] power - power in each time bin Vector_Type origin - position vector to the signal's source (default: (0,0,0)) Any length is considered to be in units of the speed of light (lt-s). If another unit is used, the internal speed of light has to be changed via the 'c' qualifier.
See also: EchoMap_makelc, EchoMap_binary, Vectory_Type
EchoMap_binary
Synopsis
reprocesses a signal on the companion of a binary
Usage
Struct_Type EchoMap_binary(
Struct_Type surface, Double_Type lum, Struct_Type signal,
Struct_Type orb, Double_Type phiorb, Double_Type mq
[, Double_Type length]
);
Qualifiers
any qualifiers are passed to 'EchoMap'
Description
Uses the 'EchoMap' function to reprocess an input 'signal' on the 'surface' of the secondary star with total luminosity 'lum'. The fields required for the 'surface' and 'signal' structure are described in the 'EchoMap' help. The 'Roche_lobe_surface' function might be useful to calculate the deformed surface of the star within the Roche potential. Note, that the velocities of the surface elements are calculated automatically if no 'v' field is specified. If not specified, the luminosity field 'L' for each surface element is calculated as well. Thereby, it is assumed that each surface element has the same area luminosity, such that the total luminosity still is 'lum'. The orbit of the binary is described by the 'orb'ital structure an has to contain the fields as described in the 'check_ephemeris' help. In addition, the inclination can be specified optionally via the 'i' field. To complete the description of the binary, the mass ratio mq = M_primary / M_companion has to be given at last. As default, the source of the signal is the position of the primary star. It can be changed by setting the 'origin' field of the signal structur to the position vector. Using the 'phiorb' parameter the orbital phase is defined, which corresponds to the viewing angle on the binary.
As a conclusion and for further details the reference
frame is shown in the following:
position of M_2 depending
omega (right handed) on orbital phase
observer (i = 0)
z 0.50
| y _ _
| / /
|/ 0.75 | x | 0.25
M_1 o ---- x _ _/
/ | observer
o M_2 (phiorb = 0) 0.00
- orbital plane = x-y-plane (fix)
- inclination and orbital phase is realized by rotating the observer accordingly
- distance M_1 to M_2 = 1 (fix)
WARNING: at the moment circular orbits (ecc = 0) are are implemented only! NOTE: for elliptical orbits, the shape of the stars depends on the orbital phase
See also: EchoMap, Roche_lobe_surface
EchoMap_makelc
Synopsis
takes reprocessing information to produce a response signal
Usage
Struct_Type EchoMap_makelc(
Struct_Type reproment, Struct_Type signal[, Double_Type length]
);
Qualifiers
noDoppler - disable Doppler Shift response - reference to a function returning the relative emissivity of a surface element as response to an incoming power peak at t = 0. For details see the description (default: delta peak at t=0) surfsigs - reference to a variable, which will get the output signal for each element (2dim-array defined as [element,power]; use the time from the returned total signal)
Description
Takes the 'reproment' structure of a former run of EchoMap and returnes the output signal produced by reprocessing the input 'signal' on the surface used to caluclate the reproment (i.e., the geometrical effect of the reprocessing). By default, the returned power vs. time is the total response to the incomming signal. However, if the optional argument 'length' is provided, the input signal is periodically repeated until the output signal has at least the given length. Furthermore, the signal is repeated backwards in time as well such that no rising or declining parts of the echo is visible in the output.
By default, it is assumed that each surface element mirrors the incoming power instantaneously. More complex behaviour can be implemented easily using the 'response' qualifier, which defines a function returning the relative 'power' emissivity over 'time' (in seconds) of a surface element. Furthermore, the response has to fullfil energy conservation such that the integrated power does not exceed 1, but may be less, i.e. power is reprocessed via non-radiative processes. The time resolution of the response has to be equal or better than the input signal!
See also: EchoMap
FriendCastor1982_dMdot_dOmega_CygX1
Synopsis
calculates the mass loss rate per solid angle for Cyg X-1
Usage
Double_Type FriendCastor1982_dMdot_dOmega_CygX1(Double_Type theta)
Description
theta
is the angle from the binary axis in degrees.
The return value is in units of 1e-6 solar masses / year / sr.
(The total mass loss rate 2 pi int dMdot/dOmega(theta) sin(theta) dtheta
is 2e-6 solar masses / year.)
See also: Friend & Castor (1982), Fig. 4
KS_test
Synopsis
computes the test statistics of a two sample Kolmogorov-Smirnov test
Usage
test_statistics = KS_test(Double_Type x1, Double_Type x2);
Description
The null hypothesis of the KS test is that two samples are distributed according to the same distribution. It is rejected if the test statistic D=max(F_1(x),F_2(x)) is greater than a certain value
MAXI_lightcurve
Usage
Struct_Type MAXI_lightcurve(String_Type source);
Description
The returned structure has the following fields:
-
time
: Modified Julian Date -
rate
,err
: 2-20 keV light curve [c/s/cm^2] -
rate_a
,err_a
: 2- 4 keV light curve [c/s/cm^2] -
rate_b
,err_b
: 4-10 keV light curve [c/s/cm^2] -
rate_c
,err_c
: 10-20 keV light curve [c/s/cm^2]
See also: http://maxi.riken.jp/top/index.php?cid=000000000036, http://ads.nao.ac.jp/abs/2009PASJ...61..999M
MJDref_satellite
Synopsis
returns the reference date (MJD) of a satellite's time
Usage
Double_Type MJDref_satellite(String_Type satellite)
MN_NFW
Synopsis
Evaluate equations of motion, total energy, or circular velocity derived from a potential with a Miyamoto & Nagai bulge and disk component and a Navarro, Frenk, & White dark matter halo
Usage
MN_NFW(Double_Types t, m[6,n]; qualifiers)
Description
Evaluate the equations of motion, the total energy, or the circular velocity at time 't' derived from a potential with a Miyamoto & Nagai bulge and disk component and a Navarro, Frenk, & White dark matter halo (see Model III in Irrgang et al., 2013, A&A, 549, A137) using either cylindrical coordinates (r [kpc], phi [rad], z [kpc]) and their canonical momenta vr [kpc/Myr], Lz [kpc^2/Myr], vz [kpc/Myr]) or cartesian coordinates (x [kpc], y [kpc], z [kpc], vx [kpc/Myr], vy [kpc/Myr], vz [kpc/Myr]), see qualifier 'coords'. Conservation of angular momentum Lz is implemented in the equations of motion for cylindrical coordinates only. The total energy E_total [kpc^2/Myr^2] is not used to integrate the equations of motion although being a conserved quantity, too. Therefore, conservation of energy, i.e., of E_total, is a measure for the precision of the numerical methods applied.
For computing orbits with n different initial conditions, the input parameter m is a [6,n]-matrix with (qualifier("coords")=="cyl") or (qualifier("coords")=="cart") m[0,*] = r; m[0,*] = x; m[1,*] = phi; m[1,*] = y; m[2,*] = z; m[2,*] = z; m[3,*] = vr; m[3,*] = vx; m[4,*] = Lz; m[4,*] = vy; m[5,*] = vz; m[5,*] = vz; If the qualifier 'eomecd' is set to "eom", the function returns a [6,n]-matrix delta with delta[0,*] = vr; delta[0,*] = vx; delta[1,*] = Lz/r^2; % = vphi delta[1,*] = vy; delta[2,*] = vz; delta[2,*] = vz; delta[3,*] = -d/dr (Potential(r,z) + Lz^2/r^2); delta[3,*] = -d/dx Potential(x,y,z); delta[4,*] = 0; % -d/dphi Potential(r,z) delta[4,*] = -d/dy Potential(x,y,z); delta[5,*] = -d/dz Potential(r,z); delta[5,*] = -d/dz Potential(x,y,z); If the qualifier 'eomecd' is set to "energy", the function returns a [1,n]-array storing the total energy for each orbit: E_total(r,z,vr,vz,Lz) = Double_Type[n] = 0.5*(vr^2+vz^2+Lz^2/r^2) + Potential(r,z) or E_total(x,y,z,vx,vy,vz) = Double_Type[n] = 0.5*(vx^2+vy^2+vz^2) + Potential(x,y,z) If the qualifier 'eomecd' is set to "circ", the function returns a [1,n]-array storing the circular velocity for each orbit: v_circ(r,z) = Double_Type[n] = sqrt( r * d/dr Potential(r,z) ) If the qualifier 'eomecd' is set to "sgcd", the function returns the Sun-Galactic center distance found to fit best to this potential.
Qualifiers
- coords [
="cyl"
]: Use cylindrical ("cyl") or cartesian ("cart") coordinates. - eomecd [
="eom"
]: Return equations of motion ("eom"), total energy ("energy"), circular velocity ("circ"), or Sun-Galactic center distance ("sgcd"). - Mb [
=439
]: Mass of bulge in Galactic mass units, see Irrgang et al. 2013. - Md [
=3096
]: Mass of disc in Galactic mass units, see Irrgang et al. 2013. - Mh [
=142200
]: Mass scale factor of halo in Galactic mass units, see Irrgang et al. 2013. - bb [
=0.236
]: Bulge scale length, see Irrgang et al. 2013. - ad [
=3.262
]: Disc scale length, see Irrgang et al. 2013. - bd [
=0.289
]: Disc scale length, see Irrgang et al. 2013. - ah [
=45.02
]: Halo scale length, see Irrgang et al. 2013.
Example
delta = MN_NFW(0, m); energy = MN_NFW(0, m; eomecd="energy"); v_circ = MN_NFW(0, m; eomecd="circ"); sgcd = MN_NFW(; eomecd="sgcd");
See also: orbit_calculator, AS, MN_TF, plummer_MW
MN_TF
Synopsis
Evaluate equations of motion, total energy, or circular velocity derived from a potential with a Miyamoto & Nagai bulge and disk component and a truncated, flat rotation curve halo model
Usage
MN_TF(Double_Types t, m[6,n]; qualifiers)
Description
Evaluate the equations of motion, the total energy, or the circular velocity at time 't' derived from a potential with a Miyamoto & Nagai bulge and disk component and a truncated, flat rotation curve dark matter halo model (see Model II in Irrgang et al., 2013, A&A, 549, A137) using either cylindrical coordinates (r [kpc], phi [rad], z [kpc]) and their canonical momenta vr [kpc/Myr], Lz [kpc^2/Myr], vz [kpc/Myr]) or cartesian coordinates (x [kpc], y [kpc], z [kpc], vx [kpc/Myr], vy [kpc/Myr], vz [kpc/Myr]), see qualifier 'coords'. Conservation of angular momentum Lz is implemented in the equations of motion for cylindrical coordinates only. The total energy E_total [kpc^2/Myr^2] is not used to integrate the equations of motion although being a conserved quantity, too. Therefore, conservation of energy, i.e., of E_total, is a measure for the precision of the numerical methods applied.
For computing orbits with n different initial conditions, the input parameter m is a [6,n]-matrix with (qualifier("coords")=="cyl") or (qualifier("coords")=="cart") m[0,*] = r; m[0,*] = x; m[1,*] = phi; m[1,*] = y; m[2,*] = z; m[2,*] = z; m[3,*] = vr; m[3,*] = vx; m[4,*] = Lz; m[4,*] = vy; m[5,*] = vz; m[5,*] = vz; If the qualifier 'eomecd' is set to "eom", the function returns a [6,n]-matrix delta with delta[0,*] = vr; delta[0,*] = vx; delta[1,*] = Lz/r^2; % = vphi delta[1,*] = vy; delta[2,*] = vz; delta[2,*] = vz; delta[3,*] = -d/dr (Potential(r,z) + Lz^2/r^2); delta[3,*] = -d/dx Potential(x,y,z); delta[4,*] = 0; % -d/dphi Potential(r,z) delta[4,*] = -d/dy Potential(x,y,z); delta[5,*] = -d/dz Potential(r,z); delta[5,*] = -d/dz Potential(x,y,z); If the qualifier 'eomecd' is set to "energy", the function returns a [1,n]-array storing the total energy for each orbit: E_total(r,z,vr,vz,Lz) = Double_Type[n] = 0.5*(vr^2+vz^2+Lz^2/r^2) + Potential(r,z) or E_total(x,y,z,vx,vy,vz) = Double_Type[n] = 0.5*(vx^2+vy^2+vz^2) + Potential(x,y,z) If the qualifier 'eomecd' is set to "circ", the function returns a [1,n]-array storing the circular velocity for each orbit: v_circ(r,z) = Double_Type[n] = sqrt( r * d/dr Potential(r,z) ) If the qualifier 'eomecd' is set to "sgcd", the function returns the Sun-Galactic center distance found to fit best to this potential.
Qualifiers
- coords [
="cyl"
]: Use cylindrical ("cyl") or cartesian ("cart") coordinates. - eomecd [
="eom"
]: Return equations of motion ("eom"), total energy ("energy"), circular velocity ("circ"), or Sun-Galactic center distance ("sgcd"). - Mb [
=175
]: Mass of bulge in Galactic mass units, see Irrgang et al. 2013. - Md [
=2829
]: Mass of disc in Galactic mass units, see Irrgang et al. 2013. - Mh [
=69725
]: Mass of halo in Galactic mass units, see Irrgang et al. 2013. - bb [
=0.184
]: Bulge scale length, see Irrgang et al. 2013. - ad [
=4.85
]: Disc scale length, see Irrgang et al. 2013. - bd [
=0.305
]: Disc scale length, see Irrgang et al. 2013. - ah [
=200
]: Halo scale length, see Irrgang et al. 2013.
Example
delta = MN_TF(0, m); energy = MN_TF(0, m; eomecd="energy"); v_circ = MN_TF(0, m; eomecd="circ"); sgcd = MN_TF(; eomecd="sgcd");
See also: orbit_calculator, AS, MN_NFW, plummer_MW
N_body_simulation_MW_kernel
Synopsis
Alternative interaction kernel for the function 'N_body_simulation'
Usage
Double_Type r[6,N] = N_body_simulation_MW_kernel(Double_Types t, m[6,N]; qualifiers)
Description
This function is an alternative interaction kernel for the function 'N_body_simulation'. It combines the mutual N-nody interactions of the standard interaction kernel 'N_body_simulation_std_kernel' with the external forces stemming from an analytical model for the gravitational potential of the Milky Way (see qualifier 'model').
Notes
Because of the unit convention used for the potentials of the Milky Way, 'psa', which is a qualifier of the function 'N_body_simulation_std_kernel' and which is assumed to be the mass of the Plummer spheres in solar masses, has to be converted to Galactic mass units and then multiplied with a constant accounting for the remaining unit conversions (see example below). The units of 'psb' have to be kpc^2.
Qualifiers
- model [
="AS"
]: Function ("AS", "MN_NFW", or "MN_TF"), which evaluates the equations of motion that result from a model for the gravitational potential of the Milky Way. - All qualifiers from the model potential function except 'coords'.
- All qualifiers from the function 'N_body_simulation_std_kernel'.
Example
% Interacting satellite galaxies in the Milky Way: s = properties_satellite_galaxies(); i = struct{ x, y, z, vx, vy, vz, psa, psb }; model = "AS"; % Milky Way mass model SunGCDist = (@(__get_reference(model)))(; eomecd="sgcd"); % Sun-GC distance of chosen mass model temp = [SunGCDist,0,0,0,0,0]; reshape(temp, [6,1]); vlsr = (@(__get_reference(model)))(0, temp; eomecd="circ")[0]; % Local standard of rest velocity of chosen mass model (i.x, i.y, i.z, i.vx, i.vy, i.vz) = cel2gal(s.ah, s.am, s.as, s.dd, s.dm, s.ds, s.dist, s.vrad, s.pma_cos_d, s.pmd; SunGCDist=SunGCDist, vlsr=vlsr); kpcmyr_to_kms = 977.7736364875057; % = conversion factor from kpc/myr to km/s = 3.0856775975*10^16 / (10^6 * 3.15582 * 10^7); i.vx /= kpcmyr_to_kms; i.vy /= kpcmyr_to_kms; i.vz /= kpcmyr_to_kms; i.psa = s.Pl_mass/2.325131802556774e+07; % conversion from solar masses to Galactic mass units Mgal to have G=1 % Mgal = 2.325131802556774e+07 = 1e8*3.0856775975*1e19/6.6742/1e-11/1.9884/1e30, see Irrgang et al., 2013, A&A, 549, A137 const = 100./kpcmyr_to_kms^2; % factor 100 because potential is given in 100 km^2/s^2, see Irrgang et al. 2013 i.psa *= const; i.psb = (s.Pl_radius)^2; r = N_body_simulation(i, -1000; kernel="N_body_simulation_MW_kernel", psa=i.psa, psb=i.psb, model=model); xrange(min_max([r.o0.x,r.o1.x])); yrange(min_max([r.o0.y,r.o1.y])); plot(r.o0.x,r.o0.y); oplot(r.o1.x,r.o1.y);
See also: N_body_simulation, N_body_simulation_std_kernel, AS, MN_NFW, MN_TF
RXTE_ASM_countrate
Synopsis
estimates the RXTE-ASM countrate from the current model
Usage
(A, B, C) = RXTE_ASM_countrate();
Description
Zdziarski et al. (http://adsabs.harvard.edu/abs/2002ApJ...578..357Z)
provide a "response matrix" for the RXTE-ASM.
The inverse of this matrix is applied to the energy flux
derived from the current fit-function and its parameters,
in order to estimate the RXTE-ASM count rates A
, B
and C
(cps).
Note that these numbers may only give a rough estimate!
See also: energyflux
RXTE_ASM_lightcurve
Synopsis
retrieves ASM lightcurves for a given source
Usage
Struct_Type RXTE_ASM_lightcurve(String_Type sourcename);
Qualifiers
- MJDmin: earliest MJD to be used
- MJDmax: latest MJD to be used
- dt: if specified, time resolution in MJD for rebinning
- no_filter_nan: do not remove empty bins after rebinning (only with
dt
) - list: lists available sources;
sourcename
may be omitted but can be a regular expression - get_list: as list, but the list of sources is returned as an array of strings
- save: saves the light curve data in a local FITS file
- verbose
RXTE_PCA_info
Synopsis
retrieves light curves for a given RXTE-PCA observation
Usage
Struct_Type RXTE_PCA_info(String_Type dirs[])
Description
The elements of dirs
may contain globbing expressions.
The returned structure has the following fields:
-
obstime = struct { start, stop}
with the times for each observation indirs
lc = struct { time, rate, error }
with the PCA light curve
obscat = struct
with information from the observation catalogue (occultation, saa, good time)
gti
xfl = struct
with information from the filter file
See also: aitlib/rxte/readxtedata.pro
RXTE_PCA_modes
Usage
RXTE_PCA_modes(String_Type obsids);
Qualifiers
- compact: only one row per ObsID, omit TSTART, TSTOP and duration
- get_struct: return information on PCA modes in a structure, too
- quiet: do not show information, implies
get_struct
Description
Reads RXTE-PCA data modes from the PCA FITS-index (FIPC) file.
The location of the RXTE data archive is determined from the
local_paths.RXTE_data
variable (defined within the isisscripts).
RXTE_filter_file_info
Synopsis
returns RXTE filter file information
Usage
Struct_Type RXTE_filter_file_info(String_Type xflfiles[])
See also: aitlib/rxte/readxfl.pro
RXTE_nr_PCUs_from_filename
Synopsis
counts how many PCUs were off from an *xyoff_excl* filename
Usage
Integer_Type RXTE_nr_PCUs_from_filename(String_Type filename)
RXTE_nr_PCUs_from_filterfile
Synopsis
returns the number of PCUs switched on over time
Usage
Integer_Type[] RXTE_nr_PCUs_from_filterfile(String_Type filterfile[, Double_Type[] time])
Description
The number of PCUs switched on during an observation may vary, which results in jumps in the lightcurve. The filter file (your_extraction/filter/*.xfl) provides time resolved information about the operating PCUs, which is read out and returned.
ATTENTION: If no time array is given the number of PCUs is returned for the full length of the observation (no GTIs applied!). In the other case the given time array HAS TO be in the SATELLITE TIME SYSTEM and in SECONDS since RXTE started operating (see 'MJDref_satellite').
See also: MJDref_satellite, RXTE_nr_PCUs_from_filename
RXTE_obscat_info
Synopsis
reads an RXTE human-readable obs(ervation)cat(alogue)
Usage
Struct_Type RXTE_obscat_info(Integer_Type obscat_days[]);
See also: aitlib/rxte/readobscat.pro
RXTE_plot_PCA_info
Synopsis
plots an overview of an RXTE observation (lc, GTI, SAA, bkg)
Usage
RXTE_plot_PCA_info(String_Type dirs[]);
or
Struct_Type RXTE_plot_PCA_info(String_Type dirs[]; get_info)
Qualifiers
- electron: set to electron ratio threshold for data extraction
- noback: set if no bkg subtraction is to be performed
- earthvle: set if EarthVLE background model is to be used
- faint: set if Faint background model is to be used
- q6: set if Q6 background model is to be used (default is to test for earthvle,faint,q6)
- skyvle: set if SkyVLE background model is to be used (default)
- exclusive: set to search for data that was extracted with the exclusive keyword to pca_standard being set.
- top: set to read top-layer data
- nopcu0: set to search for data that was extracted ignoring PCU0
- fivepcu: plot count-rates wrt to whole PCA, i.e., normalizing to 5 PCUs. Default is to plot the average countrate per PCU.
- charsize_obsid
- get_info: returns the info structure
Description
The elements of dirs
may contain globbing expressions.
See also: aitlib/rxte/rxtescreen.pro
SRT_get_mw_spectra
Synopsis
computes Milky Way spectrum from SRT data
Usage
Struct_Type[] spectra = SRT_get_mw_spectra(Struct_Type data, Integer_Type[] chunks
Description
This function requires as input the SRT data in the format given by SRT_read() and an array of integer numbers which are the chunks with the spectrum data. If an array of chunk numbers is given then all corresponding spectra are calculated and returned. The spectra are determined with SRT_spectrum and all qualifiers of this function can be used.
Example
variable data = SRT_read("milkyway.rad"); variable spectra = SRT_get_mw_spectra(data,[1,4,7]);
See also: SRT_spectrum, SRT_read
SRT_image
Synopsis
extracts an image from a SRT data structure containing an npoinscan
Usage
Float_Type img[] = SRT_image(Struct_Type data);
Qualifiers
- plot: plots the image
See also: SRT_read
SRT_plot_mw_spectra
Synopsis
plots Milky Way spectra from SRT data
Usage
SRT_plot_mw_spectra(Struct_Type[] spectrum);
or
```c
SRT_plot_mw_spectra(Struct_Type[] data, Integer_Type[] chunks);
##### Qualifiers
* pad: additional space between plot frame and data limits
* fMax: give the frequency of the cloud to be plotted, has to be of same length as spectra/chunks
* fConfMax: give the upper confidence leve of the maximal frequency of the cloud to be plotted, has to be of same length as spectra/chunks
* fConfMin: give the lower confidence leve of the maximal frequency of the cloud to be plotted, has to be of same length as spectra/chunks
* fMaxLin: type of line used to plot fMax
* fConfLine: type of line used to plot fConfMin and fConfMax
* xlabel: specify label for x-axis
* ylabel: specify label for y-axis
* title: give title for plot either as single title or for each plot individually as an array of strings
* name: file name of plot, either a single one or for each plot as a an array of string, specify file format, i.e. mw_scan.pdf
* return_xfig: return xfig-data instead of plotting
##### Description
The function automatically plots Milky Way spectra obtained with the SRT.
It can be given a single or an array of spectra or the data and an array of chunk numbers. The required data format is based on the output of SRT_get_mw_spectrum and all of its qualifiers can be used.
It is possible to give the previously determined maximal frequency of the cloud and its upper and lower confidence interval as qualifier. This will create lines in the plot at the provided frequency.
__See also__: SRT_get_mw_spectra, SRT_read, SRT_spectrum
#### SRT_read
##### Synopsis
reads SRT data structures from a SRT .rad file
##### Usage
```c
Struct_Type data[] = SRT_read(String_Type filename);
Qualifiers
- verbose
- onechunk: read all data into one structure instead of an array (default if file has no comments)
- bins_to_cut[=8]: number of bad bins at high and low frequencies
- position: [lat,longw] position of the SRT (default: lat=49.90 and longw=349.10)
See also: SRT_spectrum, SRT_image
SRT_show_chunks
Synopsis
prints information for each chunk in an array of SRT data structures
Usage
SRT_show_chunks(Struct_Type data[]);
See also: SRT_read
SRT_spectrum
Synopsis
extracts a spectrum from a SRT data structure
Usage
Struct_Type spec = SRT_spectrum(Struct_Type data);
Qualifiers
- verbose
- bins_to_cut [
=8
]: number of bins at both sides of the spectrum that are discarded - bins_cut_low [
=0
]: number of bins to additionaly cut from lower side - bins_cut_high [
=0
]: number of bins to additionaly cut from upper side - normalize:
spec.value
is normalized between 0 and 1. - vLSR: transform the frequency grid to a velocity grid
- filename=name: plots the spectrum as name.ps in the current directory
Description
See also: SRT_read
XSTAR_read_pops
Synopsis
reads the ionization balances from an XSTAR population file
Usage
Struct_Type XSTAR_read_pops(String_Type filename)
Qualifiers
- verbose
- Z: array of Z values of elements to include
- ions: array of ions to include (overrides Z qualifier)
aStar
Synopsis
pathfinding algorithm A*
Usage
Struct_Type aStar(Double_Type[][] graph, Integer_Type startX, startY, endX, endY);
Qualifiers
max - find the most expensive path instead of the cheapest one meanCost - costs are normalized by the best path length to any point in the graph, resulting in a mean cost on a path estimate - reference to a function, which estimates the cost between two nodes in the graph (parameters: x0, y0, x1, y1). By default, the distance between both nodes is used plot - plots the graph and the working progress at each step. Nodes in the open list are shown in red, in the closed list in blue, nodes of infinite cost in green and the best path, finally, in black. The function sleeps 0.01 seconds after each step or the value given with this qualifier
Description
Pathfinding algorithm with high performance and accuracy by Hart, Nilsson & Raphael, "A Formal Basis for the Heuristic Determination of Minimum Cost Paths", IEEE Transactions on Systems Science and Cybernetics SSC4 (2), 100–107, 1968
The A* algorithm (called A Star) finds the shortest (or in general cheapest) path between to nodes in a graph. The algorithm looks for the best way following the lowest known heuristic cost. This cost has to be estimated at each new discovered node and is, by default, the straight-line distance between both points. The runtime and accuracy of the search strongly depends on the choice of this estimation.
The 'graph' has to be given as 2d-array defining the cost at each point (node). Infinite costs (_Inf) are allowed and corresponds to insuperable borders. The starting node and goal node habe to be given as array indices 'startX', 'startY' and 'endX', 'endY', respectively.
The returned structure contains x - x-indices of the best path y - y-indices of the best path cost - summed costs along the best path
aglc
Synopsis
Bin a light curve from ACIS grating events.
Usage
lc = aglc(tgevt, expno_ref, tbin, wmin, wmax, tg, orders[, bkg]);
or
lc = aglc( evt_struct, tbin, wmin, wmax, tg, orders[, bkg]);
Description
Given an event file name (tgevt
) and an exposure reference file
(unfiltered events or ``stat1'') name (expno_ref
), bin a light curve
for the specified time bin in seconds (tbin
), wavelength ranges
specified by wmin
, wmax
, in Angstroms. If wmin
and wmax
are arrays,
then they must be the same length and represent low-high pairs of
wavelength regions. The argument, tg
, is a scalar or array of
grating names, and should be one or more of "HEG"
and "MEG"
, or "LEG"
(case-insensitive); only the first character is necessary. The
argument orders
should be an array of integers specifying the
diffraction orders to bin (excluding zero). If the last argument is
present, then events are binned from the background region instead of
the source region.
In the second form, an event structure as returned by
aglc_read_events
is used instead of two file names. This is
more efficient for multiple calls, since it avoids multiple file
reads.
The two files required in the first form are the grating coordinates event file, probably filtered of bad events; i.e., typically the file from which a spectrum would be binned ("Level 1.5" or "Level 2"). The second file, the exposure reference file, should be unfiltered events or the ``stat1'' file. It is used to count the exposed frames to determine the exposure; any event - cosmic ray, bad pixel, photon
- will suffice to mark a frame. Use of the ``stat1'' file is more efficient, since it has only one entry for each unique frame.
Standard binning regions are applied for source and background
events. They may be changed with aglc_set_regs
.
The return value is a structure of the form
time = Double_Type[]
% event time, bin center, since timezero.
time_min = Double_Type[]
% lower edge of time bin, seconds since timezero.
time_max = Double_Type[]
% upper edge of time bin, seconds since timezero.
counts = UInteger_Type[]
% counts per time bin.
count_rate = Double_Type[]
% Counts per second ( counts / exposure ).
stat_err = Double_Type[]
% Statistical error ( sqrt(counts) ).
exposure = Double_Type[]
% Exposure per bin, averaged over all chips included, in seconds.
timezero = Double_Type
% reference time, in seconds since MJDREF.
mjdref = Integer_Type
% Modified Julian Day reference point of Chandra data (from event header)
revidx1a = Array_Type[]
% reverse index array for the grating events.
revidx = Array_Type[]
% reverse index array for the exposure reference events.
The times are in seconds since timezero. Timezero is the number of seconds since 1998.0, which is also the reference MJD.
The count_rate
is per second, and is equivalent to counts/exposure.
The exposure is computed properly for event selections spanning
multiple CCDs with possibly different GTI tables.
The reverse indices, revidx1a
, are returned by the histogram
function; for each bin of the light curve they point to
the items in the event list which have been binned. This facilitates
subsequent operations on events in specific bins, especially when
selected on non-time criteria. (see aglc_filter
).
revidx
is like revidx1a
, but for the exposure record event list.
Example
Compute the light curve for the sum of Ne X 12A and O VIII 19A:
wlo = [12.1, 18.8];
whi = [12.2, 19.1];
g = ["H", "M"];
o = [-1,1];
c = aglc( "evt2.fits", "evt1.fits", 1000, wlo, whi, g, o );
hplot(c.time_min, c.time_max, c.counts);
See also: aglc_read_events, aglc_filter, aglc_set_regs
aglc_get_ephem
Synopsis
Retrieve a stored ephemeris, the last one used.
Usage
(jd0, pd) = aglc_get_ephem;
Description
jd0 will be the Julian day of zero phase. pd will be the period in days.
See also: aglc_pr_ephem, aglc_phased
aglc_get_regs
Synopsis
Retrieve stored source and background regions (cross-dispersion limits), and associated BACKSCAL values.
Usage
r = aglc_get_regs;
Description
Retrieves a structure containing the stored values used for binning light curves for source or background.
See also: aglc_set_regs, aglc_reset_regs
aglc_get_trange
Synopsis
Retrieve the time filter which will be applied in phase-binning.
Usage
(tmin, tmax) = aglc_get_trange;
Description
tmin,tmax are relative times in seconds, since the start of the observation.
See also: aglc_set_trange, aglc_pr_trange;
aglc_jd_to_rotnum
Synopsis
Convert Julian days to rotation number.
Usage
rotation_number = aglc_jd_to_rotnum( t_jd, hjd_epoch, period_days );
Description
rotation_number
= number of full cycles plus fractional phase, given an ephemeris.
t_jd
= times, in Julian day numbers.
hjd_epoch
= epoch of zero-phase, in HJD
period_days
= period, in days.
See also: aglc_tobs_to_jd, aglc_jd_to_phase, aglc_tobs_to_phase, aglc_phased, aglc
aglc_phased
Synopsis
Bin a phased light curve.
Usage
pc = aglc_phased(tgevt, expno_ref, phase_info, ephem, wmin, wmax, tg, orders[, bkg]);
or
pc = aglc_phased( evt_struct , phase_info, ephem, wmin, wmax, tg, orders[, bkg]);
Description
Arguments are similar to those in aglc
, with the exception of
phase_info
and ephem
. For other arguments, see aglc
.
phase_info
= a 3-element array giving the phase minimum, maximum,
and bin size desired for the resultant phased light
curve.
ephem
= a two-element array giving the JD of zero phase (NOT MJD!)
and the period in days.
Definitions for Phase Curve Structure Fields The phase curve fields are the same as for the light curve, with the exception of phase replacing time and the addition of the ephemeris. The new fields are:
-
phase
: Value of phase at the bin center. -
phase_min
: Value of phase at the low edge of the bin. -
phase_max
: Value of phase at the high edge of the bin. -
epoch
: Epoch of the ephemeris, in Julian Days (NOT MJD!). -
period: Period of the ephemeris, in days.
See also: aglc, aglc_filter, aglc_read_events
aglc_pr_ephem
Synopsis
Print the stored ephemeris used in the last phase-curve binning.
Usage
aglc_pr_ephem
Description
Retrieves the ephemeris stored by aglc_phased and prints to terminal.
See also: aglc_get_ephem, aglc_phased
aglc_pr_trange
Synopsis
Print to the terminal the current time filter.
Usage
aglc_pr_trange;
Description
Prints the min and max time filter used in binning phase curves.
See also: aglc_set_trange, aglc_get_trange;
aglc_read_events
Synopsis
Read subset of events from grating file and from exposure file into a structure.
Usage
Struct_Type s = aglc_read_events(String_Type tgevt, expno_ref);
Description
tgevt
should be a Level 1.5 (or 2) grating events file name.
expno_ref
should be a Level 1 ``stat'' or events (unfiltered) file name.
The returned value is a structure. The detailed contents of the
structure are given below in the example. While the exposure can
be computed from the EXPNO column of the unfiltered event file,
use of the `stat1' file is more efficient.
Definitions for Event Structure Fields:
-
fevt_1a
: File name for grating events. -
fevt_exp
: File name for exposure reference (unfiltered events or ``stat1'' file) -
expno[]
: Exposure number column from exposure reference file. -
ccd
: CCD_ID column from exposure reference file. -
expno_1a
: Exposure number column from grating event file. -
ccd_1a
: CCD_ID column from grating event file. -
tgpart
: TG_PART column from grating event file. (value is 1 for HEG, 2 for MEG, 0 for zero order, 99 for background) -
order
: TG_M column from grating events file, giving the diffraction order. -
tgd
: TG_D column from grating events file. This is the cross-dispersion coordinate in degrees. -
wave
: TG_LAM column from grating events file, giving the wavelength of each event. -
timedel
: TIMEDEL keyword from grating event file, needed for scaling exposure numbers to time. (See the ahelp on time.) -
timepixr
: TIMEPIXR keyword from event file. (See the ahelp on time.) -
frame_exp
: the exposure time per frame (EXPNO). (See the ahelp on time.) -
mjdref
: Modified Julian Day reference time for Chandra observations, from the event file header. -
tstart
: TSTART, from the event file header, is the observaton start time in seconds since MJDREF. -
status
: Event status column, from the grating events file. (See the definition of ACIS status bits.) -
time
: The TIME column from the exposure reference event file. -
cycle
: unused.
See also: aglc_filter, aglc, aglc_phased
aglc_set_regs
Synopsis
Set the spatial binning regions (cross-dispersion) for the source and backgrounds on both sides, and compute the background scaling factors.
Usage
aglc_set_regs( s_min, s_max, b_lo_min, b_lo_max, b_hi_min, b_hi_max );
Description
Limits are given in cross-dispersion coordinates, tg_d, in units of degrees.
s_min
= source minumum tg_d.
The center of the source region is defined as tg_d = 0.
s_max
= source region maximum tg_d.
b_lo_min
= background minimum, on the "low" side (tg_d < 0).
b_lo_max
= background maxmimum, on the "low" side (tg_d < 0).
b_hi_min
= background minimum, on the "high" side (tg_d < 0).
b_hi_max
= background maxmimum, on the "high" side (tg_d < 0).
The limits follow the constraint that
b_lo_min < b_lo_max <= s_min < s_max <= b_hi_min < b_hi_max
.
Default values are the same as the default spectral extraction regions of tgextract:
s_min = -6.6e-04
[ deg ]
s_max = 6.6e-04
[ deg ]
b_lo_min = -6.0e-03
[ deg ]
b_lo_max = s_min
[ deg ]
b_hi_min = s_max
[ deg ]
b_hi_max = 6.0e-03
[ deg ]
See also: aglc_reset_regs, aglc_get_regs
aglc_set_status_filter
Synopsis
Set a flag to specify if the grating events should be filtered on STATUS=0.
Usage
aglc_set_status_filter( flag );
Description
By default, the status filter is zero, meaning all good events. Different status bit fields qualify attributes of "bad" pixels. If, for some reason, the grating events have not been filtered on status, but should be, then setting this flag causes all events with non-zero status to be ignored.
If flag=1, then ignore events with non-zero status. If flag=0, then accept all events.
See also: aglc, aglc_phased
aglc_set_trange
Synopsis
Set a simple time range filter to be applied to a light curve before extracting a phase-binned curve. This is useful for excluding flares, which seem to predominantly occur at the beginning or end of an observation.
Usage
aglc_set_trange( tmin, tmax );
Description
tmin, tmax
are relative times in seconds, since the start of the
observation. The interval specifies the times of interest to be included
in phase binning. tmax==NULL
means to the end of the observation.
Example
aglc_set_trange( 0, 2000 ); %
first 2ks will be included in phase binning.
aglc_set_trange( 5000, NULL); %
from 5ks to the end will be selected.
See also: aglc_get_trange, aglc_pr_trange
aglc_tobs_to_jd
Synopsis
Convert Chandra obervation times, given in seconds since MJDREF, to Julian day.
Usage
jd = aglc_tobs_to_jd( t_obs, mjdref );
Description
t_obs
is a scalar or array of Chandra times in seconds since mjdref
.
mjdref
is the Modified Julian Day reference, found in Chandra
headers, or in the aglc event structure.
See also: aglc_jd_to_rotnum, aglc_jd_to_phase, aglc_tobs_to_phase, aglc_read_events
aglc_write_curve
Synopsis
Write a light or phase curve to a FITS bintable file
Usage
aglc_write_curve( outfile, hdr_ref_file, curve_struct [, history_string ] );" );
Description
outfile
= name of output file ( String_Type
).
hdr_ref_file
= Reference file for copying of header info ( String_Type
).
curve_struct
= The light or phase curve structure, as created by aglc
or aglc_phased
.
history_string
= arbitrary note (String_Type
).
See also: aglc, aglc_phased
atom_name
Synopsis
returns the symbol for atoms with proton number Z
Usage
String_Type atom_name(Integer_Type Z)
See also: element_name
bayesian_blocks
Synopsis
find bayesian blocks in the given data
Usage
Struct_Type bayesian_blocks(Struct_Type[] data);
Qualifiers
fp_rate - value of the false positive rate, which is the probability that a found change point with the current value of ncp_prior is actually not significant (default: .01) WARNING: This has no influence on data of type 3, where the default ncp_prior, fp_rate is 0.05, and it is not possible to change this currently. ncp_prior - controls the tolerance level used to find the change points. The value has to be an array of estimates for the order of number of blocks for each given dataset. The probability P that the final number of blocks N_blocks is then given by log(P) ~ ncp_prior * N_blocks Thus, ncp_prior has to be decreased to increase the final number of blocks. The default value is determined by a formula derived from simulations, given by ncp_prior = 4 - log(73.53*fp_rate*N^-0.478) for data modes 1 and 2 and ncp_prior = 1.32 + 0.577*log(N) for data mode 3, where N is the number of data points. do_iter - number of iterations on ncp_prior (default: 0) dt_min - events within the given time range are considered to be simultaneous and a new combined event at the average time is added (default: 0) num_skip - maximal number of consecutive simultaneous events (time difference < dt_min) gti - structure containing the GTIs for the given events (data mode 1) with the fields 'start' and 'stop'. The blocks are searched for each GTI individually, thus the function returns an array structures! If GTIs are given only ONE input dataset is allowed (at the moment). gti_mindt - all GTIs shorter than the given duration are ignored (i.e. no blocks are searched) (default: 0)
Description
This function is an S-lang implementation of the 'find_blocks_mult' MatLab-program described in Scargle et al., 2013, ApJ 764, 167 which source code is available on the arXiv.
WARNING: The full functionality has not been tested yet, thus please use with caution!
Depending on the fields of the input 'data' structure, the best Bayesian block representation is determined. The available data modes are: 1: Event data Double_Type[] tt - time tags of events 2: Binned data (data without specific uncertainties, e.g., counts vs. time) Double_Type[] tt - time tags Double_Type[] nn_vec - binned data at tt WARNING! Mode 2 possibly has bugs! Use with extreme caution! 3: Point measurements (data with known uncertainties, e.g., radio flux vs. time) Double_Type[] tt - time tags List_Type[] cell_data - cell data at tt defined as { Double_Type[] data, uncertainties }
The output structure contains change_points - time INDICES of the change points tt_all - combined time tags of all datasets ncp_prior_vec - ncp_prior value of each dataset data_matrix - see Fig. 2 of Scargle et al. 2013 last - index of the last change point over tt_all best - fitness (=statistics) over tt_all index_vec - dataset index over tt_all data_mode_vec - data mode of each dataset
The algorithm can handle multiple datasets, which are combined internally to one dataset (tt_all). Zero blocks are handled as well.
NOTE: in case of events (data mode 1) it is important to provide the GTIs via the corresponing qualifier as well. Otherwise the block statistics are probably wrong and thus the total block represenation!
To finally retrieve the block represenatation of the input data, the 'get_blocks_data' function can be used.
See also: get_blocks_data
bayesian_blocks_cperror
Synopsis
calculates the probability distribution of the change point positions
Usage
Struct_Type[] bayesian_blocks(Struct_Type blocks);
Description
After Scargle et al., 2013, ApJ 764, 167 the uncertainty of the position of a change point in a bayesian block representation can be estimated by calculating the fitness of a block with variing change point position, while all other change points are kept fix. Although inter-change-point dependencies are neglected in this approach, it can be used to derive proper uncertainties as long as the resulting distributions do not overlap.
This function is an implementation of this procedure based on the MatLab-programs 'figure_cp_error' and 'cp_prob' by the authors above.
The input 'blocks' have to be calculated previously by the 'bayesian_blocks' function. The returned structure array contains the probability distribution 'prob' (integral = 1.) over the time-tags 'tt' for each change point (i.e., the array index corresponds to the change point index). From each individual distribution an uncertainty for the location of the change point can be estimated via, e.g., its width.
See also: bayesian_blocks
bbootstrap
Synopsis
1-parameter bootstrap calculation of function 'funct'
Usage
Array_Type bbootstrap(Array_Type data, Integer_Type num_samp
[, funct = &mean, datatype = Double_Type]);
Qualifiers
slow - the balanced bootstrap is performed. This is more accurate, especially for long-tailed distributions, but is also much slower and requires more memory.
Description
This function has been ported from the IDL-routine bbootstrap.pro by H.T. Freudenreich, HSTX, 2/95
This program randomly selects values to fill sets of the same size as 'data'. Then calculates the 'funct' of each, which is an optional function reference with the &mean as default. It does this 'num_samp' times. The 'datatype' returned by 'funct' has to be provided to initialize the array returned by this function. If the slow-qualifier is set, the balanced bootstrap method is used to obtain the samples, hence the name 'bbootstrap'. This method requires, however, more virtual memory, and patience.
The user should choose 'num_samp' large enough to get a good distribution. The sigma of that distribution is then the standard deviation of the mean of the returned 'funct'-vector. For example, if input 'data' is normally distributed with a standard deviation of 1.0, the standard deviation of the returned 'funct'-vector (by default the mean) will be ~1.0/sqrt(N-1), where N is the number of values in 'data'.
WARNING: at least 5 points must be input. The more, the better.
bgsubHEXTElc
Synopsis
calculates a background subtracted HEXTE lightcurve
Usage
Struct_Type lc = bgsubHEXTElc(Struct_Type srclc, Struct_Type bglc);
or
Struct_Type lc = bgsubHEXTElc(Struct_Type srclc, Struct_Type bglc1, Struct_Type bglc2);
Description
HEXTE rocking -> Background is not measured simultaneously with source data.
circle_from_points
Synopsis
find the best fitting circle for a set of points
Usage
pars=circle_from_points(Array_Type x,Array_Type y);
Description
Given arrays of x/y coordinates of points which lie approximately on a circle, this function returns the center coordinate and radius of the circle. This is a numerically unstable procedure, but the routine generally returns parameters which are good enough to start a more advanced fitting routine. The function is an implementation of ideas described by L. Maisonobe in a document entitled "Finding the circle that best fits a set of points" (see http://www.spaceroots.org/downloads.html): For each tuple of three points of coordinates it finds the center and radius, which is possible in an exact manner, assuming that the points are not on a straight line. This is done for all or a subset of the points, and the resulting radius and center coordinates are then averaged. The function returns a structure with the tags xc and yc (x- and y- coordinate of the center of the circle) and radius (radius of the circle), or NULL if no solution could be found (points are on a straight line).
Qualifiers
- eps: if the determinant of the three point solution is less than than this number, the points lie on a straight line
- maxiter: average at most maxiter tuples of three points
compare_modelfits
Synopsis
to compare modelfit parameters of different epochs
Usage
compare_modelfits([array of <code>epoch fitsfiles</code>]);
Qualifiers
- plottype [="overlay"]: specifies the layout of the comparison plot, by default it shows the overlay of all component positions, alternatives are "distancevsmjd","fluxvsdistance","TBvsdistance"
- fileextent [=.eps]: output is given by default as .eps-file in working directory
- outputfile : give outputfile name and directory by hand
- sourcename[=default]: set the name of the source, by default the name is read from the .fits file,
set to NULL for not plotting a source name
- ra_mas [=[20,-20]] specifiy ra range (for "overlay")
- dec_mas [=[-20,20]] specifiy dec range (for "overlay")
- mjd_min [=54101.0]: to specify lower limit of time axis (default: 1/1/2007)
- mjd_max [=55927.0]: to specify lower limit of time axis (default: 1/1/2012)
- distance_min [=0]: lower limit of distance axis in mas
- distance_max [=50]: upper limit of distance axis in mas
- flux_min [=1e-5]: lower limit of flux axis in Jy
- flux_max [=2]: upper limit of flux axis in Jy
- TB_min [=10^5]: lower limit of brightness temperature axis in K
- TB_max [=10^15]: upper limit of brightness temperature axis in K
- linestyle [=default]: set to 0 if the flux values should not be conntected (or to another value for another line style
- ex_comps [=0]: set to 1 if components should be excluded and define the corresponding quadrant of the plot via the ex_comps-coordinates
- ex_comps_ra [=NULL]: define quadrant coordinates in mas where components should be excluded
- ex_comps_dec [=NULL]: define quadrant coordinates in mas where components should be excluded
Description
This function creates an overlay of a VLBI-clean image and the corresponding
model of Gaussian components which are over-plotted as ellipses.
The required input format are fits-file for both the clean and the modelfit images. The format of the output file depends on the suffix of the given
filename
. Possible formats of the output file are PDF, EPS,
PNG, GIF, etc. If labelling = 1 is set, the components are labeled with J0,J1,J2,... depending on their distance to [0,0]. The counterjet components (and core) can be excluded by defining the corresponding quadrant of the the plot via ex_comps_ra and _dec.
difmap_restore
Synopsis
DIFMAP is used to restore a radio image with a new beam
Usage
String_Type outname = difmap_restore(String_Type <code>fitsfile</code>, Double_Type <code>smajor</code>, <code>sminor</code>, <code>pos_angle</code>)
Qualifiers
- chan [=i]: choose channel (i,q,u)
- xsize [= key "NAXIS1"]: number of RA pixels
- xsetp [= key "CDELT1"]: step size of each RA pixel in mas
- ysize [= key "NAXIS2"]: number of DEC pixels
- ysetp [= key "CDELT2"]: step size of each DEC pixel in mas
- uvtaper [= " "]: apply a uvtaper before restoring
- outname [= "xxxxMHz_restore.fits"]: filename of the output, be default the frequency is read from the input file
- overwrite: overwrite file with outname if it already exists
Description
This function creates a radio image with a given beam using DIFMAP. The new beam has to be defined by its semimajor and semiminor axis and its position angle. The new images is saved in the current working directory. It is assumed that the UVF and MOD files have the same basename as the provided file name if the FITS image.
See also: make_spix
element_name
Synopsis
returns the name of the element with proton number Z
Usage
String_Type element_name(Integer_Type Z)
See also: atom_name
enlarge_image
Synopsis
enlarge an image to subpixel resolution by 2d interpolation
Usage
Double_Type IMG = enlarge_image(img, Integer_Type n);
or
Double_Type IMG = enlarge_image(img, Integer_Type nx, ny);
Qualifiers
- interp: reference to interpolation function (see below)
- y_first: interpolate first in y-, then in x-direction
Description
The pixels of img
are mapped to every nx
-th pixel
in x-direction and every ny
-th pixel in y-direction
(in the first usage, nx = ny = n) of IMG
:
IMG[ [::ny], [::nx] ] = img;
As no extrapolation is performed at the boundary,
width and height of IMG
are smaller than nx*w
and ny*h
,
where w
and h
are width and height of img
.
Intermediate pixels are interpolated using the function
specified by the interp
qualifier, which defaults to
gsl->interpol_cspline
if the gsl module is available,
and otherwise to ISIS' interpol
function. In general,
it can be a reference to any function of the form
new_y[] = interpol (new_x[], old_x[], old_y[]);
All qualifiers of enlarge_image
are passed to @interp
.
The interpolation is first performed in x-direction
and then in y-direction -- unless the y_first
qualifier
is specified. For linear and cubic spline interpolation,
the final result is independent of the order.
ext_info_string
Synopsis
converts the ext_line_info into a string
Usage
String_Type ext_info_string(Struct_Type info[, Integer_Type format])
Description
format=0 => string [default]
format=1 => PGPLOT
format=2 => TeX
See also: ext_line_info, ext_info_string
ext_info_string_PGPLOT
Synopsis
converts the ext_line_info into a PGPLOT string
Usage
String_Type ext_info_string_PGPLOT(Struct_Type info)
See also: ext_line_info, ext_info_string
ext_line_info
Usage
Struct_Type info = ext_line_info(Integer_Type id);
or
Struct_Type info = ext_line_info(Integer_Type Z, Integer_Type ion, Integer_Type nr);
See also: line_info
fit_gauss_to_img_noise
Synopsis
fits a gaussian profile to the distribution of pixel values
Usage
mu,sigma = fit_gauss_to_img_noise(<code>img</code>)
Qualifiers
- grid_scale [="log"]: change between fitting on "lin" or "log" grid
- cut_nsig [=NULL]: set to
N
in order to ignore values aboveN
sigma after the first iteration - keep_data: do not delete loaded data after fitting
Description
This function fits a Gaussian profile to the distribution of
pixel values in an image (arrays of different dimensions can be
given) and returns the center mu
and the width sigma
of the
best fit profile.
If the majority of pixels in the image include only noise values,
the obtained mu
and sigma
values characterize the mean value
(base level) of the noise and its amplitude.
NOTE: If the data is kept and further fitting is performed, the following
relation has to be considered mu=get_par("gauss(1).center")+min(img)-1e-5
(necessary to provide a proper grid for fitting).
Example
img = grand(500*500); % image with random numbers around 0 with sigma=1 reshape(img,[500,500]); (mu,sigma) = fit_gauss_to_img_noise (img);
See also: plot_vlbi_map
galactic_rotation_velocity
Usage
Double_Type galactic_rotation_velocity(Double_Type r; model=m)
Qualifiers
- model:
m=1
selects Model I (B&T, Fig. 2.20),m=2
selects Model II (B&T, Fig. 2.22) If neitherm=1
norm=2
is specified, the constant 220 is returned.
Description
r
is the distance from the Galactic Center in kpc.
(0.18 <= r <= 11.9
is covered by a lookup table.)
The return value is the Galactic rotation velocity in km/s.
See also: J. Binney & S. Tremaine: Galactic Dynamics (2nd ed.), Sect. 2.7 (pp. 110ff)
get_blocks_data
Synopsis
returnes the data of a Bayesian block representation
Usage
Struct_Type blockdata = get_blocks_data(Struct_Type blocks[, Integer_Type dataindex]);
or ... = get_blocks_data(Struct_Type[] blocks[, ...]);
Qualifiers
datafun - reference to a function to be called to derive the block data (default: see text) errfun - same as 'datafun', but for the block uncertainties bydt - divides the output data and uncertainties by the length of the corresponding block, can be assigned to a number used to divide the data to for further normalization
Description
Returnes the data of all block over time (e.g., a lightcurve) with variable binsize from a previous run of 'bayesian_blocks'. The default output structure has the fields: time - lower time bin of a block dt - length of each block data - the data in the blocks and error - and the corresponding uncertainties The data and their uncertainties are calculated from the original input data by functions (if multiple dataset are present use the optional 'dataindex' parameter), which will be called for each block and can be specified via qualifiers. Their calling sequence is Double_Type function(time, dt, data[, error]) where error is given in case of data mode 3 (point measurements) only. The default functions depend on the data mode used previously: 1: Event data datafun = sum(events) errfun = sqrt(sum(events)) 2: Binned data datafun = sum(nn_vec) errfun = sqrt(sum(nn_vec)) 3: Point measurements datafun = weighted_mean(cell_data[0], cell_data[1]; err) errfun = sqrt(sumsq(datafun - cell_data[0]) / (length(cell_data[0]) - 1) / length(cell_data[0]))
Example
% let 'evts' be a list of time-tagged events, % find bayesian blocks blocks = bayesian_blocks(struct { tt = evts });
% create a lightcurve of the block representation, % 'bydt' will convert counts to rate lc = get_blocks_data(blocks; bydt);
% plot the blocks hplot(lc.time, lc.time+lc.dt, lc.data);
See also: bayesian_blocks
get_component
Synopsis
to extract the component from the epoch files
Usage
get_component([array of position number in the epoch file], [array of epoch fits files]);
Qualifiers
- quadrant [=1]: Quadrant which is defined as a positive distance ( RA, DEC > 0 == 1 || RA < 0, DEC > 0 == 2 || RA < 0, DEC < 0 == 3 || RA > 0, DEC < 0 == 4)
- zero: shift all distances, so that the last component is at 0,0 ("core" component)
Description
This function returns a structure of a component. The required input is the position number starting by 1 in each epoch, and a corresponding list of epoch fits files. If a component is not existent in an epoch, the position number of the component must be 0.
get_source_counts
Synopsis
calculates source counts, background counts and background-subtracted source counts from given spectrum, with errors
Usage
Double_Type get_source_counts(Integer_Type id, Integer_Type backid, Double_Type emin, Double_Type emax);
Description
- spectrum and background have to be already loaded (id, backid)
- arf and rmf have to be already loaded and assigned
See also: load_data,assign_rmf,assign_arf
hardnessratio
Synopsis
calculates various X-ray hardness ratio from given count rates, including uncertainty if requested
Usage
Double_Type hardnessratio(Double_Type soft_count, Double_Type hard_count);
Qualifiers
- color: calculate the hardness ratio according to color=log10(S/H)
- hardness: calculate the hardness ratio according to hardness=(H-S)/(H+S)
- ratio: calculate the hardness ratio according to ratio=S/H (the default)
- back_s: background counts in the soft band
- back_h: background counts in the hard band
- backscale: ratio between the extraction regions for the source and the background
- exposure: if given, interpret soft_count, hard_count and the backgrounds as rates. Multiply the rates with exposure to get the counts.
- backexposure: if given, the background exposure. Only taken into account if exposure is also set. If not given, it is assumed that the source and background exposures are identical.
- ratio_type: (DEPRECATED) Integer_Type, for hr=s/h choose 1, for hr=(h-s)/(h+s) choose 2, for hr=log(s/h) choose 3, default = 1
- err_s: Array containing the uncertainties of the soft light curve
- err_h: Array containing the uncertainties of the hard light curve
Description
This function calculates the so-called hardness ratio or color according to the three common prescriptions used in X-ray astronomy: S/H, (H-S)/(H+S), and log10(S/H). If given, background counts are subtracted before the calculation, taking into account different source and background extraction regions and different source and background exposure times if necessary (if counts are given but the background exposure was different, set the source exposure to 1 and the background exposure to the ratio of the source and background exposure times. If rates are given, give the exposure times.
Error bars are calculated using Gaussian error propagation. Contrary to earlier versions of this routine they are ALWAYS calculated!
hardnessratio_from_dataset
Synopsis
calculates hardnessratios of the current model
Usage
Double_Type hardnessratio_from_dataset(
Integer_Type data-id, Integer_Type[2] soft_ch, hard_ch
);
Qualifiers
- soft_en: Double_Type[2]
- hard_en: Double_Type[2]
- get_counts: set to '&get_model_counts' if needed, default: get_data_counts
- subtract_background: subtract background from data. Only possible if get_counts is not set. additional qualifiers are passed to 'hardnessratio'
Description
- at least one dataset has to be defined (load_data)
- if energy bands are given, channels have to be set to [0,0]
See also: hardnessratio
hardnessratio_simulate_grid
Synopsis
calculates hardnessratios of the current model
Usage
Struct_Type hardnessratio_simulate_grid(String_Type/Integer_Type par1name, par2name,
Integer_Type[2] soft_ch1, hard_ch1, soft_ch2, hard_ch2);
or Struct_Type hardnessratio_simulate_grid(
String_Type/Integer_Type par1name, Double_Type par1min, par1max, step1,
String_Type/Integer_Type par2name, Double_Type par2min, par2max, step2,
Integer_Type soft_ch1, hard_ch1, soft_ch2, hard_ch2
);
Qualifiers
- dataindex: Integer_Type, dataset index to use, default = 1
- grid1scale: 0 = linear (default), 1 = logarithmic
- grid2scale: 0 = linear (default), 1 = logarithmic
- par1grid: Double_Type[], override value grid of paramter 1
- par2grid: Double_Type[], override value grid of paramter 2
- sample1: Integer_Type, sampling-factor along each track, default 10
- sample2: Integer_Type, sampling-factor along each track, default 10
- exposure: Double_Type, default: 1e4
- arf: String_Type, file path of arf
- rmf: String_Type, file path of rmf
- rsp: String_Type, file path of rsp (arf and rmf combined)
- soft_en1: Double_Type[2], additional qualifier, soft band for all tracks in one direction
- hard_en1: Double_Type[2], additional qualifier, hard band for all tracks in one direction
- soft_en2: Double_Type[2], additional qualifier, soft band for all tracks in other direction
- hard_en2: Double_Type[2], additional qualifier, hard band for all tracks in other direction additional qualifiers are passed to 'hardnessratio_from_dataset'
Description
- at least one dataset has to be defined (load_data)
- a model must exist (fit_fun)
- soft_ch1, hard_ch1, soft_ch2, hard_ch2 are passed to 'hardnessratio_from_dataset'; if energy-bands (soft_en1...)are given, the channels have to be set to [0,0] each.
- if logarithmic gridscale is chosen, be careful 'parmin' is not <= zero (be aware of the default values of your model)
- the output consists of a structure of tracks, where each track contains the values of par1, par2 and the corresponding hardnessratios
- either give arf AND rmf; or ONLY give rsp
See also: hardnessratio, hardnessratio_from_dataset, xfigplot_hardnessratio_grid
hardnessratio_simulate_grid_load
Synopsis
load tracks saved with hardnessratio_simulate_grid_save
Usage
Struct_Type = hardnessratio_simulate_grid_save(String_Type filename);
See also: hardnessratio, hardnessratio_from_dataset, xfigplot_hardnessratio_grid, xfigplot_hardnessratio_grid_save
hardnessratio_simulate_grid_save
Synopsis
saves the output of hardnessratio_simulate_grid into a fits file
Usage
hardnessratio_simulate_grid_save(Struct_Type tracks,
String_Type filename);
Description
- tracks is the output strcutute of hardnessratio_simulate_grid
See also: hardnessratio, hardnessratio_from_dataset, xfigplot_hardnessratio_grid, xfigplot_hardnessratio_grid_load
history
Synopsis
shows the history of commands on the interactive ISIS-shell
Usage
history();
See also: save_input
init_jet_fit
Synopsis
inititialize the fit function for the speed of jet components
Usage
init_jet_fit (Ref_Type jet_component_structure);
Description
This function initializes the fit function jet_speed
, which
can fit a linear function of the form:
dist (t) = (t - t_0) * v
to each component, where t_0 is the ejection time of the component in
MJD and v is its speed in mas/year. For each epoch an offset (in mas)
can be fitted, allowing to correct for shifts of the image (along the
line in which the distance is measured).
For the initialization a structure with the name "jet_speed_struct"
is REQUIRED, which has to contain the fields time
(date of epoch in MJD),
distance
(in mas), derr
(uncertainty of the distance), and
component
(integer values identifying the components).
The function extends this structure, and sets the initialized fit function.
Example
variable jet_speed_struct = struct { mjd = [55500, 55550, 55600, 55650, 55500, 55600, 55650], deltax = [ 0.2, 0.4, 0.5, 0.7, 0.6, 1.2, 1.5], deltay = [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], derr = [ 0.01, 0.05, 0.01, 0.1, 0.02, 0.1, 0.1], component = [ 1, 1, 1, 1, 2, 2, 2], }; init_jet_fit(&jet_speed_struct; pa = 90); ()=fit_counts; list_par; plot_jet_speed ( jet_speed_struct );
See also: plot_jet_speed
integral_cat2ds9
Synopsis
a cat2ds9 like function for INTEGRAL catalogs
Usage
()=integral_cat2ds9(String_Type inputfile,
String_Type outputfile)
Qualifiers
- source_name: name of the column with source names (default: NAME)
- ra_name: name of the column with right ascension of the sources (default: RA_OBJ)
- dec_name:name of the column with declination of the sources (default: DEC_OBJ)
- color: color to be used for the ds9 region (default: black)
- symbol: form of the ds9 region (default:box)
- noname: only boxes, no text with source names
ionized_phabs_sigma
Synopsis
calculates the photoabsorption cross section for ionized material
Usage
Double_Type sigma = ionized_phabs_sigma(Double_Type E);
Description
E
is the energy in keV (may be a scalar or an array),
sigma
is the corresponding cross section in cm^2.
The ion of interest is specified by the following
Qualifiers
- Z: atomic number (1 <=
Z
<= 30) - Nel: number of electrons (
Z
>=Nel
> 0) - n: principal quantum number for the selection of a subshell (together with l)
- l: orbital angular momentum quantum number for the selection of a subshell
See also: Verner & Yakovlev (1995)
lorentz_trafo
Synopsis
lorentz transformation for a boost in any direction
Usage
lorentz_trafo( Double_Type ct,
Vector_Type x,
Vector_Type B
);
Description
This function performs a Lorentz transformation: ( ct, x ) -> ( ct', x') ct' = gamma ( ct - B*x ) x' = x + ( (gamma-1)/B^2 * B*x - gamma * ct ) * B where x is the spartial vector, gamma the Lorentz-factor and B=v/c the velocity vector.
See also: Vector_Type
make_spix
Synopsis
creates a spectral index map using two .fits files (provided by DIFMAP)
Usage
make_spix(String_Type <code>name_lo</code>, String_Type <code>name_hi</code>);
Qualifiers
- iterations [=1]: set number of translations to average over
- cc [=cc.fits]: name of the cross correlation map. Will generate one if the file does not exist.
- lothreshold [=NULL]: if image has a pixel below this value it's set to this value. also criterion for showing spectral index if req_both_bands is set to 1
- hithreshold [=NULL]: if image has a pixel below this value it's set to this value. also criterion for showing spectral index if req_both_bands is set to 1
- shift [={NULL,NULL}]: define the shift of the second image wrt the first in format {x,y} in mas
- crpix_shift [=0.]: Extra shift for converting mas in pixel. In the past, 0.5 was used, but default is now 0., as only a relative is needed
- dec_step [=NULL]: set the pixel size in declination (value>0). Default value is the largest pixel size of the two images
- dec_size [=NULL]: set the mapsize size in declination. Default value is the largest map size of the two images
- ra_step [=NULL]: set the pixel size in right ascension (value>0). Default value is the largest pixel size of the two images
- ra_size [=NULL]: set the mapsize size in right ascension Default value is the largest map size of the two images
- beam [=[NULL,NULL,NULL]]: set a beam for restoring the two images [semi-major axis, semi-minor axis, position angle] in mas (major and minor axis) and degree (position angle)
- req_both_bands [=0]: only make calculations if both bands are above the noise level
- fit_noise [=1]: set to 0 to use fits header for noise information
- lo_nsigma [=3]: number of standard deviations from mean to define noise limit. Only matters if fit_noise==1, overrides lothreshold
- hi_nsigma [=3]: number of standard deviations from mean to define noise limit. Only matters if fit_noise==1, overrides hithreshold
- n_beams [=3]: size of beam to exclude from core on correlation calculation
- excl_core [=1]: set to 0 if you don't want the program to automatically exclude the core from correlation calculations
- overwrite: overwrite the restore fits files (if already existing)
Description
This function creates a spectral index (F=v^a) map from two fits files provided by DIFMAP. A cross correlation fits image is generated and the images are shifted by the best (iterations) translations and the spectral index, etc. are calculated for each translations. The values are averaged with the weights given by the corresponding value in the cc image. If the provided images do not have the same size, resolution and beam, the functions difmap_restore and enclosing_ellipse are used to obtain images with these properties. It is also possible to specify a desired map size and pixel size for both axes individually, otherwise the default value from the comparison of both images will be used. Note, that changing the pixel size may require changing the maps size as well.
It returns a structure that holds the information for the following values
Use the structure as an input to
struct_name.spec_map -> matrix of calculated spectral index values
struct_name.stdev -> matrix of weighted stdevs of calculated values
struct_name.avg_lum -> matrix of weighted average of the summed BRIGHTNESS of the images
struct_name.weights -> matrix of sum of weights used for calculations in each pixel
struct_name.avg_shift -> average shifts in a 2-element array [x_shift,y_shift]
struct_name.shift_weight -> weights for calculating avg shift [x_weight,y_weight]
struct_name.avg_shift_pixel -> average shift in pixel [x_shift_pixel, y_shift_pixel]
struct_name.ra_px_center -> x index of center pixel
struct_name.ra_steps -> x mas per pixel
struct_name.dec_px_center -> y index of center pixel
struct_name.dec_steps -> y mas per pixel
struct_name.major -> semi major axis of beam in mas struct_name.minor -> semi minor axis of beam in mas struct_name.source -> source name struct_name.date -> dates of images struct_name.posang -> position angle of beam in degrees
See also: plot_spix, write_spix, read_spix, difmap_restore, enclosing_ellipse
mc_sig
Synopsis
calculates the significance of a spectral component doing a Monte Carlo (MC) simulation
Usage
Struct_Type mc_sig(String_Type withComponent.fits, String_Type withoutComponent.fits);
or Struct_Type mc_sig(String_Type torqueDir);
Qualifiers
mcruns - number of MC loops (default: 10) beforeData - script to be called before any data is loaded (default: NULL) afterData - script to be called after the data have been loaded (default: NULL) beforeFit - array of scripts to be called right before a fit of faked data to the model without [0] and with [1] the component (default: NULL) id - override the dataset ID(s) after fits_load_fit with the given one(s) ignbinning - do not apply the same binning after faking the data ignnotice - do not use the same energy ranges after faking the data torque - calculate significance using given number of torque jobs, i.e. total number of runs = mcruns*torque (default: 0, i.e. don't use torque) walltime - walltime of each torque job in minutes (default: 3 minutes times number of MC runs per job) dontSubmit - don't submit the torque job-file (implies qualifier 'dontwaint') (default: submit) dontWait - don't wait on the the torque jobs to complete (implies qualifier 'dontClean') (default: wait) torqueDir - temporary directory to save necessary torque files NOTE: this directory will be deleted afterwards completely! (default: $HOME/.isis_mc_sig/) dontClean - don't delete temporary torque directory afterwards chatty - a number >0 means more chatty (default: 0)
additional qualifiers are passed to fits_load_fit
Description
This function calculates the significance of a spectral component found in real data. During each Monte Carlo loop, spectra data without the component are simulated (for each detector) and this data are then fitted with a model containing the component (that needs to be tested for) and separately fitted with a model without it. The resulting simulated differences in chi square between these fits are returned and compared to the measured difference: the number of simulated chi squares below the measured one corresponds to the significance, that the spectral component is real (i.e. in 80 cases out of 100 runs the simulated chi square difference is below the measured chi square difference, the significance is 80%).
If two FITS-files created with fits_save_fit are provided, the first includes the model "with" the component to be tested and the second one "without" it.
If only one argument is given, this has to be a temporary directory with results of a previous torque run (e.g. if "dontWait" was set). The function tries to collect and return the results as usual.
The returned structure is defined as follows: readdchisqr - the measured difference in chi square fakedchisqr - an array of simulated differences in chi square nfalse - the number of detected false positives significance - the resulting significance as defined above
Example
% FITS-files created by fitting a cutoffpl and iron line to % RXTE-PCA, -HEXTE and Swift-XRT data (Rmf_OGIP_Compliance = 0 to % load XRT data, see help of fits_load_fit) sig = mc_sig("rxte_swift_cutoffpl_ironline.fits", "rxte_swift_cutoffpl.fits"; mcruns = 1, ROC = [2,2,0], beforeData = "defineMyModels.sl", afterData = "setDataHooks.sl", chatty = 2);
See also: fits_save_fit, fits_load_fit, fakeit
mpi_master_only
Synopsis
Use function with on one host only in an MPI job
Usage
mpi_master_only (Reference_Type function, arg1, arg2, ...)
Qualifiers
- verbose: show name of host which acts as MPI master
- veryverbose: show name of host which acts as MPI master, PID and parent PID (overwrites verbose)
Description
Sometimes, for example by saving files, it is required that a function is only called once in an MPI job. mpi_master_only does this for a function with an arbitrary number of arguments. The function has to be passed by reference, the arguments simply as arguments.
Example
mpi_master_only(&sprintf, "Pi is approximately %f and Eulers number is %f", PI, E); output: Pi is approximately 3.141593 and Eulers number is 2.718281828459045 But this only once
mpi_mc_sig
Synopsis
calculates the significance of a spectral component doing a Monte Carlo (MC) simulation
Usage
Struct_Type mpi_mc_sig(String_Type withComponent.fits, String_Type withoutComponent.fits [, String_Type results.fits]);
Qualifiers
mcruns - number of MC loops (default: 10) beforeData - script to be called before any data is loaded (default: NULL) afterData - script to be called after the data have been loaded (default: NULL) beforeFit - array of scripts to be called right before a fit of faked data to the model without [0] and with [1] the component (default: NULL) id - override the dataset ID(s) after fits_load_fit with the given one(s) ignbinning - do not apply the same binning after faking the data ignnotice - do not use the same energy ranges after faking the data chatty - a number >0 means more chatty (default: 0) seed - seed for random number generator
additional qualifiers are passed to fits_load_fit
Description
This function calculates the significance of a spectral component found in real data. During each Monte Carlo loop, spectra data without the component are simulated (for each detector) and these data are then fitted with a model containing the component (that needs to be tested for) and separately fitted with a model without it. The resulting simulated differences in chi-square between these fits are returned and compared to the measured difference: the number of simulated chi squares below the measured one corresponds to the significance, that the spectral component is real (i.e. in 80 cases out of 100 runs the simulated chi square difference is below the measured chi square difference, the significance is 80%).
Two FITS-files created with fits_save_fit must be provided, the first includes the model "with" the component to be tested and the second one "without" it. An optional third filename may be given to store the results to. Otherwise, a default file will be created.
This function is supposed to be used for parallel computation with MPI. The mcruns qualifier specifies the total number of MC runs which is split accordingly among the available nodes.
The function will create a FITS file containing the following extensions: realdeltachisquare - the measured difference in chi square fakedeltachisquare - an array of simulated differences in chi square falsepositives - the number of detected false positives significance - the resulting significance as defined above
Example
% FITS-files created by fitting a cutoffpl and iron line to % RXTE-PCA, -HEXTE and Swift-XRT data (Rmf_OGIP_Compliance = 0 to % load XRT data, see help of fits_load_fit)
mpi_mc_sig("rxte_swift_cutoffpl_ironline.fits", "rxte_swift_cutoffpl.fits", "significance.fits"; mcruns = 100, ROC = [2,2,0], beforeData = "defineMyModels.sl", afterData = "setDataHooks.sl", chatty = 2); % A typical call from the command line could look like this
mpiexec -n 4 isis
to use 4 cores for the MC run. However, it is highly recommended to use a job scheduling system (e.g., SLURM) for these kind of simulations.
See also: fits_save_fit, fits_load_fit, fakeit, mc_sig
orbit_calculator
Synopsis
Calculate orbits of test particles in a Galactic gravitational potential
Usage
orbit_calculator(Double_Types ah[], am[], as[], dd[], dm[], ds[], dist[], vrad[], pma_cos_d[], pmd[], t_end; qualifiers)
% or (for error propagation)
orbit_calculator(Double_Types ah, am, as, dd, dm, ds, dist, dist_sigma[], vrad, vrad_sigma[], pma_cos_d, pma_cos_d_sigma[], pmd, pmd_sigma[], t_end; qualifiers)
or
```c
orbit_calculator(Double_Types x[], y[], z[], vx[], vy[], vz[], t_end; qualifiers)
##### Description
Calculate orbits of test particles in a Galactic gravitational potential from given
initial conditions. The latter can be given either in celestial or Galactic cartesian
coordinates (see the help on the function 'cel2gal' for format and unit conventions).
Integration starts at time t=0 and ends at t=t_end [Myr], which is the last input
parameter. A negative t_end implies backward integration, which is e.g. useful to
trace back orbits in time. The potential and its equations of motion are outsourced
to a function specified by the qualifier 'model'. In this way, switching from one model
to another one is simply done by changing the before mentioned qualifier. Make always
use of cylindrical coordinates (r [kpc], phi [rad], z [kpc]) and their canonical momenta
(vr [kpc/Myr], Lz [kpc^2/Myr], vz [kpc/Myr]) as well as of cartesian coordinates (x [kpc],
y [kpc], z [kpc], vx [kpc/Myr], vy [kpc/Myr], vz [kpc/Myr]) to set up the equations of
motion (see the qualifier 'coords'). To numerically integrate the coupled differential
equations, an adaptive Runge-Kutta method of fourth/fifth order is used (see qualifier
'ODE_solver'). The stepsize is hereby controlled such that for each step an absolute
accuracy in coordinates (in units of kpc) and velocity components (in units of km/s) is
achieved that is smaller than given by the qualifier 'tolerance'. The function is optimized
for multi-orbit calculations, i.e., when the input parameters are either arrays or 1-sigma
uncertainties are given in addition. The latter are used to create Gaussian distributed
initial conditions to perform error propagation based on a number of Monte Carlo runs as
specified by the qualifier 'MC_runs'. To assign asymmetric uncertainties, use an array
of length two, i.e., [sigma_plus, sigma_minus], instead of a single number. The Gaussian
distribution is then split up into two Gaussian distributions, one with standard deviation
sigma_plus (for values larger than the respective input parameter) and one with sigma_minus
(else). To account for correlations between the distance (or alternatively the parallax,
see the qualifier 'parallax' in the function 'cel2gal'), proper motions in right ascension,
and proper motion in declination, make use of the qualifiers 'parallax_pma_corr',
'parallax_pmd_corr', and 'pma_pmd_corr'. Note that asymmetric uncertainties and
parameter correlations are mutually exclusive for individual parameters. All orbits are
computed simultaneously and on the same time grid. Stepsize control is hereby based on
the worst-offender principle. The function returns one structure with the two fields
"i" (initial) and "f" (final) (both again structures) containing the initial and final
Galactic cartesian coordinates (in kpc) and velocities (in km/s) as well as the z-component
of angular momentum Lz (in kpc^2/Myr) and total energy E_total = E_kin + E_pot (in kpc^2/Myr^2).
To see whether conservation of energy or conservation of angular momentum is implemented
in the equations of motion or not, have a look at the help of the outsourced potential
function. The final structure contains also the minimum and maximum Galactocentric distances
Rmin = min(R(t)) and Rmax = max(R(t)) with R(t) = sqrt(x(t)^2+y(t)^2+z(t)^2) and - if the
qualifier 'coords' is set to "cyl" - the number of revolutions about the Galactic z-axis
(negative values imply a motion in direction of Galactic rotation, i.e. prograde orbits,
while positive ones imply retrograde orbits). To save the initial and final structures to
fits files, set the qualifier 'stff' to the desired prefix of the filename. If not only
the initial and final situation is of interest, but rather the entire trajectories, use
the qualifier 'set'. In this case, the additional field "tr" (trajectory) (again a structure,
one field per orbit ("o1", "o2", ...)) is added to the returned structure containing all
orbits. For tracing back orbits to the Galactic disk, use the qualifier 'r_disk' to stop
individual calculations at the moment a trajectory crosses the xy-plane inside a circle
of radius 'r_disk'.
##### Qualifiers
* coords [<code>="cyl"</code>]: Use cylindrical ("cyl") or cartesian ("cart") coordinates for the internal
computations. Note that circular orbits like that of the Sun are computed faster in cylindrical
coordinates whereas cartesian coordinates are much more efficient for straight-line trajectories
or those that come very close to the z-axis where angular momentum terms slow down cylindrical
calculations.
* dt: Deactivate adaptive mode and use this fixed stepsize instead.
* model [<code>="AS"</code>]: Function ("AS", "MN_NFW", "MN_TF", or "plummer_MW") evaluating the equations
of motion, circular velocity, and energy of the model potential.
* MC_runs [<code>=10000</code>]: Number of Monte Carlo realizations in the case that 1-sigma errors are given.
* ODE_solver [<code>="RKCK"</code>]: Choose among three different Runge-Kutta (RK) integration methods:
"RKF": RK-Fehlberg, "RKCK": RK-Cash-Karp, "RKDP": RK-Dormand-Prince.
* parallax_pma_corr [<code>=0</code>]: Correlation between parallax and proper motion in right ascension.
* parallax_pmd_corr [<code>=0</code>]: Correlation between parallax and proper motion in declination.
* pma_pmd_corr [<code>=0</code>]: Correlation between proper motion in right ascension and in declination.
* r_disk [<code>=0</code>]: If not zero, integration is stopped at the moment of crossing the xy-plane
inside a circle of radius r_disk [kpc].
* seed [<code>=_time()</code>]: Seed the random number generator via the function 'seed_random'.
* set [<code>=0</code>]: If present, trajectories will be saved ("Save entire trajectories"). As
this can be very memory-consuming for a large number of orbits, one can additionally
use the value of this qualifier to specify a lower limit on the time difference of
two consecutive moments of time that will be saved.
* stff: "Save to fits files": Prefix of fits files to which initial and final structures are written.
* SunGCDist: By default, the Sun-Galactic center distance is taken from the current model.
Use this qualifier to explicitly set a distance in kpc.
* tolerance [<code>=1e-10</code>]: Absolut error control tolerance; lower limit: 1e-15.
* verbose: Show intermediate times t.
* Any qualifiers from 'cel2gal' and the model potential function except 'vlsr' and 'eomecd'.
Important note: for consistency, the local standard of rest velocity vlsr is calculated from the
circular velocity of the current model potential evaluated at (r=SunGCDist, z=0).
##### Example
s = orbit_calculator(12,22,29.6,40,49,36,3.078,262,-13.52,16.34,-1000; r_disk=50);
s = orbit_calculator(12,22,29.6,40,49,36,3.078,[0.6,0.3],262,5,-13.52,1.31,16.34,1.37,-1000; r_disk=50, stff="HIP60350", MC_runs=100);
s = orbit_calculator(-8.4,0,0,0,242,0,250; set, model="MN_TF");
plot(s.tr.o0.x, s.tr.o0.y);
s = orbit_calculator([8:10:#500], [0:1:#500], [4:5:#500], [0:0:#500], [210:230:#500], [-10:-50:#500], 1000);
% Example using parallax, proper motions, and correlation parameters from Gaia:
s = orbit_calculator(12,22,29.6,40,49,36,0.325,0.3,262,5,-13.52,1.31,16.34,1.37,-1; parallax, parallax_pma_corr=0.1, parallax_pmd_corr=0.2, pma_pmd_corr=0.75);
% Example using the spectroscopic distance and only proper motions from Gaia:
s = orbit_calculator(12,22,29.6,40,49,36,3.077,[0.6,0.3],262,5,-13.52,1.31,16.34,1.37,-1; pma_pmd_corr=0.75);
__See also__: cel2gal, AS, MN_NFW, MN_TF, plummer_MW, xfig_3d_orbit_on_cube
#### overplot_clean_model
##### Synopsis
creates an overlay of an clean VLBI image and the modelfit components
##### Usage
```c
overplot_clean_model(String_Type <code>cleanfile</code>, String_Type <code>modfile</code>, String_Type <code>outputfile</code>);
Qualifiers
- ra_mas [={20,-20}]: specifies {left,right} limits of image in mas
- dec_mas [={-20,20}]: specifies {top,bottom} limits of image in mas
- plot_size [=15]: size of plot
- n_sigma [=3.0]: lowest contour of clean image
- sourcename [=default]: set the name of the source, by default the name is read from the .fits file,
set to NULL for not plotting a source name
- obs_date [=default]: set the observation date, by default the observation date is read from the .fits file,
set to NULL for not plotting a observation date, set to "mjd" when MJD format required
- cont_color [="gray"]: color of contours of clean image
- cont_scl [="2"]: set factor to change the separation between contour levels
- cont_lvl [=[c1,c2,..]: set contour levels [Jy] manually, overwrites other contour parameters
- model_color [="black"]: color of Gaussian ellipses
- center_symbol [="+"]: symbol stating the ellipse center, set to NULL for no symbol
- comp_label [=0]: set to 1 when labelling of components depending on distance from [0,0] is requested
- ind_inverse [=0]: if the labels should be in inverse order set to 1
- label_alt [=NULL]: give list of alternative labels
- label_color [="red"]: color of component labels
- ex_counterjet [=0]: set to 1 if counterjet components (and core) should be excluded and define the corresponding quadrant of the plot via the CJ-coordinates
- counterjet_ra [=NULL]: define counterjet coordinates in mas
- counterjet_dec [=NULL]: define counterjet coordinates in mas
Description
This function creates an overlay of a VLBI-clean image and the corresponding
model of Gaussian components which are over-plotted as ellipses.
The required input format are fits-file for both the clean and the modelfit images. The format of the output file depends on the suffix of the given
filename
. Possible formats of the output file are PDF, EPS,
PNG, GIF, etc. If labelling = 1 is set, the components are labeled with J0,J1,J2,... depending on their distance to [0,0]. The counterjet components (and core) can be excluded by defining the corresponding quadrant of the the plot via counterjet_ra and _dec.
plot_component
Synopsis
to plot modelfit parameters of different components
Usage
plot_component([array of <code>struct comp</code>]);
Qualifiers
- outputfile : give outputfile name and directory by hand
- mjd_min [=53736.0]: to specify lower limit of time axis in mjd (default: 1/1/2006)
- mjd_max [=55927.0]: to specify lower limit of time axis in mjd (default: 1/1/2012)
- date_min [=2006.0]: to specify lower limit of time axis (default: 1/1/2006)
- date_max [=2012.0]: to specify lower limit of time axis (default: 1/1/2012)
- distance_min [=-2]: lower limit of distance axis in mas
- distance_max [=20]: upper limit of distance axis in mas
- flux_min [=1e-5]: lower limit of flux axis in Jy
- flux_max [=5]: upper limit of flux axis in Jy
- TB_min [=10^5]: lower limit of brightness temperature axis in K
- TB_max [=10^15]: upper limit of brightness temperature axis in K
- size_min [=1e-3]: lower limit of size axis in mas^2
- size_max [=100]: upper limit of size axis in mas^2
- pos_min [=-180]: lower limit of pos angle in degrees
- pos_max [=180]: upper limit of pos angle in degrees
- labels [=default]: provide array of custom labels for the components
- sym [=default]: provide array of custom symbols for the components
- symcolor [=default]: provide array of custom colors for the components
- symsize [=1]: provide array of custom symbol sizes or one value to apply size to all
- legend: create a legend
- lpos_x: set the x position of the legend in world0 system (0 = left, 1=right)
- lpos_y: set the y position of the legend in world0 system (0 = bottom, 1=top)
- linestyle [=default]: set to 0 if the flux values should not be conntected (or to another value for another line style
- nocomp: provide a comp structure (or array of structures) to plot not identified components in black
- plots [=[9,10,3,4,11,12]]:provide an array of numbers to generated specified plots: 1: distance vs mjd 2: flux vs mjd 3: T_B vs distance 4: flux vs distance 5: size vs distance 6: pos angle vs distance 7: size vs mjd 8: pos angle vs mjd 9: distance vs time 10: flux vs time 11: size vs time 12: pos_angle vs time
Description
This functions creates overview plots for modelfit compoents. The required input format is the comp structure which can be obtained with get_component or an array of the comp structure.
plot_jet_speed
Synopsis
visualizes the jet_speed_struct used by init_jet_fit
Usage
plot_jet_speed (Struct_Type jet_speed_struct);
Qualifiers
- xrange: set the xrange (provide a list or array)
- yrange: set the yrange (provide a list or array)
See also: init_jet_fit
plot_spix
Synopsis
plots a spectral index map using the struct provided by make_spix
Usage
plot_spix(Struct_Type values, String_Type outputname);
Qualifiers
- inv_color [=0]: inverts color scale (i.e. roles of red and blue switch)
- min_spix [=min(values.spec_map)]: minimum spectral index to be plotted
- max_spix [=max(values.spec_map)]: maximum spectral index to be plotted
- lum_scale [=1]: use the sum of the intensity of both images to scale the brightness of the spectral index map, set to 0 in order to show the significant regions at uniform brightness
- nsigma [=0]: spectral index values are only shown for pixels where the emission (in both bands) is above the specified significance (e.g., nsigma=3 means that the emission has to be 3*sigma above the noise level)
- nsigma_lo [=0]: specify emission limit for low frequency band separately
- nsigma_hi [=0]: specify emission limit for high frequency band separately
- low_lum_limit [=0]: number of stdevs above the median that luminosity data begins to be displayed
- hi_lum_limit [=0]: number of (e^n) stdevs above the median at which color is fully saturated (0 -> no limit)
- ra_frac [={0.0,1.0}]: specifies fraction of RA to display with {left,right} limits, will keep scaling
- dec_frac [={0.0,1.0}]: specifies fraction of DEC to display with {left,right} limits, will keep scaling
- ra_mas [={NULL,NULL}]: specifies {left,right} limits of image in mas
- dec_mas [={NULL,NULL}]: specifies {top,bottom} limits of image in mas
- size [=14]: size of plot (maximum dimension)
- print_mode [=0]: set =1 for white background
- blue_bkgr [=0]: for a dark blue background set to 1
- quadratic: create a quadratic image, by default the scaling is such that scaling in RA and DEC is equal
- no_labels : plot no labels
- pmodecolor [=1]: color saturation decreases wth this root of luminosity (i.e. 2=sqrt) in print mode only
- beam_color [="gray"]: beam color
- src_name [=values.source]: manually set source name in plot
- obs_date [=values.date]: manually set date in plot
- xyunit [="mas"]: unit of the x/y labels, can be changed between "mas", "arcsec", "arcmin" and "deg"; assumption: FITS header sets unit "mas"
- xfig_tmp_dir [=xfig_get_tmp_dir]: set the path for the tmp directory used by xfig
- xfig_autoeps_dir [=xfig_get_autoeps_dir]: set the path for the autoeps directory
- color_scheme [="hot"]: setting another colorscheme (lum_scale=0 is automatically set)
- h_scalebar: if given, a horizontal scalebar is plotted above the spixmap
- no_scalebar: if given, no scalebar is plotted
Description
This function creates a spectral image plot with xfig from the structure generated by make_spix. The log of the average brightness is displayed as color saturation and the spectral index as hue so that both spectral index and brightness are apparent in the image. It will ouput a .pdf file of the plot with a color scale bar on the side. For the luminosity scaling, the data below the median value is cut out and the data above the median is fit to a gaussian. The hi_lum_limit and low_lum_limit are based on the stdev from this fit.
See also: make_spix, read_spix, write_spix
plot_vlbi_map
Synopsis
creates an image of the VLBI map with transparent background using a .fits file (provided by DIFMAP)
Usage
plot_vlbi_map(String_Type <code>fitsfile</code>, [String_Type <code>fitsfile (polflux)</code>], String_Type <code>filename</code>)
or
```c
plot_vlbi_map(Struct_Type <code>img_struct</code>, String_Type <code>filename</code>)
##### Qualifiers
* color_scheme [="ds9b"]: select a color scheme for the image
* dec_frac [={0.0,1.0}]: declination range of the image by
selecting a fraction of the file's range
* ra_frac [={0.0,1.0}]: right ascension range of the image by
selecting a fraction of the file's range
* dec_mas [={min,max}]: set declination range of the image directly in mas.
overwrites <code>dec_frac</code>
* ra_mas [={min,max}]: set right ascension range of the image directly in mas
(analog to <code>dec_mas</code>)
* fit_noise [=1]: fit the noise of a map
* n_sigma [=3.0]: start color scaling at <code>n_sigma\*sigma</code>
* cont_scl [=2.0]: set factor to change the separation between
contour levels, the levels are placed at:
<code>cont_scl^[0,1,...]\*n_sigma\*sigma</code>
* cont_lvl [=[c1,c2,..]: set contour levels [Jy] manually, overwrites
other contour parameters
* cont_width [=2]: set width of contour lines
* plot_cont [=1]: set to 0 in order to plot without contour lines
* cont_depth [=2]: set contour-line depth
* plot_vec [=0]: set to 1 in order to plot with vectors
(requires second fits file with polarized flux density)
* vec_density [=5]: density of EVPA vectors; possible values: [1..10]
* vec_width [=1]: width of EVPA vectors
* vec_color [="black"]: color of EVPA vectors
* plot_clr_img [=1]: set to 0 in order to plot without colored image
* plot_clr_key [=1]: set to 0 in order to plot without key for color scale
* plot_scale_arrows [=1]: set to 1 in order to plot with arrows denoting the size scalesin the map
* clrcut [=0]: set to 1 in order to start the color scale at the
significant emission (i.e., no color scaling below)
* clrmin [=min(image)]: set minimum value [Jy] for color scale (see png_gray_to_rgb)
* clrmax [=max(image)]: set maximum value [Jy] for color scale (see png_gray_to_rgb)
* clrmu [from fit_gauss_to_img_noise]: set noise level of the image [Jy] for color scale
(mean value of regions without significant emission)
* clrsig [from fit_gauss_to_img_noise]: set width of noise (1-sigma) [Jy] for color scale,
for identical color scale (e.g., in order to compare
different images) the following qualifiers have to be set:
clrmin, clrmax, clrmu, clrsig
* colmap_depth: set depth of png map
* bkg_white: set white map below clrmin, if clrcut is set, clrmin equals
n_sigma above the mean background.
* quadratic: create a quadratic image, by default the scaling is such that
scaling in RA and DEC is equal
* plotsize [=14]: size of the plot (scaling of image relative to font size)
* src_name [=default]: set the name of the source, by default the name
is read from the .fits file
set to NULL for not plotting a source name
* obs_date [=default]: set the observation date, by default the observation
date is read from the .fits file
set to NULL for not plotting a observation date
* axis_color [="gray"]: set axis color
* source_color [="gray"]: set color of source name
* date_color [="gray"]: set color of observation date
* date_depth [=1]: set depth of observation date
* source_depth [=1]: set depth of source name
* cont_color [="gray"]: set color of contour lines
* neg_cont_color [="gray"]: set color of contour lines below noise level
* date_xy [=[0.95,0.95]]: set world0 coordinates of date string
* plot_beam [=1]: set to 1 in order to plot the beam
* plot_components: set to plot model components on top of image
* model_circs: set to plot model component circles on top of image
* pos_comp_color [="seagreen"]: set for color of model components with positive flux
* neg_comp_color [="seagreen"]: set for color of model components with negative flux
* beam_color [="gray"]: set color of the beam
* no_labels: plot no labels
* colmap_label: set label of colormap
* neglog [=0]: use log scale for negative flux values
* frac [=0.1]: fraction to which the linear region below n_sigma\*sigma
is scaled in the color scale (if clrcut != 0)
* funit [="mJy"]: unit of the color scale bar, can be changed from "mJy" to "Jy"
* xyunit [="mas"]: unit of the x/y labels, can be changed between "mas", "arcsec", "arcmin" and "deg";
assumption: FITS header sets unit "mas"
* return_xfig: set qualifier to return the xfig object(s) instead than
rendering directly
##### Description
This function creates an image of a jet using the .fits file provided
by DIFMAP.
Alternatively a structure with the required fields (as obtained
with the function <code>read_difmap_img</code>) can be given directly to
<code>plot_vlbi_map</code> instead of the name of a fits file.
The color scheme can be selected and contour lines are calculated.
The format of the output file depends on the suffix of the given
<code>filename</code>. Possible formats of the output file are PDF, EPS,
PNG, GIF, etc.
If plot_vlbi_map is called with three arguments, the first two of them beiing maps,
the code will determine the noise level of the second map and cut the first
map at the given sigma level. A typical application would be to plot the
distribution of electric vectors or the degree of polarization only where
the polarized flux is significant.
__See also__: read_difmap_fits,fit_gauss_to_img_noise
#### plot_vlbi_pol_map
##### Synopsis
creates a polarimetric image of a VLBI map using according .fits files (provided by DIFMAP)
##### Usage
```c
(Struct_Type xfig-map, Struct_Type xfig-colormap) = plot_vlbi_pol_map(String_Type <code>fitsfile flux-map-1</code>,
String_Type <code>filename flux-map-2</code>, String_Type <code>filename plot-product</code>)
or
```c
(Struct_Type xfig-map, Struct_Type xfig-colormap, Struct_Type xfig-vectormap) =
plot_vlbi_pol_map(String_Type <code>fitsfile flux-map-1</code>, String_Type <code>fitsfile EVPA-map</code>,
String_Type <code>filename flux-map-2</code>, String_Type <code>filename plot-product</code>)
##### Qualifiers
* n_sigma_pol [=5]: start color scaling of polarized flux at <code>n_sigma\*sigma</code>
* n_sigma_vec [=5]: start plotting of EVPA vectors at <code>n_sigma\*sigma</code> in reference to the polarized flux
* colmap_pol: if qualifier exists, plot most significant polarized flux (depending on n_sigma_pol),
recommended: color map "iceandfire"
* colmap_back: if qualifier exists, plot polarized flux or total intensity flux, recommended: color map "ds9b"
* pol_cont: set to 1 in order to plot contours of flux-map-1
* pol_cont_col [="white"]: set color of contours of flux-map-1
* pol_cont_depth [=0]: set depth of contours of flux-map-1
* pol_translate [=0.4]: if plot_pol_colmap=0, the contours of polarized flux will be translated
by <code>pol_translate\*range_in_declination</code>
* render_object [=1]: set to 1 in order to render the plot object
##### Description
This function uses the existing maps flux-map-1 (polarized flux),
(the EVPA_map for the electric vector position angle) and
flux_map_2 (polarized flux or total intensity flux) to plot the polarized flux
as color-coded map (and the EVPA as vectors) - cutted at the most significant contours
depending on n_sigma_pol - on top of the total intensity contours or as contours
vertically translated to them.
One has the flexibility to plot the polarized flux distribution instead of the total
intensity flux color-coded depending on the map given for flux-map-2.
All qualifiers of <code>plot_vlbi_map</code> are forwarded to this function.
__See also__: plot_vlbi_map
#### plummer_MW
##### Synopsis
Alternative model potential for the function 'orbit_calculator'
##### Usage
```c
plummer_MW(Double_Types t, m[6,n]; qualifiers)
Description
This function provides an alternative model for the gravitational potential of the Milky Way which can be used by the function 'orbit_calculator'. The gravitational forces of a standard Milky Way potential (see qualifier 'model') are combined with those arising from the interaction with a number of moving Plummer spheres (see function 'plummer_interaction_kernel' and qualifier 'plummer_spheres') to determine the acceleration of n independent test particles at time 't'.
Notes
Because of the unit convention used for the potentials of the Milky Way, the field 'psa' of the qualifier 'plummer_spheres', which is assumed to be the mass of the respective Plummer sphere in solar masses, has to be converted to Galactic mass units and then multiplied with a constant accounting for the remaining unit conversions (see example below). The units of 'psb' have to be kpc^2.
Qualifiers
- plummer_spheres: Structure whose fields are again structures with fields 't' [Myr], 'x' [kpc], 'y' [kpc], 'z' [kpc], 'psa' [Msun/Mgal*constant] (see notes), and 'psb' [kpc^2] describing the orbits and shapes of the Plummer spheres. Always make sure that the tabulated times 't[*]' are in monotonic increasing order if 't[-1]' is positive or in monotonic decreasing order if 't[-1]' is negative.
- MW_potential [
="AS"
]: Function ("AS", "MN_NFW", or "MN_TF"), which evaluates the equations of motion that result from a model for the gravitational potential of the Milky Way. - All qualifiers from the Milky Way model potential (see qualifier 'MW_potential').
- All qualifiers from the function 'plummer_interaction_kernel'.
Example
% Test particle affected by the Milky Way and satellite galaxies: t_end = -100; % integration time in Myr model = "AS"; % Milky Way mass model s = properties_satellite_galaxies(); i = struct{ x, y, z, vx, vy, vz, psa, psb }; SunGCDist = (@(__get_reference(model)))(; eomecd="sgcd"); % Sun-GC distance of chosen mass model temp = [SunGCDist,0,0,0,0,0]; reshape(temp, [6,1]); vlsr = (@(__get_reference(model)))(0, temp; eomecd="circ")[0]; % Local standard of rest velocity of chosen mass model (i.x, i.y, i.z, i.vx, i.vy, i.vz) = cel2gal(s.ah, s.am, s.as, s.dd, s.dm, s.ds, s.dist, s.vrad, s.pma_cos_d, s.pmd; SunGCDist=SunGCDist, vlsr=vlsr); kpcmyr_to_kms = 977.7736364875057; % = conversion factor from kpc/myr to km/s = 3.0856775975*10^16 / (10^6 * 3.15582 * 10^7); i.vx /= kpcmyr_to_kms; i.vy /= kpcmyr_to_kms; i.vz /= kpcmyr_to_kms; i.psa = s.Pl_mass/2.325131802556774e+07; % conversion from solar masses to Galactic mass units Mgal to have G=1 % Mgal = 2.325131802556774e+07 = 1e8*3.0856775975*1e19/6.6742/1e-11/1.9884/1e30, see Irrgang et al., 2013, A&A, 549, A137 const = 100./kpcmyr_to_kms^2; % factor 100 because potential is given in 100 km^2/s^2, see Irrgang et al. 2013 i.psa *= const; i.psb = (s.Pl_radius)^2; ps = N_body_simulation(i, t_end; kernel="N_body_simulation_MW_kernel", psa=i.psa, psb=i.psb, model=model); % ps contains time-dependent coordinates of Plummer spheres % add information about shape of the Plummer spheres: j = 0; foreach field (get_struct_field_names(ps)) { temp = struct{ psa=i.psa[j], psb=i.psb[j] }; set_struct_field(ps, field, struct_combine(get_struct_field(ps, field), temp)); j++; }; % compute trajectories of test particles: s = orbit_calculator(4,38,12.8,-54,33,12,61,723,0.86,0.57,t_end; set, model="plummer_MW", MW_potential=model, plummer_spheres=ps); plot(s.tr.o0.x,s.tr.o0.y); % without satellite galaxies: s = orbit_calculator(4,38,12.8,-54,33,12,61,723,0.86,0.57,t_end; set, model=model); oplot(s.tr.o0.x,s.tr.o0.y);
See also: orbit_calculator, plummer_interaction_kernel, AS, MN_NFW, MN_TF
plummer_interaction_kernel
Synopsis
Evaluate equations of motion or potential energy from a number of Plummer spheres
Usage
Double_Type eom[3,n] = plummer_interaction_kernel(Double_Types t, m[6,n], Struct_Type ps; qualifiers)
or
```c
Double_Type energy[n] = plummer_interaction_kernel(Double_Types t, m[6,n], Struct_Type ps; qualifiers)
##### Description
This function computes the equations of motion or potential energy of n test particles
at time 't' caused by the interaction with a number of moving Plummer spheres. Depending
on whether cylindrical coordinates (r,phi,z) and their canonical momenta (vr,Lz,vz) or
cartesian coordinates (x,y,z,vx,vy,vz; see qualifier 'coords') are used, the second
input parameter 'm' is a [6,n]-matrix with
(qualifier("coords")=="cyl") or (qualifier("coords")=="cart")
m[0,\*] = r; m[0,\*] = x;
m[1,\*] = phi; m[1,\*] = y;
m[2,\*] = z; m[2,\*] = z;
m[3,\*] = vr; m[3,\*] = vx;
m[4,\*] = Lz; m[4,\*] = vy;
m[5,\*] = vz; m[5,\*] = vz;
If the qualifier 'eomecd' is set to "eom", the function returns a [3,n]-matrix 'delta' with
delta[0,\*] = -d/dr Phi(r,phi,z); delta[0,\*] = -d/dx Phi(x,y,z);
delta[1,\*] = -d/dphi Phi(r,phi,z); delta[1,\*] = -d/dy Phi(x,y,z);
delta[2,\*] = -d/dz Phi(r,phi,z); delta[2,\*] = -d/dz Phi(x,y,z);
whereby Phi is the sum over all Plummer potentials.
If the qualifier 'eomecd' is set to "energy", the function returns an array of length n
storing the potential energy for each orbit:
E(r,phi,z) = Double_Type[n] = Phi(r,phi,z)
or
E(x,y,z) = Double_Type[n] = Phi(x,y,z)
For each Plummer sphere, the third input parameter 'ps', which is a structure, contains
a field which is again a structure with fields 't', 'x', 'y', and 'z' (all are arrays of
the same length) that list the time-dependent positions of the sphere. Additionally, the
shape of the respective sphere (see notes below) is given by the two fields 'psa' and
'psb' (both scalars).
##### Notes
The potential of a Plummer sphere at distance r is
Phi(r) = -psa\*(psb+r^2)^(-1/2)
The resulting acceleration in radial direction is
d^2/dt^2 r = -d/dr Phi = -psa\*(psb+r^2)^(-3/2)\*r
Cubic spline interpolation is used to determine the positions of the Plummer spheres
at time 't' if the GSL-module is available. Otherwise, linear interpolation is applied.
Extrapolation is not allowed. Always make sure that the tabulated times '.t[\*]' are in
monotonic increasing order if '.t[-1]' is positive or in monotonic decreasing order if
'.t[-1]' is negative.
##### Qualifiers
* coords [<code>="cyl"</code>]: Use cylindrical ("cyl") or cartesian ("cart") coordinates.
* eomecd [<code>="eom"</code>]: Return equations of motion ("eom") or potential energy ("energy").
##### Example
ps = struct{ o0, o1 };
ps.o0 = struct{ t=[0,1,2], x=[-1,0,1], y=[0,1,2], z=[0,1,2], psa=1, psb=4 };
ps.o1 = struct{ t=[0,1,2], x=[1,0,-1], y=[0,1,2], z=[0,1,2], psa=1, psb=4 };
m = Double_Type[6,1];
m[0,0] = 0; m[1,0] = 2; m[2,0] = 0; m[3,0] = 0; m[4,0] = 0; m[5,0] = 0;
r = plummer_interaction_kernel(0,m,ps; coords="cart");
r = plummer_interaction_kernel(0,m,ps; eomecd="energy");
__See also__: N_body_simulation_std_kernel, plummer_MW
#### powerlaw_noise
##### Synopsis
simulates a light curve with a power law distributed power spectrum
##### Usage
```c
Double_Type rate[] = powerlaw_noise(Integer_Type n);
Qualifiers
- beta [= 1.5]: power law index of the power spectrum
- mean [= 0]: mean rate of the simulated light curve
- sigma [= 1]: standard deviation of the simulated light curve
Description
n
is the number of bins of the simulated lightcurve.
It should be a power of two best performance of the fast Fourier transform.
See also: Timmer & Koenig (1995): "On generating power law noise", A&A 300, 707-710
principal_components
Synopsis
performs a principal components analysis
Usage
Struct_Type PCA = principal_components(Struct_Type s);
Description
The normalized components (which are stored in PCA.components.c
#i)
are calculated from the fields of the structure s
such that they have a mean of 0 and a variance of 1.
From them, the covariance matrix PCA.cov_matrix
is calculated,
which is diagonalized (see PCA.eigenvalues
and PCA.eigenvectors
).
The principal components are stored in PCA.components.pc
#i
in ascending order of their contribution to the total variance.
Qualifiers
- table[="tab"]: name of the structure
s
See also: cov_matrix
print_statistics
Synopsis
shows statistical information on an array of numbers
Usage
print_statistics(Double_Type a[]);
or
print_statistics(Struct_Type s);
Description
If the argument of print_statistics
is an array a
,
its minimum, average, standard deviation and maximum
are shown. If a
contains irregular numbers (nan
or +/-inf
),
the same quantities are also shown for the regular numbers only.
If print_statistics
is called with a structure argument,
the above mentioned task is performed on all array fields.
See also: moment
psd_lc
Synopsis
simulate a random light curve that follows a given PSD using the algorithm of Timmer and Koenig
Usage
Double_Type rate[] = psd_lc(Integer_Type n, Double_Type dt, Ref_Type PSD);
or
(rate1, rate2) = psd_lc(Integer_Type n, Double_Type dt, Ref_Type PSD; time_lag_spectrum=...);
Qualifiers
- mean [= 100]: mean count rate of the simulated lightcurve
- sigma [= 20]: standard deviation of the simulated lightcurve
- poisson: if poisson noise should be applied on the lightcrurve
- nr_PCUs [= 1]: number of PCUs
- time_lag_spectrum: spectrum of time lags (two lightcurves will be returned)
- mean_2 [= 100]: mean count rate of the second lightcurve
- sigma_2 [= 20]: standard deviation of the second lightcurve
Description
n
is the number of bins of the simulated lightcurve.
It should be a power of two for best performance of the FFT.
dt
is the time resolution.
PSD
is a reference to a function which takes one argument
-- the frequency -- and calculates the corresponding PSD value.
The number of PCUs is needed for the calculation of Poisson noise
if a mean RXTE count rate is given in counts/PCU.
If time_lag_spectrum
is reference to a function of one argument
-- the frequency -- that calculates the time lag spectrum,
an additional second lightcurve is returned that has
the corresponding time lag with respect to the first one.
See also: Timmer & Koenig (1995): "On generating power law noise", A&A 300, 707-710
radio_mod2img
Synopsis
creates an image using a model and a beam
Usage
Struct_Type img = radio_mod2img( Struct_Type <code>mdl</code>, Double_Type <code>beam</code>);
Qualifiers
- src_name [=NULL]: set the name of the source
- date_mjd [=_NaN]: set the observation date (MJD)
- obs_date [=NULL]: set the observation date (string). if only date_mjd or only obs_date are set, the other one is calculated.
- nrpix_ra [=512]: number of pixels for right ascension
- nrpix_dec: number of pixels in declination (by default calculated from RA-DEC-range)
- ra [=[ra_min,ra_max]]: RA range in mas
- dec [=[dec_min,dec_max]]: DEC range in mas, by default the RA-DEC-range is calculated from the distribution of model components and the beam size
- delt [=[ra_delt,dec_delt]]: resolution in mas/pixel, overwrites nrpix qualifiers
- sigma [=1e-3*max(img)]: set the value for 1 sigma (used by plot_vlbi_map)
Description
This function uses a model to generate an image and convolves it with the beam.
Here the "beam" is simply a Gaussian profile, which is given as an array:
beam = [smajor_axis,sminor_axis,position_angle];
The model has to have the fields:
flux
flux of each model component [Jy]
ra
relative RA of component (change to delta_x?)
dec
relative RA of component (change to delta_y?)
smajor
component size (smajor axis) [mas] (0 for point source)
sminor
component size (sminor axis) [mas]
pang
position angle of component's smajor axis
Currently the last two fields (sminor
, and pang
) are
ignored and only circular components (smajor
) are used.
Example
variable mdl = struct { flux = [0.6 , 0.9, 0.2, 1.2, 3], ra = [2 , 1.4, 0.7, 0.25, 0], dec = [1 , 0.6, 0.35, 0.1, 0], smajor = [0.3 , 0.2, 0, 0, 0], sminor = [0.3 , 0.2, 0, 0, 0], pang = [0 , 0, 0, 0, 0] }; variable beam = [0.2, 0.07, 0.4]; variable img = radio_mod2img (mdl, beam); plot_vlbi_map (img, "test.pdf");
See also: plot_vlbi_map, read_difmap_fits
read_difmap_fits
Synopsis
read a fits image provided by DIFMAP
Usage
Struct_Type img_struct = read_difmap_fits(String_Type <code>fitsfile_name</code>)
Qualifiers
- fit_noise: fit properties of the image noise
Description
This function reads a .fits file provided by DIFMAP and returns
a structure containing the image and keywords.
The returned structure can be used as input for the function
plot_vlbi_map
.
See also: plot_vlbi_map
read_ep
Synopsis
read in the epoch files
Usage
read_ep([array of epoch fits files]);
Qualifiers
- quadrant [=1]: Quadrant which is defined as a positive distance ( RA, DEC > 0 == 1 || RA < 0, DEC > 0 == 2 || RA < 0, DEC < 0 == 3 || RA > 0, DEC < 0 == 4)
Description
This functions returns a structure of an Epoch . The require input is an array of epoch fits files.
read_spix
Synopsis
reads the FITS file created with write_spix
Usage
Stuct_Type <code>spix</code> = read_spix( String_Type <code>filename</code>);
See also: make_spix, write_spix, plot_spix
require_atoms
Usage
require_atoms();
Description
This function loads the atomic data from the Astrophysical Plasma Emission Database via
atoms(aped);
unless _isis->Dbase
is already initialized.
See also: atoms
ridge_line
Synopsis
calculates the ridge line
Usage
Struct_Type ridge_line = ridge_line (Array_Type <code>img</code>)
Qualifiers
- ref_pix [=[y,x]]: starting pixel of the ridge line (index convention: img[y,x]) brightest pixel used by default
- dr [=0.8]: step size in pixels
- steps: number of steps (default is maximal length of dimensions)
Description
This function calculates the ridge line in a (jet) image.
Starting from the brightest point (assumed to be the core),
or a point specified by the ref_pix
qualifier, the ridge
line is calculated. The ridge line is given in different ways:
The brightest points at each distance (in steps of dr
pixels)
from the reference pixel, are given by the fields peak_x
and
peak_y
of the returned structure. The field peak_flux
contains the corresponding pixel value.
As an alternative to the peak ridge line, a flux-centered ridge
line is provided. The fields flux_cent_x
and flux_cent_y
specify the points, at which the integrated flux (along the same
distance from the refence pixel) to the left is equal to that on
the right side.
Notes
- elliptical beams will create artifacts (wiggles in the ridge line), in order to avoid this effect, obtain the ridge line from an image convolved with a circular beam
- valid points can be obtained, e.g., by filtering the radii at which the peak flux exceeds the map's noise level (fit_gauss_to_img_noise)
- currently no smoothing of the points is done (e.g., fit spline to ridge line)
- there can be problems if the radius (used for ridge line calculation)
lies completely within the jet/beam,
ref_pix
qualifier can be used to select another starting point (on the "jet axis")
Example
variable mdl = struct { flux = [0.3, 0.5, 0.1, 0.2, 0.3, 0.9, 3], ra = [3.7, 2.8, 2.2, 1.4, 0.9, 0.5, 0], dec = [1.2, 0.8, 0.8, 0.7, 0.6, 0.4, 0], smajor = [0.4, 0.2, 0.1, 0.03, 0.1, 0.01, 0.01], sminor = [0.4, 0.2, 0.1, 0.03, 0.1, 0.01, 0.01], pang = [0, 0, 0, 0, 0, 0, 0] }; variable beam = [0.3, 0.3, 0.7]; variable img = radio_mod2img (mdl, beam); variable r = ridge_line (img.img;); struct_filter(r, where(r.peak_value>1e-4)); plot_image (log(img.img+1e-4)); color(4); oplot(r.peak_x,r.peak_y); color(13); oplot(r.flux_cent_x,r.flux_cent_y);
See also: fit_gauss_to_img_noise, radio_mod2img, plot_vlbi_map
rndbknpwrlc
Synopsis
simulates a light curve with a broken power law distributed power spectrum
Usage
(t, rate) = rndbknpwrlc(Integer_Type nt);
Qualifiers
- beta1[=1.5]: first power law index of the power spectrum
- beta1[=1.0]: second power law index of the power spectrum
- perc[]
- mean[=0]: mean rate of the simulated light curve
- sigma[=1]: standard deviation of the simulated light curve
- dt[=1]: time resolution of the simulated light curve
Description
nt
is the number of bins of the simulated lightcurve.
It should be a power of two best performance of the fast Fourier transform.
rndpwrlc
Synopsis
simulates a light curve with a power law distributed power spectrum
Usage
(t, rate) = rndpwrlc(Integer_Type nt);
Qualifiers
- beta [=1.5]: power law index of the power spectrum
- mean [=0]: mean rate of the simulated light curve
- sigma [=1]: standard deviation of the simulated light curve
- dt [=1]: time resolution of the simulated light curve
Description
nt
is the number of bins of the simulated lightcurve.
It should be a power of two best performance of the fast Fourier transform.
See also: Timmer & Koenig (1995): "On generating power law noise", A&A 300, 707-710
sltable
Synopsis
Generate a LaTeX table
Usage
String_Type table = sltable(par1[,par2,par3,...])
Description
sltable() is intended to streamline the production of tables of, e.g., spectral parameters with error bars, although it should be flexible enough to produce most simply-structured tables. There is a set of examples on the Remeis wiki: http://www.sternwarte.uni-erlangen.de/wiki/doku.php?id=isis:sltable
By default this builds tables using the "deluxetable" package; this can be disabled to produce standard "tabular" tables by setting the "deluxe" qualifier to zero. The deluxetable environment has been tested using the latest aastex61.cls stylefile. Tabular tables will require the "booktabs" package, as they provide the \toprule, \midrule, and \bottomrule commands used here.
Each argument given to sltable() provides the values of one column of the table (or one row, if the "horiz" qualifier is nonzero). The length of each argument (or the length of the "value" field of a Struct_Type argument) defines the number of rows (or columns, if horiz=1). The arguments can be String_Type (in which case they are printed literally), Double_Type (or some other numerical type), in which case they are printed to three decimal places (TODO: figure out a better, more flexible way of doing this), or Struct_Type, as defined below.
Struct_Type arguments should at least have a "value" field, which should be an array of some type. If only the "value" field is present, the value printed to the table follows the rules for String_Type and Double_Type arguments above. Further fields which can be provided are: min: the minimum values of this parameter max: the maximum values of this parameter freeze: if 1, the parameter is frozen, if 0, it is thawed limit: if 1, this can be an upper/lower limit nodata: if 1, there is no data here and "..." should be printed These should all be arrays of the same length as the "value" field, and allow for the printing of error bars and indicating upper/lower limits and frozen parameters.
If the min and max fields are provided and the parameter is not frozen and not an upper or lower limit, the confidence intervals are printed using the TeX_value_pm_error() function. If it is an upper or lower limit, it is rounded to three decimal places and indicated accordingly with "<" or ">" symbols. If frozen, it is rounded to three decimal places and surrounded by parentheses (these delimiters can be customized by changing the "frozendelim" qualifier).
If the "nodata" field is nonzero, then a "no data" indicator will be printed. By default, for a deluxetable, this is "\nodata" and for a tabular table it is "\ldots". This is useful if you are putting multiple models which do not share the same parameters into the same table.
This is a lot of words, so here is an
EXAMPLE
This will produce a simple 3x3 table with values with error bars, some
frozen parameters, some upper limits, and some cells in the table filled
with literal strings.
variable par1 = struct{value=[1,2,3],min=[0.9,1.9,0.0],max=[1.1,2.3,3.2],
freeze=[0,1,0],limit=[1,1,1]};
variable par2 = struct{value=[1,2,3],min=[0.9,1.9,2.9],max=[1.1,2.3,3.2],
freeze=[1,0,0],limit=[1,1,1]};
variable par3 = struct{value=["foo","bar","baz"]};
variable parnames = ["par1","par2","par3"];
variable parunits = ["unit1","unit2","unit3"];
variable obsNames = ["obs1","obs2","obs3"];
variable notes = ["note1","note2"];
variable noteSym = ["1","2"];
variable t = sltable(par1,par2,par3;
colnames={parnames,parunits},rownames=obsNames,
notes={notes,noteSym},label="tab:example",
caption="this is the table's caption");
The variable "t" now contains the full LaTeX table - drop it into a LaTeX
document and see what it ends up looking like. Note the use of the
"colnames" and "rownames" qualifiers to define the names and units for the
columns and rows, and how List_Type and Array_Type values are interpreted
differently there. Footnotes are also possible via the "notes" qualifier
(although currently we can't place the footnote markers into the table
automatically - you'll have to handle that yourself).
This script is under development. Contact Paul Hemphill (pbh@space.mit.edu) if you find any bugs.
Qualifiers
- deluxe: Integer_Type. If nonzero, we'll produce a deluxetable. Otherwise, tabular. Tabular tables aren't as fancy and aren't implemented as well at the moment. In theory we could add other table types. Default 1.
- horiz: Integer_Type. If nonzero, arguments to sltable indicate rows of the table. If 0, they indicate columns.
- rownames: String_Type[] or List_Type[]. Names for the rows of the table. Multiple columns for the row names (e.g., parameter names and units in separate columns) can be handled by providing a List_Type containing multiple String_Type arrays.
- colnames: String_Type[] or List_Type[]. Names for the columns of the table (i.e., the column headers). Similar to rownames above, multi-line headers can be produced by passing a List_Type to this qualifier.
- caption: String_Type. LaTeX caption for the table. Default is blank. Don't include a label statement here - that's handled by the "label" qualifier below.
- label: String_Type. LaTeX label for the table. Default is "placeholder," so make sure to set this.
- nodata: String_Type, set to what should be printed if there is no data for a particular field (e.g., if you are putting multiple models into the same table and they do not share the same parameters).
- frozendelim: String_Type[], two-element string array containing opening and closing delimiters for a frozen value. By default this is ["(",")"], i.e., the frozen value is surrounded by parentheses.
- star: Integer_Type. If nonzero, the table will be a table* or deluxetable*. If 0, the table will be a table or deluxetable.
- tabletypesize: String_Type. Sets the font size of the table. Should be a LaTeX font size without the backslash. Default is "footnotesize".
- longtable: Integer_Type. If nonzero (and we're making a deluxetable), this is a long table, so a startlongtable statement will be included in the table header. I don't think tabular tables do anything with this right now. Default 0.
- landscape: Integer_Type. If nonzero, this is a landscape-orientation table. Includes a \rotate command if deluxetable, or wraps everything in a begin/end landscape block if tabular. Default 0.
- notes: List_Type[2]. Element 0 should be a String_Type[] array of the text of any footnotes; element 1 should be a String_Type[] array of the symbols associated with those notes.
See also: TeX_value_pm_error
sltableBuildBody
Synopsis
Build a 2D array of table cells
Usage
String_Type table[] = sltableBuildBody(Array_Type par1[, Array_Type par2,...])
Description
Builds the "body" of a table (in the form of a 2D array of strings, containing the fields of the table).
Arguments should be either 1D arrays or structs, the same as one would give to sltable(). An array of strings will be used as-is. Arrays of numbers will be rounded to 3 significant figures. Structs allow the handling of error bars and offer considerably more flexibility, as defined below:
Struct arguments should, at minimum, have a "value" field, which must be an array containing the value of the parameter for each row (or column, if the table is horizontally-aligned). Further fields can be used to modify the output. All should be arrays of the same length as the value field, and they have the following names and uses:
min: minimum value of the parameter. max: maximum value of the parameter. Min & max determine rounding & errors. freeze: if 1, parameter is frozen, if 0, parameter is thawed. limit: if 0, parameter can be zero without being an upper/lower limit.
Note that most of these are the same as those in the struct returned by get_par_info(). The exception is the "limit" field, which determines how a value of 0 for the parameter is handled. If "limit" is 1, then a confidence interval containing zero is interpreted as an upper or lower limit. If "limit" is 0, the value and error bars are typeset as a standard confidence interval. E.g., one would set limit to 1 for a power-law normalization and to 0 for a power-law index (the system is not flexible enough at the moment to handle pathological cases like upper/lower limits on logarithmic parameters like photon indices).
The number of "rows" in the table is determined by the length of the first array or the first struct's value field. All later arrays and struct fields should have the same number of elements.
Qualifiers
- horiz: Integer_Type. If zero, each argument defines one column of the table. If nonzero, each argument defines one row. Default 0.
- frozendelim: String_Type[2]. Defines strings to use before and after a frozen parameter value. Default ["(",")"].
- nodata: String_Type. Defines string to use when no data is present in a cell. Default is "\nodata" for deluxetables and "\ldots" for tabular tables.
- rownames: String_Type[] or List_Type. Names for rows of table. If given as a List_Type of String_Type arrays, each element of the list will be included in its own column (this is useful for, e.g., putting the parameter name and the units in separate columns).
See also: sltable
solarAbund
Synopsis
returns the solar abundance of an element
Usage
Double_Type solarAbund(Integer_Type Z)
Description
Z
is the element's proton number.
See also: Grevesse et al., 1996: Standard Abundances
stacking_vlbi_images
Synopsis
stack VLBI images of same size
Usage
stacking_vlbi_images(Array_Type <code>string(cleanfiles)</code>, String_Type <code>outputfile</code>);
Qualifiers
- ra_mas [={20,-20}]: specifies {left,right} limits of image in mas
- dec_mas [={-20,20}]: specifies {top,bottom} limits of image in mas
- plot_size [=15]: size of plot
- n_sigma [=3.0]: lowest contour of clean image
- sourcename [=default]: set the name of the source, by default the name is read from the .fits file, set to NULL for not plotting a source name
- obs_date1 [=default]: set the observation date, by default the date of the first observation is used
- obs_date2 [=default]: set the observation date, by default the date of the last observation is used
- cont_scl [="2"]: set factor to change the separation between contour levels
- cont_lvl [=[c1,c2,..]: set contour levels [Jy] manually, overwrites other contour parameters
- major: specify a vector for the major tic marks of the plot
- minor: specify a vector for the minor tic marks of the plot
Description
This function creates a stacked image of individual VLBI-clean images.
Important note: all images need to be of the same size! Different beam sizes are NOT recognized or corrected.
The required input format are fits-file for both the clean and the modelfit images. The format of the output file depends on the suffix of the given
filename
. Possible formats of the output file are PDF, EPS,
PNG, GIF, etc.
struct_print_field_statistics
Usage
struct_print_field_statistics(Struct_Type s);
See also: moment
tcp_client
Synopsis
connects to a TCP server
Usage
Struct_Type tcp_client(String_Type host[, Integer_Type port]);
or
Integer_Type tcp_client(String_Type host[, Integer_Type port];
msg_handler = Ref_Type, obj_handler = Ref_Type);
Qualifiers
nogreet - disable greet messages (see below) chatty - chattiness (default: 1)
Further qualifiers are listed at HOOKS below.
Description
Uses the socket' module to connect to a TCP server. After a connection attempt, a greet message is exchanged with the server by default unless the
nogreet' qualifier is set. In the latter
case the tcp_server' has be set up with the same qualifier. After the connection has been established, a structure with functions for communicating with the server is returned. Alternatively, if the
msg_handlerand
obj_handler` hooks are
set (see the HOOK section below) the client remains in a main
loop until the connection to the server has been closed. In this
case, the function returns 0 if any error occured during the
lifetime of the connection or 1 otherwise. The client structure,
which is passed to all hooks and returned if these hooks are no
set, has the following fields:
Integer_Type &send(String_Type msg or Any_Type obj)
Function for sending a string message or an SLang object
to the server (see pack_obj' for list of supported types). Returns the number of bytes sent. Warning: Currently, sending objects is slow (~40 kB/s, i.e., ~20 s for 100,000 doubles). Qualifiers: receipt - wait for a confirmation by the client for the receipt (the server's
receive' function needs the
same qualifier to be set).
BString_Type or Any_Type &receive()
Wait for and receive a message from the server. Returns
either the message itself (BString_Type) or in case of a
received object the object itself.
Qualifiers:
dontwait - do not wait for a message and return NULL in
case no message or object is pending.
receipt - send a receipt to the client after having
received a message or an object
Integer_Type &disconnect()
Disconnect from the server, i.e., close the communication
socket. Returns 1 in success or 0 otherwise.
Integer_Type connected - boolean value indicating whether
the connection is still established (1) or not (0)
Struct_Type config - internal configuration
Struct_Type hook - defined hooks
HOOKS The following functions are called via qualifiers of the same name, which need to be set to references (Ref_Type) to user- defined functions (see below for an example). The passed `client' is the client's structure as defined above.
String_Type or Integer_Type &greet_hook(client, greetmsg)
Called after having received the greet message from the server
(if not disabled by the nogreet' qualifier). Should return the greet message to be sent back to the server. Return 0 in case the greet message from the server should be rejected, which will cancel the connection. Void_Type &connect_hook(client) Called after a connection to the server is established. Can be used to start the communication with the server in combination with the
msg_handler' hook.
Void_Type &closed_hook(client)
Called after the connection has been closed from the server's
side.
Void_Type &msg_handler(client, msg)
Called after a message (BString_Type) has been received from
the server.
Void_Type &obj_handler(client, obj)
Called after an object has been received from the server.
Example
% simple functions for receiving messages and objects define getmsg(c, msg) { vmessage("message from server: %s", msg); }; define getobj(c, obj) { message("received an object:"); print(obj); }; % start the client variable stat = tcp_client("localhost"; msg_handler = &getmsg, obj_handler = &getobj); % check its exit status if (stat == 0) { message("lost connection to server"); }
% a more sophisticated example can be found on the Remeis wiki % www.sternwarte.uni-erlangen.de/wiki/doku.php?id=isis:socket
See also: tcp_server, socket, unpack_obj
tcp_get_server_handle
Synopsis
get a handle to a running TCP server
Usage
Struct_Type tcp_get_server_handle();
Qualifiers
kill - shut down / kill the server
Description
In case the handle, i.e., the structure of a running
tcp_server' got lost by accident, this function returns this structure again. Alternatively, the
kill' qualifier tries to shut down the server. If
this fails all communication sockets are closed,
which basically means to kill the server.
See also: tcp_server
tcp_server
Synopsis
implements a basic TCP server
Usage
Integer_Type tcp_server();
or Struct_Type tcp_server(; background);
Qualifiers
port - bind port (default: 1149) maxclients - maximum allowed number of clients connected simultaneously (default: 10) nogreet - disable greet messages (see below) background - run server in the background (see below) ip - bind address (default: "0.0.0.0", i.e., accept connections from outside; only change if you know what you are doing) chatty - chattiness (default: 1)
Further qualifiers are listed at HOOKS below.
Description
Uses the socket' module to implement a basic TCP server. It listens for clients trying to connect and accepts a connection after a greet message has been exchanged successfully. The greet message can be disabled using the
nogreet' qualifier,
which in this case needs to be set for a `tcp_client' as well.
Furthermore, any message or object received from a client is
read automatically. After any of these actions taken by the
server, a user-defined function may be called in order to handle
these events further (see the HOOK section below).
The background' qualifier allows the server to be started in background mode, such that commands can still be entered into the ISIS prompt. This background mode is realized by calling the main server function, which handles client events, each second scheduling an
alarmsignal (SIGALRM). NOTE: the alarm signals are triggered in the main program loop of ISIS, i.e., while it's "working". That means the main server function is not executed while ISIS awaits input from the prompt (ISIS "sleeps" here). Thus, you might want to do at least something, e.g, press enter. WARNING: further alarms should not be scheduled using
alarm'
or `signal'. Otherwise, the server might crash.
If not in background mode, `tcp_server' returns either 1 after the server has been shut down successfully or 0 otherwise. Else the structure of the following form is returned. Note that this structure is also passed to each hook (see below):
Struct_Type[] client - list of clients Integer_Type &shutdown() Function to shut down the server. After all clients are requested to disconnect, all sockets are closed. Returns whether the shut down was successful (1) or not (0). Qualifiers: trigger - shut down the server at the end of the main loop, which does not interrupt ongoing communication Void_Type &cleanclients() Function to remove clients no longer connected from the list. Integer_Type running - boolean value indicating the server's state (1: running; 0: shut down) Struct_Type config - internal configuration Struct_Type hook - defined hooks
Each `client' structure has the following fields:
Integer_Type &send(String_Type msg or Any_Type obj)
Function for sending a string message or an SLang object
to a client (see pack_obj' for list of supported types). Returns the number of bytes sent. Warning: Currently, sending objects is slow (~40 kB/s, i.e., ~20 s for 100,000 doubles). Qualifiers: receipt - wait for a confirmation by the client for the receipt (the client's
receive' function needs the
same qualifier to be set).
BString_Type or Any_Type &receive()
Wait for and receive a message from the client. Returns
either the message itself (BString_Type) or in case of a
received object the object itself.
Qualifiers:
dontwait - do not wait for a message and return NULL in
case no message or object is pending.
receipt - send a receipt to the client after having
received a message or an object
Integer_Type &kick()
Kick the client, i.e., close the communication socket.
Returns 1 in success or 0 otherwise.
Integer_Type connected - boolean value indicating whether
the client is still connected (1) or not (0)
Integer_Type number - the client's number, which is a serial
increasing identifier starting with zero
Struct_Type config - internal configuration
BString_Type[] unhmsg - in case the msg_handler'-hook is not set, messages of the client, which have been received without a call to
receiveduring the server's main loop are appended here List_Type unhobj - in case the
obj_handler'-hook is not set,
objects without a call to `receive' are appended here
HOOKS
The following functions are called via qualifiers of the same
name, which need to be set to references (Ref_Type) to user-
defined functions (see below for an example). The passed
server' and
client' are the server's and corresponding
client's structures as defined above.
Integer_Type &connect_hook(server, client)
Called after a client tries to connect to the server. Should
return 1 or 0 for accepting or rejecting the connection,
respectively.
String_Type &greet_hook(server, client)
or Integer_Type &greet_hook(server, client, greetmsg)
Called after accepting a client's connection, which should
return the greet message to be sent to the client. Is called
again after the re-greet message has been received from the
client, which should return 1 or 0 for accepting or rejecting
the greet, i.e, the connection, respectively. The hook is not
called if the server was set up with the nogreet' qualifier. Void_Type &established_hook(server, client) After a connection has been accepted including the greet, this function allows to start the communication between the client and the server by, e.g., sending messages to the client in combination with the
msg_handler' hook.
Void_Type &disconnect_hook(server, client)
Called after a client disconnected from the server by itself,
i.e., the client has not been kicked.
Void_Type &msg_handler(server, client, msg)
Called after a message (BString_Type) has been received from
the client. If this hook is not set, then the message is
appended to the unhmsg
field of the client's structure.
Void_Type &obj_handler(server, client, obj)
Called after an object has been received from the client. If
this hook is not set, then the object is appended to the
unhobj
field of the client's structure.
Example
% after a client's connection has been accepted, % send a data structure and disconnect client define do_client(s,c) { vmessage("sending data to client #%d", c.number); ()=c.send("welcome client! sending data..."); ()=c.send(struct { data = [1:5]*10, err = [1:5] }); ()=c.kick(); % remove kicked client from the list s.cleanclients(); } % start the server variable stat = tcp_server(; established_hook = &do_client); % check its exit status if (stat == 0) { message("server exited unexpectedly"); }
% a more sophisticated example can be found on the Remeis wiki % www.sternwarte.uni-erlangen.de/wiki/doku.php?id=isis:socket
See also: tcp_client, socket, alarm, pack_obj
transpose_table
Synopsis
transforms a structure of arrays into an array of structures and vice versa
Usage
Array_Type table_by_rows = transpose_table(Struct_Type table_by_columns);
or
Struct_Type table_by_columns = transpose_table(Array_Type table_by_rows);
Description
table_by_columns.field[i] = table_by_rows[i].field;
% for each field
and index i
voigt_profile
Synopsis
computes a normalized Voigt profile
Usage
Double_Type voigt_profile(Double_Type x, [[N,] x0,] sigma, gamma)
Description
The Voigt profile V is the convolution
of a Gaussian profile (with width sigma>0
)
G(x; sigma) = exp{-x^2/(2 sigma^2)} / [sqrt(2 pi) sigma]
and a Lorentzian profile (with width gamma>0)
L(x; gamma) = gamma/pi / (x^2 + gamma^2)
i.e.:
V(x; sigma, gamma) = int G(t; sigma) * L(x-t; gamma) dt
The voigt_profile
function computes
N * V(x-x0; sigma, gamma)
i.e., x0
is the center of the Voigt profile (default: 0)
and N
is its normalization (default: 1), i.e.
int voigt_profile(x, N, x0, sigma, gamma) dx = N
.
Do not confuse this "Voigt profile" with the "Voigt function H(a, u)".
See also: voigt, get_cfun2
wcs_axes_plot
Synopsis
gets WCS coordinates of axes and returns xfig plot object
Usage
xfig_plot_object = wcs_axes_plot (String_Type image_filename);
or
```c
xfig_plot_object = wcs_axes_plot (image, wcs);
##### Qualifiers
* width: set plot width
* height [=14]: set plot height, unless width and height are both
set, the other one is calculated based on the
format of the image (currently CDELT1 != CDELT2
is not taken into account, but only #pixels)
* x1label [=WCS of major]: set ticlabels of x1-axis
* x2label [=0]: set ticlabels of x2-axis
* y1label [=WCS of major]: set ticlabels of y1-axis
* y2label [=0]: set ticlabels of y2-axis
* axis_color [="white"]: set the color of the axes
* x1major: set WCS values for major tics of x1-axis.
similar qualifiers exist for x2, y1, and y2-axes.
unless x1 and x2 are both set, the qualifier affects
both x-axes. the same holds for the y-axes and the
minor qualifiers.
* x1minor: set WCS values for minor tics of x1-axis.
similar qualifiers exist for x2, y1, and y2 axes
* return_axes: set qualifier to additionally return
data (Assoc_Type) including axes information.
the returned object <code>wcsaxes</code> can be accessed
with <code>wcsaxes["x1"]</code> (keys are ["x1","x2","y1","y2"])
* info: print the WCS (RA-DEC) ranges of the axes
(helpful for setting minor and major by hand)
##### Description
Creates an xfig plot object based on a sky image with WCS coordinates
and determines the major and minor tics as well as the tic labels.
The underlying coordinates in the xfig plot are pixels (with integers
corresponding to the pixel centers).
The tic marks of each axis correspond to the "leading" (most varying)
coordinate along this axis (see wcs.xtype). There can/will be problems
when poles (and similar things) are within the field of view.
It is possible (and recommended) to specify major and minor tics and
the ticlabels (e.g., in order to use sexagesimal labels) of each axis
using the available qualifiers.
The alternative usage allows to specify the image and the wcs
information separately. They can be given either as a string of the
filename or directly as the image (array) and the WCS struct. In the
latter case they can be modified before.
WARNINGS:
- check WCS coordinates (the wcsfuns module uses less keywords
than, e.g., ds9)
- when plotting sub-images (i.e., cutting the image before
calling this function), make sure that the WCS structure is
consistent (e.g., modified center pixel)
##### Example
variable filename = "img.fits";
variable p = wcs_axes_plot (filename);
p.plot_png (log(img); cmap="ds9b");
p.render("img.pdf");
% if, e.g., WCS has to be modified/corrected:
variable wcs = fitswcs_get_img_wcs(filename);
wcs.cdelt[1] \*= -1; % can be necessary if a keyword was not read
variable p = wcs_axes_plot ( filename, wcs );
% use sexagesimal axis labels (example for declination from 44° 8' - 44° 20')
variable dec_major = [8:20:2]; % arcmin
variable y1major = 44 + dec_major / 60.0;
variable y1minor = 44 + [8:20:1] / 60.0; % minor in 1' steps
variable dec_ticlabels = array_map (String_Type, &sprintf, "$44^\circ\,%02d'$", dec_major);
variable p = wcs_axes_plot (filename; y1major = y1major, y1minor = y1minor, y1label = dec_ticlabels);
#### write_spix
##### Synopsis
writes the structure provided by make_spix in a FITS file
##### Usage
```c
write_spix(Stuct_Type <code>spix</code>, String_Type <code>filename</code>);
See also: make_spix, read_spix, plot_spix
xfigplot_hardnessratio_grid
Synopsis
plots hardness-grid and/or hardness-datapoints
Usage
xfigplot_hardnessratio_grid(Struct_Type tracks);
or xfigplot_hardnessratio_grid(Struct_Type tracks, String_Type fits_hr_table
);
Qualifiers
- outputdir: String_Type, name of output directory, default: 'testplot.pdf'
- W: Integer_Type, width of new xfig-plot
- H: Integer_Type, height of new xfig-plot
- xrange: Double_Type[2], range on x-axis
- yrange: Double_Type[2], range on y-axis
- padx: Double_Type, percentage of padding of whole plot in x direction, default: 0.05
- pady: Double_Type, percentage of padding of whole plot in y direction, default: 0.05
- xlabel: String_Type, label on x-axis
- ylabel: String_Type, label on y-axis
- par1label: String_Type, label within plot of par1
- par1label_x: Double_Type, x-position of 'par1label' as percentage (needs to be set, if 'par1label' is given)
- par1label_y: Double_Type, y-position of 'par1label' as percentage (needs to be set, if 'par1label' is given)
- par1label_c: Integer_Type/String_Type, color of 'par1label' as RGB, default: 0x000000 (black)
- par2label: String_Type, label within graph of par2
- par2label_x: Double_Type, x-position of 'par2label' as percentage (needs to be set, if 'par2label' is given)
- par2label_y: Double_Type, y-position of 'par2label' as percentage (needs to be set, if 'par2label' is given)
- par2label_c: Integer_Type/String_Type, color of 'par2label' as RGB, default: 0x000000 (black)
- extralabel: String_Type, extra label within graph
- extralabel_x: Double_Type, x-position of 'extralabel' as percentage (needs to be set, if 'extralabel' is given)
- extralabel_y: Double_Type, y-position of 'extralabel' as percentage (needs to be set, if 'extralabel' is given)
- extralabel_c: Integer_Type/String_Type, color of 'extralabel' as RGB, default: 0x000000 (black)
- extralabel2: String_Type, extra label within graph
- extralabel2_x: Double_Type, x-position of 'extralabel2' as percentage (needs to be set, if 'extralabel' is given)
- extralabel2_y: Double_Type, y-position of 'extralabel2' as percentage (needs to be set, if 'extralabel' is given)
- extralabel2_c: Integer_Type/String_Type, color of 'extralabel2' as RGB, default: 0x000000 (black)
- cstart: Integer_Type/String_Type, start color of grid as RGB, default: 0xC1CDCD (gray)
- cstop: Integer_Type/String_Type, stop color of grid as RGB default: 0xC1CDCD (gray)
- lend1: if exists, switch end (or begin) of tracks1 (where each tracks' value will be assigned)
- lend2: if exists, switch end (or begin) of tracks2 (where each tracks' value will be assigned)
* shift1x: Double_Type, x-position of par1-values at end of track, default: 0.
* shift1y: Double_Type, y-position of par1-values at end of track, default: 0.
* shift2x: Double_Type, x-position of par2-values at end of track, default: 0.
* shift2y: Double_Type, y-position of par2-values at end of track, default: 0.
* hr_x_col: String_Type, name of hr_x-column ('fits_hr_table') for plotting datapoints
* hr_y_col: String_Type, name of hr_y-column ('fits_hr_table') for plotting datapoints
* err_x_col: String_Type, name of error_x-column ('fits_hr_table') for plotting datapoints
* err_y_col: String_Type, name of error_y-column ('fits_hr_table)' for plotting datapoints
* err_c: Integer_Type/String_Type, color of errorbars as RGB
default: 0x000000 (black)
* symbol_shape: String_Type, shape of datapoints
* symbol_size: Double_Type, size of datapoints, default: 1.
* symbol_cstart: Integer_Type/String_Type, (start) color of datapoints as RGB
* symbol_cstop: Integer_Type/String_Type, stop color of datapoints as RGB,
if not given 'symbol_cstart' is color of all
* source_col: String_Type, name of source-column ('fits_hr_table') for plotting datapoints,
needs to be given only if different colors for different sources
* extralabel_sources_x: Double_Type, x-position of source_names as percentage
(if given, source names are displayed in their respective color)
* extralabel_sources_y: Double_Type, y-position of source_names as percentage
* tracks1: Integer_Type, if given, only tracks1 are plottet
* tracks2: Integer_Type, if given, only tracks2 are plottet
* notracks: Integer_Type, if given, no tracks are plottet, 'fits_hr_table' must be
given for datapoints
##### Description
- for only grid: give only tracks
- for grid with datapoints: give tracks and FITS-table with hardnessratios (and errors)
- for only datapoints: give tracks and FITS-table, but use 'notracks'-qualifier
__See also__: hardnessratio_from_dataset