Reference for slang functions
Aitoff_projection
Synopsis
computes an Aitoff projection
Usage
(Double_Type x, y) = Aitoff_projection(Double_Type l, b);
Qualifiers
- deg:
l
andb
are in degrees, not in radian - normalized:
x
andy
are normalized (byPI/2
) such thatabs(x) <= 2
andabs(y) <= 1
.
Description
NOTE: for astronomical purposes you probably want to use the equal area Hammer-Aitoff projection rather than an Aitoff-projection. Erroneously many astronomers call the Hammer-Aitoff projection an Aitoff projection
x = 2 * cos(b) * sin(l/2) / sinc(alpha);
y = sin(b) / sinc(alpha);
where
cos(alpha) = cos(b) * cos(l/2)
and
sinc(alpha) = sin(alpha)/alpha
See also: Hammer_projection, Lambert_Equal_Area_projection
AzEl_from_RAdec
Synopsis
computes horizontal (azimut, elevation) coordinates from a point's (RA, dec) at a time MJD
Usage
(az, el) = AzEl_from_RAdec(RA, dec, MJD, longw, lat);
Description
This function is deprecated. Please use equatorial2horizon instead.
See also: equatorial2horizon,horizon2equatorial
B_lambda
Synopsis
Calculates the spectral density for black body radiation (wavelength space)
Usage
Double_Type Blambda = B_lambda(lambda,T);
Qualifiers
- keV: T argument is kT in keV
- A: Wavelength is given in Angstroms
- Angstrom: Wavelength is given in Angstroms
Description
This function returns the spectral energy density of a black body, in units erg/cm^2/s/cm/sr Either lambda or T can be an array.
B_nu
Synopsis
Calculates the spectral density for black body radiation (frequency space)
Usage
Double_Type Bnu = B_nu(nu,T);
Qualifiers
-
keV: 1st argument is energy in keV, T argument is kT in keV
-
MHz: Frequency is given in MHz
-
GHz: Frequency is given in GHz
Description
This function returns the spectral energy density of a black body, in units erg/cm^2/s/Hz/sr Either nu or T can be an array.
BinaryCor
Synopsis
Removes the influence of the doublestar motion for circular or eliptical orbits.
Usage
Array_Type t = BinaryCor(Array_Type OR Double_Type time), time in MJD
Qualifiers
- asini: Projected semi-major axis [lt-secs], Mandatory
- porb: Orbital period at the epoch [days], Mandatory
- eccentricity: Eccentricity, (0<=e<1)
- omega: Longitude of periastron [degrees], mandatory
- t0: epoch for mean longitude of 0 degrees (periastron, MJD)
- t90: epoch for mean longitude of 90 degrees (MJD)
- pporb: rate of change of the orbital period (s/s) (default 0)
- limit: absolute precision of the correction of the computation (days, default: 1D-6)
- maxiter: stop Newton-Raphson iteration after maxiter steps if limit is not reached (default: 20)
Description
(transcribed from IDL-programm BinaryCor.pro)
For each time, the z-position of the emitting object is computed and the time is adjusted accordingly. This is iterated until convergence is reached (usually only one iteration is necessary, even in high elliptic cases).
Follows equations from Hilditch's book and has also been checked against fasebin/axBary. All codes give identical results, (to better than 1d-7s) as checked by a Monte Carlo search using 1d7 different orbits.
qualifiers t90 and t0 have to be in days and in the same time system as time (e.g. JD or MJD)
Circular orbits: * if time of lower conjunction Tlow is known, set t0=Tlow and omega=0 * if time of ascending node is known, Tasc, set t90=Tasc and omega=0 * if time of mid eclipse is known, Tecl, set t0=Tecl-0.25*porb and omega=0
BinaryPos
Synopsis
computes the position of a star in a binary system
Usage
Double_Type[] BinaryPos(Double_Type[] time)
Qualifiers
- porb: Orbital period at the epoch [same units as time argument], mandatory
- eccentricity: Eccentricity, (0<=e<1)
- omega: Longitude of periastron [degrees], mandatory
- asini: Projected semi-major axis (a sin i), mandatory for spectroscopic binaries, give it in light seconds (a sini/c) if you want to compute a time correction
- semi: semimajor axis, mandatory for visual binaries
- bigomega: Longitude of the ascending node [degrees], mandatory for visual binaries
- incl: inclination [degrees], mandatory for visual binaries
- t0: epoch for mean longitude of 0 degrees (periastron, MJD)
- t90: epoch for mean longitude of 90 degrees (MJD)
- pporb: rate of change of the orbital period (s/s) (default 0)
Description
(transcribed from IDL-programm BinaryPos.pro)
For spectroscopic binaries, the z-position of the emitting object is computed (negative: closer to Earth). In case of visual binaries the output is a 3d array containing the xyz position of the object in a right handed coordinate system where the z axis points away from Earth, the x axis points towards the North Celestial Pole, and where the y-axis increases towards East.
Follows equations from Hilditch's book.
Qualifiers t90 and t0 have to be in days and in the same time system as time (e.g. JD or MJD)
Circular orbits: * if time of lower conjunction Tlow is known, set t0=Tlow and omega=0 * if time of ascending node is known, Tasc, set t90=Tasc and omega=0 * if time of mid eclipse is known, Tecl, set t0=Tecl-0.25*porb and omega=0
CCF_1d
Synopsis
calculates the cross correlation function of two (1d) arrays
Usage
Double_Type CCF = CCF_1d(a1, a2, Integer_Type dx);
or
Array_Type CCF = CCF_1d(a1, a2);
Qualifiers
- notperiodic: arrays are not extended periodically.
Description
The arrays a1
and a2
have to have the same length n
.
They can be passed directly or as references to arrays.
By default (unless the notperiodic
qualifier is given), both arrays
are considered to be extended with periodic boundary conditions.
When using notperiodic
, chose |dx| < n-1
, since the larger |dx|
,
the smaller are the parts of the arrays that can actually be compared.
This function computes the correlation of a1[x]
and a2[x-dx]
:
CCF = sum_x { (a1[x]-)/|a1| * (a2[x-dx]-)/|a2| }
For periodic boundary conditions, the entire array is considered,
else the index x
takes only values such that x
and x-dx
are contained.
The norm |a| = sqrt{ sum_x a[x]^2 }
ensures that for any array a
,
CCF_1d(a, a, 0) = +1
and CCF_1d(a, -a, 0) = -1
.
If only the two arrays a1
and a2
are given, the function returns
an array with all possible correlations between them, i.e.:
[ CCF_1d(a1, a2, 0), CCF_1d(a1, a2, 1), ..., CCF_1d(a1, a2, n-1) ]
This is most reasonable for periodic boundary conditions.
Notes
The function is vectorized on dx
, i.e., can handle an array dx
.
Example
variable a1 = Double_Type[n]; % some data variable a2 = Double_Type[n]; % some other data
variable dx = [-n/2 : n/2]; variable CCF_np = CCF_1d(a1, a2, dx; notperiodic);
variable CCF_p = CCF_1d(a1, a2); % This is equivalent to: variable CCF_p = array_map(Double_Type, &CCF_1d, &a1, &a2, [0:n-1]);
See also: CCF_2d
CCF_2d
Synopsis
calculates the cross correlation function of two 2d-arrays (e.g., images)
Usage
Double_Type CCF = CCF_2d(img1, img2, Integer_Type dx, Integer_Type dy);
Description
The 2d-arrays img1 and img2 from which the CCF is to be computed can
either be passed directly or by a reference to the arrays.
This function computes the correlation of img1[y, x]
and img2[y-dy, x-dx]
,
where the indices (x,y) take all values such that (x,y) and (x-dx,y-dy) are contained:
CCF = sum_{x,y} { (img1[y,x]-)/|img1| * (img2[y-dy,x-dx]-)/|img2| }
The norm |img| = sqrt{ sum_{x,y} img[y,x]^2 }
ensures that CCF_2d(img, +-img, 0, 0) = +-1
.
for any image img
.
If dx
or dy
are arrays, CCF
will be a 2d array.
Qualifiers
- verbose: turns on verbosity
- savebest=&best: if
dx
ordy
are arrays, the best combination will be saved inbest
.
See also: CCF_1d
COPY
Synopsis
makes a copy of a nested Struct/Array/List_Type
Usage
Struct_Type copy = COPY(Struct_Type copyme);
or
```c
Array_Type copy = COPY(Array_Type copyme);
or
List_Type copy = COPY(List_Type copyme);
##### Description
<code>COPY</code> returns a properly dereferenced copy of an
arbitrarily nested <code>Struct/Array or List_Type</code> with
all its entries.
Depending on the Data_Type of the initial given
copyme and the Data_Type of its entries, <code>COPY</code>
iteratively calls the according sub-function
(<code>struct_copy</code>, <code>array_copy</code> or <code>list_copy</code>).
If the given Data_Type or one of its entries is not
one of privouse mentioned ones, <code>COPY</code> returns
@(Data_Type) or respectively @(entry).
NOTE: \* Ref_Type entries will not be dereferenced, which
allows to copy XFIG-OBJECTS !!!
\* Double_Type[] is an Array_Type
##### Example
Array_Type:
s = Struct_Type[1]; s[0]=struct{ a=Array_Type[1,2] };
s[0].a[[0]]=[0:9];
copy = COPY(s); copy[0].a[[0]] = ["modified"];
print(s[0].a);
List_Type:
s = Struct_Type[1]; s[0]=struct{ a=Array_Type[1,2] };
s[0].a[[0]]=[0:9];
l = [{ array_copy(s), 1., [1:10] }, {"abs"}];
copy = COPY(l); copy[0][0][0].a[[0]] = ["modified"];
print(l[0][0][0].a);
Struct_Type:
s = Struct_Type[1]; s[0]=struct{ a=Array_Type[1,2] };
s[0].a[[0]]=[0:9]; S = struct{ f=s };
C = struct_copy(S); C.f[0].a[[0]] = ["modified"];
print(S.f[0].a);
__See also__: struct_copy, array_copy, list_copy
#### Compton_crosssection
##### Usage
```c
Double_Type sigma = Compton_crosssection(Double_Type E);
Description
E
is the energy in keV.
sigma
is the total Klein-Nishina crosssection in cm^2:
sigma = 3/4 sigma[Th] * [ (1+y)/y^3 * { 2y*(1+y)/(1+2y) - log(1+2*y) } + log(1+2*y)/(2*y) + (1+3y)/(1+2y)^2 ]
where y = E/(511 keV)
.
See also: Rybicki & Lightman (1979)
DateOfJD
Synopsis
Calculates the date corresponding to a JD value
Usage
Struct_Type date = DateOfJD(JD);
Description
The return value is a structure of the following form:
date = struct { year, month, day, hour, minute, second };
The algorithm used is that by Meeus, Astronomical Formulae for
Calculators. This routine is array safe.
Qualifiers
- julianswitch: For JD smaller than this date, return the date in the Julian calendar rather than in the Gregorian calendar. The default is JD2299161 (15 October 1582), which is appropriate for most of catholic Europe. Note, however, that there are countries that switched only much later (see https://en.wikipedia.org/wiki/Julian_calendar), for example, Russia only switched on 1 February 1918 (Gregorian).
}
See also: JDofDate, MJDofDate
DateOfMJD
Synopsis
calculates the date from its MJD value
Usage
Struct_Type date = dateOfMJD(MJD);
Description
The return value is a structure of the following form:
date = struct { year, month, day, hour, minute, second }
and always in the Gregorian Calendar. See DateOfJD if you want more
functionality. This function is array safe.
See also: DateOfJD, MJDofDate
EpochofJD
Synopsis
Returns the Julian or Besselian epoch of a (M)JD
Usage
epoch=EpochofJD(jd;qualifiers)
Qualifiers
- mjd: argument is a modified julian date, not a JD
- besselian: return the Besselian epoch
- julian: return the Julian epoch (the default)
Description
This function converts a (modified) julian date into the corresponding Julian or Besselian epoch (e.g., J2000, B1950.0 etc.) using the formulae given by Lieske (1979, A&A 73, 282) for Besselian dates and the standard definition of the Julian Date.
If high precision is relevant, note that JD is assumed to be in TT, note that the epoch will slowly drift with respect to Gregorian year fractions such as those calculated by jd2year, but contrary to that function the Julian/Besselian epoch linearly maps to TT.
This routine is array safe.
See also: JDofEpoch,jd2year
Faddeeva
Synopsis
Compute w(z) = exp((-iz)^2 erfc(-iz) for complex z
Usage
Complex_Type[] Faddeeva(Complex_Type[]);
Description
This function uses continued fraction expansion and the algorithms described by Humlicek () and Hui () to compute an approximation to the Faddeeva function
The algorithm used only allows to compute w(z) for Im(z)>=0, but using the relation w(-z) = 2*exp(-z^2)-w(-z) gives the remaining half.
The derivative is given via dw/dz = 2i/sqrt(pi) - 2*z*w(z)
See also: Faddeeva_dz
It is claimed that the algorithm has an accuracy of <4e-4.
Faddeeva_dz
Synopsis
Compute derivative of Faddeeva function
Usage
Complex_Type[] Faddeeva_dz (Complex_Type[]);
Description
Computes complex derivative of Faddeeva function
See also: Faddeeva
GMST
Synopsis
Calculates the Greenwich sidereal time and the local siderial time
Usage
Double_Type GST = GMST(Double_Type JD);
Qualifiers
- hour: Return the GMST in hours, not in seconds.
- dut1: The difference UT1-UTC in seconds (see IERS Bull. A).
- mjd: Argument is in MJD, not in JD.
- era: Return the Earth rotation angle (in rotations).
- lst: Return the local mean siderial time LST.
- lon: longitude for the calculation of the LST (east is positive; in geodetic [ITRS] coordinates, for practical purposes these are identical to the WGS-84 coordinates given by GPS).
- lat: latitude for the calculation of the LST (north is positive, only needed if xp and yp are given)
- deg: lon and lat are in degrees, not in radian.
- xp: x-coordinate of the Celestial Intermediate Pole (CIP) ALWAYS in arcseconds, as given by the IERS bulletin.
- yp: y-coordinate of the CIP, ALWAYS in arcseconds.
Description
This code calculates the Greenwich Mean Siderial Time in seconds.
The argument JD
is the Julian Date in the UT1 time system.
For most applications, one can take UT1=UTC. For the most precise
applications, a correction term DUT1=UT1-UTC can be looked up in
IERS Bulletin A
(see https://www.iers.org/IERS/EN/Publications/Bulletins/bulletins.html).
The algorithms used are described by Kaplan (2005, USNO Circular 179, Sect. 2.6.2).
This routine is array safe.
Note that contrary to the coordinate transform routines the deg qualifier only affects the interpretation of the lon, lat qualifiers and does NOT affect the other arguments or return arguments!
Gauss_complex
Synopsis
Compute complex Gauss profile
Usage
Complex_Type = Gauss_complex(z, z0, sigma);
GhoshLamb79
Synopsis
calculates the spin-up of a neutron star after Ghosh & Lamb (1979)
Usage
Double_Type GhoshLamb79(Double_Type pulse_period; disk);
or Double_Type GhoshLamb79(Double_Type pulse_period; wind);
Qualifiers
neutron star/binary parameters L, M, R, B, mu, V0, cs, Vorb, a, xi (see the description for details) fastness - hook into the issue of a fastness parameter >=0.9. If this qualifier is set to "limit" then the fastness is set to 0.9 in this case. If set to a reference function (Ref_Type) then a hook function of the form Double_Type fastness_hook(fastness, P, L, M, R, B) is called, which is given the computed 'fastness' and all neutrons star parameters. It has to return the new value for the fastness.
Description
This implements the accretion torque theory by Ghosh & Lamb, 1979, ApJ 234, 296 (here after GL79) which calculates the spin-up of a neutron star spinning with a period of 'pulse_period' (given in seconds). The returned \dot{P} will be in seconds per seconds.
The theory depends on the physical parameters of the neutron star and the accretion mode. All of these are controlled by qualifiers.
Two accretion modes are available (set as corresponding qualifier): disk - accretion from an accretion disk (Eq. 15 of GL79) wind - accretion from a stellar wind (Eq. 24 of GL79)
Three neutron star parameters are used independently from the mode: L - the X-ray luminosity (in 10^37 erg/s; default: 1) M - the neutron star mass (in solar masses; default: 1.4) R - the neutron star radius (in km; default: 11.5)
For the disk-mode the magnetic field strength is needed. One of the following quantities has to be given: mu - magnetic moment (in 10^30 G cm^3; default: 1) B - surface magnetic field strength (in 10^12 G; overwrites mu) For the wind-mode several parameters are needed: V0 - the capture velocity (see Eq. 20 of GL79) (in km/s; default: 1000) cs - the sound speed at the capture radius (in km/s; default: 0) Vorb - the orbital velocity (in km/s; default related to: V0^2 = Vorb^2 + Vwind^2 + cs^2 assuming Vwind = Vorb) a - the binary separation (in lt-s; default: 100) xi - "on the order of unity" (see Eq. 22 of GL79; default: 1)
The function should be vectorized, i.e., an array of pulse periods can be given or the qualifiers can be arrays.
In case of disk accretion, the so-called fastness parameter is computed after Eq. 16 of GL79. This equation is accurate by 5% for fastnesses <0.9. For larger values a warning message is displayed. For fastnesses >=1.0 the solution is a complex number and therefore a _NaN value is returned and an error message is displayed. The 'fastness' qualifier can be used to hook into this issue.
See also: pulsarGL79
GreenwichSiderealTime_from_MJD
Synopsis
computes the Greenwich sidereal time
Usage
Double_Type GST = GreenwichSiderealTime_from_MJD(Double_Type MJD);
Description
MJD
is the Modified Julian Date (*in UT*).
This routine is a wrapper around GMST(JD), and only
kept for compatibility reasons.
This function is deprecated. Please do not use for future work.
See also: GMST
Hammer_projection
Synopsis
Computes the Hammer-Aitoff projection
Usage
(Double_Type x, y) = Hammer_projection(Double_Type l, b);
Qualifiers
- deg:
l
andb
are in degrees, not in radian - normalized:
x
andy
are normalized (bysqrt(2)
) such thatabs(x) <= 2
andabs(y) <= 1
. - astronomical: flip x-axis for astronomical maps, where east is to the left
- inverse: calculate the inverse projection, interpreting l as the x- and b as the y-coordinate; all other qualifiers are also interpreted as expected.
Description
This is the projection erroneously called the Aitoff projection by many astronomers. It is an equal area projection. The projection equations are
x = 2 * sqrt(2) * cos(b) * sin(l/2) / sqrt(1 + cos(b)*cos(l/2));
y = sqrt(2) * sin(b) / sqrt(1 + cos(b)*cos(l/2));
The inverse function returns nan if the arguments given are not possible.
See also: Aitoff_projection, Lambert_Equal_Area_projection
JD2MJD
Synopsis
Convert Julian Dates to Modified Julian Dates.
Usage
JD2MJD(JD);
Description
Helper function to convert Julian Dates to MJD by subtracting off 2400000.5.
See also: j2dmjd, mjd2jd
JDofDate
Synopsis
Calculates the Julian date
Usage
Double_Type JD = JDofDate( struct { year, month, day, hour, minute, second } );
or
Double_Type JD = JDofDate(year, month, day[, hour[, minute[, second]]]);
Description
This routine calculates the Julian Date for a given Gregorian date following the routine by Fliegel and van Flandern (1968, Comm. ACM 11, 657) and further optimized as discussed by Pulkkinen and van Flandern (1979, ApJ Suppl 41: 391). The date is either given by a structure or by at least 3 integer arguments. Fractional parts are ignored. The default values for hour, minute and second are 0, fractional hours etc. are properly taken into account. This function is array safe.
See also: MJDofDate,dateOfMJD
JDofEpoch
Synopsis
Returns the (M)JD of a Julian or Besselian epoch such as J2000.0
Usage
jd=JDofEpoch(epochstring;qualifiers)
Qualifiers
- mjd: Return the modified julian date (default: JD)
Description
This function returns the julian date corresponding to a Julian or Besselian epoch (e.g., J2000, B1950.0 etc.). For highest precision, the most commonly used values are hardcoded, other values are calculated using the formulae given by Lieske (1979, A&A 73, 282) for Besselian dates and the standard definition of the Julian epoch (based on Julian years with 365.25 days).
The epoch string must start with a "B" for a Besselian epoch and "J" for a Julian epoch, followed by a number that designates the year fractional values such as 2017.25 are possible.
Note that the time system in which the dates are defined is TT, for most applications that do not require microsecond accuracy for the conversion of positions the precise time system won't matter.
If the routine is called with a floating point number, it is assumed that this number is the epoch and this number is returned. This greatly simplifies functions in which the epoch is either a string or a number.
See also: EpochofJD,tt2tai,tt2tdb
KeplerEquation
Synopsis
solves Kepler's Equation E - e sin E = M
Usage
Double_Type E = KeplerEquation(Double_Type M, e)
Description
Solve Kepler's Equation by using the result of the method by S. Mikkola (1987) Celestial Mechanics, 40, 329-334 as starting value for a Newton-Raphson iteration to extend the applicability of this function to higher eccentricities.
Qualifiers
- threshold: stopping criterion for the Newton Raphson iteration, default = 1e-5
Lorentz_complex
Synopsis
Compute complex Lorentz profile
Usage
Complex_Type[] = Lorentz_complex(z, z0, gamma);
MJD2JD
Synopsis
Convert Modified Julian Dates to Julian Dates.
Usage
MJD2JD(MJD);
Description
Helper function to convert MJD to Julian Dates by adding 2400000.5.
See also: jd2mjd, JD2MJD, MJD2JD
MJD2UNIXtime
Synopsis
converts modified Julian date to UNIX time
Usage
Long_Type MJD2UNIXtime(Double_Type MJD)
Description
Time formats and conversion functions: Long_Type secs = MJD2UNIXtime(Double_Type MJD); Double_Type MJD = UNIXtime2MJD(Long_Type secs);
Struct_Type tm = localtime(Long_Type secs); % or gmtime Long_Type secs = mktime(Struct_Type tm);
String_Type str = strftime(String_Type fmt, Struct_Type tm);
Current time and date: Long_Type secs = _time(); String_Type str = time();
See also: UNIXtime2MJD, localtime, gmtime, mktime
MJD2fermi
Synopsis
calculate Fermi MET used for LC and spectral analysis
Usage
Double_Type = MJD2fermi (Double_Type);
Description
Calculates the Fermi Mission Elapsed Time (MET) in seconds for a given MJD(UTC).
Example
variable m = 54900.345; variable fermi_s = MJD2fermi(m);
See also: fermi2MJD, MJDref_satellite
MJDofDate
Synopsis
calculates the modified Julian date
Usage
Double_Type MJD = MJDofDate( struct { year, month, day, hour, minute, second } );
or
Double_Type MJD = MJDofDate(year, month, day[, hour[, minute[, second]]]);
Description
The Gregorian date is either given by a structure or by at least 3 integer arguments. Fractional parts are ignored. The default values for hour, minute and second are 0, fractional hours etc. are properly taken into account.
See also: JDofDate,dateOfMJD
MJDofDateString
Synopsis
calculates the MJD from a string including year, month, day [, hour [, minute [, second]]]
Usage
Double_Type <code>mjd</code> = MJDofDateString(String_Type <code>date_string</code>);
Qualifiers
- verbose: prints the found values for YYYY-MM-DD hh:mm:ss.sss
- no_2digit_correction: see
parseDateString
Description
This function calculates the Modified Julian Date from a String
using the parseDateString
function.
Example
MJDofDateString("69/07/21"); MJDofDateString("1969-07-21 02:56:30.44"); MJDofDateString("sdf1969--$07..21QQ_02/56((30,,44xyz");
See also: parseDateString
MJDstrftime
Synopsis
formats a modified Julian date
Usage
String_Type MJDstrftime(String_Type fmt, Double_Type MJD);
See also: strftime, MJD2UNIXtime, localtime
N_body_simulation
Synopsis
Compute orbits for N interacting particles
Usage
Struct_Type ret = N_body_simulation(Struct_Type in, Double_Type t_end; qualifiers)
Description
By numerical integration and without any simplifying approximations, this function directly solves the equations of motion of a system of N particles under the influence of their mutual forces (The interaction is specified in a separate function, see qualifier 'kernel'.) from time t=0 to time t=t_end. Negative values for t_end imply backward integration. Cartesian coordinates are used throughout.
The input structure 'in' has to contain the fields "x", "y", "z", "vx", "vy", "vz". Each of these fields has to be an array of length N. The corresponding index gives the particle id, i.e., particle 0 refers to x[0], y[0], z[0], vx[0], vy[0], vz[0]. The return structure 'ret' contains for each particle a field, e.g., for particle 0 a field named "o0". This field is again a structure with fields "t", "x", "y", "z", "vx", "vy", and "vz" giving the time-dependent coordinates of the respective particle.
Notes
An adaptive Runge-Kutta-Fehlberg method of fourth/fifth order is applied to solve the coupled system of first-order differential equations. The stepsize is hereby controlled such that for each step an absolute accuracy in coordinates and velocity components is achieved that is smaller than given by the qualifier 'tolerance'.
Qualifiers
- kernel [
="N_body_simulation_std_kernel"
]: Name of the function which describes the mutual interaction. Note that all qualifiers are passed to this function as well. - threshold [
=0
]: Lower limit on the time difference of two consecutive moments of time that will be saved. - tolerance [
=1e-10
]: Absolut error control tolerance; lower limit: 1e-15. - verbose: Show intermediate times t.
Example
% Four interacting particles: s = struct{x, y, z, vx, vy, vz}; s.x = [-10,0,10,0]; s.y = [0,10,0,-10]; s.z = [0,0,0,0]; s.vx = [0,-1,0,1]; s.vy = [-1,0,1,0]; s.vz = [0,0,0,0]; r = N_body_simulation(s, 50; kernel="N_body_simulation_std_kernel", psa=10*[1,1,1,1], psb=0.1*[1,1,1,1]); xrange(min_max([r.o0.x,r.o1.x,r.o2.x,r.o3.x])); yrange(min_max([r.o0.y,r.o1.y,r.o2.y,r.o3.y])); plot(r.o0.x,r.o0.y); oplot(r.o1.x,r.o1.y); oplot(r.o2.x,r.o2.y); oplot(r.o3.x,r.o3.y);
See also: N_body_simulation_std_kernel, N_body_simulation_MW_kernel
N_body_simulation_std_kernel
Synopsis
Default interaction kernel for the function 'N_body_simulation'
Usage
Double_Type r[6,N] = N_body_simulation_std_kernel(Double_Types t, m[6,N]; qualifiers)
Description
This function is the default interaction kernel of the function 'N_body_simulation'. It computes the equations of motion of N interacting particles at time t. The input parameter 'm' is a [6,N]-matrix with m[0,j] = x_j; m[1,j] = y_j; m[2,j] = z_j; m[3,j] = vx_j; m[4,j] = vy_j; m[5,j] = vz_j;
For each particle, the potential of particle i (at position (x_i,y_i,z_i)) exerted on particle j (at position (x_j,y_j,z_j)) is a Plummer sphere of the form Phi(x_i,y_i,z_i,x_j,y_j,z_j) = -psa_i*(psb_i+dis^2)^(-1/2) with dis^2 = (x_i-x_j)^2+(y_i-y_j)^2+(z_i-z_j)^2 and psa_i and psb_i model constants (see qualifiers). The resulting acceleration of particle j is then d^2/dt^2 x_j = -d/dx_j Phi = sum( psa_i*(psb_i+dis^2)^(-3/2)*(x_i-x_j), i!=j) d^2/dt^2 y_j = -d/dy_j Phi = sum( psa_i*(psb_i+dis^2)^(-3/2)*(y_i-y_j), i!=j) d^2/dt^2 z_j = -d/dz_j Phi = sum( psa_i*(psb_i+dis^2)^(-3/2)*(z_i-z_j), i!=j) yielding the following system of first-order differential equations, which are the equations of motion of this system and which are returned as a [6,N]-matrix 'r' d/dt x_j = r[0,j] = vx_j d/dt y_j = r[1,j] = vy_j d/dt z_j = r[2,j] = vz_j d/dt vx_j = r[3,j] = -d/dx_j Phi = sum( psa_i*(psb_i+dis^2)^(-3/2)*(x_i-x_j), i!=j) d/dt vy_j = r[4,j] = -d/dy_j Phi = sum( psa_i*(psb_i+dis^2)^(-3/2)*(y_i-y_j), i!=j) d/dt vz_j = r[5,j] = -d/dz_j Phi = sum( psa_i*(psb_i+dis^2)^(-3/2)*(z_i-z_j), i!=j)
Qualifiers
- psa [
=Double_Type[N]+1
]: Parameter used to parametrize the interaction potential. - psb [
=Double_Type[N]+0
]: Parameter used to parametrize the interaction potential.
Example
N = 4; m = Double_Type[6,N]; r = N_body_simulation_std_kernel(0,m; psa=Double_Type[N]+1, psb=Double_Type[N]+1);
See also: N_body_simulation
RAdec_from_AzEl
Synopsis
computes equatorial (RA, dec) coordinates from a point's (azimut, elevation) at a time MJD
Usage
(RA, dec) = RAdec_from_AzEl(az, el, MJD, longw, lat);
Description
This function is deprecated. Please use horizon2equatorial instead.
See also: horizon2equatorial,equatorial2horizon
RAdec_from_galLB
Synopsis
convert galactic (l, b) coordinates to equatorial (RA, dec) coordinates
Usage
(RA, dec) = RAdec_from_galLB(l, b);
Qualifiers
- l_unit [
="deg"
]: set the unit of the galactic longitudel_unit="rad"
=>
l in radl_unit="hms"
=>
l in hours as a scalar or an array of the form[H, M]
or[H, M, S]
- b_unit [
="deg"
]: set the unit of the galactic lattitude
Description
This function is deprecated. Please use galactic2equatorial instead.
See also: galLB_from_RAdec, AzEl_from_RAdec, RAdec_from_AzEl, angle_to_rad
RD2rad
Synopsis
Convert right ascension (h, m, s) and declination (deg, arcmin, arcsec) to radians
Usage
RD2rad(Double_Types ah[], am[], as[], dd[], dm[], ds[])
Description
Given right ascension in hours, minutes, seconds and declination in degrees, minutes of arc, seconds of arc, the function converts them to radians. Always provide six numbers, i.e., fill up with zeros where necessary. Right ascension is assumed to be positive. If the degrees argument of declination is negative the minus sign is automatically applied to the minutes and seconds parameters. If the degrees argument of declination is zero, the sign of the minutes argument (if nonzero) is applied to the seconds parameter.
Example
(raInRad, declInRad) = RD2rad(01, 45, 12.54, 87, 55, 45.34); (raInRad, declInRad) = RD2rad(01, 45.209, 0, -87, 55, 45.34); (raInRad, declInRad) = RD2rad(01.753472, 0, 0, 0, -55, 45.34); (raInRad, declInRad) = RD2rad(01, 45, 12.54, 0, 0, 45.34); (raInRad, declInRad) = RD2rad([01, 01], [45,45.209], [12.54,0], [87,-87], [55,55], [45.34,45]); (raInRad, declInRad) = RD2rad(rad2RD(0,-PI/4.));
See also: hms2deg,dms2deg,rad2RD
Roche_critical
Synopsis
calculates the critical potential at L1 and the corresponding radius in z-direction
Usage
Struct_Type crit = Roche_critical(Double_Type mq);
Description
The Roche potential describes the motion of a binary consis- ting of two objects, M1 and M2. The mass ratio mq = M1/M2 must be given. From the position of the first Lagrange point 'xL1', which is a saddle point in the Roche potential, the value of the critical potential 'critpot' at this point is calculated. The corresponding equipotential surface describes the inner most closed surface, which is also known as 'Roche lobe'. Since the Roche potential does not depend on the binary rotation, the radius 'rcrit' of the star in z-direction (perpendicular to the orbital plane) is often used as the star's radius.
All three values are returned by the following structure: xL1 - distance of M1 to L1 in units of the binary displacement critpot - value of the critical potential at L1 rcrit - radius of M1 in z-direction in units of the binary displacement
See also: Roche_potential
Roche_equipotential_volume
Synopsis
computes the volume inside an Roche equipotential surface
Usage
Double_Type Roche_equipotential_volume
Qualifiers
- phi [= critical potential]: the surface's Roche potential
- N [= 100]: number of subdivisions in x and y for numerical integration
Roche_lobe_dimension
Synopsis
determines the size of the Roche lobe
Usage
Struct_Type Roche_lobe_dimension(Double_Type q)
Description
The structure has the following fields:
potential
: the critical Roche potential of the Roche lobe
{x
/y
/z
}{min
/max
}: the coordinates of the envelopping box
Roche_lobe_surface
Synopsis
calculates surface elements of a binary star
Usage
Struct_Type = Roche_lobe_surface(Double_Type mq, Double_Type fill, Integer_Type ntheta);
Qualifiers
noRoche - neglect the Roche potential -> spherical star noNorm - do not return the normal vectors noArea - do not return the areas poly - return the polygons defining the elements rotate - rotate the star around its (z-)axis by the given degrees (default: 0) eps - numerical precision (default: 1e-8) chatty - be chatty
Description
This function calculates the surface of the primary star in a binary, which is deformed by the Roche potential. Since this potential can not be solved analytically the surface must be divided into elements. The number of elements into theta direction can be specified by the 'ntheta' argument. The total number of surface elements taking into account the singularity at the poles and the symmetry in theta is here 2*ntheta * (ntheta-2) + 2 The Roche potential for a spherical orbit depends on the mass ratio mq = M_primary / M_secondary To determine the size of the star, the undeformed radius Rz in direction of the angular velocity vector (omega) can be used via the so-called fill factor: fill = Rz / R_crit where R_crit is the critical undeformed radius for a star, which fills its Roche lobe completely. In that case, the potential value at the surface is equal to that at the first Lagrange point, which might lead to mass transfer onto the companion (Roche lobe overflow).
The returned structure contains the following fields: Vector_Type[] r - position vector to each element starting at the stars center Vector_Type[] n - normal vector of each element Vector_Type[] A - area of each element Vector_Type[] p - if the 'poly' qualifier is given, it contains the polygons defining each surface element
NOTE: The origin of the reference frame is in the center of the star and the first Lagrange point is on the x-axis, while the rotation axis of the star is along the z-axis! If you want to rotate the star around its rotation axis (z-axis), use the 'rotate'-qualifier!
See also: Roche_critical, Roche_potential, radius_to_unit, Vector_Type
Roche_potential
Usage
Double_Type Roche_potential(Double_Type x[, y[, z]]);
Qualifiers
- q: mass ratio
- sph: coordinates are considered as spherical coordinates r, phi, theta
Description
The default values for y and z is 0.
The Roche potential has the following form:
-1/sqrt(x^2 + y^2 + z^2 ) - q/sqrt( (x-1.)^2 + y^2 + z^2 ) - (1+q)/2.*( (x-q/(1.+q))^2 + y^2 )
Roman
Synopsis
replaces an integer in the with capital Roman numeral strings
Usage
String_Type rom = Roman(Integer_Type);
See also: roman
SIprefix
Usage
Double_Type SIprefix(String_Type prefix)
SaturationVaporPressure
Synopsis
calculate the saturation vapor pressure of water in air
Usage
Double_Type[] SaturationWaterPressure(T;qualifiers)
Qualifiers
- centigrade: temperature argument is in centigrade (default: K)
- kPa : Return saturation vapor pressure in kilopascals
- hPa : Return saturation vapor pressure in hectopascals (=mbar)
- water : Calculate saturation vapor pressure over water
- ice : Calculate saturation vapor pressure over ice
- iapws : Use the more precise IAPWS prescription (see below)
Description
This function calculates the saturation vapor pressure (svp) of water over ice and water for a given temperature (default: K, but see the centigrade qualifier), the routine returns the svp in Pa unless the kPa or hPa keywords are given. By default, for T>0C (273.15K) the svp over water is returned, for smaller temperatures that for ice. Use the "water" and "ice" qualifiers if this is not what you want.
As a default, for the svp over water the equation of Davis (1992), Metrologia 29, 67-70 is used, while for ice the prescription given by Marti and Mauersberger (1993), Geophys. Res. Lett. 20, 363-366 (1993) is applied. For the default settings, the relative deviation between these equations and the IAPWS recommended procedure is very small. It does not exceed 0.03% in the 0-100C range, and is <1.5% in the -50-0C range and less than 2.5% between -100 and -50C, while the absolute difference in the range below 50C does not exceed 1 Pa.
If the iapws qualifier is given, the routine uses the recommendations of the International Association for the Properties of Water and Steam (IAPWS), which was originally given by Peter H. Huang, "New equations for water vapor pressure in the temperature range -100 deg C to 100 deg C for use with the 1997 NIST/ASME steam tables," in: Papers and abstracts from the third international symposium on humidity and moisture, Vol. 1, p. 69-76, National Physical Laboratory, Teddington, Middlesex, UK, April 1998. The algorithm used here is as described by http://emtoolbox.nist.gov/Wavelength/Documentation.asp The uncertainty is 20kPa at 100 C, less than 2 kPa at 45 C and 0.7 kPa at 20C.
Note that formally the calculation is only valid in the range from -100 to +100 centigrade. The function returns "NaN" for T>100C and extrapolates formalism smoothly to lower temperatures (the values are very small in this regime anyway).
This function is array safe.
TeX_value_pm_error
Synopsis
formats a (value, min, max)
triple as TeX string
Usage
String_Type TeX_value_pm_error(Double_Type value, Double_Type min, Double_Type max)
or
```c
String_Type TeX_value_pm_error(Struct_Type output_of__round_conf)
##### Qualifiers
* factor [=<code>1</code>]: factor to multiply <code>value, min, max</code>, will not be shown,
works only if the function is called with three arguments
* neg: factorize overall minus sign
* sci [=<code>2</code>]: exponent to switch to scienficic notation
* scismall [=<code>-sci</code>]: exponent for small numbers
* scilarge [=<code>+sci</code>]: exponent for large numbers
##### Description
<code>TeX_value_pm_error</code> tries to produce a reasonable TeX output
from a <code>(value, min, max)</code> tripel or the output of the function
<code>round_conf</code>.
Rounding and obtaining the number of digits is done with the
function <code>round_conf</code>.
__See also__: round_conf,round_err
#### UTC2UNIXtime
##### Synopsis
converts Universal Time (GMT) to seconds since 1970-01-01, 0:00 UTC
##### Usage
```c
Long_Type secs = UTC2UNIXtime(Struct_Type tm);
or
Long_Type secs = UTC2UNIXtime(Integer_Type Y[, m[, d[, H[, M[, S]]]]]);
Description
UTC2UNIXtime
can be seen as the inverse of the gmtime
function,
operating on Coordinated Universal Time instead of dates in a local time zone,
as mktime
does, which is the inverse of localtime
.
UTC2UNIXtime
can also take the components of the date Y
-m
-d
, H
:M
:S
UTC
as separate arguments. The year Y
is then given directly (unlike tm.tm_year = Y - 1900
),
and likewise the month m
(unlike tm.tm_mon = m-1
).
For omitted arguments, m=1
, d=1
, H=0
, M=0
, and S=0
are assumed.
See also: mktime, gmtime, localtime
Voigt_complex
Synopsis
Compute complex Voigt profile
Usage
Complex_Type[] = Voigt_complex(z, z0, sigma, gamma);
WorldCoastLine
Usage
Struct_Type WCL[] = WorldCoastLine();
Description
WCL
is an array of structures with x
and y
fields
which specify the coordinates (longitude east and lattitude north)
in degrees of the World's coast lines.
Example
xrange(-2, 2); yrange(-1, 1); variable cl; foreach cl (WorldCoastLine()) oplot( Hammer_projection(cl.x, cl.y; deg, normalized) );
X11color2rgb
Synopsis
return rgb value for x11's rgb.txt
Usage
variable rgb=X11color2rgb(colorname)
Qualifiers
- rgbfileFile defining the colors, default: /usr/share/X11/rgb.txt
Description
This function returns the integer specifying the rgb code of a color named in the file rgb.txt accompanying the X11 distribution.
See also: xfig_lookup_w3c_color
__push_array
Synopsis
pushes the values of an array onto the stack
Usage
(Any_Type v1, v2, ...) = __push_array(Any_Type a[]);
Description
NOTE that the Stack is limted to 2499 Elements, if the length of a exceeds this value a Stack_Overflow Error is thrown!
_maximum
Synopsis
computes the maximum of >=2 values
Usage
Array_Type _maximum(Array_Type x1, x2, ...)
Description
x1
, x2
, ... have to be arrays of the same length, or scalars.
The _maximum
function returns an array of the maximum values
by successively calling the _max
function.
If all x1
, x2
, ... are scalars, _maximum
returns a scalar value,
in this case, _maximum(x1, x2, ...)
is equivalent to max([x1, x2, ...])
.
See also: _max, max, _minimum, _min, min
_minimum
Synopsis
computes the minimum of >=2 values
Usage
Array_Type _minimum(Array_Type x1, x2, ...)
Description
x1
, x2
, ... have to be arrays of the same length, or scalars.
The _minimum
function returns an array of the minimum values
by successively calling the _min
function.
If all x1
, x2
, ... are scalars, _minimum
returns a scalar value, too,
in this case, _minimum(x1, x2, ...)
is equivalent to min([x1, x2, ...])
.
See also: _min, min, _maximum, _max, max
air_to_vacuum
Synopsis
Convert air wavelengths to vacuum wavelengths
Usage
Double_Type l_v[] = air_to_vacuum(Double_Type l_a[])
Description
Convert air wavelengths l_a (in Angstroem) to vacuum wavelengths l_v (in Angstroem) for dry air at 15 degrees Celsius, 101.325 Pa, and with 450 parts per million CO2 content. Valid in the wavelength range between 2300 to 17000 Angstroem.
Notes
According to Equation 1 in Ciddor 1996, Applied Optics, 35, 1566: 1e8*(l_v-l_a)/l_a = n-1 = k1/(k0-l_v^(-2)) + k3/(k2-l_v^(-2)) => l_a = l_v / ( k1*1e-8/(k0-l_v^(-2)) + k3*1e-8/(k2-l_v^(-2)) + 1 ) with k0 = 238.0185 * 1e-8 / Angstroem^2, k1 = 5792105 * 1e-8 / Angstroem^2, k2 = 57.362 * 1e-8 / Angstroem^2, k3 = 167917 * 1e-8 / Angstroem^2. This function is not analytically solvable for l_v. To zeroth order, however, l_v ~ l_a so that l_v ~ l_a * ( k1*1e-8/(k0-1/l_a^2) + k3*1e-8/(k2-1/l_a^2) + 1 ).
Qualifiers
- verbose [=1]: Set to 0 to suppress warnings.
Example
l_v = [2000,4000,8000,16000]; print(l_v - air_to_vacuum(vacuum_to_air(l_v)));
See also: vacuum_to_air
alpha_ff
Synopsis
Calculate the absorption coefficient for free-free absorption
Usage
Double_Type alpha = alpha_ff(nu,T);
Qualifiers
- Z: average nuclear charge (default=1)
- ne: electron particle density (cm^-3; default: 1e10)
- np: proton particle density (cm^-3; default: 1e10)
Description
This function returns the absorption coefficient for free-free radiation (bremsstrahlung). The frequency, nu, is measured in Hz and can be an array, the temperature T is measured in K.
For the moment this function assumes the Gauntt factor equals unity.
Note: per Kirchhoff's law the total bremsstrahlung spectrum from a slab of size R is B_nu(nu,T)*(1.-exp(-R*alpha_ff(nu,T)))
See also: j_ff
angle2string
Synopsis
convert an angle to a string for pretty printing
Usage
string=angle2string(angle;qualifiers)
Qualifiers
- hoursign: angle is a hour angle, always include sign
- declination: angle is a declination, always include sign
- latitude: angle is a latitude, always include sign
- hours: display hours, not degrees
- deg: input angle is in degrees (default: radians)
- separator: string to insert between the degrees/hours, minutes, and after the seconds. If a scalar: insert between degrees and minutes only [e.g., ":"]. If an array: insert the three array elements between the output numbers (e.g., ["h","m","s"]). Default is ["d","m","s"] and ["h","m","s"] depending on whether the hours-qualifier is set or not.
- blank: do not display leading zeros
- secfmt: sprintf format for the seconds (default: %04.1f, i.e., a precision of 0.1 sec)
- mas: display to a precision of milliseconds corresponds to secfmt="%06.3f"
- muas: display to a precision of microseconds corresponds to secfmt="%08.5f"
Description
This function is used to produce well formatted strings out of angular quantities which are given in radians or degrees. (the default is radian).
Angles can be displayed either in degrees or in hours, various separators can be used to separate the hms-terms, and the routine can be told to always display a sign (e.g., for declination- or latitude-like quantities), or not.
This routine is array safe.
See also: hms2deg,dms2deg,generate_iauname
angle_to_rad
Synopsis
converts an angle in degrees or h:m:s format into radian
Usage
Double_Type angle_to_rad(Double_Type x);
Qualifiers
- unit [
="deg"
]: unit ofx
("deg"
/"hms"
/"rad"
)
Description
x
can be a scalar value or (unless unit="rad"
) an array of the form
-
[deg, arcmin]
(respectively[hour, min]
forunit="hms"
) -
[deg, arcmin, arcsec]
(respectively[hour, min, sec]
forunit="hms"
) -
[sign, deg, arcmin, arcsec]
(respectively[sign, hour, min, sec]
).
The latter representation is needed for -1 < deg
/hour < 0
.
See also: hms2deg,dms2deg
angular_separation
Usage
sep=angular_separation(ra1,dec1,ra2,dec2);
Synopsis
calculates the angular distance between two points on a sphere
Qualifiers
- deg: if set, the input coordinates and output are in degrees (default: radian)
- radian: if set, the input coordinates and output are in radian (the default))
Description
This routine calculates the angular separation between points on the sky with coordinates ra1/dec1 and ra2/dec2.
The function is equivalent to greatcircle_distance but is compatible in terms of its qualifiers with the other coordinate system routines. It also returns the angular separation in the units of the input parameters rather than always in rad.
This function is array safe (for either ra1/dec1 or ra2/dec2).
See also: hms2deg,dms2deg,angle2string,greatcircle_distance
ansi_escape_code
Synopsis
transforms a text format, e.g., a color into its corresponding ANSI code
Usage
String_Type ansi_escape_code(String_Type format[, String_Type text]);
Qualifiers
list - print a list of available keywords
Description
Transforms the given (human readable) format string into a sequence of ANSI escape codes. The format string consists of keywords, e.g., colors separated by semicolons. A list of supported keywords is printed setting the 'list' qualifiers.
The cursor movement keywords 'up', 'down', 'forwards', and 'back' allow an optional preceding number, which specifies the amount of the movement.
In order to reset a previous format you may use the 'reset' keyword. In case an optional text is provided, this text and the reset ANSI code is appended to the returned string automatically.
Example
ansi_escape_code("red;blink", "ALERT"); ansi_escape_code("3up"); % move three lines up ansi_escape_code("blue"); % blue color *from now on* ansi_escape_code("reset"); % turn off the format again
append_struct_arrays
Synopsis
appends a structure's fields to another structures's fields
Usage
append_struct_arrays(Struct_Type &s, Struct_Type additional_s[]);
Description
s
and additional_s
have to be structures with the same fields.
s
, which is changed by this function, has to be passed by reference.
For every field
,
s.field = [s.field, additional_s.field];
The following two statements are equivalent:
append_struct_arrays(&s, additional_s);
s = merge_struct_arrays( [s, additional_s] );
See also: merge_struct_arrays
aproposer
Synopsis
recall object names and the documentation satisfying a regular expression
Usage
aproposer("s")
Description
This is an extended version of S-Lang's 'apropos' function, which may be used to get a list of all defined objects whose name matches the regular expression "s". In addition, the SYNOPSIS of all functions described in the documentation files is checked on matches as well. Finally, the output is formatted such that the matching substring and the object's name are printed in different colors (see below for format options).
In order to use this function instead of the intrinsic 'apropos' function, you can put alias("aproposer", "apropos"); into, e.g., your .isisrc file.
Further variables can be defined within the .isisrc file for
more control options:
APROPOS_ENABLE_SYNOPSIS - 0 = turn off the synopsis search
APROPOS_FORMAT_NAME - format string for the object's name
APROPOS_FORMAT_MATCH - format string for the substring match
APROPOS_LINE_WIDTH - number of columns of the terminal,
which is needed to format the output (default: output of
stty size
or 80 if the first attempt fails)
Read the documentation of 'ansi_escape_code' for details about
the format string.
See also: apropos, .apropos, help, who, get_doc_files, ansi_escape_code
array2matrix
Synopsis
Transforms an Array_Type into a matrix
Usage
Any_Type[l,...] = array2matrix(Array_Type[l] array);
Qualifiers
- filler[=-_NaN]:Filler for empty array entries
Description
Similar to array_flatten, this function flattens the array into a single one. In contrast to array_flatten, the array_shape of the individual entries is conserved. If the array_shape of the individual array entries is not the same, the maximal dimension is determined and empty entries are being filled with the value given with the filler qualifier.
Example
variable arr = Array_Type[3]; arr[0] = _reshape( [1:3*4*5], [3,4,5] ); arr[1] = _reshape( [1:3*4], [3,4] ); arr[2] = _reshape( [1:2*3*4], [2,3,4] );
variable mat = array2matrix(arr); print(mat);
See also: array_flatten
array_copy
Synopsis
makes a copy of a nested array
Usage
Array_Type copy = array_copy(Array_Type[] array);
Description
array_copy
copies an Array_Type[]
with all its
entries, which can be of any Type! If an entry is a Struct_Type
array_copy
calls struct_copy
and if an entrie is
an Array_Type[]
, array_copy
calls itself.
Example
See also: struct_copy
array_density
Synopsis
calculates the 'intensity' of an array, i.e. the invers of the spacing
Usage
array_density(bin_lo, bin_hi);
array_flatten
Synopsis
flattens an array of Array_Type[] into a single array
Usage
Any_Type array_flatten(Array_Type[] array);
Description
The given 'array' itself contains arrays of any data type. This function flattens the 'array' such that all internal arrays are merged.
Also the array can have entries of any data type, they must be of the same type! Otherwise the function will throw an error. To avoid this split the array, e.g., with
arr1 = arr[where(_typeof(arr)==Double_Type)];
or if there are NULL entries,
arr1 = arr[wherenot( arr == NULL )];
and use arr1 as argument for array_flatten!
Example
variable arr = Array_Type[3]; arr[0] = [1,2,3,4]; arr[1] = [5,6]; arr[2] = [7,8,9];
variable arrF = array_flatten(arr); % returns an Integer_Type[9] containing % [1,2,3,4,5,6,7,8,9]
array_insert
Synopsis
Inserts an array into another array
Usage
Any_Type[] = array_insert(Any_Type[] A, Any_Type[] IA, Integer_Type i);
Description
Inserts the array or single value IA into the array A at the index position i. A and IA have to be of the same DataType_Type! Basically it is:
Example
variable A = [1:10]; variable IA = [2,1]; variable NA = array_insert( A, IA, 3 ); print(NA);
array_map_index
Synopsis
Call a function on a subset of arrays
Usage
Array_Type array_map_index (Int_Type[] start, Int_Type[] stop,
DataType_Type type, Ref_Type func, args, ...);
or
```c
Array_Type array_map_index (Int_Type[] index, DataType_Type type,
Ref_Type func, args, ...);
or
Array_Type array_map_index (Ref_Type start, Ref_Type stop,
DataType_Type type, Ref_Type func, args, ...);
or
Array_Type array_map_index (Ref_Type index, DataType_Type type,
Ref_Type func, args, ...);
##### Description
The <code>array_map_index</code> function is somewhat the opposite of the
<code>array_map</code>. As <code>array_map</code> allows it to apply a function to every
array element individually, <code>array_map_index</code> can be used to apply a
function (taking an array) on one or more subsamples of the target array.
Case 1: If called with <code>start</code> and <code>stop</code> the function <code>func</code> is
applied to <code>args[start[i]:stop[i]]</code> for each <code>i</code>.
Case 2: If called with only one index array <code>index</code> the function is
applied to args[index]. This is the same as <code>func(args[index]);</code>.
IMPORTANT: array_map_index applys <code>func</code> to an array!
Case 3: If called with <code>&start</code> and <code>&stop</code> the function behaves the
same as for case 1, but this is usefull when you want to nest <code>array_map</code>
calls.
Case 4: If called with <code>index</code> the it behaves the same as for case 2 but,
again, this might be usefull for <code>array_map</code> nesting.
NOTE: This function can hide a lot of important operations so it should be
avoided for implementing algorithms. However, it is very useful for fast
scripting and interactive (shell) operations.
##### Example
% Case 1: apply mean() to a subset of data
data = rand_gauss(0.1, 100) % generate 100 normal distributed random numbers
means = array_map_index([0,5,25], [99,94,74], &mean, data);
print(means);
% this outputs the same as
print(mean(data)); print(mean(data[[5:94]])); print(mean(data[[25:74]]));
% Case 2: apply mean() to specified subset
data = rand_gauss(0.1, 100) % 100 random numbers
idx = where(data < 0) % select only negative data
neg_mean = array_map_index(idx, Double_Type, &mean, data);
print(neg_mean);
% this results in the same as
print(mean(data[n]));
% Case 3: apply mean() to array of arrays
data = Array_Type[3];
data[0] = rand_gauss(0.1, 100);
data[1] = rand_gauss(0.01, 100);
data[2] = rand_gauss(1, 100);
start = [0,5,25];
stop = [99,94,74];
means = array_map(Array_Type, &array_map_index, &start, &stop,
Double_Type, &mean, data);
% returns the mean as in example for case 1 but applied to all sub arrays
% note that you have to declare the start and stop arrays explicitly because
% &[5:94] is not a valid expression
% Case 4: Same as case 3 but with only one index array
% Use array_map_index to call a function n times
rn = array_map_index([1:5], Array_Type, &rand_gauss, 0.1, 100);
% gives 5 arrays of 100 random numbers each
__See also__: array_map, array_map_qualifiers
#### array_map_qualifiers
##### Synopsis
Apply a function to each element of an array (also if given as qualifier)
##### Usage
```c
Array_Type array_map (DataType_Type type, Ref_Type func, args, ... );
or
```c
(Array_Type, ...) array_map (DataType_Type type, ..., Ref_Type func, args, ... );
##### Description
The 'array_map_qualifiers' function is an extension of the 'array_map'
function and supports all its features.
This function also works for qualifiers! First of all qualifiers given to
the 'array_map_qualifiers' are passed onto the called function 'func'.
Furthermore if the qualifier(s) itself is an array, it will be mapped too!
This is especially helpfull when used for calls of integration function
such as 'qromb' as these require parameters to be given as qualifiers!
This functionality can be turned off with the "noqualmap"
qualifier, which is necessary for, e.g., epochfolding routines (and much
faster!)
NOTE that this function works as 'array_map', i.e., the mapping will apply
for the first argument which is an array! As the qualifiers are passed after
the arguments, an argument array will supersede a qualifier array as long as
they do not have the same length.
##### Example
define t(x){ return qualifier("a",0)\*x;};
%% Mapping function using an argument array:
print(array_map_qualifiers( Double_Type, &t, [1,2,3,4] ; a=1 ));
%% Mapping function using a qualifier array:
print(array_map_qualifiers( Double_Type, &t, 1 ; a=[1,2,3,4] ));
%% Mapping function argument array superseding qualifier array:
print(array_map_qualifiers( Double_Type, &t, [3,4,5,6] ; a=[1,2] ));
%% Mapping function using argument array AND qualifier array:
print(array_map_qualifiers( Double_Type, &t, [3,4,5,6] ; a=[1,2,3,4] ));
%% Map function using a qualifier array application on 'qromb' integrater
print(array_map_qualifiers( Double_Type, &qromb, &t, 0, 1 ; a=[0:2] ));
__See also__: array_map, struct_of_arrays_2_struct_array
#### array_permutation_matrix
##### Synopsis
returns all possible permutations of the given arrays
##### Usage
```c
Array_Type array_permutation_matrix(Array_Type[] or List_Type arrays);
Description
Calculates the index arrays of all possible permutations of the elements in the given arrays (see the below example). The input has to be either an array of arrays (Array_Type) or a list of arrays (List_Type). The returned matrix (J,K) is an array of integer arrays. The element (J,K) is the index of the K^th input array, which corresponds to the J^th permutation.
Example
variable arr = Array_Type[3]; arr[0] = [10,20,30,40]; arr[1] = [50,60]; arr[2] = [70,80,90];
variable matrix = array_permutation_matrix(arr); % Returns an Array_Type[24] corresponding to 4*2*3 = 24 possible % permutations. Each element is an Integer_Type[3] containing the % indices for each of the 3 input arrays: % matrix[0] = [0,0,0] -> [10,50,70] % matrix[1] = [1,0,0] -> [20,50,70] % matrix[2] = [2,0,0] -> [30,50,70] % matrix[3] = [3,0,0] -> [40,50,70] % matrix[4] = [0,1,0] -> [10,60,70] % matrix[5] = [1,1,0] -> [20,60,70] % ... % matrix[23] = [3,1,2] -> [40,60,90]
array_permute
Synopsis
generates a random permutation of [0 : n-1]
Usage
Integer_Type[] random_array(Array_Type a)
or
Integer_Type[] random_array(Integer_Type n)
Description
The return value is an array of indices.
To randomize an array a
, use a[array_permute(a)]
.
See also: array_sort
array_remove
Synopsis
removes one element from an array
Usage
Array_Type a_ = array_remove(Array_Type a, Integer_Type i);
Description
a_
contains all elements of a
excecpt of the one indexed by i
.
array_unique
Synopsis
takes elements of an array only once
Usage
Array_Type array_unique(Array_Type a)
Description
Duplicated array elements are removed from the array. This function is basically a[unique(a)] where unique gives the indicies of the unique array elements.
See also: unique
array_zip
Synopsis
"zip" two or more arrays together
Usage
Any_Type array_zip(Any_Type[] a, b, ...);
Description
Returns an array of tuples (or higher orders), where each tuple i contains the i-th element from each of the given arrays. The returned array is truncated to the length of the shortest array. Thus, the number of tuples is equal to this smallest length.
Example
a = [1:10]; b = [11:20]; c = array_zip(a,b); % c = [1,11,2,12,3,13,4,14,...,10,20];
assoc_struct_combine
Synopsis
associate structure fields, combine structures and fill entries
Usage
Struct_Type <code>comb</code> = assoc_struct_combine(Struct_Type <code>s1</code>, Struct_Type <code>s2</code>, String_Type <code>f1</code>, String_Type <code>f2</code>);
Qualifiers
- double_fill [=_NaN]: filler for missing entries in Double_Type fields of
s2
- float_fill [=_NaN]: filler for missing entries in Float_Type fields of
s2
- integer_fill [=_0]: filler for missing entries in Integer_Type fields of
s2
- string_fill [=""]: filler for missing entries in String_Type fields of
s2
- flag_field_name [="both_struct_flag"]: additional field name for
comb
Description
This function combines two structures. It uses the field with name f1
of the
structure s1
and that with name f2
of s2
to associate the entries with the function
get_intersection. The fields of the returned structure comb
have the length of the
field f1
. The entries in these fields are filled according to the qualifiers, if
there is no association found for this entry.
Fields of a type for which no fill qualifier is given are filled with NULL and
returned as a list. Fillers can be defined via qualifiers for all variable types.
The suffix _fill is required, i.e., for example the qualifier to define a filler
for Complex_Type is complex_fill
.
WARNING: This function fails for more dimensional fields. If there are fields with
the same name in both structures, the field of s2
is used, as it is done by
struct_combine. The field f2
is not included.
The function adds a field (name can be indentified via qualifier flag_field_name
),
which indicates if the entry was found in both structures (1) or if it exists
only in f1
(0).
Examples
a = struct{source = ["1", "2", "3"], flux=[0.1, 2.5, 0.7]}; b = struct{source = ["3", "2"], redshift=[0.65, 0.03], opt_id = ["Quasar","BL Lac"]}; c = assoc_struct_combine (a,b,"source", "source"); print_struct(c);
% without filling fields, but using just common fields this corresponds to: (i1,i2) = get_intersection(a.source,b.source); struct_filter(a,i1); struct_filter(b,i2); d = struct_combine(a,b); print_struct(d);
See also: get_intersection, struct_combine, struct_filter
associate_arrays
Synopsis
computes the smallest pairwise difference between elements of two arrays and returns the indices
Usage
Struct_Type result = associate_arrays(Double_Type[n1] A, Double_Type[n2] B);
Description
For each element A[i]
this function searches the element B[j]
with the smallest difference to A[i]
. It returns the indices j
as result.index[i]
and the corresponding difference as result.diff[i]
.
Qualifiers
- nr_closest [=1]: set to N, to obtain the N closest elements.
In this case additionally the index of the second, ..., N-th closest element is
returned as
result.index2[*,0], ..., result.index2[*,N-2]
, the same holds forresult.diff2
.
See also: array_sort
associate_coords
Synopsis
computes the smallest pairwise angular separation between arrays of positions and returns the indices
Usage
Struct_Type result = associate_coords(Double_Type[n1] RA1, Double_Type[n1] DEC1),Double_Type[n2] RA2, Double_Type[n2] DEC2);
Qualifiers
- nr_closest [=1]: set to N, to obtain the N closest elements.
In this case additionally the index of the second, ..., N-th closest element is
returned as
result.index2[*,0], ..., result.index2[*,N-2]
, the same holds forresult.sep2
. - unit [
="deg"
]: unit of the angular coordinates - alpha1_unit [
="deg"
]: unit of the RA1 - alpha2_unit [
="deg"
]: unit of the RA2 - delta1_unit [
="deg"
]: unit of the DEC1 - delta2_unit [
="deg"
]: unit of the DEC2
Description
For each position (RA1[i], DEC1[i])
this function searches the element (RA2[j],DEC2[j])
with the smallest difference to (RA1[i],DEC[i])
. It returns the indices j
as result.index[i]
and the corresponding separation in rad as result.sep[i]
.
The units of each coordinate can be set with the qualifier unit
, but also independently,
as for the function greatcircle_distance
.
See also: greatcircle_distance, associate_arrays
bezNp
Synopsis
calulates a Bezier curve between vectors
Usage
Vector_Type bez = bezNp(Vector_Type vectors);
Description
The Bezier curve between each pair of the arrays in 'vectors' will be calculated. The 'vectors' arrays needs to be at least two entries. Returned value 'bez' will be one entry shorter than 'vectors'. Main working horse for the bezier function.
See also: bezier
bezier
Synopsis
calulates a Bezier curve from given supporting points
Usage
Vector_Type bez = bezier(List_Type points, Integer_Type nbins);
Description
Bezier calculates the Bezier curve between the first and last [x,y] entry of 'points', using the other entries of 'points' as supporting points. 'points' must be a List_Type, with each element containing the [x,y] values of the corresponding point. 'nbins' gives the numbers of points used to calculate each part of the curve, as well as the length of the returned curve.
Qualifiers
- allbez : if set, all intermediate Bezier curves are returned as a List, containing all the Vector_Type Bezier curves.
- conline : if set, the initial lines connecting the points are also returned, i.e., the return value becomes (Vector_Type conlines, Vector_Type bez).
See also: xfig_plot_allbezier, bezNp
binarydigits
Synopsis
retrieves a number's dual representaion
Usage
String_Type binarydigits(Integer_Type x)
Qualifiers
- n [
=16
]: number of bits
blackbody_to_rgb
Synopsis
computes RGB colors of hot blackbody radiation.
Usage
Integer_Type blackbody_to_rgb(Double_Type T)
Description
T
is the temperature of the blackbody in Kelvin.
The return value is its 24 bit RGB value.
It is computed after the code "RGB VALUES FOR HOT OBJECTS" by Dan Bruton (astro@sfasu.edu), which can be found at http://www.physics.sfasu.edu/astro/color.html.
blocks_between_gaps
Synopsis
retrieves blocks in an array whose elements may contain gaps
Usage
Integer_Type blocks[] = blocks_between_gaps(Double_Type a[], Double_Type gap);
Description
The array a has to contain monotonically increasing values. The return value will be an integer array of the same length, whose elements start with 0 and increase by 1 whenever the difference of sequential elements of the array a is larger than gap.
See also: split_struct, split_lc_at_gaps
cart2cyl
Synopsis
Convert Cartesian to cylindrical coordinates
Usage
cart2cyl(Double_Types x[], y[], z[], vx[], vy[], vz[]);
Description
Convert Cartesian (x,y,z,vx,vy,vz) to cylindrical (r,phi,z,vr,vphi,vz) coordinates with phi rotating counter-clockwise from the positive x-axis.
Qualifiers
Example
(r,phi,z,vr,vphi,vz) = cart2cyl(1,1,1,2,3,4); (r,phi,z,vr,vphi,vz) = cart2cyl([1,0],[1,0],[1,0],[2,0],[3,0],[4,0]); (r,phi,z,vr,vphi,vz) = cart2cyl( cyl2cart(1,PI/2,1,2,3,4) );
See also: cart2sphere, cyl2cart
cart2sphere
Synopsis
Convert Cartesian to spherical coordinates
Usage
cart2sphere(Double_Types x[], y[], z[], vx[], vy[], vz[]);
Description
Convert Cartesian (x,y,z,vx,vy,vz) to spherical (r,phi,theta,vr,vphi,vtheta) coordinates with polar angle phi rotating counter-clockwise from the positive x-axis and azimuth angle theta being zero for the positive z-axis.
Qualifiers
Example
(r,phi,theta,vr,vphi,vtheta) = cart2sphere(0,0,0,1,2,3); (r,phi,theta,vr,vphi,vtheta) = cart2sphere(0,0,1,1,2,3); (r,phi,theta,vr,vphi,vtheta) = cart2sphere(1,1,1,3,3,3); (r,phi,theta,vr,vphi,vtheta) = cart2sphere(1,0,0,0,2,1); (r,phi,theta,vr,vphi,vtheta) = cart2sphere([0,1],[0,0],[0,0],[1,0],[2,2],[3,1]);
See also: cart2cyl
cel2gal
Synopsis
Transform celestial to cartesian Galactic coordinates
Usage
cel2gal(Double_Types ah[], am[], as[], dd[], dm[], ds[], dist[], vrad[], pma_cos_d[], pmd[]; qualifiers)
Description
Transform celestial coordinates (right ascension [h, m, s], declination [deg, arcmin, arcsec], distance [kpc], radial velocity [km/s], proper motion in right ascension times cosine of declination [mas/yr], proper motion in declination [mas/yr]) to right-handed, cartesian Galactic coordinates (x [kpc], y [kpc], z [kpc], vx [km/s], vy [km/s], vz [km/s]) with the Galactic center at the origin, the Sun on the negative x-axis and the z-axis pointing to the north Galactic pole implying clockwise Galactic rotation when seen from the half space with positive z. See function 'RD2rad' for detailed information on how to properly enter right ascension and declination.
Qualifiers
- parallax: Set this qualifier to use parallaxes [mas] instead of distances [kpc].
- SunGCDist: Sun-Galactic center distance [kpc]; default: 8.4 (see Model I in Irrgang et al., 2013, A&A, 549, A137)
- vxs: Sun's x-velocity component [km/s] relative to the local standard of rest; default: 11.1 (Schoenrich, Binney & Dehnen, 2010: MNRAS 403, 1829)
- vys: Sun's y-velocity component [km/s] relative to the local standard of rest; default: 12.24 (Schoenrich, Binney & Dehnen, 2010: MNRAS 403, 1829)
- vzs: Sun's z-velocity component [km/s] relative to the local standard of rest; default: 7.25 (Schoenrich, Binney & Dehnen, 2010: MNRAS 403, 1829)
- vlsr: Local standard of rest velocity [km/s]; default: 242 (see Model I in Irrgang et al., 2013, A&A, 549, A137)
- GC_NGP: Celestial coordinates of the Galactic center and the north Galactic pole in the format [ah_GC, am_GC, as_GC, dd_GC, dm_GC, ds_GC, ah_NGP, am_NGP, as_NGP, dd_NGP, dm_NGP, ds_NGP]; default: [17, 45, 37.224, -28, 56, 10.23, 12, 51, 26.282, 27, 07, 42.01] (J2000, Reid & Brunthaler, 2004, ApJ, 616, 872)
Example
(x, y, z, vx, vy, vz) = cel2gal(0, 0, 0, 0, 0, 0, 0, 0, 0, 0); % -> Sun's position and velocity (x, y, z, vx, vy, vz) = cel2gal(17, 45, 37.224, -28, 56, 10.23, 8.4, -11.1, -3.171, -5.544); % -> Galactic center's position and velocity (x, y, z, vx, vy, vz) = cel2gal([0,17], [0,45], [0,37.224], [0,-28], [0,56], [0,10.23], [0,8.4], [0,-11.1], [0,-3.171], [0,-5.544]); % -> position and velocity of Sun and Galactic center stored in the same array (x, y, z, vx, vy, vz) = cel2gal(17, 45, 37.224, -28, 56, 10.23, 8.4, -11.1, -3.171, -5.544; GC_NGP=[17, 45, 37.1991, -28, 56, 10.221, 12, 51, 26.2755, 27, 7, 41.704]); % -> use non-standard coordinates for Galactic center and north Galactic pole (x, y, z, vx, vy, vz) = cel2gal( gal2cel(0, 0, 0, 0, 0, 0) ); % compute vr and vphi: vr = (x*vx+y*vy)/sqrt(x^2+y^2); vphi = (y*vx-x*vy)/sqrt(x^2+y^2);
See also: RD2rad, gal2cel
cerf
Synopsis
Complex error function
Usage
Complex_Type[] = cerf(Complex_Type[]);
Description
Compute complex error function. Only useful in a region with abs(z)<10.
See also: Faddeeva, cerfc
cerfc
Synopsis
Complex error function complement
Usage
Complex_Type[] = cerfc(Complex_Type[]);
Description
Compute complex error function complement. Only useful in a region with abs(z)<10.
See also: Faddeeva, cerf
cholesky_decomposition
Synopsis
Decompose a matrix into the product of a lower triangular matrix and its transpose
Usage
Double_Type cd[n,n] = cholesky_decomposition(Double_Type m[n,n])
Description
The Cholesky decomposition is a decomposition of a symmetric, positive-definite matrix into the (matrix) product of a lower triangular matrix and its transpose: m = cd # transpose(cd). The Cholesky–Banachiewicz algorithm is used here to carry out the decomposition.
Notes
The Cholesky decomposition can be used to generate correlated variables that obey a given covariance matrix. See the example section for details.
Example
m = Double_Type[3,3]; m[0,*] = [ 9,-3,-6]; m[1,*] = [-3,10, 5]; m[2,*] = [-6, 5, 6]; print(m); % symmetric, positive-definite matrix cd = cholesky_decomposition(m); print(cd); % lower triangular matrix print(cd#transpose(cd)); % matrix product of 'cd' and its transpose
% From three uncorrelated normal variables 'g', generate three correlated normal % variables 'x' that obey the covariance matrix 'm':
g = Double_Type[3,100000]; g[0,*] = grand(100000); g[1,*] = grand(100000); g[2,*] = grand(100000); (cov_mat, cor_mat) = covariance_correlation_matrix(g); print(cov_mat); % the covariance matrix of 'g' shows that there are no correlations
x = cd#g; (cov_mat, cor_mat) = covariance_correlation_matrix(x); print(cov_mat); % the covariance matrix of 'x' is indeed (almost) 'm'
See also: covariance_correlation_matrix
complementary_array
Synopsis
returns the "difference" of arrays
Usage
Array_Type C = complementary_array(Array_Type B, A)
Description
C = B
\ A
complex2rgb
Synopsis
converts a complex number to a RGB color
Usage
Integer_Type complex2rgb(Complex_Type z);
Description
arg(code{z}) determines the hue,
|z
| determines the lightness:
0<=|z
|<=1 is mapped from black to color,
1<=|z
|<inf is mapped from color to white.
The return value is a 24-bit RGB value.
Qualifiers
- sat[=1]: saturation (from 0 to 1) of the color
See also: hsl2rgb, png_write
continued_fraction
Synopsis
computes a continued fraction
Usage
Double_Type continued_fraction(UInteger_Type a[])
Description
1 | 1 |
continued_fraction(a) = a[0] + ------ + ------ + ...
| a[1] | a[2]
continued_fraction_expansion
Synopsis
expands a floating point number as a continued fraction
Usage
UInteger_Type[] continued_fraction_expansion(Double_Type x)
Qualifiers
- verbose
- maxlen [
=64
]
convert_units
Synopsis
convers units of a physical quantity
Usage
convert_units(PhysicalQuantity_Type q);
or
PhysicalQuantity_Type convert_units(PhysicalQuantity_Type q; copy)
Qualifiers
- leng: unit of length
- time: unit of time
- mass: unit of mass
- curr: unit of electrical current
- temp: unit of temperature
- SI: sets leng="m", time="s", mass="kg", curr="A", temp="K"
- copy: does not change unit of q, but returns a copy of q
Description
The allowed units are specified by the associative arrays baseunits_length_in_m, baseunits_time_in_s, baseunits_mass_in_kg, baseunits_current_in_A, and baseunits_temperature_in_K which specify the factor to the corresponding SI units.
correlation_coefficient
Synopsis
calculates the [weighted] linear correlation coefficient between two arrays
Usage
Double_Type rho = correlation_coefficient(Double_Type x[], Double_Type y[] [ , Double_Type w[] ]);
Description
The function of the unweighted correlation coefficient reads:
n * sum(x*y) - sum(x) * sum(y)
rho = -------------------------------------------------------
sqrt(n*sum(x^2)-sum(x)^2) * sqrt(n*sum(y^2)-sum(y)^2)
w defines the weight of each point. By default it is set to w = 1.
See also: find_correlations
cosmo_param
Synopsis
Derive full set of cosmological density parameters from a partial set
Usage
Struct_Type P = cosmo_param();
Qualifiers
- omega_m: set value of the normalized matter density (Omega_M)
- omega_lambda: set value of the normalized cosmological constant (Omega_Lambda)
- omega_k: set value of the the normalized curvature term (Omega_k)
- q0: set value of the deceleration parameter (q0)
Description
This procedure is called by lumdist to allow the user a choice in defining any two of four cosmological density parameters. Given any two of the four input parameters this program will derive the remaining two. Omega_M - Normalized matter energy density, non-negative numeric scalar Omega_Lambda - Normalized cosmological constant, numeric scalar Omega_k - Normalized curvature parameter, numeric scalar. This is zero for a flat universe q0 - Deceleration parameter, numeric scalar = -R*(R'')/(R')^2 = 0.5*Omega_m - Omega_lambda Here "normalized" means divided by the closure density so that Omega_m + Omega_lambda + Omega_k = 1. For a more precise definition see Carroll, Press, & Turner (1992, ArAA, 30, 499).
If less than two parameters are defined, this procedure sets default values of Omega_k=0 (flat space), Omega_lambda = 0.7, Omega_m = 0.3 and q0 = -0.55
If more than two parameters are defined upon input (overspecification), then the first two defined parameters in the ordered list Omega_m, Omega_lambda, Omega_k, q0 are used to define the cosmology.
(This procedure is translated from the IDL COSMO_PARAM function.)
Example
variable default_set = cosmo_param(); print (default_set);
variable my_set = cosmo_param(; omega_k = 0.1, omega_lambda = 0.65); print (my_set);
See also: lumdist
cov_matrix
Synopsis
computes a covariance matrix
Usage
COV = cov_matrix(Struct_Type s);
or
COV = cov_matrix(Double_Type X[][]);
or
COV = cov_matrix(Double_Type x0[], x1[], ...);
Description
COV[i,j] =
cov(x_i, x_j) = < (x_i - <x_i>) * (x_j - <x_j>) >
covariance_correlation_matrix
Synopsis
Compute the covariance and correlation matrices of a set of variables
Usage
(Double_Types cov_mat[n,n], cor_mat[n,n]) = covariance_correlation_matrix(Double_Type a[n,N])
Description
This function estimates the covariance and correlation matrices of a set of random variables. Consider a sample 'a' of 'N' column vectors of length 'n': a = Double_Type[n,N]. Each of the 'n' rows of 'a' may represent the realization of a random variable. The 'n' times 'n' covariance matrix 'cov_mat' between these 'n' variables is cov_mat[i,j] = 1/(N-1) * sum( (a[i,*]-mean(a[i,*]))*(a[j,*]-mean(a[j,*])). The 'n' times 'n' correlation matrix 'cor_mat' is cor_mat[i,j] = cov_mat[i,j]/sigma_i/sigma_j where 'sigma_i' is the square root of the variance of 'a_i': sigma_i = sqrt( 1/(N-1) * sum( (a[i,*]-mean(a[i,*]))^2 ) ).
Notes
The Cholesky decomposition of the covariance matrix can be used to generate correlated variables that obey a given covariance matrix. See the example section for details.
Example
% From three uncorrelated normal variables 'g', generate three correlated normal % variables 'x' that obey the covariance matrix 'm':
g = Double_Type[3,100000]; g[0,*] = grand(100000); g[1,*] = grand(100000); g[2,*] = grand(100000); (cov_mat, cor_mat) = covariance_correlation_matrix(g); print(cov_mat); % the covariance matrix of 'g' shows that there are no correlations print(cor_mat);
m = Double_Type[3,3]; m[0,*] = [ 9,-3,-6]; m[1,*] = [-3,10, 5]; m[2,*] = [-6, 5, 6]; cd = cholesky_decomposition(m); x = cd#g; % matrix product of 'cd' and 'g' (cov_mat, cor_mat) = covariance_correlation_matrix(x); print(cov_mat); % the covariance matrix of 'x' is indeed (almost) 'm' connect_points(0); plot(x[0,*],x[2,*]);
See also: cholesky_decomposition
create_struct_field
Synopsis
create and set the value associated with a structure field
Usage
Struct_Type create_struct_field(s, field_name, field_value);
Qualifiers
skip - do not update an existing field
Description
Does the same as the `set_struct_field' function except that the given field is created first if it does not exist.
See also: set_struct_field, struct_combine
crossing_number_polygon
Usage
cn=crossing_number_polygon(P,V)
Synopsis
return the crossing number for a point P in a polygon
Description
Return the crossing number for a point P in a (potentially complex polygon defined by the vertex points V). The polygon has to be closed, i.e. V.x[n]==V.x[0] and V.y[n]==V.y[0] where n is the number of polygon points.
The crossing number is the number of times that a ray starting from point P crosses the polygon. The point is outside of the polygon if the crossing number is even, and inside if the crossing number is odd (even-odd rule).
The point P is a struct{x,y}, while the polygon V is defined by its vertex points, which are stored in a struct{x[],y[]} where the arrays are the x- and y-coordinates.
See the URL below for more explanations.
Based on code by Dan Sunday, http://geomalgorithms.com/a03-_inclusion.html and patterned after Randolph Franklin (2000), http://www.ecse.rpi.edu/Homepages/wrf/research/geom/pnpoly.html
See also: winding_number_polygon,point_in_polygon,simplify_polygon
cyl2cart
Synopsis
Convert cylindrical to Cartesian coordinates
Usage
cyl2cart(Double_Types r[], phi[], z[], vr[], vphi[], vz[]);
Description
Convert cylindrical (r,phi,z,vr,vphi,vz) to Cartesian (x,y,z,vx,vy,vz) coordinates with phi rotating counter-clockwise from the positive x-axis.
Qualifiers
Example
(x,y,z,vx,vy,vz) = cyl2cart(1,PI/2,1,2,3,4); (x,y,z,vx,vy,vz) = cyl2cart([1,0],[PI/2,0],[1,0],[2,0],[3,0],[4,0]); (x,y,z,vx,vy,vz) = cyl2cart( cart2cyl(1,1,1,2,3,4) );
See also: cart2cyl
dateOfMJD
Synopsis
calculates the date from its MJD value
Usage
Struct_Type date = dateOfMJD(MJD);
Description
The return value is a structure of the following form:
date = struct { year, month, day, hour, minute, second }
and always in the Gregorian Calendar. See DateOfJD if you want more
functionality. This function is array safe.
See also: DateOfJD, MJDofDate
dayOfWeek
Synopsis
returns the day of the week of the given date
Usage
Integer_Type dayOfWeek(Integer_Type year, month, day);
Description
This function uses the shell command `date' in order to return the day of the week from 1 (Monday) to 7 (Sunday) of the given date.
Example
dayOfWeek(1837, 9, 9);
daysPerMonth
Synopsis
returns the number of days of a month
Usage
Char_Type daysPerMonth(m[, y])
Description
m
= 1, 2, ... 12 specifies the month Jan, Feb, ..., Dec.
If y
is not specified, a non-leapyear is assumed.
See also: leapyear
deg2dms
Synopsis
Convert a floating point angle in degrees to degrees, minutes, seconds
Usage
(d,m,s) = deg2dms(degree)
Qualifiers
%
- hours: If set, return the coordinate in hours, minutes, sec (RA)
- radian: If set, the angle is in radians, not degrees
Description
This is a convenience routine to convert astronomical coordinates given as a floating point number into a number in degrees, minutes, seconds.
For negative angles, the sign of the coordinate is given to the highest non-zero integer.
Use angle2string to generate strings from angles, and generate_iauname if you want to build strings to use in source names.
See also: dms2deg, hms2deg, angle2string, generate_iauname
diff
Synopsis
returns the differences between adjacent elements in an array
Usage
Double_Type diff(Double_Type[] array);
Description
a = sqr([1:4]); % [1,4,9,16] b = diff(a); % [3,5,7]
dms2deg
Synopsis
Convert angle in degree, minute, seconds to floating points degrees or radian
Usage
degree = dms2deg(d,m,s);
Qualifiers
- hours: If set, the coordinate is given in hour, minutes, sec (RA)
- radian: If set, return angle in radians, not degrees
Description
This is a convenience routine to convert astronomical coordinates given in degrees, minutes, and seconds into a floating point number. The sign of the highest non-zero argument decides the sign of the returned floating point number. In other words: The routine tries to be intelligent for negative declinations in the 0>=dec>-1 strip. For example, if the declination is -0d 14' 30.2", call this routine as dms2deg(0,-14,30.2).
Example
% Convert the coordinates % alpha=19h 49m 35.49s and dec=-30d 12' 31.8" variable ra=hms2deg(19,49,35.49); variable dec=dms2deg(-30,12,31.8);
See also: hms2deg
douglas_peucker
Synopsis
Simplify polygons using the Douglas Peucker Algorithm
Usage
(xx,yy)=douglas_peucker(x,y,d2);
Description
This function implements an algorithm to simplify complex polygons (Douglas & Peucker, 1973, Cartographica 10(2), 112). A recursive version presented by Hershberger & Snoeyink (1992, Proc. 5th Intl. Symp. on Spatial Data Handling, 134-143) is used. For a polygon P defined by positions in given by the arrays (x,y), the algorithm returns the polygon P2 defined by coordinates (xx,yy) with the property that the maximum distance between the lines segments of P and P2 is smaller than sqrt(d2). Note 1: The algorithm only yields a good result if P is not self-intersecting. Note 2: As discussed by Hershberger & Snoeyink, the worst performance of the Douglas-Peucker-algorithm is O(N^2). Hershberger & Snoeyink (1998, Computational Geometry 11(3-4), 175-185) present a O(N log N) algorithm. Unfortunately for me (J. Wilms), this algorithm is too complex to be implementable on the ICE between Hamburg and Bamberg... Note 3: This function is called recursively and does not perform sanity checks on x,y,d2. If you do not know about the properties of your data, use simplify_polygon.
See also: simplify_polygon
draw_progress_bar
Synopsis
Draw a progress bar
Usage
draw_progress_bar( position , maximum )
Qualifiers
- tip: String to draw tip of the arrow (default: ">")
- bar: String to draw bar of the arrow (default: "=")
- front: String to draw space in front of arrow (default: ".")
- append: String to append to each drawn progress bar (default: "")
- columns: Columns to use to write progress bar (default: Terminal width)
- fmt: Format to use for writing the percentage info in front of the bar. The printing function gets passed the percentage of the current state to the function as the second argument after the format (default: "%.1f%%", for example 01.3%)
Description
Draw a progress bar across the prompt. The width of the bar (arrow) is calculated from the fraction of position/maximum. In case position is greater than maximum nothing happens. Note that the function automatically determines the width of the current terminal. This procedure takes some tens of milliseconds. The function can be sped by setting the "columns" qualifier manually.
ecliptical2equatorial
Synopsis
convert ecliptical (lambda,beta) to equatorial coordinates
Usage
(alpha,delta)=ecliptical2equatorial(lambda,beta;qualifiers)
or
```c
Vector_Type eqp=equatorial2ecliptical(ecl;qualifiers)
##### Qualifiers
* equinox: equinox of the transformation. Float: JD, string: equinox
designation (e.g., "J2000.0" or "B1950.0";
default: J2000.0)
* deg: interpret angular arguments in degrees (default is radians!)
applies also to the return value.
* mjd: interpret equinox as a MJD (default: JD)
##### Description
The routine takes coordinates in the ecliptical (lambda, beta)
system and converts them into equatorial coordinates (right ascension,
declination), using the obliquity of the ecliptic for the equinox
and vice versa.
Alternatively, the routine also accepts a Vector_Type and then
returns an Vector_Type vector in equatorial coordinates (or
vice versa if the inverse qualifier is given).
The default equinox is J2000.0.
This routine just calls equatorial2ecliptical with all arguments and the
inverse qualifier, but is a simpler interface.
__See also__: Vector_Type, equatorial2ecliptical, JDofEpoch,dms2deg,hms2deg
#### edit_var
##### Synopsis
allows to edit S-Lang variables in an editor
##### Usage
```c
edit_var(&x);
or
Any_Type y = edit_var(Any_Type x);
Qualifiers
- tmpfile: temporary file [default: /tmp/edit_var_$UID_$PID]
Description
edit_var supports the following data types:
- Undefined_Type (=Void_Type), Null_Type,
- Integer_Type, Double_Type, Complex_Type
- String_Type, BString_Type as well as
- Array_Type
- Assoc_Type
- Struct_Type
- List_Type in a recursive way. (Circular linked list are currently not supported.)
S-Lang code defining the variable x is shown in the editor specified by the EDITOR environ- ment variable or jed, if EDITOR is undefined. edit_var uses jed's folding mode (one should run 'Buffers' => 'Folding' => 'Enable Folding') for nested data structures, which can hence be very easily investigated.
The S-Lang code for x can be modified. After saving the temporary file and closing the editor, the file is evaluated, which should return an S-Lang object that is either stored in x (if passed by reference) or returned.
Example
i = edit_var(struct { uc='A', s=1H, us=1UH, l=1L, ul=1UL, ll=1LL, ull=1ULL });
ellipse
Synopsis
calculates points on an ellipse centered at the origin
Usage
(Double_Type X,Y) = ellipse(Double_Type smaj, smin, posang, phi)
Qualifiers
- x [
=1
]: compression/streching factor along the x-axis - y [
=1
]: compression/streching factor along the y-axis
Description
This function calculates the coordinates X
, Y
of an ellipse with
semimajor axis smaj
, semiminor axis smin
, and position angle posang
(measured in rad with respect to the x-axis) which is parameterized with phi
.
The complete ellipse is covered when phi
covers the values from 0 to 2*PI.
Examples
variable phi=[0:2*PI:#200]; plot ( ellipse(5,3,0.2*PI,phi) ); oplot( ellipse(5,3,0.2*PI,phi ; x=0.5) );
See also: enclosing_ellipse
empty_struct
Synopsis
Returns a struct with 0 fields
Usage
Struct_Type = empty_struct();
Description
This function returns a struct with 0 fields, which is usable with other funtions, e.g., struct_field_exists.
enclosing_ellipse
Synopsis
calculates the smallest ellipse enclosing two ellipses centered at the origin
Usage
(Double_Type smaj,smin,posang) = enclosing_ellipse(Double_Type smaj1, smin1, posang1, smaj2, smin2, posang2)
Description
This function calculates the semimajor axis smaj
, the semiminor axis smin
,
and the position angle posang
of the smallest ellipse enclosing
two ellipses centered at the origin.
Examples
smaj1 = 5.; smin1 = 2.; posang1 = 0.3; smaj2 = 4.; smin2 = 1.; posang2 = 0.7*PI; variable phi=[0:2*PI:#200];
plot ( ellipse( enclosing_ellipse(smaj1, smin1, posang1, smaj2, smin2, posang2) ,phi) ); oplot( ellipse(smaj1, smin1, posang1 ,phi) ); oplot( ellipse(smaj2, smin2, posang2 ,phi) );
See also: ellipse
equatorial2ecliptical
Synopsis
convert equatorial (alpha,delta) to ecliptical (lambda,beta) coordinates
Usage
(lambda,beta)=equatorial2ecliptical(alpha,delta;qualifiers)
or
```c
Vector_Type ecl=equatorial2ecliptical(eqp;qualifiers)
##### Qualifiers
* equinox: equinox of the transformation. Float: JD, string: equinox
designation (e.g., "J2000.0" or "B1950.0";
default: J2000.0)
* deg: interpret angular arguments in degrees (default is radians!)
applies also to the return value.
* mjd: interpret equinox as a MJD (default: JD)
* inverse: perform the conversion from ecliptical to equatorial
coordinates (see ecliptical2equatorial() for an
equivalent interface)
##### Description
The routine takes coordinates in the equatorial (right ascension,
declination) system and converts them into ecliptical coordinates
using the obliquity of the ecliptic for the equinox and vice versa..
Alternatively, the routine also accepts an Vector_Type and then
returns a vector in ecliptical coordinates (or vice versa if the inverse
qualifier is given).
The default equinox is J2000.0.
__See also__: Vector_Type, ecliptical2equatorial, JDofEpoch,dms2deg,hms2deg
#### equatorial2galactic
##### Synopsis
convert J2000.0 equatorial (alpha,delta) to galactic (l,b) coordinates
##### Usage
```c
(l,b)=equatorial2galactic(alpha,delta;qualifiers)
or
```c
Vector_Type gal=equatorial2galactic(eqp;qualifiers)
##### Qualifiers
* deg: : interpret angular arguments in degrees (default is radians!)
applies also to the return value.
* inverse: : perform the conversion from galactic to equatorial
coordinates (see galactic2equatorial() for an
equivalent interface)
* blaauw: : use the original Blaauw definition of the IAU galactic
coordinate system; see below for caveats
##### Description
The function takes coordinates in the equatorial (right ascension,
declination) system and converts them into galactic
coordinates (or vice versa, if the inverse qualifier is given) for
the IAU 1958 Galactic coordinate system (Blaauw et al., 1960,
MNRAS 121, 123).
Alternatively, the routine also accepts a Vector_Type and then
returns a vector in galactic coordinates (or vice versa if the
inverse qualifier is given).
The function is for J2000.0 coordinates only, use the
precess function to convert to that equinox if needed.
See Lane (1979, PASP 91, 405) for a discussion of the subtleties of
this conversion and Johnson & Soderblom (1987, AJ 93, 864) for a
non-standard derivation of the conversion matrix (the one used here
is basically identical, but uses proper Euler angles). By default,
the conversion is done using the precessed pole and inclination for
the FK5 system (J2000.0) given by Liu et al (2011, A&A 526, A16),
which is a thorough derivation of values that are essentially
identical to those also given in the appendix of Reid & Brunthaler
(2004, ApJ 616, 872).
If the qualifier Blaauw is given, then the coordinates are first
precessed to B1950.0 and then transformed. Note that this approach
is not good for the highest precision, since no conversion from,
e.g., FK5 to FK4 coordinates is performed as this would result in a
non-orthogonal coordinate system. Please use this qualifier only
if you know what you are doing and read Liu et al. (2011) before
using the qualifier.
__See also__: Vector_Type, galactic2equatorial,precess,dms2deg,hms2deg
#### equatorial2horizon
##### Synopsis
convert equatorial (alpha,delta) to horizon (elevation, azimuth) coordinates
##### Usage
```c
(azi,ele)=equatorial2horizon(alpha,delta;qualifiers)
or
```c
Vector_Type hor=equatorial2horizon(eqp;qualifiers)
##### Qualifiers
* JD: JD for which the calculation is to be performed. Mandatory.
* lon: geographic longitude of the observer, positive towards the east. Mandatory.
* lat: geographic latitude of the observer, positive towards the north. Mandatory.
* deg: interpret all angular arguments in degrees (default is radians!)
applies also to the return value.
* mjd: interpret JD argument as a MJD (default: JD)
* inverse: perform the conversion from horizontal to equatorial
coordinates (see horizontal2equatorial() for an
equivalent interface)
##### Description
The routine takes coordinates in the equatorial (right ascension,
declination) system for the equinox and ecliptic of the date and converts
them into the elevation and azimuth for a given observer's position and
time. Note that the JD, lon, and lat qualifiers are mandatory!
Formally the coordinates have to be topocentric coordinates in
order to include parallax effects. For most applications this is not
important as long as a precision of a few arcsec or worse is needed.
Eventually a routine applying all of these effects will be provided.
Alternatively, the routine also accepts an Vector_Type and then
returns a vector in horizon coordinates (or vice versa if the inverse
qualifier is given).
__See also__: Vector_Type, horizontal2equatorial, JDofEpoch,dms2deg,hms2deg
#### erfinv
##### Synopsis
Computes the inverse error function in abs(z)<1
##### Usage
```c
Double_Type erfinv(Double_Type z);
Description
Can be used to calculate a confidence level by
sqrt(2)*erfinv(fraction)
The function is implemented according to Sergei Winitzki's 2008
publication: "A handy approximation for the error function and its
inverse". The approximation has a relative error <1.3E-4.
See also: cerf,cerfc
eval_array
Synopsis
evaluate an expression for an array of values
Usage
Array_Type eval_array(Type_Type type, String_Type expression, Array_Type values)
or
eval_array(String_Type expression, Array_Type values);
Qualifiers
- marker [
="*"
]: place holder in theexpression
string, where thevalues
are to be inserted - n_markers [=all]: number of markers in
expression
to be replaced; if positive, the firstn_markers
are replaced, if negative, the last|n_markers|
are replaced.
Description
The marker
string in expression
is sequentially replaced with
each value of the array values
(note that the values are converted
through the string
function), and the resulting string is evaluated.
Unless type==Void_Type
, each evaluation is expected to result in one value
of this type
, and eval_array
returns an array of all these values.
If type==Void_Type
(which can be omitted), no return value is expected.
Examples
%
set fields in an array of structures
public variable s = Struct_Type[8]; _for $1 (0, 7, 1) s[$1] = struct { a, b };
eval_array("s[*].a = 1", [0:3]); % => s[[0:3]].a = 1;
eval_array("s[*].a = 2", [4:7]); % => s[[4:7]].a = 2;
%
using the array several times
eval_array("s[*].b = *", [0:3]); % => s[[0:3]].b = [0:3];
eval_array("s[*].b = * * 2", [4:7]; n_markers=2); % => s[[4:7]].b = [4:7] * 2;
%
using a different marker to clarify the last example
eval_array("s[#].b = 2 * #", [4:7]; marker="#"); % => s[[4:7]].b = 2 * [4:7];
See also: array_map, array_struct_field, eval
factorial
Synopsis
calculates n!
Usage
Double_Type factorial(Integer_Type n);
Description
n
! = n * (n-1) * ... * 2 * 1
Note that n
is always converted to an integer (without rounding);
for fractional n
use, e.g., the GSL module's Gamma function.
See also: gsl->gamma
fermi2MJD
Synopsis
calculate MJD from Fermi seconds
Usage
Double_Type = fermi2MJD (Double_Type);
Description
Calculates MJD(UTC) for a given Fermi Mission Elapsed Time (MET) in seconds.
Example
variable m = 239557420.0; variable fermi_s = fermi2MJD(m);
See also: MJD2fermi, MJDref_satellite
find_correlations
Synopsis
calculates all correlations between columns of a table
Usage
Struct_Type info = find_correlatins(Struct_Type s);
Description
info.corr[i] = correlation_coefficient(s.<info.x[i]>, s.<info.y[i]>);
See also: correlation_coefficient
find_function_maximum
Synopsis
looks for the position of a function's maximum value
Usage
Double_Type x0 = find_function_maximum(Ref_Type &f, Double_Type x1, Double_Type x2);
or
Double_Type x0 = find_function_maximum(Ref_Type &f, Double_Type x1, Double_Type x2, &f_x0);
Qualifiers
- qualifiers: structure of qualifiers to be passed to f
- eps [=1e-12]
Description
f
has to be a real function with one argument.
A binary search is performed to find x0
such that
f(x0) = max( f([x1:x2]) )
.
If f
is not convex in [x1:x2]
, the algorithm does not need
to succeed. Otherwise, the accuracy of x0
is (x2-x1)*eps
.
See also: find_function_value
find_function_value
Synopsis
computes an inverse function
Usage
Double_Type x0 = find_function_value(Ref_Type &f, Double_Type val, x1, x2);
Qualifiers
- qualifiers: structure of qualifiers to be passed to f
- eps [=1e-12]
- quiet: do not show error message
Description
f
has to be a real function with one argument.
A binary search is performed to find x0
such that
f(x0) = val
.
If f
is not strictly monotonic in [x1:x2]
, the algorithm does not
need to succeed. Otherwise, the accuracy of x0
is (x2-x1)*eps
.
See also: find_multiargumentfunction_value, find_function_maximum
find_multiargumentfunction_value
Synopsis
computes an inverse function
Usage
Double_Type xi0 = find_multiargumentfunction_value(Ref_Type &f, Double_Type val, x1, x2, ..., xn);
Qualifiers
- qualifiers: structure of qualifiers to be passed to f
- eps [=1e-12]
Description
f
has to be a real function with n arguments.
While xi = [xi_min, xi_max]
is an array, xj
(for j!=i
) are constants.
A binary search is performed to find xi0
such that
f(x1, x2, ..., xn) = val
for xi = xi0
.
If f
is not strictly monotonic for xi_min < xi < xi_max
, the algorithm does not
need to succeed. Otherwise, the accuracy of xi0
is (xi_max-xi_min)/1e12
.
See also: find_function_value
first_valid_digit
Synopsis
gives the position of the first valid digit
Usage
Integer_Type fvd = first_valid_digit( Integer/Double_Type x )
;
fread_struct
Synopsis
Read binary data from a file into a pre-defined structure
Usage
fread_struct(Struct_Type s, File_Type fp);
Qualifiers
char_to_string - convert fields of an array of chars (Char_Type[]) with length greater one into strings chatty - be verbose
Description
This function uses fread' to reads binary data from a file into the fields of a structure. All fields have to have a defined data type! That is at least each field value has to be set to a certain
DataType_Type'. In case field values are arrays, the
corresponding amount of objects of the array's data type are read
from the file. The fields have to be defined in the same order as
their corresponding objects should be read from the file. Note
that the function does not return anything, but the fields of the
given structure are updated instead! Finally, make sure to use
the correct number of bits to be read for each field. In doubt
always use, e.g., Int32_Type instead of Integer_Type as the
latter might depend on how S-lang was compiled.
Example
% define the structure to be read variable s = struct { count = Int16_Type, % first, read one 16-bit integer list = Float32_Type[10], % secondly, read 10x 32-bit floats msg = Char_Type[64] % finally, read 64x characters }; % read the file variable fp = fopen("mybinaryformat.file", "r"); fread_struct(s, fp; char_to_string); ()=fclose(fp); % print the message, which will be a string due to the % char_to_string-qualifier message(s.msg);
See also: fread
frequency2velocity
Synopsis
Calculate the velocity to a given frequency.
Usage
frequency2velocity(freq,restfreq);
Qualifiers
veldef - The velocity definition which must be one of OPTICAL, RADIO, or TRUE. Defaults to RADIO.
Description
Convert frequency to velocity (m/s) using the given rest frequency and velocity definition. The units (Hz, MHz, GHz, etc) of the frequencies to convert must match that of the rest frequency argument.
Adopted from the gbtidl script freqtovelo.pro
ftest_xspec
Synopsis
calculates the F-test probability as in xspec
Usage
Double_Type ftest_xspec(chisq2, dof2, chisq1, dof1)
Description
The new chi-square and DOF, chisq2 and dof2, should come from adding an extra model component to (or thawing a frozen parameter of) the model which gave chisq1 and dof1. If the F-test probability is low then it is reasonable to add the extra model component.
WARNING: It is not correct to use the F-test statistic to test for the presence of a line (see Protassov et al astro-ph/0201457).
g_earth
Synopsis
Calculates the acceleration at the Earth's surface as a function of latitude
Usage
Double_Type g = g_earth(Double_Type phi);
Qualifiers
- deg: phi is in degrees, not in radian.
- altitude: altitude of the observer above the geoid [m]
- mks: return g in m/s^2, not in cm/s^2
- allen: use the equation of Allen (1973, Astrophys. Quantities)
- almanac: Urban & Seidelman (2013; Expl. Suppl. Astron. Almanac)
Description
This code calculates the surface acceleration on the earth as a function of the latitude (in rad). The equations used are from a compilation of different approximations to g by Mangum & Wallace (2015), PASP 127, 74.
The default returned is for the WGS84 geoid.
This routine is array safe in phi or altitude.
See also: dms2deg
gal2cel
Synopsis
Transform cartesian Galactic to celestial coordinates
Usage
gal2cel(Double_Types x[], y[], z[], vx[], vy[], vz[]; qualifiers)
Description
Transform right-handed, cartesian Galactic coordinates (x [kpc], y [kpc], z [kpc], vx [km/s], vy [km/s], vz [km/s]) with the Galactic center at the origin, the Sun on the negative x-axis and the z-axis pointing to the north Galactic pole implying clockwise Galactic rotation when seen from the half space with positive z to celestial coordinates (right ascension [h, m, s], declination [deg, arcmin, arcsec], distance [kpc], radial velocity [km/s], proper motion in right ascension times cosine of declination [mas/yr], proper motion in declination [mas/yr]).
Qualifiers
- SunGCDist: Sun-Galactic center distance [kpc]; default: 8.4 (see Model I in Irrgang et al., 2013, A&A, 549, A137)
- vxs: Sun's x-velocity component [km/s] relative to the local standard of rest; default: 11.1 (Schoenrich, Binney & Dehnen, 2010: MNRAS 403, 1829)
- vys: Sun's y-velocity component [km/s] relative to the local standard of rest; default: 12.24 (Schoenrich, Binney & Dehnen, 2010: MNRAS 403, 1829)
- vzs: Sun's z-velocity component [km/s] relative to the local standard of rest; default: 7.25 (Schoenrich, Binney & Dehnen, 2010: MNRAS 403, 1829)
- vlsr: Local standard of rest velocity [km/s]; default: 242 (see Model I in Irrgang et al., 2013, A&A, 549, A137)
- GC_NGP: Celestial coordinates of the Galactic center and the north Galactic pole in the format [ah_GC, am_GC, as_GC, dd_GC, dm_GC, ds_GC, ah_NGP, am_NGP, as_NGP, dd_NGP, dm_NGP, ds_NGP]; default: [17, 45, 37.224, -28, 56, 10.23, 12, 51, 26.282, 27, 07, 42.01] (Reid & Brunthaler, 2004, ApJ, 616, 872)
Example
(ah, am, as, dd, dm, ds, dis, vrad, pma, pmd) = gal2cel(-8.4, 0, 0, 11.1, 254.24, 7.25); % -> Sun's position and velocity (ah, am, as, dd, dm, ds, dis, vrad, pma, pmd) = gal2cel(0, 0, 0, 0, 0, 0); % -> Galactic center's position and velocity (ah, am, as, dd, dm, ds, dis, vrad, pma, pmd) = gal2cel([-8.4,0], [0,0], [0,0], [11.1,0], [254.24,0], [7.25,0]); % -> position and velocity of Sun and Galactic center stored in the same array (ah, am, as, dd, dm, ds, dis, vrad, pma, pmd) = gal2cel(-8.4, 0, 0, 11.1, 254.24, 7.25; GC_NGP=[17, 45, 37.1991, -28, 56, 10.221, 12, 51, 26.2755, 27, 7, 41.704]); % -> use non-standard coordinates for Galactic center and north Galactic pole (ah, am, as, dd, dm, ds, dis, vrad, pma, pmd) = gal2cel( cel2gal(17, 45, 37.224, -28, 56, 10.23, 8.33, -11.1, -2.91, -5.12) );
See also: cel2gal, rad2RD
galLB_from_RAdec
Synopsis
convert equatorial (RA, dec) coordinates to galactic (l, b) coordinates
Usage
(l, b) = galLB_from_RAdec(RA, dec);
Qualifiers
- ra_unit [
="deg"
]: set the unit of the right ascensionRA
l_unit="rad"
=>
RA
in radl_unit="hms"
=>
RA
in hours as a scalar or an array of the form[H, M]
or[H, M, S]
- dec_unit [
="deg"
]: set the unit of the declinationdec
Description
This function is deprecated. Please use equatorial2galactic instead.
See also: equatorial2galactic,galactic2equatorial
galactic2equatorial
Synopsis
convert galactic (l,b) coordinates to J2000.0 equatorial (alpha,delta)
Usage
(alpha,delta)=galactic2equatorial(l,b;qualifiers)
or
```c
Vector_Type eqp=equatorial2galactic(gal;qualifiers)
##### Qualifiers
* deg: interpret angular arguments in degrees (default is radians!)
applies also to the return value.
##### Description
The function takes coordinates in the galactic (l,b) system
and converts them into the equatorial J2000.0 system. Alternatively
the routine also accepts a Vector_Type and then returns a
vector in equatorial coordinates for J2000.0.
The function is for J2000.0 coordinates only, use the
precess function to convert the result to another equinox
if needed.
See the description of the function equatorial2galactic for
references and further information. This routine is just a
wrapper around that function, with the qualifier "inverse" set.
__See also__: Vector_Type, equatorial2galactic,precess,dms2deg,hms2deg
#### gcrs2j2000_matrix
##### Synopsis
return the frame bias matrix to convert coordinates from the GCRS to J2000.0
##### Usage
```c
Matrix33_Type mat=gcrs2j2000_matrix();
Description
This function returns the frame bias matrix needed to convert a vector in the Geocentric Celestial Reference System (GCRS) into the J2000.0 dynamical reference system (see astronomical almanac, section B)
See also: precess,precession_matrix
generate_iauname
Synopsis
for a given position, generate a coordinate string obeying the IAU convention
Usage
string=generate_iauname(ra,dec)
Qualifiers
%
- prefix: prefix for the string (e.g., "XMMU_")
- radian: if set, the angles are in radians, not degrees
Description
This is a convenience routine to produce strings of the type XMMU Jhhmmss.s+ddmmss from J2000.0 positions. The routine obeys the IAU convention of truncating (not rounding!) the coordinate to the digits shown.
The default-string is "Jhhmmss.s+ddmmss", use the prefix qualifier to prepend the mission name.
Use angle2string to format coordinates with appropriate rounding.
See also: dms2deg, hms2deg, angle2string, deg2dms
geographic2vector
Synopsis
Calculate a geocentric vector from geographic or geodetic coordinates.
Usage
Vector_Type=geographic2vector(lambda,phi,hei;qualifiers)
Qualifiers
- deg: if set, then lambda and phi are measured in degrees, (default: radian).
- datum: a string defining the geodetic datum. The default is WGS 84, which corresponds to GPS coordinates. Other geodetic data are IERS 2010, GRS 80, IAU 1976, GRS 67, IAU 1964, Hayford (1924), Clarke (1866), and Airy (1830), per the parameters given in the Astronomical Almanac.
- geocentric: The coordinates are not geodetic, but geocentric. In this case hei is interpreted as the distance from the geocenter (otherwise it is the height with respect to the reference ellipsoid). Note: it is very rare for coordinates to be geocentric! GPS is geodetic.
Description
For a given set of geodetic coordinates, i.e., geographic longitude (positive towards the east!), latitude (positive on the Northern hemisphere), and a height above the reference spheroid (close to height above mean sea level), calculate the x,y,z-position in a geocentric coordinate system. Here, the xy-plane is the Earth's equator, the xz-plane the ITRS meridian (close to Greenwich), and the z-axis points towards the Earth's north pole. The xyz- values are in meters.
Note that the longitude and latitude are assumed to be in rad, unless the deg qualifier is given. The height is in meters.
The routine returns a Vector_Type. If called with arrays, then an array of Vector_Type is returned.
See also: dms2deg
get1match
Synopsis
returns one (the first) matching substring of a regular expression matching
Usage
String_Type match_str = get1match(String_Type str, String_Type regexp);
Description
get1match is deprecated. Use S-Lang's string_matches instead.
See also: string_matches
get_arg_struct
Synopsis
obtains the arguments/options the current isis instances was called with.
Usage
Struct_Type = get_arg_struct( );
or
```c
Struct_Type = get_arg_struct( String_Type[] name, type);
##### Qualifiers
* delim: [="="] Delimeter between argument name and value (e.g., also, "=.,").
* prefix: [="--"] Argument prefix. Arguments missing this prefix are ignored.
##### Description
This function returns a Struct_Type with fields corresponding to
the names of those arguments/options the current isis instance
was called with, which have the given 'prefix'. Thereby the 'name'
is defined as the string enclosed by the 'prefix' and the 'delimeter'.
Expected SYNTAX: <prefix>"argname"<delimeter>"value"
It is possible to give several 'delimeter' as a string-chain (e.g.,
"=,.-"). If an argument/option contains several delimeters, only the
string between the 1st and 2nd delimeter is taken as value for the
corresponding name.
This function can be given an string array 'name' and 'type', assigning
a Data_Type to the argument 'name'. Thereby 'type' can be:
"d" : Integer_Type
"f" : Double_Type
Otherwise, if used without arguments ('name' & 'type'), the
values are returned as String_Type.
If an argument has no value given, i.e. does not contain a delimeter
the corresponding field of the returned structure is equal NULL.
##### Example
/> isis -g --arg --arg1=1e-3 --arg2=1.01 --arg3=abcd --arg4=1=bla
isis> print( get_arg_struct );
{arg=NULL,
arg1="1e-3",
arg2="1.01",
arg3="abcd",
arg4="1"}
isis> print( get_arg_struct("arg1","f") );
{arg=NULL,
arg1=0.001,
arg2="1.01",
arg3="abcd",
arg4="1"}
isis> print( get_arg_struct(["arg1","arg2","arg4"],["f","f","d"]) );
{arg=NULL,
arg1=0.001,
arg2=1.01,
arg3="abcd",
arg4=1}
__See also__: atof, atoi, strreplace, strtok, set_struct_field, array_map, is_substr
#### get_coordinatearrays_of_image
##### Synopsis
returns two 2d arrays of image coordinates
##### Usage
```c
(XX, YY) = get_coordinatearrays_of_image(Any_Type img[,]);
Description
XX[y, x] = x
and YY[y, x] = y
for every coordinate (x,y
) of img
.
get_ephemeris
Synopsis
returns the epoch and period of an ephemeris
Usage
Struct_Type get_ephemeris(String_Type ephemeris)
Description
The list of known ephemerides is stored internally
and is shown if an unknown ephemeris is requested.
Ephemerides are returned as a structure, containing
at least the internal name
of the ephemeris (which
usually contains the source name and a reference),
the epoch T0
, and the period P
.
More information may be stored in further reasonably named
fields, e.g., Pdot
for the change of the period.
See also: orbitalphase
get_grid
Usage
(XX, YY) = get_grid(Double_Type X[], Double_Type Y[]);
Description
XX
and YY
are two-dimensional arrays
which span the grid defined by the arrays X
and Y
.
Example
variable XX, YY; (XX, YY) = get_grid([-2:2:0.01], [-2:2:0.01]); plot_image( atan2(YY, XX) );
get_intersection
Synopsis
finds elements occuring in two arrays
Usage
(Integer_Type i1[], i2[]) = get_intersection(array1[], array2[]);
Description
array1[ i1 ] == array2[ i2 ]
Qualifiers
- ordered: assumes that
array1
andarray2
are ordered increasingly
get_par_combinations
Synopsis
calculates all possible parameter combinations
Usage
Struct_Type[] parcomb = get_par_combinations( Struct_Type par )
Description
All possible parameter combinations are calculated for the given parameter structure. That means if there are multiple parameters, each with there own grid all possible combinations are returned.
'par' must be a Struct_Type, but there are no restriction otherwise. Each field is representing a parameter and its value its grid, which can be any Data_Type (e.g., Double_Type or String_Type)!
NOTE that the combination is done, s.t. the last parameter in the struct is varied first and the first parameter last, i.e., the first parameter in the parameter combination will stay constant longest!
Examples
variable par = struct{ a = [0:1:#3 (closed)], b = ["hello","world"], c = [ 2001, 2012] }; variable parcomb = get_par_combinations( par ); print(parcomb);
get_struct_fields
Synopsis
returns several fields of a structure
Usage
(Any_Type val1, val2, ...) = get_struct_fields(Struct_Type s, String_Type fieldname1, fieldname2, ...);
Qualifiers
- i: indices used for filtering array fields
Description
Each value corresponds to the according field of the structure.
If an i
qualifier is given, array-typed field values
are filtered with these indices, i.e.,
val = get_struct_feld(s, fieldname)[i];
It also possible to create lists of structure field values, by just passing lists (or, equivalently, arrays) as arguments.
Examples
variable table = struct { x=[1:5], y=[1:5]^2, err1=[1:5], err2=[1:5]\\*2 }; plot_with_err( get_struct_fields(table, "x", "y", "err1"; i=[1,0,4,2,3]) ; connect_points); % equivalent to: % plot_with_err( table.x[1,0,4,2,3], table.y[1,0,4,2,3], table.err1[1,0,4,2,3] ; connect_points); plot_with_err( get_struct_fields(table, "x", "y", {"err1", "err2"}) ); % or (even less redundant): plot_with_err( get_struct_fields(table, "x", "y", "err"+["1", "2"]) ); % equivalent to: % plot_with_err( table.x, table.y, {table.err1, table.err2} );
See also: get_struct_field
get_variable_name
Synopsis
returns the namespace and name of a given reference
Usage
(namespace, name) = get_variable_name(&var);
Description
The name of the variable or function the given reference points to will be returned as a string. The namespace where it is defined is returned as well. In case of a private namespace, however, NULL will be returned instead.
See also: current_namespace
greatcircle_coordinates
Synopsis
calculates the coordinates of the greatcircle between two points on a sphere
Usage
(Double_Type lambda[], phi[]) = greatcircle_coordinates(lambda1, phi1, lambda2, phi2);
Qualifiers
- unit [
="deg"
]: unit of the angular coordinates - delta [
=0.5
]: angular step in degrees
Description
(lambda
i, phi
i) are the spherical coordinates of point i.
unit="deg"
: lambda
i and phi
i are in degrees.
They can be scalar values or arrays of the form
[deg, arcmin]
, [deg, arcmin, arcsec]
or [sign, deg, arcmin, arcsec]
.
unit="rad"
: lambda
i and phi
i are scalars in radian.
unit="hms"
: code{lambda}i are scalars hour angles (24h = 360deg)
or arrays in h:m:s format, i.e., [h, m]
or [h, m, s]
.
The phi
i are nevertheless in degrees as above.
See also: greatcircle_distance
greatcircle_distance
Synopsis
calculates the angular distance between two points on a sphere in radians
Usage
Double_Type greatcircle_distance(alpha1, delta1, alpha2, delta2)
Qualifiers
- unit [
="deg"
]: unit of the input angular coordinates - alpha1_unit [
="deg"
]: unit of the alpha1 - alpha2_unit [
="deg"
]: unit of the alpha2 - delta1_unit [
="deg"
]: unit of the delta1 - delta2_unit [
="deg"
]: unit of the delta2
Description
(alpha
i, delta
i) are the spherical coordinates of point i.
unit="deg"
: alpha
i and delta
i are in degrees.
They can be scalar values or arrays of the form
[deg, arcmin]
, [deg, arcmin, arcsec]
or [sign, deg, arcmin, arcsec]
.
unit="rad"
: alpha
i and delta
i are scalars in radian.
unit="hms"
: alpha
i are scalars hour angles (24h = 360deg)
or arrays in h:m:s format, i.e., [h, m]
or [h, m, s]
.
The delta
i are in degrees as above.
The units of each coordinate can be set independently.
Note that independent of the unit setting, the returned great circle
distance will always be in radian.
See also: greatcircle_coordinates
gridmapping
Synopsis
computes a table of a 2d mapping
Usage
Struct_Type data = gridmapping(&getxy, X, Xfine, Y, Yfine);
hardnessratio_error_prop
Synopsis
calculates a ratio and error propagation
Usage
(hr, err) = hardnessratio_error_prop(h, h_err, s, s_err);
Description
hr = (h-s) / (h+s)
err = sqrt[ (2s/(h+s)^2 * h_err)^2 + (2h/(h+s)^2 * s_err)^2 ]
See also: ratio_error_prop
hms2deg
Synopsis
Convert angle in hour, minute, seconds to degrees or radian
Usage
degree = hms2deg(h,m,s);
Qualifiers
- radian: If set, return angle in radians, not degrees
Description
This is a convenience routine to convert astronomical coordinates given in hours, minutes, and seconds into a floating point number.
The routine is equivalent to calling dms2deg with the hours qualifier.
See also: dms2deg
horizon2equatorial
Synopsis
convert horizon (azi,ele) coordinates to equatorial (alpha,delta)
Usage
(alpha,delta)=horizon2equatorial(azi,ele;qualifiers)
or
```c
Vector_Type eqp=horizon2equatorial(hor;qualifiers)
##### Qualifiers
* JD: JD for which the calculation is to be performed. Mandatory.
* lon: geographic longitude of the observer, positive towards the east. Mandatory.
* lat: geographic latitude of the observer, positive towards the north. Mandatory.
* deg: interpret all angular arguments in degrees (default is radians!)
applies also to the return value.
* mjd: interpret date as a MJD (default: JD)
##### Description
The function takes coordinates in the horizon (azimuth,elevation) system
and converts them into the equatorial system for the ecliptic and
elevation of the date. Alternatively the routine also accepts a
Vector_Type and then returns avector in equatorial coordinates.
See the description of the function equatorial2horizon for
references and further information. This routine is just a
wrapper around that function, with the qualifier "inverse" set.
__See also__: Vector_Type, equatorial2horizon,dms2deg,hms2deg
#### hsl2rgb
##### Synopsis
converts a (hue, saturation, lightness) color to (red, green, blue)
##### Usage
```c
Integer_Type rgb = hsl2rgb(Double_Type h, s, l);
Description
The hue h
is an angle in color-space. Here, h
is between 0 and 1.
h
of red is 0, h
of green is 1/3, and h
of blue is 2/3.
The lightness l
of black is 0, and l
of white is 1.
Like the hue, the saturation 0<=s<=1
matters only for 0<l<1
.
The return value is a(n array of) 24 bit RGB value(s).
See also: rgb2hsl
init_time_structure
Synopsis
creates a time structure
Usage
Struct_Type init_time_structure(Integer_Type Y[, m[, d[, H[, M[, S]]]]])
Description
The tm_year
field of the returned structure
is Y-1900
, and tm_mon
is m-1
.
See also: localtime, mktime, strftime
inner_product
Synopsis
Usage
M[n,...,j] = inner_product( A[n,...,m], B[m,...,j] );
Description
In ISIS the operator # calculates the inner product of two arrays A, B, i.e., contracts the last and first dimension, respectively. A and B can have an arbitrary number of dimensions, it is just required that the last dimension of A and first dimension of B is of the same size!
Example
A = _reshape( [1:2*3], [2,3] ); B = _reshape( [1:3*4], [3,4] ); M = inner_product( A, B ); N = A # B; vmessage("M
"); print(M); vmessage("N
"); print(N);
See also: operator #
inner_products
Synopsis
Usage
M[n,...,l] = inner_products( A[n,...,m], B[m,...,j] [, ..., Z[k,...,l] ] );
Description
In ISIS the operator # calculates the inner product of two or more arrays. This function is basically nesting the inner_product function!
See also: inner_product, operator #
integrate2d
Synopsis
numerical computation of an 2-dim. integral
Usage
Double_Type I = integrate2d( Ref_Type func, y1, y2,
Double_Type x1, x2
);
Qualifiers
- intfun[&qromb]:Reference to the integrator NOTE: Qualifiers are passed to all sub-functions!
Description
This function numerically computes the 2-dimensional integral 'I' over the integrant-function 'func' of the form:
I = int_x1^x2 int_y1(x)^y2(x) func(x,y) dx dy
The given limits 'x1' and 'x2' are the lower and upper integration limits for the 'x' variable, respectively. 'y1' and 'y2', on the other hand, are references to functions, which calculates the lower and upper limit y1(x) and y2(x) for the 'y' variable depending on 'x'.
The implemented solution is adopted from the 'Numerical Recipes in C' Chap. 4.6. !
The integrator can be changed by the qualifier 'intfun', which is set to 'qromb' as default (see 'help qromb). Qualifiers given to 'integrate2d' are passed to the 'intfun' function!
Example
% Calculation of the area of an rectangle: % I = int_0^1 dx int_0^2 dy % = [ x ]_0^1 * [ y ]_0^2 % = 2
variable x1 = 0., x2 = 1.; define y1(x){ return 0.; }; define y2(x){ return 2.; }; define func( x, y ){ return 1.; }; variable I = integrate2d( &func, &y1, &y2, x1, x2 ); vmessage("Area of rectangle is A = %g",I);
% Calculation of the area of the unit circle in cartesian coordinates: % I = int_-1^1 int_-sqrt(1-sqr(x))^sqrt(1-sqr(x)) dx dy % = int_-1^1 2 * sqrt(1-sqr(x)) dx
define tfun( x, y ){ return 1; } define y1(x){ return -sqrt(1-sqr(x)); } define y2(x){ return sqrt(1-sqr(x)); } variable I = integrate2d( &tfun, &y1, &y2, -1., 1. ; qromb_max_iter=24, qromb_eps=1e-3 ); vmessage("Area of unite circle is A = %g, abs(A-PI)/PI = %g", I, abs(I-PI)/PI );
%
See also: qromb, qsimp, integrate2d_test
integrateAB5
Synopsis
integrates an ordinary differential equation with the 5th order Adams-Bashforth algorithm
Usage
x = integrateAB5(&f, t1, t2, dt[, x0]);
Description
The integral of the ordinary differential equation
dx/dt = f(x(t), t) with x(t0) = x0
reads
x(t) = x0 + int_{t0}^{t} f(x(t'), t') dt'
&f
is a reference to a function with two arguments:
define f(x, t)
{
return ...;
}
x
may be be a scalar as well as an array.
See also: integrateRK4
integrateRK4
Synopsis
integrates an ordinary differential equation with the 4th order Runge-Kutta algorithm
Usage
x = integrateRK4(&f, t1, t2, dt[, x0]);
Description
The integral of the ordinary differential equation
dx/dt = f(x(t), t) with x(t0) = x0
reads
x(t) = x0 + int_{t0}^{t} f(x(t'), t') dt'
&f
is a reference to a function with two arguments:
define f(x, t)
{
return ...;
}
x
may be be a scalar as well as an array.
See also: integrateAB5
integrateRK45_adaptive
Synopsis
Integrate an ODE with an adaptive 4th/5th order Runge-Kutta algorithm
Usage
(t, x(t)) = integrateRK45_adaptive(&f, t0, t, [, x0]);
Description
Given an ordinary differential equation of the form
dx/dt = f(t, x(t))
(with x
either a scalar or an array)
this function numerically computes the solution
x(t) = int_{t0}^{t} f(t', x(t')) dt' [ + x0 ]
by means of a 4th/5th order Runge-Kutta (RK) algorithm.
&f
is a reference to a function with two arguments:
define f(t, x)
{
return ...;
}
If x is an array, the adaptive stepsize is chosen according
to the needs of the "worst-offender" equation implying that
the desired accuracy is reached by each individual component.
Note: All qualifiers are also passed to the function f.
Qualifiers
- eps [
=1e-12
]: absolut error control tolerance; lower limit: 1e-15 - method [
="RKCK"
]: choose among three different RK45 methods: "RKF": RK-Fehlberg, "RKCK": RK-Cash-Karp, "RKDP": RK-Dormand-Prince - path: return the entire path and not only the final result
- verbose: show intermediate "times"
t
- plus all qualifiers of the function f
Example
% d/dt x(t) = -sin (t^2)*2t; -> analytical solution: x(t) = cos(t^2) for x(0)=1 define f(t,x) { return -sin(t^2)*2*t; }; (t, x) = integrateRK45_adaptive(&f, 0, 10, 1; path); plot(t,x-cos(t^2));
% d/dt [x1(t),x2(t)] = [x2(t),-x1(t)]; -> analytical solution: x1(t) = cos(t) for x1(0)=1, x2(0)=0 define f(t,x) { [x[1],-x[0]]; } (t, x) = integrateRK45_adaptive(&f, 0, PI, [1,0]; path); plot(t,x[*,0]-cos(t));
See also: integrateRKF45, integrateRK4, integrateAB5
integrateRKF45
Synopsis
integrates an ODE with the adaptive 4th/5th order Runge-Kutta-Fehlberg algorithm
Usage
x = integrateRKF45(&f, t1, t2, dt[, x0]);
Qualifiers
- eps [
=1e-12
]: error control tolerance - verbose: show intermediate "times"
t
Description
This implementation of the adaptive Runge-Kutta-Fehlberg algorithm is still considered EXPERIMENTAL and MAY REQUIRE FURTHER TESTING!
The integral of the ordinary differential equation
dx/dt = f(x(t), t) with x(t0) = x0
reads
x(t) = x0 + int_{t0}^{t} f(x(t'), t') dt'
&f
is a reference to a function with two arguments:
define f(x, t)
{
return ...;
}
x
may be be a scalar as well as an array.
See also: integrateRK4, integrateAB5
integrate_trapez
Synopsis
numerical integration with the trapezoidal rule
Usage
Double_Type int = integrate_trapez(Double_Type x[], Double_Type y[]);
alternative:
Double_Type int = integrate_trapez(Ref_Type function, Double_Type min, Double_Type max);
Description
If two arguments (x[]
and y[] = function(x[])
) are given the integral is calculated via:
int = sum_i (y[i-1] + y[i])/2 * (x[i]-x[i-1])
An alternative usage with three arguments requires a reference to
the function and the integration limits. The function is evaluated
on an equidistant grid. The number of steps can be set with a qualifier.
The second case allows a forth argument specifying a previous result
with a smaller number of steps, which is meant for iterative usage by
the function qsimp
.
Qualifiers
- steps [=100]: number of equidistant grid points, on which the function is evaluated
Example
% Integration of the sin function from 0 to PI/3 % The result should be 0.5 variable x = [0:PI/3:#100];
% The first method: providing grid and corresponding function values integrate_trapez (x, sin(x));
% Alternative method using reference to the function and integration limits integrate_trapez (&sin, 0, PI/3); integrate_trapez (&sin, 0, PI/3; steps = 5000);
See also: integrateRK4, qsimp
interpol2D
Synopsis
Usage
img2 = interpol2D(X2, Y2, X, Y, img);
Description
img = Double_Type[ny, nx];
X = Double_Type[nx];
Y = Double_Type[ny];
min(X) <= X2 <= max(X)
min(Y) <= Y2 <= max(Y)
interpol_2d_table
Usage
Double_Type interpol_2d_table(Double_Type x, y, Struct_Type table)
interpol_polynomial
Synopsis
interpolates a polynomial function between data points
Usage
Double_Type y = interpol_polynomial(Double_Type x, Double_Type X[], Double_Type Y[]);
Description
X
and Y
are arrays of same length n
.
y
is the Lagrange polynomial (of degree n-1
) interpolating
the data points (X[i]
, Y[i]
), evaluated at x
:
n-1 n-1 x - X[j]
y = sum Y[i] * product -----------
i=0 j=0,j!=i X[i]-X[j]
If x
is an array (which doesn't need to be ordered), y
will also be an array.
See also: interpol_points
inversemapping
Synopsis
computes the inverse of a 2d mapping numerically
Usage
(x, y) = inversemapping(&getxy, cache, x_, y_);
Description
(cache.x_, cache.y_) = getxy(cache.x, cache.y); (x_, y_) = getxy(x, y);
j_ff
Synopsis
Calculates the emission coefficient for free-free emission (bremsstrahlung)
Usage
Double_Type j = j_ff(nu,T);
Qualifiers
- Z: average nuclear charge (default=1)
- ne: electron particle density (cm^-3; default: 1e10)
- np: proton particle density (cm^-3; default: 1e10)
Description
This function returns the emission coefficient for free-free radiation (bremsstrahlung). The frequency, nu, is measured in Hz and can be an array, the temperature T is measured in K.
Note: per Kirchhoff's law the total bremsstrahlung spectrum from a slab of size R is B_nu(nu,T)*(1.-exp(-R*alpha_ff(nu,T)))
See also: alpha_ff
jd2mjd
Synopsis
Convert Julian Dates to Modified Julian Dates.
Usage
jd2mjd(JD);
Description
Helper function to convert Julian Dates to MJD by subtracting off 2400000.5.
See also: mjd2jd
jd2year
Synopsis
transform a JD into a fractional year
Usage
Double_Type year = jd2year(Double_Type JD);
Qualifiers
- mjd: The argument is given in MJD, not in JD.
Description
NOTE: This routine takes into account leap years. This is good for plotting, but this means that there are slight phase shifts that could be problems when doing time series analysis. For such purposes best use EpochofJD!
See also: EpochofJD
kendall_tau_censored
Synopsis
calculation of Kendall rank correlation coefficient
Usage
Struct_Type cc = kendall_tau(Double_Type x[], Double_Type y[]);
or
```c
Struct_Type cc = kendall_tau(Double_Type x[], Double_Type y[], Char_Type y_is_upper_limit[]);
or
Struct_Type cc = kendall_tau(Double_Type x[], Char_Type x_is_upper_limit[], Double_Type y[], Char_Type y_is_upper_limit[]);
##### Description
This function calculates the Kendall rank correlation coefficient <code>tau</code>,
as defined by M. G. Kendall (1938, Biometrika, 30, 81), between two
arrays of equal length.
The result <code>tau</code> is within the interval [-1,1], where 1 means a perfect
agreement between both rankings, -1 means perfect disagreement, and values
around 0 indicate that both arrays are independent.
The value of <code>tau</code> is the difference of concordant and discordant pairs
divided by the total number of pairs.
In the definition of Helser (Dennis R. Helsel, "Nondetects and Data Analysis",
Wiley, 2005) only the number of determined (valid) pairs is used instead of
the total number. The structure returned by this function includes the values
of <code>tau</code>, <code>sigma</code>, and <code>pval</code> in both ways.
The value <code>sigma</code> characterizes the expected width of the distribution of <code>tau</code>
for random data. (If the "gsl" module is not available, the p-value is not calculated.
It can be obtained from <code>tau</code> and <code>sigma</code>.)
See statistic module for further functions: require("stats");
##### Examples
variable x = [1:2:#100];
variable y = grand(100)+x;
variable cc_no_limits = kendall_tau_censored(x,y);
variable cc_incl_limits = kendall_tau_censored(x,y,nint(urand(100)));
__See also__: kendall_tau, pearson_r, spearman_r, spearmanrho, correlation_coefficient
#### kerr_rms
##### Synopsis
calculates the radius of marginal stability of a black hole
in units of GM/c^2
##### Usage
```c
kerr_rms(a)
kerr_rplus
Synopsis
calculates the event horizon of a black hole in units of GM/c^2
Usage
kerr_rplus(a)
keyinput
Synopsis
reads input from the keyboard
Usage
String_Type = keyinput([String_Type message]);
Qualifiers
silent - characters are not echoed, e.g. useful for password input time - read timout after given seconds. Returns an empty string nchr - returns automatically after given number of chars have been read prompt - output the given string before reading default - initial text to be read
Description
Prompts the user to input a line, which is ended
with return. If a limited number of chars is given
using the nchr' qualifier, the line is automatically ended after reaching this number without pressing return. It is also possbile to limit the time for an input with the
time' qualifier and to disable echoing the input using
the silent' qualifier. The latter is useful, e.g., for password inputs or for single key events. If
nchr==1', the returned keycode is
checked against special keys such as the arrow keys.
If one of these is detected, the returned string contains
the name if the pressed key.
leapyear
Synopsis
tells whether a year is a leapyear in the Gregorian calendar
Usage
Char_Type leapyear(Integer_Type y)
Description
For a leapyear in the Gregorian calendar, y
is either
divisible by 4 but not by 100,
or y
is divisible by 400.
See also: daysPerMonth
linear_regression
Synopsis
computes a linear regression fit
Usage
Struct_Type = linear_regression(Double_Type[] x, y[, err]);
Description
err
is the uncertainty of the y
values.
If err
is not specified, all data points are weighted equally.
The parameters a
and b
of the best fit y = a + b*x
are a = (Sy*Sxx - Sx*Sxy)/D
and b = (S *Sxy - Sx*Sy )/D
where S = sum(1/err^2)
Sx = sum(x/err^2)
, Sxx = sum(x*x/err^2)
Sy = sum(y/err^2)
, Sxy = sum(x*y/err^2)
and D = S*Sxx - Sx^2
.
The returned structure contains the fields a
and b
,
and the errors da
and db
, respectively,
as well as the reduced chi square chisq
of the fit.
See also: linear_fit_xerr_yerr_data
list_copy
Synopsis
makes a copy of a nested list
Usage
List_Type copy = list_copy(List_Type list);
Description
list_copy
returns a properly dereferenced copy of an
arbitrarily nested List_Type
with all its entries.
Depending on the Data_Type of the entries, list_copy
iteratively calls the according sub-function
(struct_copy
, array_copy
or list_copy
).
If the Data_Type of one of its entries is not
one of privouse mentioned ones, list_copy
returns
@(Data_Type) or respectively @(entry).
Example
See also: COPY, struct_copy, array_copy
lumdist
Synopsis
Calculate luminosity distance (in Mpc) of an object given its redshift
Usage
result = lumdist(z);
Qualifiers
- silent: If set, the program will not display adopted cosmological parameters at the terminal.
- h0 [=70]: Hubble parameter in km/s/Mpc
- omega_k [=0]: curvature constant, normalized to the closure density. Default is 0, indicating a flat universe.
- omega_m [=0.3]: Matter density, normalized to the closure density, default is 0.3. Must be non-negative.
- omega_lambda [=0.7]: Cosmological constant, normalized to the closure density.
- q0 [=-0.55]: Deceleration parameter, numeric scalar = -R*(R'')/(R')^2
Description
INPUTS: z = redshift (positive scalar or vector) OUTPUTS: The result of the function is the luminosity distance (in Mpc) for each input value of z. The luminosity distance in the Friedmann-Robertson-Walker model is taken from Caroll, Press, and Turner (1992, ARAA, 30, 499), p. 511. Uses a closed form (Mattig equation) to compute the distance when the cosmological constant is zero. Otherwise integrates the function using QSIMP. See help for cosmo_param for a description of the cosmological parameters, note that only two out of the four parameters should be given!
The routine can fail to converge at high redshift for closed universes with non-zero lambda. This can presumably be fixed by replacing QSIMP with an integrator that can handle a singularity.
Example
% Plot the distance of a galaxy in Mpc as a function of redshift out % to z = 5.0, assuming the default cosmology (Omega_m=0.3, Lambda = 0.7, % H0 = 70 km/s/Mpc) variable z = [0:5:0.1]; plot (z, lumdist(z));
% Now overplot the relation for zero cosmological constant and % Omega_m=0.3 oplot (z, lumdist(z ; omega_lambda=0 , omega_m=0.3) );
See also: cosmo_param, qsimp
mass_function
Synopsis
Calculates the mass function of a binary.
Usage
Double_Type mass = mass_function(Double_Type porb, asini[, error_porb, error_asini]);
or Double_Type[] i = mass_function(Double_Type mass, Double_Type[] Mopt[, error_mass, error_Mopt]; i);
Qualifiers
i - calculate and return the inclination (case 2a) Mx - neutron star mass in units of the solar mass (default: 1.4) rad - return the inclination angle in radiant (default: degrees) chatty - set to zero to suppress output messages
Description
The mass function of a binary is given by (see, e.g., Hilditch, 2001)
f(M) = (M_opt sin i)^3 / (M_x + M_opt)^2 = (4 Pi^2 (a sin i)^3) / (G P_orb^2)
with the masses of the neutron star, M_x, and its optical ompanion star, M_opt, the inclination, i, the semi-major axis, a, the orbital period, P_orb, and the gravitational constant, G.
There are two ways this function can be used:
- The orbital period (in days) and the semi-major axis (in lt-s) are provided, which then is used to calculate the right hand side of the above equation. The returned mass is given in units of the solar mass.
- The value of the mass function, mass, and the companion mass(es), Mopt, are provided and the i-qualifier is set. The returned value(s) is (are) the orbital inclination(s), i. In case of 2) the neutron star mass is assumed to be Mx = 1.4 solar masses by default, but can be changed by the Mx-qualifier. The thirs possible case, that the companion mass is returned, is formally solving a third order polynomial, which is yet not implemented here.
If uncertainties for the given parameters are provided the error propagated values are returned as well for the specific case via (value,error) = mass_funtion(..., error(s));
matrix2array
Synopsis
Transforms a matrix to an Array
Usage
Array_Type[] matrix2array( Any_Type[...,d,...] M, d );
Description
From the given dimension d of a matrix M an array is created, i.e., the given dimension of the matrix is shifted into the array which then contains matrices without this dimension.
Example
variable M = _reshape( [1:2*3*4], [2,3,4] ); variable A0 = matrix2array(M,0); variable A2 = matrix2array(M,2); vmessage("A0");print(A0); vmessage("A2");print(A2);
matrix33_as_array
Synopsis
return the contents of the matrix as a 3x3 matrix
Usage
arr=matrix33_as_array(m);
Description
Return the contents of the matrix as a 3x3 array
See also: vector, Vector_Type, matrix33_new
matrix33_determinant
Synopsis
return the determinant of a 3x3 matrix
Usage
det=matrix33_determinant(m);
Description
Calculates the determinant of a 3x3 matrix.
See also: vector, Vector_Type, Matrix33_Type
matrix33_diag
Synopsis
return a diagonal-matrix
Usage
Matrix33_Type=matrix33_diag(m11,m22,m33);
Description
Returns a 3x3 diagonal matrix. If m11,m22,m33 are given, then the three diagonal values are initialized to these three values. If only m11 is given, then a scalar matrix where all three elements are equal to m11 is returned.
See also: Vector_Type, Matrix33_Type
matrix33_identity
Synopsis
return a identity-matrix
Usage
Matrix33_Type=matrix33_identity();
Description
Returns a 3x3 identity matrix.
See also: Matrix33_Type
matrix33_new
Synopsis
instantiate a 3x3 matrix
Usage
Matrix33_Type=matrix33_new();
Description
Instantiates a new 3x3 matrix object. The following initializers are available: * No arguments: a zero-Matrix is returned * One argument: If the argument is of type Matrix33_Type: a copy of the argument is returned If the argument is a scalar: all matrix elements are initialized to this scalar (use matrix33_diag to initialize a diagonal matrix!) If the argument is an array with 9 elements: the matrix is initialized to these elements * Nine arguments: matrix elements m11,m12,m13,m21,m22,m23,m31,m32,m33, i.e., the elements are in row order and the matrix is:
m11 m12 m13 M = m21 m22 m23 m31 m32 m33
See also: Matrix33_Type, Vector_Type
matrix33_null
Synopsis
return a null-matrix
Usage
Matrix33_Type=matrix33_null();
Description
Returns a 3x3 null-matrix
See also: Matrix33_Type
matrix33_reflect
Synopsis
return a 3x3 reflection matrix to change the handedness of the ith axis
Usage
Matrix33_Type=matrix33_reflect(i);
Description
Returns a reflection matrix, i.e., a matrix that changes the handedness of the ith axis (where i=1: x-axis, i=2: y-axis, and i=3: z-axis) when multipliying it with a vector.
See also: vector, Vector_Type, matrix33_new
matrix33_rot
Synopsis
return a standard rotation matrix about the x-, y-, or z-axis
Usage
Matrix33_Type=matrix33_rot(i,angle;qualifiers);
Qualifiers
- deg: angle is given in deg [default: radians]
Description
Returns a rotation matrix to transform a column-3 vector from one cartesian coordinate system to another. The new coordinate system is given by rotating the original system in a counter clockwise way around the ith axis (where i=1: x-axis, i=2: y-axis, and i=3: z-axis)
The nomenclature follows Kaplan et al, The IAU Resolutions on Astronomical Reference Systems, Time Scales, and Earth Rotation Models, USNO circular 179, 2005.
See also: vector, Vector_Type, matrix33_new
matrix33_scalar
Synopsis
return a scalar-matrix
Usage
Matrix33_Type=matrix33_scala(m11);
Description
Returns a 3x3 scalar matrix (a matrix where all diagonal elements have the same value and where all other elements are zero)
See also: Vector_Type, Matrix33_Type
matrix33_transpose
Synopsis
return the transpose of a 3x3 matrix
Usage
Matrix33_Type=matrix33_transpose(m);
Description
Returns the transpose of a 3x3 matrix
See also: vector, Vector_Type, matrix33_new
\datatype{Matrix33_Type}
Synopsis
3x3 matrix type
Description
3x3 Matrix type that is compatible with Vector_Type
The following operations are defined for the Matrix33_Type: * M1+M2: Addition and subtraction of matrices * -M1: Changing sign of a matrix * M1*M2: Matrix - Matrix multiplication * M1*V: Matrix - Vector multiplication (V is a Vector_Type) * M1*f and f*M1: Matrix - real multiplication * M1/M2: Multiply M1 with the inverse of M2. Do NOT use this...
Objects of type Matrix33_Type can be instantiated with the following functions (see there for detailed descriptions): * matrix33_new: general initialization * matrix33_diag: return a diagonal matrix * matrix33_scalar: return a scalar matrix * matrix33_null: return a null matrix * matrix33_identity: return an identity matrix * matrix33_rot: return a rotation matrix for the x-, y-, or z-axis * matrix33_reflect: return a reflection matrix for the x-, y-, or z-axis
Other functions operating on matrices (functions marked with + are also available through accessor functions): * matrix33_determinant: calculate the determinant of the matrix (+) * matrix33_get_diag: return the diagonal elements of the matrix (+) * matrix33_get_trace: return the sum of the diagonal elements (+) * matrix33_transpose: return the transpose of the matrix (+) * matrix33_adjoint: return the adjoint of the matrix * matrix33_cofactors: return the matrix of cofactors
The Matrix33_Type object has the following accessors: * m.determinant(): return the determinant of the matrix * m.as_array(): return the elements of the matrix as a 3x3 array * m.transpose(): return the transpose of the matrix (not in place!) * m.diag(): return the diagonal elements as a vector * m.trace(): return the trace of the matrix * m.inverse(): return the inverse of the matrix (not in place!)
matrixmul
Synopsis
Matrix multiplication
Usage
M[n,m] = matrixmul( A[n,j], B[j,m] );
Description
(Efficient) Calculation of the product of two matrices A and B. If A is a (n x j)-Matrix and B a (j x m)-Matrix the resulting Matrix has the dimensions (n x m).
If A or B is one dimensional, the neccessary second dimension is assumed to be of length 1, i.e., A[n]*B[m] = A[n,1]*B[1,m] = M[n,m]
This function is identical to the ISIS intrinsic operator # or the inner product and only available for backwards compatibility!
See also: transpose, inner_product
max_time
Synopsis
finds the latest time in an array of structures with a time field
Usage
Double_Type tmax = max_time(Struct_Type structs[]);
Description
structs
is an array of structures containing the field time
.
tmax
is computed in the following way:
tmax = max( array_map(Double_Type, &max, array_map(Array_Type, &get_struct_field, [structs], "time")) );
See also: min_time, min_struct_field, max_struct_field
median_2d
Synopsis
computes the median image of a list of images (2d arrays)
Usage
med_img = median_2d(img1, ...);
Description
merge_struct_arrays
Synopsis
creates a structure whose fields are merged from all structures' fields
Usage
Struct_Type merged_s = merge_struct_arrays(Struct_Type s[]);
Description
All elements of s have to be structures with the same fields. The return value merged_s is another structure of this kind, and merged_s.field = [s[0].field, s[1].field, ..., s[-1].field]; holds for every field. (If s[i].field is NULL, it is skipped.)
Qualifiers
- remove_excess_fields: remove fields not present in all structures.
- reshape = Integer_Type dim: reshape the merged fields to the dimensions of the original fields, using the sum of dimensiom 'dim' to account for the increased array. Care has to be taken that the other dimension need to have the same length in all structures.
See also: append_struct_arrays, get_intersection, reshape
message_system
Usage
Integer_Type message_system(String_Type cmd)
Description
message(cmd);
system(cmd);
midas_bary_helio_corr
Usage
Struct_Type midas_bary_helio_corr(Y, m, d, H, M, S, RAh, RAm, RAs, decdeg, decmin, decsec)
Description
Y, m, d, H, M, S specify the date in UT.
The fields of the returned structure have the following meaning:
bary_corr_time: Barycentric correction time, in days
helio_corr_time: Heliocentric correction time, in days
bary_RV_corr: Total barycentric RV correction, in km/s
helio_RV_corr: Total heliocentric RV correction, in km/s
diurnal_RV_corr: (incl.) diurnal RV correction, in km/s
min_max
Synopsis
returns the minimum and maximum value of an array
Usage
(mn, mx) = min_max(Double_Type a[]);
or
(mn, mx) = min_max(Double_Type a1[], a2[], ...);
Qualifiers
- pad [=0]: additive fraction to be padded on both sides:
-dmin = dmax = (max-min) * pad
This qualifier may be overwritten by the following ones:
-
dmin [=0]: difference added to the minimum
-
dmax [=0]: difference added to the maximum
-
logpad [=0]: multiplicative fraction to be padded on both sides:
fmin^-1 = fmax = (max/min)^logpad
All array elements must be either positive or negative.
(For negative arrays -A
, [min_max(-A; logpad=lp)]
is the same as array_reverse([min_max(A; logpad=lp))
.)
This qualifier may be overwritten by the following ones:
- fmin [=1]: factor to multiply the minimum
- fmax [=1]: factor to multiply the maximum
Description
mn = min(a) * fmin + dmin;
mx = max(a) * fmax + dmax;
See also: min, max, moment
min_struct_field
Synopsis
returns the minimal field value of an array of structures
Usage
Double_Type mn = min_struct_field(Struct_Type structs[], String_Type fieldname);
Description
structs
is an array of structures containing the field called fieldname
.
mn
can be computed in the following way:
mn = min( array_map(Double_Type, &min, array_map(Array_Type, &get_struct_field, [structs], fieldname)) );
See also: max_struct_field, min_time, max_time
min_time
Synopsis
finds the earliest time in an array of structures with a time field
Usage
Double_Type tmin = min_time(Struct_Type structs[]);
Description
structs
is an array of structures containing the field time
.
tmin
is computed in the following way:
tmin = min( array_map(Double_Type, &min, array_map(Array_Type, &get_struct_field, [structs], "time")) );
See also: max_time, min_struct_field, max_struct_field
mix_rgb_colors
Synopsis
mix two rgb colors
Usage
rgb=mix_rgb_colors(rgb1,rgb2,fraction)
Description
This function mixes two colors in rgb space. Given the rgb values in the usual encoding in one 24bit integer, the color is mixed according to new rgb=rgb1*fraction+rgb2*(1-fraction) The operations are analoguous to the color mixing performed by the xcolor package of LaTeX (the operation is similar to LaTeX's color1!fraction!color2 syntax). The function returns the integer specifying the new rgb values.
See also: xfig_mix_color
mjd2jd
Synopsis
Convert Modified Julian Dates to Julian Dates.
Usage
mjd2jd(MJD);
Description
Helper function to convert MJD to Julian Dates by adding 2400000.5.
See also: jd2mjd
mkdir_rec
Synopsis
Create a new directory (recursively)
Usage
Int_Type mkdir_rec (String_Type dir [,Int_Type mode])
Description
Does the same as 'mkdir' with the exception that this function creates also subdirectories. See the help of 'mkdir' for a detailed description.
See also: mkdir
normalized_in_0_1
Synopsis
computes a normalized array with min=0, max=1
Usage
Double_Type a_norm[] = normalized_in_0_1(Double_Type a[]);
Description
a_norm = (a - min(a)) / (max(a)-min(a));
ntree
Synopsis
creates a new node of a tree with n children
Usage
Struct_Type ntree(DataType_Type objectType[, Struct_Type ntreeType]);
or Struct_Type ntree(DataType_Type objectType[, Integer_Type numChilds]);
Description
An n-tree is a special kind of a tree data structure. Each node has exactly n children, which are n-tree nodes themselves. The number of children is specified, e.g., by the 'numChilds' parameter, which is 2 by default. Initially, all children are null pointers. Most importantly, a node contains a list of objects, which is an array of the type specified by the 'objectType' parameter. Using the functions provided in the n-tree structure, objects can be inserted or removed.
The n-tree structure contains the following functions: insert - appends an object to the list. Usage: ntree.insert(objectType object); remove - removes an object from the list and returns 1 on success, 0 otherwise. Usage: Integer_Typer ntree.remove(objectType object); get - returns all objects stored in the node. Usage: objectType[] ntree.get(); getChilds - returns all objects stored in the children. Usage: objectType[] ntree.getChilds(); getAll - iterates over all nodes recursively to return all objects stored in the tree. Usage: objectType[] ntree.getAll(); childN - returns the N's child of the actual node. If the child does not exist, NULL will be returned. The function names can be overwritten by the type definition (see below). Usage: Struct_Type ntree.childN(); Qualifiers: new - creates the child if it does not exist renew - creates a new child no matter if it exists already delete - deletes the child
In addition, the following fields exist: objects - array of objects stored in the node childs - array of length N containing the children of the node type - structure describing the type of the n-tree (see below)
The type of the n-tree can be specified by the optional 'ntreeType' parameter. This structure defines additional behaviours of the n-tree, and the number of children and their the return function names eventually. The default type sets these names to the ones described above (childN). Additional informations about the usage of n-tree types can be found in the help of the 'ntree_type' function.
Example
% creates an 4-tree (called quad tree) with one % initial node and an object array of Integer_Type tree = ntree(Integer_Type, 4);
% insert '5' into the initial node tree.insert(5);
% create and insert '10' into the third child tree.child3(; new).insert(10);
% checks if child2 exists if (tree.child2() == NULL) message("child2 does not exist");
% insert '8' into child3, but it is not % created again since it exists already tree.child3(; new).insert(8);
% get all objects of the children % (will return [10,8]) print(tree.getChilds());
See also: ntree_type
ntree_type
Synopsis
returns a specific stucture defining an n-tree type
Usage
Struct_Type ntree_type(String_Type typeName);
Description
This functions returns the structure defining the specific n-tree type named 'typeName'. A list of all available types is shown if no name is passed. The 'help' qualifier shows additional information about the given type.
The type of an n-tree is specified and can be extended by a structure with the following fields:
numChilds Defines the number of children of each node. node. It has the same effect as the 'numChilds' parameter of the 'ntree' function.
childNames (optional) Array of strings defining the return function names for each children. The default names are "child1" to "childN".
new (optional) Reference to a callback function, which gets called right before a new node created by 'ntree()' is returned. Parameters passed are (1) the new node as a structur, (2) as which child number it is created (starting at 1, NULL if it is the root) and (3) the parent node it belongs to (NULL if there is none). The function has to return the new node.
insert (optional) Reference to a callback function, which gets called before an object is inserted into a node by 'node.insert()'. Parameters passed are (1) the node as a structure and (2) the object, which should be inserted. The function has to return the node. If NULL is returned, the object will not be inserted into the node.
remove (optional) Reference to a callback function, which gets called before an object is removed from a node by 'node.remove()'. Parameters passed are (1) the node as a structure and (2) the object, which should be removed. The function has to return the node. If NULL is returned, the object will not be removed from the node.
Example
% first define a callback function, which % prevents inserting negative numbers define ntree_myType_insert(node, object) { if (object < 0) return NULL; else return node; }
% structure defining a 4-tree (called quadtree) % and the children return functions are renamed % to match directions. The insert callback % function defined above is assigned also. variable myType = struct { numChilds = 4, childNames = ["left","right","top","bottom"], insert = &ntree_myType_insert };
% create a new n-tree of myType variable tree = ntree(Double_Type, myType);
See also: ntree
nutation_angles
Synopsis
calculate the nutation in longitude and obliquity for a given date
Usage
(dpsi,deps)=nutation_angles(JD;qualifiers)
Qualifiers
- JD_frac: fractional part of the Julian Date (for highest precision)
- theory: nutation theory to use. Possible string values are IAU2000A - full IAU2000A theory, i.e., an implementation of the theory of Matthews et al. (2002, J. Geophys. Res, 107, 2068) IAU2000B - truncated IAU2000A theory per McCarthy & Luzum (2003, Cel. Mech. Dyn. Astron, 85, 37); accuracy is 1 mas in 1995-2050 NU2000K - truncated IAU2000A theory from Kaplan (2009, USNO circular 181); accuracy is 0.1mas in 1700-2300 [the default]
- mjd: the time argument is in MJD, not in JD
- JD_frac: fractional part of the time argument. Use this if highest precision is needed. The nutation angles are then returned for the date JD+JD_frac.
- deg: return nutation angles in degrees
- arcsec: return nutation angles in arcseconds
- mas: return nutation angles in milliarcseconds
- iaufile: FITS file containing the theory coefficients (usually found in ISISSCRIPTS_REFPATH)
Description
This routine calculates the nutation angles for the forced nutation of the non-rigid earth according to one of the nutation theories described above for date given by the JD-argument (which is interpreted as TT). For most applications, usage of the NU2000K or IAU2000B routines is fully sufficient.
This routine is array safe (but JD_frac MUST be a scalar in this case!), the angles returned are in radians unless one of the deg, arcsec, or mas keywords is given.
The routine is based on the NOVAS 3.1 C-code (Bangert et al., 2011, see http://aa.usno.navy.mil/software/novas/novas_info.php )
See also: precess
nutation_matrix
Synopsis
Calculate the nutation matrix for the mean equinox JD
Usage
Matrix33_Type N=nutation_matrix(JD)
Qualifiers
- theory: :nutation theory to use, see help for function nutation_angles for details. Default: NU2000K
Description
This function returns the nutation matrix to convert coordinates for the mean equinox of JD to the true equinox of JD.
Note: This is mainly an internal function, it therefore does NOT have the usual qualifiers such as mjd, however, for some applications it is useful to have direct access to the nutation matrix
object_visibility_satellite
Synopsis
Calculate the time intervals when an object is observable with an astronomical satellite
Usage
visibility=object_visibility_satellite(ra,dec,mjdstart,mjdstop);
Qualifiers
- dt: Time resolution of the visibility calculation (days, default: 1)
- numerical: see below
Qualifiers to select the satellite. Constraints are ANDed if multiple satellite qualifiers are given, i.e., the joint visibility is returned
- swift : 46-180 deg from Sun; so far ignores lunar constraint (>23deg)
- suzaku : 70-110 deg from Sun
- xmm : 70-110 deg from Sun
- chandra: >46deg from Sun; ignores lunar constraint (>6deg)
- hst : >50deg from Sun
- nustar : no solar avoidance constraints(!)
Description
This routine calculates the time intervals during which an astronomical source with position ra, dec (both given in degrees) is visible to a given astronomical satellite. Most popular facilities are supported (see above). For each time between mjdstart and mjdstop the routine calculates the angular separation between the source and the Sun. By default it returns an array of structures giving the time intervals during which a source is visible. The tags of the structure describing each interval are MJDstart and MJDstop: start and stop MJD of the visibility interval MJDstring : human readable string of the above An empty array is returned if the object is not visible during the time interval specified by mjdstart and mjdstop.
If the qualifier "numerical" is set, the function returns a structure with the following tags mjd: MJD theta: separation betwee source and Sun (deg) visible: boolean visibility (0: no, 1: visible) instead of the array of time intervals
Note: if the visibility interval starts before mjdstart or extends beyond mjdend, the start and stop times of the search interval are shown instead. Note 2: several satellites also have moon avoidance angles, or avoidance angles for the Earth. Both are NOT (yet) taken into account.
Example
% convert source position into degrees variable ra=hms2deg(19,49,35.49); variable dec=dms2deg(+30,12,31.8);
% visibility during the year 2013 variable tstart=MJDofDate(2013,1,1); variable tstop=tstart+365;
variable vis=object_visibility_satellite(ra,dec,tstart,tstop;suzaku,dt=0.1);
variable i; for(i=0;i<length(vis);i++) { printf("%s
",vis[i].MJDstring); }
obliquity
Synopsis
calculate the obliquity of the ecliptic
Usage
eps=obliquity(JD;mjd);
Qualifiers
- mjd: if set, the argument is in MJD, not in JD
- laskar: use the polynomial expression by Laskar (1986, Astron. Astrophys. 157, 59), which is valid for a time span of around 10000 years around J2000, rather than the IAU 2006 expression.
- true: return the true obliquity for the date by applying nutation
- theory: nutation theory to use, see help for function nutation_angles; default: NU2000K
- deg: return the obliquity in degrees
- arcsec: return the obliquity in arcseconds
Description
This routine calculates the obliquity of the ecliptic for the Julian Date(s) JD using the IAU 2006 resolutions, which are based on Hilton et al., (2006, Cel. Mech. Dyn. Astron. 94, 351). JD can be an array, it is formally measured in TT.
The IAU 2006 expression is good to better than 1 arcsecond for +/-500 years around J2000. Outside of this interval use the polynomial expression by Laskar (laskar keyword).
The angles returned are in rad, use the deg or arcsec qualifiers if you want them in degrees or arc seconds.
See also: nutation_angles
orbitalphase
Synopsis
calculates an orbital phase at a given MJD time from an ephemeris
Usage
Double_Type phi = orbitalphase(Double_Type MJD, T0, P0, [Pdot, [Pddot ]]);
or
Double_Type phi = orbitalphase(Double_Type MJD, String_Type ephemeris);
Description
phi = (MJD-T0)/P - floor( (MJD-T0)/P );
where P = P0 + ( Pdot/2. + Pddot*MJD/6. )*MJD ;
The string ephemeris
can be used to refer to an internally
stored (T0, P) pair or, if available, (T0, P, Pdot, Pddot) quadruple.
If orbitalphase
is called with an unknown ephemeris
,
the allowed values are shown by get_ephemeris
.
Qualifiers
- numpulse: return pulse number instead of phase
See also: get_ephemeris
outer_product
Synopsis
computes the outer product of a number of vectors
Usage
Double_Type P[] = outer_product(Double_Type f1[], f2[], ..., fn[]);
Description
The n
arguments f1
, f2
, ..., fn
have to be one-dimensional arrays.
The return value P
is an n
-dimensional array with
P[i1, i2, ..., in] = f1[i1] * f2[i2] * ... fn[in]
.
outer_products
Synopsis
computes multi-dimensional arrays from one-dimensional ones
Usage
(Double_Type P1[], P2[], ...) = outer_products(Double_Type f1[], f2[], ..., fn[]);
Description
P1[i1, i2, ..., in] = f1[i1];
P2[i1, i2, ..., in] = f2[i2];
...
Pn[i1, i2, ..., in] = fn[i1];
The return value Pj
is calculated as outer product of x1, x2, ..., xn
,
where xj
= fj
and
xi = [1, 1, ..., 1]
(same length as fi
) for i!=j
.
See also: outer_product
pack_obj
Synopsis
converts an SLang object into a binary string
Usage
BString_Type[] pack_obj(Any_Type obj);
Description
Uses the `pack' function to convert a single SLang object into an array of binary strings (BString_Type). The format specifier is included as first character of the strings.
If the object is a single number or string then the returned array consists of a single item only. In case the given object is an array or a structure, the returned array starts with the type of the input array and its length or the field structure definition, respectively. The remaining items are the input array items or the structure items. For structures, the contained objects are converted recursively.
The following data-types are supported: Int16_Type, Int32_Type, Int64_Type, Float_Type, Double_Type, Char_Type, String_Type, Array_Type, Struct_Type and its aliases.
Example
s = pack_obj(PI); % s[0] = "d\030-DT\373!\011@"
s = pack_obj([1:3] + 4); % s[0] = "ak\003\000\000\000" % s[1] = "\005\000\000\000" % s[2] = "\006\000\000\000" % s[3] = "\007\000\000\000"
variable obj = struct { number = 598105, float = 341.12e8, more = struct { msg = "hello" } }; s = pack_obj(obj); print(s);
See also: unpack_obj, pack, unpack
parseDateString
Synopsis
parses a string including year, month, day [, hour [, minute [, second]]]
Usage
Struct_Type parseDateString(String_Type date_string);
Qualifiers
- no_2digit_correction: prevents the correction of two digit year specification, e.g., 83/12/13 instead of 1983/12/13. If this qualifier is not set, 2000 is added to years from 00-19 and 1900 to years 20-99.
Description
This function tries to interprete the given string as a data and returns
a structure consistent with the return value of localtime'. It is assumed that the first number occuring in the string specifies the year, the second the month, and the third the day. If further numbers are found, they are interpreted as hour, minute, second and fraction of seconds. Note that the seconds contained in the structure as defined in
localtime' has
to be an integer, while here a double is used in case of fraction of
seconds. There can be any separators between the numbers (except numbers).
Example
parseDateString("69/07/21"); parseDateString("1969-07-21 02:56:30.44"); parseDateString("sdf1969--$07..21QQ_02/56((30,,44xyz");
See also: localtime
period_search_via_string_length
Synopsis
A period-finding method for sparse randomly spaced observations
Usage
(Double_Type p[], SL[]) = period_search_via_string_length(Double_Type t[], s[], pmin, pmax; qualifiers)
Description
This function may be useful to determine the period of a variable signal 's' observed at a relatively small number 'N' of randomly spaced observations at times 't' possibly taken many periods apart.
To identify candidate periods in the given range between 'pmin' and 'pmax', string-lengths 'SL' for trial periods 'p' (sampled in equal frequency steps, see qualifier "delta_phi_max") are computed and tabulated. The trial period giving the lowest string-length results in the most smooth-looking curve and is thus the most likely period.
Notes
The underlying idea is to minimize the so-called "string-length" in a phase diagram. The phase diagram for a trial period 'p' is created by computing the phase 'phi' according to phi = (t-min(t))/p mod 1. The string-length 'SL' is then simply the sum of the lengths of line segments joining successive points (phi_i, s_i) in a phase diagram: SL = sum ( ((s_i-s_i-1)^2+(phi_i-phi_i-1)^2)^(1/2) )
- ((s_1-s_N)^2+(phi_1-phi_N+1)^2)^(1/2). To get rid of different scales caused by different units in ordinate (e.g., km/s or magnitude) and abscissa (unit = trial period), the observations are scaled so that s' = 0.5*(s-min(s))/(max(s)-min(s))-0.25 in order to give equal weight to measures and phases (for details, see Dworetsky 1983, MNRAS, 203, 917).
For reference: the string-length of a sinusoidal string is 1.4637.
Qualifiers
- delta_phi_max [=0.01]: Trial periods 'p' are sampled in equal frequency steps of delta_freq = delta_phi_max/(max(t)-min(t)). This ensures that phase errors are always less than 'delta_phi_max' since phi = (t-min(t))/p mod 1 = (t-min(t))*freq mod 1 -> delta_phi = (t-min(t))*delta_freq -> delta_phi_max = (max(t)-min(t))*delta_freq.
Example
t = 1000*urand(50); % time of observations s = sin(t)+0.05*grand(length(t)); % signal of observations; true period is 2*PI connect_points(0); point_style(2); plot(t, s); (p, SL) = period_search_via_string_length(t, s, 1, 10; delta_phi_max=0.01); connect_points(-1); point_style(-1); plot(p, SL); ind = where(SL==min(SL)); SL[ind[0]]; p[ind[0]]; % index, value, and period of minimum string-length connect_points(0); point_style(2); plot((t-min(t))/p[ind[0]] mod 1, s); % phase diagram
See also:
physical_quantity
Synopsis
initializes a physical quantity with number and units
Usage
PhysicalQuantity_Type physical_quantity(Double_Type number)
Qualifiers
- leng: length unit and dimension
- time: time unit and dimension
- mass: mass unit and dimension
- curr: electrical current unit and dimension
- temp: temperature unit and dimension
- unit: array of compound units (see
physical_quantity_unit
)
Description
The units are specified as strings with their dimensionality indicated by "^" as powers, i.e., "[unit]^[dim]". [dim] can be a negative number as well. For dim==1, "^1" can be omitted.
Example
variable c = physical_quantity(299792.458; length="km", time="s^-1");
See also: physical_quantity_unit
physical_quantity_current_in_A
Synopsis
defines a new current unit or gets the conversion for an existing one
Usage
physical_quantity_current_in_A(String_Type name, Double_Type value);
or
Double_Type value = physical_quantity_current_in_A(String_Type name);
Description
The unit name
= value
A can be used in
physical_quantity(x; curr=name); % = x * value
A
See also: physical_quantity
physical_quantity_length_in_m
Synopsis
defines a new length unit or gets the conversion for an existing one
Usage
physical_quantity_length_in_m(String_Type name, Double_Type value);
or
Double_Type value = physical_quantity_length_in_m(String_Type name);
Description
The unit name
= value
m can be used in
physical_quantity(x; leng=name); % = x * value
m
See also: physical_quantity
physical_quantity_mass_in_kg
Synopsis
defines a new mass unit or gets the conversion for an existing one
Usage
physical_quantity_mass_in_kg(String_Type name, Double_Type value);
or
Double_Type value = physical_quantity_mass_in_kg(String_Type name);
Description
The unit name
= value
kg can be used in
physical_quantity(x; mass=name); % = x * value
kg
See also: physical_quantity
physical_quantity_temperature_in_K
Synopsis
defines a new temperature unit or gets the conversion for an existing one
Usage
physical_quantity_temperature_in_K(String_Type name, Double_Type value);
or
Double_Type value = physical_quantity_temperature_in_K(String_Type name);
Description
The unit name
= value
K can be used in
physical_quantity(x; curr=name); % = x * value
K
See also: physical_quantity
physical_quantity_time_in_s
Synopsis
defines a new time unit or gets the conversion for an existing one
Usage
physical_quantity_time_in_s(String_Type name, Double_Type value);
or
Double_Type value = physical_quantity_time_in_s(String_Type name);
Description
The unit name
= value
s can be used in
physical_quantity(x; time=name); % = x * value
s
See also: physical_quantity
physical_quantity_unit
Synopsis
defines a new compound unit or gets the value of an existing one
Usage
physical_quantity_unit(String_Type name, Double_Type value;; qualifiers);
or
physical_quantity_unit(String_Type name, PhysicalQuantity_Type value);
or
PhysicalQuantity_Type value = physical_quantity_unit(String_Type name);
Description
If the second argument of the physical_quantity_unit
function
is a double value, the value of the unit is
constructed with the physical_quantity
function,
i.e., all its qualifiers can be applied.
The unit name
= value
can be used in
physical_quantity(x; unit=name); % = x * value
See also: physical_quantity
planetpos
Synopsis
Calculate the apparent position of the Sun and the planets (Mercury,..., Neptune).
Usage
(r,alp,del) = planetpos(JD;qualifiers);
or
```c
vector = planetpos(JD;qualifiers);
##### Qualifiers
* object: a string (valid values: Sun, Mercury, Venus,
Mars, Jupiter, Saturn, Uranus, Neptune; default:
Sun).
* of_date: if set, return the apparent position in the geocentric
celestial reference system for the ecliptic and
equinox of date (this is what is listed in the
astronomical almanac).
* vector: return a geocentric vector for the position.
* deg: return the right ascencion alp and the declination del
in degrees (default: radian).
* mjd: the given date is the MJD
##### Description
Returns the distance (in AU), and apparent right ascension and declination
of a planet in the geocentric celestial reference system for the
date JD. The default coordinate systems is the J2000.0 system, if you
want to include the effects of precession and nutation, use the
of_date qualifier.
The JD is to be given in TT.
#### png_read_curve
##### Synopsis
reads a curve from a x-y plot in an image
##### Usage
```c
(X, Y) = png_read_curve(filename, xpix1, xval1, xpix2, xval2, ypix1, yval1, ypix2, yval2);
Description
filename
is a png file containing a color-defined curve.
Its x
and y
coordinates are calibrated by pix
/val
pairs
of pixel coordinate and corresponding value assuming a linear scale.
The pixel coordinates start with (0, 0) in the upper left corner.
The return values are calibrated x-values and y-values
averaged over all pixels of the curve in the corresponding column.
Qualifiers
- color [
=0x000000
(black)]: considered color of the curve - wherenot: consider any color except the one specified above
point_in_polygon
Usage
ret=point_in_polygon(p0,V);
or
```c
ret=point_in_polygon(x,y,Vx,Vy);
##### Synopsis
determine whether a point is in a polygon
##### Qualifiers
* evenodd: use the crossing number method
* crossing: use the crossing number method
* winding: use the winding number method (the default)
##### Description
The function returns 1 if the point p0 is located inside the
polygon defined by the vertices V, and 0 if it is located
outside of the polygon.
The polygon has to be closed, i.e. V.x[n]==V.x[0] and
V.y[n]==V.y[0] where n is the number of polygon points.
Either the winding number (default) or the even-odd-rule
define what is meant by inside.
The point is either defined by a struct{x,y} and the vertices
by a struct{x[],y[]}, or the coordinates can be directly
given in the respective arrays.
See the URL below for more explanations.
Based on code by Dan Sunday,
http://geomalgorithms.com/a03-_inclusion.html
__See also__: crossing_number_polygon,point_in_polygon,simplify_polygon
#### point_is_left_of_line
##### Usage
```c
ret=point_is_left_of_line(p0,p1,p2);
or
```c
ret=point_is_left_of_line(p0x,p0y,p1x,p1y,p2x,p2y);
##### Synopsis
determines whether point p2 is left of a line through p0 and p1
##### Description
This function tests if point p2 is to the left of an infinite line defined
by points p0 and p1. The points are defined by structs p=struct{x,y},
where x is the x-coordinate and y is the y-coordinate.
The function returns
+1 if P2 is left of the line
0 if P2 is on the line
-1 -1 if P2 is right of the line
Based on code by Dan Sunday, http://geomalgorithms.com/a03-_inclusion.html
__See also__: crossing_number_polygon,winding_number_polygon,point_in_polygon,simplify_polygon
#### polint
##### Synopsis
Polynomial Interpolation and Extrapolation with error estimation
##### Usage
```c
(Double_Type y, dy ) = polint (Double_Type[] xa, ya, Double_Type x);
Description
Returns an interpolated value y
at the given point x
and an
error estimate dy
for the given arrays xa
and ya
.
If P(x) is the polynomial of degree n-1 such that P(xa_i) = ya_i, i=0,n-1,
then the retruned value y=P(x).
Adapted from algorithm in Numerical Recipes,
by Press et al. (1992, 2nd edition), Section 3.1.
See also: qromb
polytropic_standard_model
Synopsis
Compute the structure of a polytropic standard model
Usage
Struct_Type s = polytropic_standard_model(Double_Type M, R, X, Y)
Description
Based on the stellar mass 'M' (in solar masses), stellar radius 'R' (in solar radii), hydrogen mass fraction 'X', and helium mass fraction 'Y', this function computes the structure of a polytropic standard model. The output structure 's' contains the fields "r" (radial coordinate in solar radii), "rho" (density in g per cm^3), "P" (pressure in dynes per cm^2), and "T" (temperature in K), i.e., the structure of the polytrope.
Notes
Polytropes are gaseous spheres in hydrostatic equilibrium which obey a polytropic equation of state P = K rho^((n+1)/n) where 'P' is the pressure, 'rho' the density, 'K' the polytropic constant, and 'n' the polytropic index. The so-called polytropic standard model is a polytropic model of index n=3. Its structure has a surprisingly good resemblence to those of main sequence stars and is, thus, used as a starting point to study the internal structure of stars. For details, see, e.g., Chapter 2.4 in the book "Principles of Stellar Evolution and Nucleosynthesis" by D.D. Clayton.
Example
s = polytropic_standard_model(1, 1, 0.73, 0.26);
See also: solve_Lane_Emden_equation, solve_Eddington_quartic_equation
pos_modulo
Synopsis
computes the modulo and ensures that it is positive
Usage
m = pos_modulo(a, b);
Description
a
and b
can either be arrays or scalars of a numerical type.
m = a mod b; %
if m is positive, otherwise:
m = (a mod b) + b;
precess
Synopsis
precess equatorial coordinates between two equinoxes
Usage
(alpha,delta)=precess(alpha,delta;qualifiers)
or
```c
AstroVector3_Type ecl=precess(eqp;qualifiers)
##### Qualifiers
* fromequinox: starting equinox of the transformation.
Float: JD, string: epoch designation (e.g.,
"J2000.0" or "B1950.0";
default: J2000.0)
* toequinox: end equinox of the transformation (same interpretation)
as fromequinox; required!
* deg: interpret angular arguments in degrees (default is radians!)
applies also to the return value.
* mjd: interpret equinoxes as a MJD (default: JD)
##### Description
Precess a coordinate from one equinox to another using
the IAU 2000A model.
The default equinox is J2000.0, formally the equinoxes are in TDB.
__See also__: Vector_Type,JDofEpoch,dms2deg,hms2deg,precession_matrix,gcrs2j2000_matrix
#### precession_matrix
##### Synopsis
return the precession matrix for the transform from J2000.0 to a given JD
##### Usage
```c
Matrix33_Type mat=precession_matrix(JD)
Description
This function calculates the precession matrix for the calculation of the precession FROM J2000.0 TO JD using the IAU 2000.0A theory (Capitaine et al., 2003, A&A 412, 567, Sect. 7)
To get the matrix for the conversion from JD to J2000.0, transpose the matrix returned by this function.
In most cases users will want to use the precess-function to precess astronomical coordinates.
Note: This is mainly an internal function, it therefore does NOT have the usual qualifiers such as mjd or any safety checks, however, for some applications it will be useful to have direct access to the precession matrix.
See also: precess,precession_matrix
precession_matrix2
Synopsis
return the precession matrix for the transform between two equinoxes
Usage
Matrix33_Type mat=precession_matrix2(fromJD,toJD)
Description
This function calculates the precession matrix for the calculation of the precession from fromJD TO toJD using the IAU 2000.0A theory (Capitaine et al., 2003, A&A 412, 567, Sect. 7)
In most cases users will want to use the precess-function to precess astronomical coordinates.
Note: This is mainly an internal function, it therefore does NOT have the usual qualifiers such as mjd or any safety checks, however, for some applications it will be useful to have direct access to the precession matrix.
See also: precess,precession_matrix
primefactors
Synopsis
factorizes an integer number into primes
Usage
Integer_Type factors[] = primefactors(Integer_Type x);
Description
The array of prime factors will be ordered:
factor[i] <= factor[j] %
for i<j
.
print_array
Synopsis
Prints the entries of an array in one line.
Usage
print_array(Array_Type);
Description
Prints the entries of an array in one line instead of different lines. Helpful for displaying dependent values of different arrays like a table. EXAMPLE UVOT-filter names and their wavelengths.
isis> filters=["UVW2","UVM2","UVW1","U","B","V"]; isis> wave=[2078.546871,2257.186625,2659.009185,3475.482080,4359.034280,5429.548507]; isis> print_array(filters); print_array(wave); UVW2 UVM2 UVW1 U B V 2078.546871 2257.186625 2659.009185 3475.482080 4359.034280 5429.548507
print_struct
Synopsis
prints the fields of a structure as columns of a table
Usage
print_struct([File_Type F,] Struct_Type s);
Description
If s
is a structure of arrays or lists,
these are displayed as columns of a table.
The output is written to stdout
unless another F
is specified, which may either be a file pointer
or a string containing the filename.
Qualifiers
- i: array of rows which are to be shown at all
- mark: array of rows which are to be marked
- fields: array of fieldnames (columns) which are to be shown
- fmt: format string(s) for the columns, see
help("sprintf");
- sep: string separator between the columns (default =
" "
) - initial: initial separator on a line (default =
""
) - final: final separator on a line (default =
""
) - html: set
sep
,initial
andfinal
such that an HTML table is produced - tex: set
sep
,initial
andfinal
such that a TeX table is produced - nohead: don't show head line with field names
Examples
print_struct( s); %
The fields of s
are printed.
print_struct(F, s); %
The fields of s
are printed to the file F
.
print_struct( , s; fmt="%.2f"); %
The fields are floats which are printed with 2 decimals.
print_struct( , s; fmt="%.2f", sep=" "); %
As before, but with smaller separation between the columns.
See also: writecol, sprintf
properties_satellite_galaxies
Synopsis
Retrieve the properties of several satellite galaxies of the Milky Way
Usage
Struct_Type s = properties_satellite_galaxies()
Description
Retrieve the properties (current kinematics, spatial extent, mass) of several satellite galaxies of the Milky Way.
The output structure contains for each object the celestial coordinates in right ascension [h, m, s] and declination [deg, arcmin, arcsec], distance [kpc], radial velocity [km/s], proper motion in right ascension times cosine of declination [mas/yr], proper motion in declination [mas/yr]), half-light radius [kpc], dynamcial mass inside half-light radius [solar masses], Plummer softening radius [kpc], and the Plummer mass [solar masses].
Notes
Positions, distances, radial velocities, half-light radii, and dynamcial masses are from McConnachie 2012 when available. Proper motions (PMs) are available only for a subset of satellites and are set to zero otherwise.
The parameters of the Plummer potentials (softening radius, mass) are determined under the assumption that the dynamcial mass inside the half-light radius is 0.85 times the total mass, which gives a softening radius of about one third of the half-light radius.
Main reference: McConnachie 2012, AJ, 144, 4. Additional references or comments:
- LMC, SMC: PMs and Plummer parameters from Kallivayalil et al. 2013, ApJ, 764, 161. Half-light radius and mass derived from Plummer parameters.
- Canis Major: PMs from Dinescu et al. 2005, AJ, 631, L49. Dynamical mass assumed to be ten times stellar mass, half-light radius set to 1 kpc.
- Sagittarius dSph: PMs from Pryor et al. 2010, AJ, 139, 839.
- Fornax: PMs from Piatek et al. 2007, AJ, 133, 818.
- Sextans: PMs from Walker et al. 2008, ApJ, 688, L75.
- Sculptor: PMs from Piatek et al. 2006, AJ, 131, 1445.
- Draco: PMs from Pryor et al. 2014, AJ, submitted [arXiv:1407.3509].
- Ursa Minor: PMs from Piatek et al. 2005, AJ, 130, 95.
- Carina: PMs from Piatek et al. 2003, AJ, 126, 2346.
- Leo I: PMs from Sohn et al. 2013, ApJ, 768, 139.
- Leo II: PMs from Lepine et al. 2011, ApJ, 741, 100.
- Bootes III: Dynamical mass assumed to be ten times stellar mass, half-light radius set to 10 pc.
Example
s = properties_satellite_galaxies(); print_struct(s);
See also: N_body_simulation_MW_kernel
pushStructFieldArray
Synopsis
appends a value to an array which is the field of structure
Usage
pushStructFieldArray(s, field, value);
Description
s.field = [s.field, value]; %
if s.field
was not NULL
qromb
Synopsis
Integrate using Rombergs's rule to specified accuracy.
Usage
Double_Type int = qromb (Ref_Type function, Double_Type min, Double_Type max);
Description
Integrate a function between the limits min
and max
to specified
accuracy using the extended trapezoidal rule. Adapted from algorithm
in Numerical Recipes, by Press et al. (1992, 2nd edition), Section 4.3.
The precision and number of maximal iterations can be set via qualifiers.
Any other keywords are passed directly to the user-supplied function.
NOTE: Romberg is more efficient then Simpson.
Qualifiers
-
qromb_eps [=1e-6]: Scalar specifying the fractional accuracy before ending the iteration.
-
qromb_max_iter [=16]: Integer specifying the total number iterations at which
qsimp
will terminate even if the specified accuracy has not yet been met. The maximum number of function evaluations is 2^(qromb_max_iter-2)-1. -
k [=5]: Integration is performed by Romberg’s method of order 2k, where, e.g., k=2 is Simpson’s rule.
Example
% Compute the integral of sin(x) from 0 to PI/3. % The value obtained should be cos(PI/3) = 0.5 variable val = qromb( &sin, 0, PI/3);
See also: integrate_trapez,polint,qsimp
qsimp
Synopsis
Integrate using Simpson's rule to specified accuracy.
Usage
Double_Type int = qsimp (Ref_Type function, Double_Type min, Double_Type max);
Description
Integrate a function between the limits min
and max
to specified
accuracy using the extended trapezoidal rule. Adapted from algorithm
in Numerical Recipes, by Press et al. (1992, 2nd edition), Section 4.2.
The precision and number of maximal iterations can be set via qualifiers.
Any other keywords are passed directly to the user-supplied function.
Qualifiers
- qsimp_eps [=1e-6]: Scalar specifying the fractional accuracy before ending the iteration.
- qsimp_max_iter [=16]: Integer specifying the total number iterations
at which
qsimp
will terminate even if the specified accuracy has not yet been met. The maximum number of function evaluations is 2^(qsimp_max_iter+5)-1.
Example
% Compute the integral of sin(x) from 0 to PI/3. % The value obtained should be cos(PI/3) = 0.5 variable val = qsimp( &sin, 0, PI/3);
See also: integrate_trapez,qromb
quantile
Synopsis
gets an arbitrary quantile of an array
Usage
Any_Type quantile(Double_Type p, Any_Type a[]);
Description
quantile(0, a) = min(a)
and quantile(1, a) = max(a)
.
For values 0 < p < 1
, intermediate values of the a
will be returned.
See also: median
rad2RD
Synopsis
Convert right ascension and declination from radians to (h, m, s) and (deg, arcmin, arcsec)
Usage
rad2RD(Double_Types RA[], D[])
Description
Given right ascension and declination in radians, the function converts them to hours, minutes, seconds and degrees, minutes of arc, seconds of arc. In case of negative declinations, the minus sign is assigned only to the highest non-vanishing term, i.e., to degrees if degrees is nonzero, to minutes if minutes is nonzero but degrees is zero, ... . Negative right ascensions are not expected.
Example
(ah, am, as, dd, dm, ds) = rad2RD(0,-PI/4.); (ah, am, as, dd, dm, ds) = rad2RD(0,-PI/181.); (ah, am, as, dd, dm, ds) = rad2RD(3/4.*PI,0); (ah, am, as, dd, dm, ds) = rad2RD([0,3/4.*PI],[-PI/4,0]); (ah, am, as, dd, dm, ds) = rad2RD(RD2rad(01, 45, 12.54, 87, 55, 45.34));
See also: RD2rad
radial_velocity_LSR
Synopsis
computes the radial velocity component of the local standard of rest
Usage
Double_Type radial_velocity_LSR(Doule_Type RA, dec, MJD)
or
Double_Type radial_velocity_LSR(Doule_Type az, el, MJD, longw, lat)
Description
right ascension and declination, Modified Julian Date azimuth, elevation, geographic longitude (west), geographic lattitude
radial_velocity_binary
Synopsis
computes the radial velocity of a binary system as a function of orbital phase or time
Usage
Double_Type v = radial_velocity_binary(Double_Type x);
Qualifiers
- v0: systemic velocity
- K: velocity semi amplitude, in km/s
- asini: value of a * sin i, in kilometers
- P: orbital period, in days
- T0: epoch of periastron, x means time in days
- T90: epoch of mean longitude 90 degr, x means time in days
- e: eccentricity (0 <= e < 1)
- omega: longitude of periastron in radian (unless degrees is set), default = 0
- degrees: omega is measured in degrees instead of radian
Description
x has the meaning of orbital phase (0 <= x < 1) unless T0 or T90 is specified.
See also: Ch. 3 of R.W. Hilditch, An Introduction to Close Binary Stars, Cambridge Univ. Press, 2001
radius_to_unit
Synopsis
converts a radius from a given unit into another
Usage
Double_Type = radius_to_unit(Double_Type radius, String_Type from_unit, to_unit);
Qualifiers
asini - projected semi major axis (lt-s) i - inclination (degrees) mq - mass ratio of star to companion
Description
Converts the 'radius' from a given unit 'from_unit' into another unit 'to_unit'. There might be additional information needed for the calculation, e.g. the projected semi major axis, which is passed via qualifiers.
The following units are supported (with additional needed qualifiers): rsun - solar radii cm - centimeters lts - light seconds disp - binary displacement (asini, i, mq) fill - fill factor = radius/critical radius (asini, i, mq)
Note: It might be possible, that less qualifiers may be passed, if the variables are canceled during the calculation. For example, converting from binary displacement to fill factor requires the mass ratio 'mq' only.
Note: Let the projected semi major axis be the distance from CM to M1. Then, mq is defined as M1/M2. That means, to convert the radius of the companion, mq has to be inverted and asini has to be the distance from CM to M2, which is simply asini*mq[not inverted]
See also: Roche_potential, Roche_critical
ratio_error_prop
Synopsis
calculates a ratio and error propagation
Usage
(rat, err) = ratio_error_prop(a, a_err, b, b_err);
or
(rat, err) = ratio_error_prop(a, b; Poisson);
Description
rat = a/b
err = sqrt[ (1/b * a_err)^2 + (a/b^2 * b_err)^2 ]
Qualifiers
- Poisson:
a_err = sqrt(a); b_err = sqrt(b);
See also: hardnessratio_error_prop
rectangle_where
Synopsis
finds the rectangle in a 2d array, where an expression is true
Usage
(Integer_Type Y[], X[]) = rectangle_where(Integer_Type expr);
Description
The rectangle in the 2d array defined by the index-arrays
Y
and X
contains all elements for which expr!=0
,
i.e., one can crop the elements for which expr==0
.
(Note that all(expr[Y,X])
is not necessarily true.)
Example
img = img[rectangle_where(img>0)];
reduce_struct
Synopsis
remove one or more fields from a strucure
Usage
Struct_Type reduced_struct = reduce_struct(Struct_Type s, String_Type fieldsnames[]);
Qualifiers
- extract:If given, the returned struct does only contain the fields 'fieldnames'!
Description
Either removing the fields 'fieldnames' from the given structure 's' (if they even exist) or if the qualifier 'extract' is given removing all other fields but those given with 'fieldnames'!
refraction
Synopsis
calculate the correction for astronomical refraction
Usage
Double_Type refraction(z0;qualifiers)
Qualifiers
- lambda: wavelength (in A, between 3000 and 17000 A)
- mum: wavelength is in microns
- nm: wavelength is in nanometers
- temperature: temperature at observer [K; default: 288.15K=15C]
- lapse_rate: temperature lapse rate [K/m; default: 0.0065K/m = 6.5K/km]
- centigrade: temperature is given in C
- pressure: ambient total pressure at observer (Pa, default: 1013.25kPa)
- kPa: ambient pressure is in kPa
- hPa: ambient pressure is in hPa (or mbar)
- rel_humidity: relative humidity at observer (between 0 and 1)
- altitude: altitude of observer above geoid (m; below 11000m)
- latitude: geographical latitude of the observer (rad; default: 0)
- deg: all angles are given in degrees, not radians
- exact: use numerical integration also for z0 below 80deg
Description
For a given zenith distance z and local observing conditions, this function calculates the refraction angle, that is the difference R=z-z0 where z is the unrefracted zenith distance that would be measured if the Earth did not have an atmosphere, and z0 is the observed zenith distance.
The most common use of this function will be to calculate z0 for a given topocentric zenith distance. For all practical purposes, R is so small and changes slowly enough, such that users can call the function with z and determine z0 from the return value. Note that the default arguments are in radians, use the deg qualifier to switch to degrees.
The refraction depends very slightly on the atmospheric properties, i.e., its temperature profile. These can be set with the respective qualifiers.
The default settings of the routine are that for zenith distances smaller than 80 degrees the approximations given by Saastamoinen (1972, Bull. Geodesique 105, 279 and 1972, Bull. Geodesique 106, 383) are used.
For larger zenith distances, or if the "exact"-qualifier is set, a numerical approach is chosen, following the approach discussed by Auer & Standish (2000, AJ 119, 2472 [first submitted in 1979!]), Hohenkerk and Sinclair (1985, HM Nautical Almanac Office, NAO Technical Note 63), and Mangum and Wallace (2015, PASP 127, 74). The function uses the atmospheric model discussed by Hohenkerk and Sinclair, but uses the exact refraction formula for air. This exact approach yields good results for zenith distances up to a few degrees larger than 90 degrees (i.e., observation of the horizon from a mountain top) and is the one on which the refraction formulae of the Astronomical Almanac are based.
For large zenith angles, there are slight differences at the arcsecond or less level between the results discussed here and the numbers listed by Hohenkerk or Auer. These are due to the different treatment of numerical instabilities and the use of numerical differentiation formulae. Given that for large zenith angles the simplified atmospheric model of Hohenkerk (constant temperature lapse rate in the troposphere, constant temperature in the stratosphere) results in larger systematic errors anyway (see Nauenberg, 2017, PASP 129, 44503), this should not be seen as an error of the function.
As a caveat, the calculations using the Saastamoinen formulae assume that all pressure terms there are in Pa, rather than in mbar. This reproduces the results with respect to the numerical simulations and values tabulated elsewhere to <1". But it is not what Saastamoinen claims. I (J. Wilms) am puzzled...
If the exact qualifier is set or for large z, all input into the function MUST be scalars - in this case this routine is NOT array safe.
For the Saastamoinen formulae, the routine is array safe in all relevant parameters.
refractive_index_air
Synopsis
calculate the refractive index of dry and moist air
Usage
Double_Type refractive_index_air(lambda;qualifiers)
Qualifiers
- mum: wavelength is in microns
- nm: wavelength is in nanometers
- temperature: temperature [K; default: 288.15K]
- centigrade: temperature is given in C
- pressure: ambient total pressure (Pa, default: 1013.25kPa)
- kPa: all pressure arguments are in kPa
- hPa: all pressure arguments are in hPa (or mbar)
- CO2_ppm: CO2 fraction in ppm (default: 450)
- water_pressure: partial water vapor pressure
- rel_humidity: relative humidity
- silent: do not emit warning messages
Description
This function calculates the phase refractive index of dry and moist air for light in the optical and IR following the standard paper by Ciddor (1996, Appl. Optics 35(9), 1566) and the discussion by Stone and Zimmerman (NIST Engineering Metrology Toolbox; http://emtoolbox.nist.gov/Wavelength/Documentation.asp) for a given wavelength (default A, but see the mum and nm qualifiers) and ambient conditions (temperature, pressure, humidity or partial vapor pressure, and CO2 concentration). The function has been verified against the values given by Ciddor and in the NIST Engineering Metrology Toolbox.
The relative uncertainty of the approximations used here is claimed to be around 2e-8. The range of validity is 3000 A<lambda<17000 A, temperatures between -40C and 100C (233-373K), and pressures between 10 and 140kPa. The CO2 fraction is allowed to vary between 0 and 2000ppm.
This routine is array safe. If lambda is an array, the other qualifiers can be either arrays of the same length as lambda or single valued.
rename_struct_fields
Usage
Struct_Type new_s = rename_struct_fields(Struct_Type s, String_Type fieldnames[]);
replicate_table
Synopsis
repeats columns of a table, possibly with a periodic shift
Usage
Struct_Type replicate_table(Struct_Type t);
Qualifiers
- P [= 1]: period
- back [= 0]: repeat periods backwards
- ahead [= 1]: repeat periods forwards
- periodic [
= ["bin_lo", "bin_hi"]
]: periodic structure fields
Description
The return value is a structure with the same fields as t
.
All array fields are repeated (back+1+ahead) times; periodic ones as
[ val-back*P, ..., val-P, val, val+P, ..., val+P*ahead ]
,
while other ones are just replicated:
[ val , ..., val , val, val , ..., val ]
.
rescale_range
Synopsis
rescales a value
Usage
Double_Type y = rescale_range(Double_Type x);
Qualifiers
- in: inputrange
- out: outputrange
Description
possibilities for in/out and the corresponding scaling function:
"0:1"
, no scaling
"-inf:inf"
, arctan
rgb2gray
Synopsis
Weighted RGB to GRAY convertion
Usage
Integer_Type gray = rgb2gray( Integer_Type[] RGB );
Usage
Integer_Type gray = rgb2gray( Integer_Type[] R, G, B );
Description
Converts a RGB color, which is either encoded in one Integer or as single channels (R,G,B), into a weighted gray scale using:
GRAY = (0.3 * R) + (0.59 * G) + (0.11 * B)
The 8 bit gray value(s) are return as Integer_Type[] in the interval [0,255].
Qualifiers
- norm: return values are in Double_Type interval [0,1].
- rgb: return values as single RGB encoding: GRAY << 16 | GRAY << 8 | GRAY
See also: rgb2hsl, hex2rgb
rgb2hex
Synopsis
converts r, g, b values between 0 and 1 to a 24-bit integer
Usage
Integer_Type rgb2hex(Double_Type r, g, b)
Qualifiers
- str: return value is a string with preceding "#".
- des: input color will be desaturated to gray-scale values. If a value is given, the saturation will be scaled accordingly.
See also: rgb2hsl
rgb2hsl
Synopsis
converts (red, green, blue) values to (hue, saturation, lightness)
Usage
(Double_Type h, s, l) = rgb2hsl(Integer_Type rgb);
or
(Double_Type h, s, l) = rgb2hsl(Double_Type r, g, b);
Description
In the first usage, rgb is an (array of) 24 bit RGB values. In the second form, r, g, and b are values between 0. and 1.
The return values are always normalized to [0, 1].
See hsl2rgb
for a definition of hue, saturation, lightness.
See also: hsl2rgb
rgb2r_g_b
Synopsis
converts a ( rgb ) color to (red, green, blue)
Usage
Integer_Type r,g,b = rgb2r_g_b( Integer_Type[] );
Description
Converts a color encoded in one Integer into its single r, g, b values. The return value is a(n array of) 24 bit RGB value(s) in Integer_Type interval [0,255].
Qualifiers
- norm: return values are in Double_Type interval [0,1].
- str: If qualifier is given ('str'=NULL) r,g,b will be retruned as strings. If 'str' != NULL, e.g., 'str'=',', than the singe color channels are combined and returned: r + str + g + str + b
See also: rgb2hsl, hex2rgb
rgb_interpol_colmap
Synopsis
interpolate between multiple rgb colors
Usage
col=rgb_interpol_colmap(rgbind, rgbcol, nintcol)
Description
This function builds a color array by interpolating between multiple colors. The resulting array has nintcol colors. In this array, the colors rgbcol[ii] will be placed at col[rgbind[ii]]. The interpolation between the individual colrs is linear. The input color indices rgbind[*] have to be sorted in ascending order. The result can be assigned as a color map directly using png_add_colormap.
See also: mix_rgb_colors
roman
Synopsis
replaces an integer with lowercase Roman numeral strings
Usage
String_Type rom = roman(Integer_Type);
See also: Roman
round2
Synopsis
Round to the nearest integral value or to given digit
Usage
Double_Type[] = round2( Double_Type[] value);
or
```c
Double_Type[] = round2( Double_Type[] value, Integer_Type digit );
##### Description
This function rounds its argument to the nearest integral value and
returns it as a floating point result. If the argument is an array,
an array of the corresponding values will be returned.
If a 2nd argument is given it is used as digit the value is supposed
to be rounded to.
__See also__: round, floor2, ceil2, nint
#### round_conf
##### Synopsis
converts confidence intervals after DIN 1333 and gives the rounded decimal place.
##### Usage
```c
Struct_Type round_conf(Double_Type conf_lo, conf_val, conf_hi);
or
```c
Struct_Type round_conf(Double_Type val, sym_err);
##### Qualifiers
##### Description
<code>conf_lo</code> is the lower confidence limit
<code>conf_val</code> is the best fit
<code>conf_hi</code> is the upper confidence limit
The alternative usage with two arguemts allows to specify a
value <code>val</code> with a symmetric uncertainty <code>sym_err</code>.
The returned structure contains the fields
- <code>err_lo</code> (the rounded lower error)
- <code>value_val</code> (the rounded best fit)
- <code>err_hi</code> (the rounded upper error)
- <code>digit</code> (the rounded decimal place of the error)
__See also__: round_err, TeX_value_pm_error
#### round_err
##### Synopsis
rounds an error after DIN 1333 and gives the rounded decimal place
##### Usage
```c
Struct_Type round_err(Double_Type x);
Qualifiers
digits: decimal place where rounding will be applied (suspends DIN 1333)
Description
x
is the error which will be rounded
The returned structure contains the fields
-
value
(the rounded error) -
digit
(the rounded decimal place of the error)
EXAMPLES round_err(0.1278) will return value=0.13 and digit=-2 round_err(1278) will return value=1300 and digit= 2 since after DIN 1333, the decimal point may only be right after the rounded digit or left of it, the error should be in LaTeX assigned as $13.0\times 10^{2}$ (? comment: in order to have the correct number of significant digits shouldn't it be: $0.13\times 10^{4}$, $1.3\times 10^{3}$, or $13\times 10^{2}$)
See also: round_conf
save_slang_variable
Synopsis
allows to save S-Lang variables into a file
Usage
save_slang_variable(file, &var1, &var2, ...);
Qualifiers
edit - open an editor to modify the variables; in that case all variables have to be passed as references and will be set to their new values after the editor is closed delete - delete the file after editing (requires the 'edit' qualifier to be set); note that the given variables are still modified!
Description
The S-Lang code defining the given variables is saved into a file, specified by either the filename or an already opened file-pointer. In order to handle arrays with a large number (>1000) of items as well as complex structures, the S-Lang code uses temporary variable names. Their values are assigned step by step to avoid a stack overflow. The file can be evaluated later to push the saved variables onto the stack (see the example).
This function allows to modify the given variables as well. In that case the 'edit'-qualifier has to be set and all passed variables have to be given as references. The file the variables are saved into is shown in the editor specified by the EDITOR environment variable or jed, if EDITOR is undefined. After saving the file and closing the editor, the file is evaluated, which should push the (modified) S-Lang objects onto the stack. These objects are finally assigned to the given variables.
NOTE: the latter feature is based on the function 'edit_var', which does the same except that the main purpose is to edit the variables using a temporary file. Here, the S-lang code is human readable as well, since no temporary variables are used to assign the values stepwise. In that case, however, stack overflow errors may occure. But who wants to edit such large variables...?
The function supports the following data-types: Integer_Type, Double_Type, Complex_Type, Char_Type, String_Type, BString_Type, Null_Tpye, Void_Type (=Undefined_Type), as well as Array_Type, Assoc_Type, Struct_Type, List_Type Vector_Type, Ref_Type (as structure fields)
Example
% define a structure variable a = struct { example = "foo" };
% save the structure into a file save_slang_variable("mystruct.sl", a);
% restore the variable into a new one variable b; (b,) = evalfile("mystruct.sl");
% edit the original variable % using a temporary filename save_slang_variable("/tmp/myedit", &a; edit, delete);
% this operation can be performed using % 'edit_var' as well (but better readable) edit_var(&a);
See also: evalfile, edit_var
show_slang_code
Usage
show_slang_code(String_Type function);
Description
The function show_slang_code
reads all .sl
files
in the directories contained in the S-Lang load path,
and *tries* to find the definition of function
by parsing the code for {}-brackets and comments.
In its current version, show_slang_code
gets confused, e.g.,
from {}-brackets and % characters in strings.
simplify_polygon
Synopsis
Simplify polygons using the Douglas Peucker Algorithm
Usage
(xx,yy)=simplify_polygon(x,y,d);
Description
Use the Douglas Peucker Algorithm to simplify the polygon P defined by the positions in the arrays x,y such that the maximum distance between all segments of the resulting polygon P2 defined by (xx,yy) is smaller than d. See the help for the function douglas_peucker for caveats Note that the argument d in simplify_polyon defines the distance, while the corresponding argument in douglas_peucker defines the distance squared!
See also: douglas_peucker
solveODEbyIntegrate
Synopsis
solves an ordinary differential equation on a grid by integration
Usage
x = solveODEbyIntegrate(&f, tgrid);
Qualifiers
rdt - integration step size relative to 'tgrid' binsize (default: 0.1) x0 - integration constant (default: 0) t0 - time reference for x0 (default: first 'tgrid' bin) method - integratipn method (default: &integrateRK4)
Description
The integral of the ordinary differential equation (ODE)
dx/dt = f(x(t), t) with x(t0) = x0
reads
x(t) = x0 + int_{t0}^{t} f(x(t'), t') dt'
&f
is a reference to a function with two arguments:
define f(x, t)
{
return ...;
}
See also: integrateRK4
solve_2d_system_of_equations
Usage
(x, y) = solve_2d_system_of_equations(ax1, ay1, c1, ax2, ay2, c2)
Description
ax1 * x + ay1 * y = c1
ax2 * x + ay2 * y = c2
solve_Eddington_quartic_equation
Synopsis
Solve the Eddington quartic equation
Usage
Double_Type beta = solve_Eddington_quartic_equation(Double_Type M, mu)
Description
For a given mass 'M' (in solar masses) and mean molecular weight 'mu', solve the Eddington quartic equation 1-beta = 0.003*M^2*mu^4*beta^4 numerically for 'beta' by using Newton's method to find the root of a function.
Notes
According to Equation (2-15) in the book "Principles of Stellar Evolution and Nucleosynthesis" by D.D. Clayton, the mean molecular weight can be approximated by mu ~ 2/(1+3*X+0.5*Y) where 'X' and 'Y' are the mass fractions of hydrogen and helium.
The Eddington quartic equation needs to be solved to compute the polytropic constant of the polytropic standard model, see the function 'polytropic_standard_model'.
Example
X = 0.73; Y = 0.26; mu = 2/(1+3*X+0.5*Y); beta = solve_Eddington_quartic_equation(1, mu);
See also: polytropic_standard_model
solve_Lane_Emden_equation
Synopsis
Solve the Lane-Emden equation numerically
Usage
Struct_Type s = solve_Lane_Emden_equation(Double_Type n)
Description
The Lane-Emden equation for index 0<'n'<5 is a second order differential equation of the form
1 d ( dtheta) ---- --- (xi^2 ------) = -theta^n xi^2 dxi ( dxi)
with the boundary conditions
dtheta theta = 1, ------ = 0 at xi = 0 . dxi
The solution 'theta' satisfying these boundary conditions is called Lane-Emden function of index 'n' and is returned by this function together with its derivative dtheta.
The Lane-Emden functions are used to study the structures of so-called polytropes, that is, gaseous spheres in hydro- static equilibrium which obey a polytropic equation of state: pressure = constant times density^((n+1)/n). For details, see, e.g., Chapter 2.4 in the book "Principles of Stellar Evolution and Nucleosynthesis" by D.D. Clayton. The output structure contains also the polytropic constants Rn = first zero of the Lane-Emden function of index n,
dtheta Mn = -xi^2 ------ evaluated at Rn. dxi
( 3 dtheta)^(-1) Dn = -(-- ------) evaluated at Rn, (xi dxi)
Rn^((n-3)/n)(3Dn)^((3-n)/3n) Bn = ---------------------------- . (n+1)Mn^((n-1)/n)
Notes
The numerical solution of the Lane-Emden equation is achieved here by rewritting it as a system of first order differential equations d --- [y_0,y_1] = [y_1, -(y_0)^n - 2/xi y_1 ] dxi which then allows an adaptive Runge-Kutta method of 4./5. order to be applied.
Qualifiers
- tol [=1e-10]: Absolute error control tolerance. Lower limit is 1e-15.
Example
s = solve_Lane_Emden_equation(3); print(s);
See also: polytropic_standard_model
sort_arrays
Synopsis
Usage
Integer_Type[] sort_arrays( Array_Type[] A1 [, A2 [, A3 ...]] );
Qualifiers
- dir [=1]:1: Ascending sorting; -1: Descending sorting
- method [='msort']: Choose sorting method, 'msort' for merge- and 'qsort' for quick-sort.
Description
Similar to 'array_sort' this funtion sorts arrays, but also sorts multiple appearing array entries based on an addition given array and so on. In other words this function allows to sort several arrays (e.g., columns in fits-file) depending on the sorting of the previouse array.
This is especially usefull if the arrays represent parameter combinations and you want to perform a resorting. See example below.
NOTE THAT all arrays must have the same length!
Example
% Example 1: a1 = [ 5,6,1,7,5,6,7 ]; a2 = [ 8,9,4,4,1,2,3 ]; i12= sort_arrays(a1,a2); array_map( Void_Type, &vmessage, "%.2d %.2d", a1[i12], a2[i12] ); i21 = sort_arrays(a2,a1); array_map( Void_Type, &vmessage, "%.2d %.2d", a1[i21], a2[i21] );
% Example 2: resorting parameter combinations: par = struct{ a = [0:1:#3 (closed)], b = ["hello","world"], c = [ 2001, 2012] }; parcomb = get_par_combinations( par ); print(parcomb);
parcomb = struct_array_2_struct_of_arrays(parcomb); struct_filter( parcomb, sort_arrays(parcomb.b,parcomb.a,parcomb.c) ); print(struct_of_arrays_2_struct_array(parcomb));
See also: array_sort, get_par_combinations
sort_struct_arrays
Usage
Struct_Type sort_struct_arrays(Struct_Type s, String_Type fieldname);
sphere_urand
Synopsis
generate spherical coordinates uniformly distributed on a sphere
Usage
(Double_Type theta, phi) = sphere_urand();
or
(Double_Type theta[], phi[]) = sphere_urand(Integer_Type n);
Description
The function generates n
(default: 1
) pairs (theta[i], phi[i])
of spherical coordinates -- such that the points
[ sin(theta)*cos(phi), sin(theta)*sin(phi), cos(theta) ]
are uniformly distributed on the surface of the unit sphere.
spherical_volume
Synopsis
computes the volume of an object parameterized in spherical coordinates
Usage
Double_Type vol = spherical_volume(Ref_Type &r);
Qualifiers
- Ntheta[=180]
- Nphi[=180]
- logfile[=NULL]: If logfile is a filehandle, the function evaluations are logged.
Description
r
has to be a real function with two arguments (theta, phi).
vol = int_0^2pi dphi int_0^pi dtheta sin(theta) int_0^r(theta, phi) dr r^2
= int_0^2pi dphi int_0^pi dtheta sin(theta) r(theta, phi)^3/3
The numerical integration of phi and theta is performed on the following grid:
theta = [0 : Pi : #Ntheta]; phi = [0 : 2*Pi : #Nphi];
split_lc_at_gaps
Synopsis
splits a light curve at gaps of a certain length
Usage
Struct_Type split_lc[] = split_lc_at_gaps(Struct_Type lc, Double_Type gap_threshold);
Qualifiers
- time [
="time"
]: time field
Description
lc
has to be a structure containing a time
field, which is an array of ascending values,
and other fields, which are arrays of the same length. As soon as the difference of
sequential times is larger than gap_threshold
, the structure is split, such that finally
an array of structures like lc
is returned.
The following constructions are equivalent:
split_lc_at_gaps(lc, gap_threshold)
split_struct(lc, blocks_between_gaps(lc.time, gap_threshold))
See also: split_struct, blocks_between_gaps
split_struct
Synopsis
splits a structure of related arrays into several structures of the same kind
Usage
Struct_Type structs[] = split_struct(Struct_Type s, Integer_Type group[]);
Description
All fields of the structure s
, as well as group
, have to be related arrays
of equal length. The array elements with the same corresponding group
value
are selected to form a new structure, unless the group
value is negative.
In the latter case (group<0
), the array elements are discarded.
The split structures are returned in an array, indexed by the group
value.
In other words:
structs[g].field = s.field[where(group==g)]
for 0 <= g <= max(group)
and any field
of the structure s
.
Examples
variable prime_struct = struct { i, p };
prime_struct.i = [ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
prime_struct.p = [ 2, 3, 5, 7, 11, 13, 17, 19, 23, 29];
variable group = [-1, 1, -1, 2, 0, 1, 2, 3, 1, 3];
% which might be obtained from the following selection:
group = Integer_Type[length(prime_struct.p)];
group[*] = -1;
group[where(prime_struct.p mod 10 == 1)] = 0;
group[where(prime_struct.p mod 10 == 3)] = 1;
group[where(prime_struct.p mod 10 == 7)] = 2;
group[where(prime_struct.p mod 10 == 9)] = 3;
variable primes = split_struct(prime_struct, group);
writecol(stdout, primes[1].i, primes[1].p); % primes with last digit 3
See also: split_lc_at_gaps
split_struct_cols
Synopsis
split structure in an array of structures with maximal num columns
Usage
Struct_Type str_array = split_struct_cols(Struct_Type str,
Integer_Type num);
Description
This function splits a structure in an array of structures, with maximal 'num' columns in each structure.
This can be useful, as fits-tables only allow 999 colums to be written in the same extension. With split_struct_cols these can be splitted an written into different extensions. Afterwards they can be easily combined with struct_combine().
See also: fits_save_fit, struct_combine
storevar
Synopsis
use `save_slang_variable' instead
strftime_MJD
Synopsis
formats an Modified Julian Date as a string
Usage
String_Type strftime_MJD([String_Type fmt,] Double_Type MJD)
Description
If fmt
is not specified, "%Y-%m-%d, %H:%M:%S"
is assumed.
See also: strftime, gmtime, MJD2UNIXtime
string_intersection
Synopsis
Returns the common substrings within all strings of the given array
Usage
String_Type string_intersection(String_Type[] strings);
Qualifiers
minlen - restrict the returned substrings to a minimum length longest - return the longest substring only
Example
a = string_intersection(["mister", "twister"]; minlen = 3); % -> ["ister", "ster", "ter", "ste", "iste", "ist"]
string_match_perl
Synopsis
matches a string against a RE and stores the results in $1, $2, ...
Usage
Integer_Type string_match(String_Type str, pat[, Integer_Type pos])
Description
The string str
is matched against the regular expression pat
,
starting at the position pos
(numbered from 1, which is the default).
The function returns the position of the start of the match in str
.
The entire matching string is stored in the global variable $0,
the substrings of patterns enclosed by "\( ... \)" pairs
are stored in $1, $2, ..., $9 (as in Perl).
string_match_perl
is deprecated. Use S-Lang's string_matches instead.
Example
if( string_match_perl("hello world", "\([a-z]+\) \([a-z]+\)"R) )
()=printf("1=%s, 2=%s
", $1, $2);
%
Equivalent way using string_matches
:
variable m = string_matches("hello world",
\([a-z]+\) \([a-z]+\)
, 1);
if(m!=NULL) vmessage("1=%s, 2=%s", m[1], m[2]);
See also: string_matches
strjoin_wrap
Synopsis
concatenates elements of an array and inserts linebreaks
Usage
Sting_Type strjoin_wrap(Array_Type x, String_Type delim);
Qualifiers
- maxlen[=72]: maximum number of characters in a line
- newline[=" "]: new line string
- initial[=""]: initial delimiter
- final[=""]: final delimiter
See also: strjoin
strreplace1
Synopsis
replaces one occurence of a substring in a string by another string
Usage
String_Type res = strreplace1(String_Type str, search_str, repl_str);
Description
replaces the first occurence of search_str
in str
by repl_str
:
(res, ) = strreplace(str, search_str, repl_str, 1);
strreplaces
Synopsis
Replace one or more substrings
Usage
String_Type strreplaces(String_Type a, String_Type[] b, String_Type[] c[, Integer_Type max_n]);
Description
This is a modificaiton to 'strreplace', where multiple, different substrings (given as arrays) are replaced at the same time.
See also: strreplace
strsplit
Synopsis
split a string at a separator and return contents as an array
Usage
Array_Type strsplit(String_Type s, String_Type separator);
Description
This routine takes a string s and splits it into parts separated by a separator character.
Example
variable ra="12:34:56.78"; variable rastr,hh,mm,ss; rastr=strsplit(ra,":"); hh=atof(rastr[0]); mm=atof(rastr[1]); ss=atof(rastr[2]); ra=hms2deg(hh,mm,ss);
struct_array
Synopsis
creates an array of structures
Usage
Struct_Type[] struct_array(Integer_Type n, Struct_Type s);
Description
An array of n structures is created, where each structure is a copy of the given one s.
Example
variable sarr = struct_array(4, struct { firstname, lastname }); sarr[2].firstname = "Karl"; sarr[2].lastname = "Remeis";
print(sarr); % will return % {firstname=NULL, lastname=NULL} % {firstname=NULL, lastname=NULL} % {firstname="Karl", lastname="Remeis"} % {firstname=NULL, lastname=NULL}
See also: table_copy
struct_array_2_struct_of_arrays
Synopsis
Converts a Struct_Type[] to Struct_Type with array fields
Usage
Struct_Type struct_array_2_struct_of_arrays(Struct_Type[] SA);
Qualifiers
- depth [0]: If a field is a Struct_Type itself, the this function calls itself iteratively 'depth' times going down the rabbit hole.
- dim [0]: If a field is multi-dimensional 'resh' is the dimension to which the array is appended to.
Description
This function converts a Struct-Array 'Struct_Type[]' to a single Struct_Type by converting the COMMON fields of all structs into arrays. The function works iteratively in case a field of the given struct array SA itself is a Struct_Type. The function will call itself maximaly 'depth' times. In case a field is multi-dimensional the array will be appended to the 'dim'-th dimension.
NOTE that this function runs orders faster than 'merge_struct_arrays' while offering the same features as, e.g., reshaping multi-dimensional arrays (the reshaping in this function also works if not every field is multidimensional!).
The only drawback is that the function fails if an entry is NULL as this case cannot be catched. In that case you may use 'merge_struct_arrays'!
Examples
%%%%% SIMPLE example variable sa = struct_array( 3, struct{ a=1, b=2 } ); sa[0] = struct_combine( struct{c}, sa[0] ); variable s = struct_array_2_struct_of_arrays( sa ); print(s); {a=Integer_Type[3], b=Integer_Type[3]}
%%%%% COMPLEX example variable sa = struct_array( 30000, struct{ a=1, b=_reshape([1:120],[1,2,60]), c=struct{ d="hello" } } ); %%% only merge the parent structure variable s = struct_array_2_struct_of_arrays( sa ; depth=0 ); isis> print(s); % {a=Integer_Type[30000], % b=Integer_Type[30000,2,60], % c=Struct_Type[30000]}
%%% merge also the first generation variable s = struct_array_2_struct_of_arrays( sa ; depth=1 ); print(s); % {a=Integer_Type[30000], % b=Integer_Type[30000,2,60], % c=Struct_Type with 1 fields} print(s.c); % {d=String_Type[3]}
%%% use the second dimension of mudlti-dimensional arrays for merging: variable s = struct_array_2_struct_of_arrays( sa ; dim=1 ); print(s); % {a=Integer_Type[30000], % b=Integer_Type[1,60000,60], % c=Struct_Type[30000]}
%%% Runtime comparision to 'merge_struct_array' tic; s = struct_array_2_struct_of_arrays( sa ); vmessage(" sa2soa took %g seconds to call",toc); print(s); tic; m = merge_struct_arrays( sa ); vmessage(" msa took %g seconds to call",toc); print(m);
See also: merge_struct_array, array_struct_field, array_flatten, unique, intersection
struct_copy
Synopsis
makes a copy of a struct with all its fields
Usage
Struct_Type copy = struct_copy(Struct_Type struct);
Description
struct_copy
copies a Struct_Type
with all its
fields, which can be of any Type! If a field is an Array_Type[]
struct_copy
calls array_copy
and if a field is
Struct_Type
, struct_copy
calls itself.
Further it copies Ref_Type, which allows
!!! copiing XFIG-OBJECTS !!!
Example
See also: array_copy
struct_field_or_default
Synopsis
return value of a structure field or a default value
Usage
value=struct_field_or_default(Struct_Type s, String_Type tag, default);
Description
return s.tag if it exists, otherwise default
struct_fields_list_to_array
Synopsis
convert all fields of a structure from lists to arrays
Usage
struct_fields_list_to_array(Struct_Type s[, DataType_Type t]);
See also: list_to_array
struct_of_arrays_2_struct_array
Synopsis
Converts a Struct_Type with array fields to a Struct_Type[]
Usage
Struct_Type[] struct_of_arrays_2_struct_array(Struct_Type S);
Description
This function converts a 'Struct_Type' containing fields which are arrays (all arrays must have the same dimension!) to a Struct-array 'Struct_Type[]' with the same fields but one dimensional. That means the returned Struct_Type[] has then the dimension (array_shape) of the fields. The returned Struct-Array will contain all fields the original 'S' also had, BUT for those fields with a different array_shape then the first field and the value will be set to the first entry of its field!
Example
variable s = struct{ a=[1:3], b=[4:6], c=[3,20] }; variable sa = struct_of_arrays_2_struct_array( s ); print(sa); {a=1, b=4, c=3} {a=2, b=5, c=3} {a=3, b=6, c=3}
See also: array_map
struct_to_assoc_array
Synopsis
converts a structure to an associative array
Usage
Assoc_Type struct_to_assoc_array(Struct_Type structure[, DataType_Type data-type]);
Description
The given structure is converted into an associative array. Thereby, the keys of the returned array will be equal to the field names of the structure.
If the second argument is provided, the values of the associative array will be of that type. Otherwise, the type of the first field of the structure will be used. Note, that all fields have to be of the same data-type!
Example
c = struct_to_assoc_array(struct { pi = 3.141592653589793, c = 2.99792458e10, h = 6.62606876e-27 }); foreach key (c) using ("keys") vmessage("c[%s] = %e", key, c[key]); ;
See also: Assoc_Type, assoc_get_keys, get_struct_field_names
sunpos
Synopsis
computes the RA and Dec of the Sun at a given date.
Usage
(Double_Type eq, Double_Type ecl) = sunpos(Double_Type JD);
Description
Computes the position of the Sun at the date JD (or dates, if JD is an array), given in Julian Date. Position calculation takes into account pertubations from the planets, aberration and nutation. Output value 'eq' is a structure in equatorial coordinates: eq = struct{ra, dec}, 'ecl' is in ecliptic plane coordinates: ecl = struct{lon, obl} (Longitude and true obliquity).
Note: Function adopted from the IDL function with the same name. See there for accuracy tests, which should also apply here.
Qualifiers
- MJD: give time in MJD instead of JD
- radian: return all variables as radians rather than degrees
- hour: return RA in hours rather than degrees (not possible when radian is set).
Example
(eq, ecl) = sunpos(2455055.5) ; ()=printf("Right ascension: %.4f deg, Declination: %.4f deg
", eq.ra, eq.dec) ;
table_add_row
Usage
table_add_row(table, ...);
table_copy
Synopsis
makes a copy of all columns of a table
Usage
Struct_Type copy = table_copy(Struct_Type table);
Description
copy.field = @(table.field); %
for every field
NOTE: This also works for XFIG-PLOT-OBJECTS
table_print_TeX
Synopsis
prints certain columns of a structure as a TeX table
Usage
table_print_TeX(Struct_Type tab[, String_Type cols[]]);
Description
If no columns are specified by the string-array cols
,
all columns of the table tab are used.
Qualifiers
- align: alignment of all columns: "l", "c", "r" [default: "l"]
- heading[=
cols
]: headings for TeX table - output[=stdout]: output can be written in a file by stating a File_Type or the filename (String_Type)
tai2tt
Synopsis
Convert a (M)JD in International Atomic Time into Terrestrial Time (TT, also known as TDT or ET)
Usage
tt=tai2tt(JD)
Description
This routine converts a Julian Date that is given in International Atomic Time (TAI) into Terrestrial Time, which is roughly corresponding to Terrestrial Dynamic Time or Ephemeris Time. This is done by adding 32.184s to the TT. This routine is array safe.
Qualifiers
- mjd: The argument is given in MJD, not in JD.
- deltat: Return the TT-TAI in seconds, rather than the corrected (M)JD.
See also: utc2tai,tai2utc,tt2tai,tt2tdb
tai2utc
Synopsis
Converts a (M)JD in International Atomic Time into civil time (UTC)
Usage
utc=tai2utc(JD)
Description
This routine converts a Julian Date that is given in international atomic time into Coordinated Universal Time (UTC) by applying the relevant correction for leap seconds. This routine is array safe.
NOTE: This function is only approximately correct. Within about a minute of a leap second this routine can be off by as much as one second.
Qualifiers
- mjd: The argument is given in MJD, not in JD.
- leapfile: path to a file containing the leap second information. This file is available, e.g., from ftp://maia.usno.navy.mil/ser7/tai-utc.dat Most astronomical software systems distribute this file. The default is to look for it in $LHEASOFT/refdata
See also: utc2tai,tt2tai,tai2tt
taylor
Synopsis
Returns the taylor series for given coefficients
Usage
Double_Type taylor(Double_Type[] x, coefficients)
Description
Returns the taylor series around x with the given 'coefficients', which lengths determines the order. The sum is evaluated using the Horner (1819) schema for numerical accuracy and computation speed:
c0 + c1*x + 1/2!*c2*x^2 + 1/3!*c3*x^3 + ... = (((... + c3)*x/3 + c2)*x/2 + c1)*x + c0
Example
x = [0, 1]; c = [1, 4e-3, 0, 3e-6]; % up to third order
taylor(x, c); % returns [1.0, 1.0040005]
taylorcoeff_from_struct
Synopsis
extracts taylor coefficients from the field names of a structure
Usage
Double_Type[] taylorcoeff_from_struct(Struct_Type struct[,
String_Type 0th_pattern, String_Type Nth_pattern]);
Description
If a structure has fields which names reflect the different orders of a taylor series, this function extracts these coefficients by matching the field names against two patterns.
'0th_pattern' is the field name of the zero order and 'Nth_pattern' this of all higher order terms. Both patterns have to be regular expressions. In particular, the expression of the Nth order must contain a number extraction two determine the order N of that field.
By default, '0th_pattern' is set to "[a-zA-Z]0" and 'Nth_pattern' to "a-zA-Zdot"R. If the latter _pattern matches, but the extracted number is an empty string N=1 is assumed.
The returned array contains the taylor coefficients in ascending order as used by the 'taylor'-function.
See also: taylor, string_matches
time_array
Synopsis
Converts a time [s] (e.g., _time) into an array [sec,min, ...]
Usage
Integer_Type[] = time_array( Integer_Type );
or
```c
String_Type = time_array( Integer_Type ; str);
##### Qualifiers
* str: Instead of a time array a String_Type is returned!
##### Description
Converts a time given in seconds into a time array, where
the first entry corresponds to the seconds, the second to
the minutes, etc. The length of the array is variable,
i.e., if the given time is of the order of hours, the array
will not include days and years.
If the qualifier 'str' is given the array will be converted
into a string (see example).
##### Example
isis> ta = format_time ( _time ); ta;
Integer_Type[5];
isis> print(ta);
54
1
14
302
43
%or
isis> ts = format_time ( _time ;str ); ts;
43y 302d 14h 1m 54s
__See also__: _time
#### tt2tai
##### Synopsis
Convert a (M)JD in Terrestrial Time (TT, aka TDT or ET) into TAI
##### Usage
```c
tai=tt2tai(JD)
Description
This routine converts a Julian Date that is given in Terrestrial Time (formerly known as Terrestrial Dynamic Time or Ephemeris Time) into TAI. This is done by subtracting 32.184s from the TT. This routine is array safe.
Qualifiers
- mjd: The argument is given in MJD, not in JD.
- deltat: Return the TAI-TT in seconds, rather than the corrected (M)JD.
See also: utc2tai,tai2utc,tai2tt,tt2tdb
tt2tdb
Synopsis
Convert a (M)JD in Terrestrial Time (TT, aka TDT or ET) into Terrestrial Barycentric Time
Usage
tdb=tt2tdb(JD)
Description
This routine converts a Julian Date that is given in Terrestrial Time into Terrestrial Dynamic Barycentric Time (the timescale of JPLs DE405 ephemeris). The routine uses equation 2.6 of Kaplan (2005, USNO Circular 179), which is better than 10mus between AD1600 and AD2200. Note that normally it is sufficient to use TT as the argument to DE405 or VSOP1987.
An inverse routine is not provided, since there is typically no need to convert back from TDB, but such a routine can be written using the delta t provided here.
This routine is array safe.
Qualifiers
- mjd: The argument is given in MJD, not in JD.
- deltat: Return the TAI-TT in seconds, rather than the corrected (M)JD.
See also: utc2tai,tai2utc,tai2tt
unique_n
Synopsis
finds the indicies of unique tupels
Usage
Integer_Type[] unique_n(Array_Type a1, ..., an)
Description
The arrays a1
, ..., an
(with n>=1
) must have the same length,
but may contain different data types.
The return value is an array of all indices i
of all unique tupels (a1[i], ..., an[i])
.
unpack_obj
Synopsis
converts a binary string back into an SLang object
Usage
Any_Type pack_obj(BString_Type obj);
Description
An input array of binary strings (BString_Type), which has
been created by the pack_obj' function, is converted back into an SLang object using the
unpack' function.
See the documentation of pack_obj
for more details.
Example
s = pack_obj(PI); % s[0] = "d\030-DT\373!\011@"
a = unpack_obj(s); % a = 3.141592653589793;
See also: pack_obj, unpack, pack
utc2tai
Synopsis
Convert a (M)JD in civil time (UTC) into International Atomic Time
Usage
tai=utc2tai(JD)
Description
This routine converts a Julian Date that is in civil time (Coordinated Universal Time, UTC) into International Atomic Time (TAI) by applying the relevant correction for leap seconds. This routine is array safe.
Qualifiers
- mjd: The argument is given in MJD, not in JD.
- deltat: Return the TAI-UTC in seconds, rather than the corrected (M)JD.
- leapfile: path to a file containing the leap second information. This file is available, e.g., from ftp://maia.usno.navy.mil/ser7/tai-utc.dat Most astronomical software systems distribute this file. The default is to look for it in $LHEASOFT/refdata
See also: utc2tai,tt2tai,tai2tt
v_relativistic_to_optical
Synopsis
Convert relativistic velocity to optical velocity.
Usage
v_relativistic_to_optical(V_rel);
Description
Calculate the optical velocity definition from a given relativistic velocity. Velocities must be given in km/s.
vacuum_to_air
Synopsis
Convert vacuum wavelengths to air wavelengths
Usage
Double_Type l_a[] = vacuum_to_air(Double_Type l_v[])
Description
Convert vacuum wavelengths l_v (in Angstroem) to air wavelengths l_a (in Angstroem) for dry air at 15 degrees Celsius, 101.325 Pa, and with 450 parts per million CO2 content. Valid in the wavelength range between 2300 to 17000 Angstroem.
Notes
According to Equation 1 in Ciddor 1996, Applied Optics, 35, 1566: 1e8*(l_v-l_a)/l_a = n-1 = k1/(k0-l_v^(-2)) + k3/(k2-l_v^(-2)) => l_a = l_v / ( k1*1e-8/(k0-l_v^(-2)) + k3*1e-8/(k2-l_v^(-2)) + 1 ) with k0 = 238.0185 * 1e-8 / Angstroem^2, k1 = 5792105 * 1e-8 / Angstroem^2, k2 = 57.362 * 1e-8 / Angstroem^2, k3 = 167917 * 1e-8 / Angstroem^2.
Qualifiers
- verbose [=1]: Set to 0 to suppress warnings.
Example
l_a = vacuum_to_air([2000,4000,8000,16000]);
See also: air_to_vacuum
vector_astro
Synopsis
vector initializer that is more appropriate for astronomy
Usage
Vector_Type vector_astro(Double_Type x, y, z)
or
```c
Vector_Type vector_astro(Double_Type r, phi, theta; sph)
or
Vector_Type vector_astro(Double_Type r, lon, lat; astro)
or
Vector_Type vector_astro(Double_Type phi, theta; sph)
or
Vector_Type vector_astro(Double_Type lon, lat; astro)
##### Qualifiers
* astro: consider the given coordinates to be astronomical (r, lon, lat)
following the astronomical definition (i.e., lon ranges from
0 to 2pi, lat from -pi/2 to pi/2; default for two arguments)
* sph: consider the given coordinates to be spherical (r, phi, theta)
following the mathematical definition (i.e., theta ranges from
0 to pi)
* deg: interpret angular arguments in degrees
##### Description
The components of the vector are returned within the Vector_Type
and are accessible like a structure with the fields x, y, and z.
If the astro qualifier is given, cartesian coordinates are
calculated as
<pre>
[x, y, z] =
r \\* [cos(lon)\\*cos(lat), sin(lon)\\*cos(lat), sin(lat)]
</pre>
If the sph qualifier is given, the cartesian coordinates are
calculated as
<pre>
[x, y, z] =
r \\* [cos(phi)\\*sin(theta), sin(phi)\\*sin(theta), cos(theta)]
</pre>
If the function is called with only two arguments, then these represent the
spherical or astronomical coordinates, and we assume that r=1.
If the initializer is called with a single argument of type Vector_Type,
the corresponding Vector_Type is returned
__See also__: Vector_Type,Matrix33_Type,vector_to_spherical,dms2deg,hms2deg
#### vector_to_spherical
##### Synopsis
Returns the spherical coordinates corresponding to a 3D Vector
##### Usage
```c
(r,lon,lat)=vector_to_spherical(Vector_Type v)
Qualifiers
- astro: return coordinates using the astronomical convention (the default)
- sph: return coordinates (r,theta,phi) using the mathematical definition (i.e., theta ranges from 0 to pi)
- deg: return angular arguments in degrees (default: radian)
Description
This function returns the spherical coordinates corresponding to a vector. See the description of function vector_astro for the relevant equations.
See also: Vector_Type,vector_astro,dms2deg,hms2deg
vsop87
Synopsis
calculates high precision planetary positions using the VSOP87 theory
Usage
Struct_Type vsop87(JD,planet;qualifiers)
Description
This routine implements the planetary theory VSOP87 of Bretagnon and Francou (1988, A&A 202, 309). The routine can return high precision planetary positions for all planets for a given (array) of times (in TT).
Depending on the qualifiers the routine can return rectangular coordinates (X,Y,Z) or ecliptical coordinates (longitude, latitude, distance) in a heliocentric coordinate system either for epoch J2000 or for the epoch of date. In addition, elliptical orbital elements and standard orbital elements can also be calculated.
By default the routine returns rectangular heliocentric coordinates in a structure with elements x, y, z, vx, vy, vz. Here x, y, z are in an ecliptical coordinate system for the epoch J2000, where the x-axis is in the direction of the vernal equinox, the z-axis points to the ecliptical North pole, and the y-axis forms a right handed coordinate system with the x- and z-axes. The default units are AU and AU/day.
For heliocentric coordinates the planet can have the values "Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", and "Neptune". For heliocentric rectangular coordinates the value "Earth-Moon" will return the position of the Earth-Moon barycenter.
In the case of solar system barycentric data, it is not possible to obtain the elements of the Earth-Moon barycenter, only the coordinates of the Earth can be queried. In addition, for this set of coordinates it is also possible to request the rectangular position of the Sun.
For the years 1900-2100 interval the precision of the positions is at the level of a few milli-arcseconds, with the exception of Saturn, for which a precision of "only" 0.1'' is reached.
The precision is better than 1" for the following time intervals around AD2000: Mercury, Venus, Earth-Moon Barycenter, Mars: 4000 years Jupiter, Saturn: 2000 years Uranus, Neptune: 6000 years
The time argument to the routine is the Julian Date in Terrestrial Time (TT), which corresponds to the international atomic time TAI+32.184s.
This is a powerful routine with a large number of options. A set of helper functions provide simpler interfaces.
Qualifiers
-
mjd: The argument is in MJD(TT), not in JD(TT)
-
barycentric: Return coordinates with respect to the solar system barycenter rather than heliocentric coordinates. This implies that coordinates are for J2000.
-
spherical: Return spherical coordinates in the ecliptical coordinate system. Only possible for heliocentric data. In this case the structure returned contains the elements lat, lon, r, latdot, londot, rdot, where lat and lon are the ecliptical latitude and longitude (in rad), r is the distance in AU, and latdot, londot, rdot are their time derivatives (in rad/day or AU/day, respectively.)
-
mean_equinox: Coordinates are given in the frame defined by the mean equinox and ecliptic of the date. In VSOP87, this system is derived from the J2000 one by pure precessional motion. The default is to return J2000 coordinates.
-
elements: Return the orbital elements. The structure contains the following tags: a: semi-major axis (in AU) l: mean longitude (rad) k=e cos(p): elliptical element k h=e sin(p): elliptical element h q=sin(g) cos(omega): elliptical element q p=sin(g) sin(omega): elliptical element p e: orbital eccentricity plon: perihelion longitude (in rad) g: semi inclination (rad) i: inclination (rad) omega: longitude of the ascending node (also called G)
-
deg: All angle quantities to be returned are in deg, or deg/day.
-
mks: All distances returned by the function are given in meters or m/s, depending on the deg qualifier, angles are returned in rad and rad/s or deg and deg/s.
-
cgs: Dito, but use centimeters instead of meters.
-
vsopfile: Path to file containing the VSOP quantities; only necessary if this is a non-standard installation of isisscripts. This qualifier is only used if the vsop87 data have not yet been read. By default, the file with the coefficients is searched for in the directory pointed to by the environment variable ISISSCRIPTS_REFPATH.
-
forcearray: Force the tags in the returned struct to be arrays. The default is that they are arrays only if JD is an array.
See also: XXXXX
wavelength_to_rgb
Synopsis
computes the color values of visual light
Usage
Integer_Type[] wavelength_to_rgb(Double_Type lambda[]);
Description
380 <= lambda <= 780
is the wavelength in nm.
The return value is its 24bit RGB color value.
The conversion is performed after the code by Dan Bruton, see http://www.physics.sfasu.edu/astro/color.html.
weighted_mean
Usage
Double_Type m = weighted_mean(x, w);
or
Double_Type m = weighted_mean(x; err=err); % or weighted_mean(x, err; err)
Description
m = sum(w*x) / sum(w);
%
for Gaussian errors:
w = 1/err^2;
where_max
Synopsis
finds the indices of an array where the values are maximal
Usage
Integer_Type i[] = where_max(Array_Type a);
Description
i = where(a==max(a)); % => a[i] == max(a)
See also: where_min
where_min
Synopsis
finds the indices of an array where the values are minimal
Usage
Integer_Type i[] = where_min(Array_Type a);
Description
i = where(a==min(a)); % => a[i] == min(a)
See also: where_max
winding_number_polygon
Usage
wn=winding_number_polygon(P,V)
##### Synopsis
return the winding number for a point P in a polygon
##### Description
return the winding number for a point P in a (potentially
complex polygon V)
The winding number counts the number of times the polygon V winds
around point P. The point is outside if the winding number is 0.
The algorithm used here is as efficient as the determination
of the crossing number.
The point P is a struct{x,y}, while the polygon is a
struct{x[],y[]} where the arrays are the x- and y-coordinates.
The polygon has to be closed, i.e. V.x[n]==V.x[0] and
V.y[n]==V.y[0] where n is the number of polygon points.
See the URL below for more explanations.
Based on code by Dan Sunday,
http://geomalgorithms.com/a03-_inclusion.html
__See also__: crossing_number_polygon,point_in_polygon,simplify_polygon
#### write_slurm_script
##### Synopsis
writes a slurm job script, which can be submitted with squeue
##### Usage
```c
write_slurm_script(file,jobname,walltime,cmds);
Qualifiers
- queue:set the queue (remeis by default)
- silent: don't show any output
- serial: instead of parallel, execute commands one after the other
- option: setting commands or environment variables before executing the actual commands
- ntaks[=1]: number of tasks per CPU
- memory[=1000]: MByte in memory allocated (--mem-per-cpu in slurm)
- output: output logfile name
- error: error logfile name
- account: if the selected queue needs a specific account set this qualifier accordingly
- dir[=getcwd()]: absolut path where the script should be run (default is cwd)
Description
Writing a slurm job script, which can be submitted with squeue. The function might yet not include all options and works best for simple cases. The cmds variable is a String_Type array, with each field being one command line. The
Example
% define your commands variable cmds = [ "echo 123" , "echo 456" ]; % call the function variable file = "slurmscript.slurm"; variable jobname = "script"; variable walltime = "00:30:00"; write_slurm_script(file,jobname,walltime,cmds);
xfig_mix_colors
Synopsis
mix two named xfig colors
Usage
xfig_mix_colors(color1,color2,fraction)
Qualifiers
- name: xfig name of the color mix
Description
This function mixes two named xfig colors in rgb space. The function looks up the rgb values of color1 and color2, and mixes their rgb values according to new color=color1*fraction+color2*(1-fraction) The operations are analoguous to the color mixing performed by the xcolor package of LaTeX (the operation is similar to LaTeX's color1!fraction!color2 syntax). The function returns the xfig name of the new color, or nothing if the name qualifier is given (recommended).
See also: mix_rgb_colors
xfig_new_X11color
Synopsis
define a new color from x11's rgb.txt for xfig
Usage
xfig_new_X11color(colorname)
Qualifiers
- x11nameName of the color in rgb.txt
- rgbfileFile defining the colors, default: /usr/share/X11/rgb.txt
Description
This function is a wrapper around X11color2rgb. It searches for the color using x11color2rgb and then defines a new xfig color with the same name. Note that some of xfig's colors can have the same name as names in rgb.txt but different rgb values. In this case the default of this function is to override xfig's definition. If this is not desired, use the x11name qualifier.
See also: x11color2rgb,xfig_lookup_w3c_color
xfig_plot_allbezier
Synopsis
adds the plot of all Bezier curves from the bezier function to a xfig plot object
Usage
Xfig_Struct pl_new = xfig_plot_allbezier(Xfig_struct pl, bez [, conlines ]);
Description
Adds the plot of all Bezier curves to a previously defined xfig-object 'pl'. 'conlines' and 'bez' are the values returned from the bezier-function called with both qualifiers "allbez" and "conlines". Colors of each Bezier line collection is changed automatically. If conlines are plotted, they will be plotted in "gray".
See also: bezier
yearOfMJD
Synopsis
transforms a date given as MJD into the year (with fractional parts)
Usage
Double_Type year = yearOfMJD(Double_Type MJD);
See also: jd2year
zams
Synopsis
return properties of the ZAMS per Tout et al., 1996, MNRAS 281, 257
Usage
(radius,luminosity,temperature)=zams(mass;z=metallicity);
Qualifiers
z - metallicity (between 0.0001 and 0.03; default: 0.02 (solar))
Description
This function uses the rational function fits of Tout et al. (1996, MNRAS 281, 257) and calculates the radius, luminosity, and temperature for zero age main sequence stars of a given mass (in solar units; can be an array, must be between 0.1 and 100 Msun). The radius and luminosity are returned in solar units, the temperature is in K. Tout et al. state that in general the fits are better than 7.5per cent in luminosity and 5 per cent in radius, for solar metallicity they are good to 3 per cent in L and 1.2 per cent in R.
zang
Synopsis
Determine the angular size of an object as a function of redshift
Usage
Double_Type angsize = zang (Double_Type size, Double_Type redshift);
Description
Requires an input size in kpc and returns an angular size
in arc seconds. (Analog to the function zsize
,
which converts angular size to projected size size given the
redshift.)
The default cosmology according is used. The function passes
qualifiers to the functions lumdist
and cosmo_param
allowing
to change the cosmology.
Example
variable size = 50.; % kpc variable z = 1.5; % redshift zang(50,1.5; omega_m = 0.3, omega_lambda = 0.0, silent); % ---> 6.58 arc seconds
See also: lumdist, cosmo_param
zsize
Synopsis
Determine the projected size of an object as a function of redshift
Usage
Double_Type proj_size = zsize (Double_Type angle, Double_Type redshift);
Description
Requires an input angular size in arc seconds and returns
the projected size in kpc. (Analog to the function zang
,
which converts projected size to angular size given the
redshift.)
The default cosmology according is used. The function passes
qualifiers to the functions lumdist
and cosmo_param
allowing
to change the cosmology.
Example
variable angsize = 1.; % arc seconds variable z = 1.5; % redshift zsize(angsize,z); % ---> 8.46 kpc
See also: zang, lumdist, cosmo_param