What to do to determine broadening functions (BF's) at DDO?

Slavek Rucinski

Version 4.3: July 3, 2007

The most recent set of IDL routines. Read in and compile as one file.


1. Always use the most recent update of this text. It is evolving over time!
2. Problems with reading in of the FITS files:
The routines BFpro*.pro utilize the routine readfits.pro from IDL-Astro. It is absolutely crucial that the FITS files are not corrupted by ascii transfer via FTP (must be binary!) or by packing/compression programs. It has been established that the use of the gzip (UNIX) and WinZip (Windows) program pair for packing & unpacking .gz files corrupts the FITS files. Such corrupted files are sometimes difficult to recognize as the garbage data may appear only in some sections of the spectra.
3. There may be problems when the FITS file headers have epoch & equinox 1900 or 1950. Make sure to use 2000 coordinates right from the start.


This description is somewhat specific to the current setup of the binary star program at DDO utilizing the grating 1800 l/mm. In what follows, we will use terms "old system" and "new system".

The old system (T or THX; files with prefix C) was based on the Thompson thick, square CCD chip with 1024 pixels of 19um (0.2A with the 1800 l/mm grating) giving about 200A of a spectral window at the central wavelength of 5184A, the magnesium triplet. This is the backup system when the JY2 system is not working. After the failure of the old, good original chip, we normally use the blue coated chip #3 (T3).

The new system (JY2), in use since October 2004 (files with K prefix), is based on a thinned Marconi chip, 512x2048, with pixels 13.5um. This system is heavier and may produce a slightly larger differential flexure. Migration to the new system appeared to increase the errors, although we were a bit too quick in proclaiming in 2004 that all was due to the increased flexure. We did experiments with changing the band, moving to the one with the centre at 6290A, right at the telluric oxygen band 6270-6330A. The dense forest of the band lines in 6287-6330A (typical depth 2-6%) was processed separately using the BF formalism against a master spectrum of Regulus, to determine the geocentric flexure shifts. The stellar BF was determined from the rest of the spectrum, spliced together into one spectrum. We removed the end vignetting problems (shutter?) as well so we used: 6165-6270 and 6330-6430A. This did not work that well: The lines in the near-IR are weak, not as deep as the MgI 5184A triplet (in majority of stars from middle-A down to middle-K). Thus, in the Spring of 2005 we returned to the MgI and this works well.

In October 2005, the JY2 system failed and we returned to the T3 system, but this is hopefully just a short glitch.


This description assumes that one follows the "alternative" route of the BF derivation, as described in the "Cookbook", i.e. first, the full BF is determined (without any singular values removed), then the result, with all pixels being independent, is locally smoothed by Gaussians (sigma typically 1.0 to 3.0) to match the natural smoothing of the spectrograph slit.. The appropriate smoothing is sigma=1.5 pix for the old system and sigma=2.0 for the new system.

All used IDL routines are kept as a file BFall_IDL.pro. It is suggested to save this file in the default library/routine directory of IDL. In Windows, the IDL interface permits to read this file in (File->Open) and then compile all its components (Run->Compile...) in one go.

Because of the inter-dependence of the routines, it may be advisable to compile (.run) the IDL procedures, as given in the collection BFall_IDL.pro, a few times in succession. Please let me know if you miss a routine. Better ask me than take from IDL-astro as they may not be compatible.

The routines discussed here make use of IDL-Windows settings for plotting (set_plot,'win'). You will have to change this for unix/linux settings for eg. X-windows: set_plot,'x'. For PostScript, use: set_plot,'ps'
which gives idl.ps (note that this device must be closed with: device,/close).


1. You must have spectra fully reduced in IRAF: wavelength calibrated (geocentric system is fine, i.e. no RV helio corrections applied), normalized to unity at maxima (rectified), etc. You must have all spectra of your program star in FITS format in one directory.

If your spectra are not rectified, use the routine:
as in:
This, after some cursor work (read the comments of the two IDL routines used in this process), writes into your directory a file
with "r" added to the name, in this case: K0003456r.fits.

2. Make sure that the FITS keywords will be recognized by the IDL programs. Only the following are:
Keyword   Typical      Meaning
CRVAL1    5079.53793d0 first wavelength
CDELT1    0.20360869d0 step in wavelength
NAXIS1    1024         length of spectrum
EQUINOX   2000.        equinox
EPOCH     2000.        treated as equinox
DATE-OBS '2002-08-08'  UT date new format
         '29/08/96'    old format
TIME-OBS '07:06:39'    UT start
UTSTART  '07:06:39'    same
ENDTIME  '07:23:30'    UT end
UTEND    '07:23:30'    same
RA       '23:11:31.00' hh:mm:ss.s
DEC      '-36:53:34.00' dd.mm.ss.s
OBJECT   'HD157856'

3. Place one good spectrum of the template in the same directory where the program stars observations are located. Your may want to have another directory for all standard stars, also with the same template in front because standards observed on individual nights should be analyzed in the same way as the program binary star to determine the nightly corrections or to be sure that the selected template observation did not have any problems. This is described below in the description of the handling of the std star observations.

It is usually the most important to check the velocity for the selected spectrum of the template. We sometimes see external random errors exceeding our nominal precision of about 1 km/s and reaching >2-3 km/s; it is therefore important to check if the selected template did not happen to be just accidentally wrongly observed (with the star guided on one side of the slit, etc.).

I suggest to spend some time selecting the best template. A good choice is usually a spectrum with the highest number of exp-meter counts.

4. List all spectra in one file with name preferably: star-name.lst.
ls *.fts > star-name.lst
In MS-DOS command window:
dir *.FTS /B > star-name.lst

This list file will be used by several routines. The file will have to be edited to move the name of the FITS file with the template spectrum to the first position. In the arrays created by IDL, the first position (actually zero-th position in IDL) will be used for results relating to the template.

Edit file star-name.lst by putting the full name of the FITS file with the template spectrum in the first line, to be used as a template. This spectrum must be in the main directory, not in /std.

5. It is practical to process all standard star spectra separately. This may be done in the sub-dir /std, usually with the same
template as this will simplify finding of the final "night corrections". It will be your choice, if you will process the std star spectra first, or after the program star. If you process the std stars first, you will have more freedom in selecting a good template. However, if you see - from the beginning - a good template spectrum, you can use it to process the program star first and then correct all the velocities later, on the basis of the subsequent analysis of the standards.

6. Processing of the template spectrum will impact on the rest of  analysis, because it will produce several auxiliary vectors. The most important is to establish the span of the re-binned spectra. Thus, you must check the selection of the parameters
w00, n, stepV (which are used to form the log-wave vector w1) in BFpro1.pro.

You can run an auxiliary routine BFpro4.pro which simply compares the ends of the new vector w1 with the ends of  individual spectra. w1 has unequal (power-law distributed) spacing between the wavelengths. Nothing is computed here; only the ends are checked.
Typically the input line for BFpro4.pro will look like:
here 5079.7 is the start of the new vector which will have 1018 elements corresponding to the velocity step of 11.8 km/s (the old system at 5184A).

Other settings:
For the JY2 at 5184A and the 1800 grating, use 8.19 km/s as the sampling interval:
(avoid a line at the left edge at 5040A).

For the JY2 at 5184A and the 2160 grating, use 6.7 km/s instead of 8.19.

JY2 at 6290A, from 6165 to 6432A::
This choice removes the short-end vignetting problem and retains the whole spectrum almost to the other end. At this point disregard the presence of the telluric band.

Note: BFpro*.pr and all other IDL routines must be always compiled first using .run or .compile.


Run the IDL program called BFpro1.pro to prepare the template spectrum for derivation of BF's from all program spectra. The template spectrum must be in the current IDL directory.

The routine now permits removal of a section of the spectrum which will not be used. This was developed for removal of the telluric lines, but may be useful for more general use.

The vector w1 keeps the wavelength log-lambda scale with any gaps required to remove some sections of the spectra. But do not remove intrinsically wide lines, like say those of hydrogen in A-type stars: they carry as good information about the BF's as the rest and this is the big advantage of the BF technique over the CCF. So the "blank" vector should be used only when you know very well that one section of all spectra is bad.
Typical use:
where the 2-element vector with zeros simply means that nothing is removed.

For the geocentric correction determination from the telluric band 6287-6330A (the strong band head within 6275-6285A is not used):
Note "t" added to names to signify "telluric".

For the stellar BF, use two spliced sections 6165-6275A (110A) and 6330-6430 (100A):
You may have an old version of this routine. The old BFPro1.pro routine was used as:
the new routine is:
where blank=[0.,0.]

Note: We have found that the 5194A feature is still the best spectral region, because (1) the JY2 chip is more sensitive towards blue, (2) the excellent new 2160K grating (in place of 1800F) does not work beyond 6500A and is already poor at 6000A. Thus blanking is usually not needed.

std_fts  =  'rC0048807.fits' = the name of the template file;
w00      = 5081.0            = start of the log-wavelength vector;
                               make sure that it is a bit (about 1-2A) longerward of the
                               lowest wavelength in the original spectrum; the
                               program will inform you about that;
n        = 1020              = number of pixels in the re-binned spectrum, must be EVEN,
         = 2000                and approximately similar to the original number;
                               check also the long end, also some 3A shorter
                               than the end of spectra - the program will inform;
m        = 121               = number of pixels in the broadening function,
         = 151                 must be ODD, a bit of balancing act here:
                               you want to have as short a BF as possible, yet you
                               don't want to have the baseline
around the BF poorly                                defined; because of that, you may have to re-run
                               Bfpro1.pro later on,
after determination of the BF's with
                               Bfpro2.pro and Bfpro3.pro. CCF may help in
                               selecting the value of m.
stepV    = 11.8              = the velocity step in km/s;
         = 6.7                 This is specific to the instrument. Select a
step a bit
                               smaller than that corresponding to
the original spacing
                               of your wavelengths. Calculate

                               basically the speed of light/resolution, then round down
                               a bit. See the routine itself.
                               stepV=11.8 km/s for the old system (19um/pix) and the
                               1800 grating at 5184A;
                               stepV=6.7 km/s for the
JY2 chip (13.5um/pix) and the
                               2160 l/mm grating at 5184A;

blank =[6272.,6330.]           added molecular feature or section removal; not needed
                               for the 5184A region.


w1                           = the log-wave vector, to be used in Bfpro2.pro.
des                          = design array, to be used in Bfpro2.pro
ww,u,v                       = results of the Singular Value Decomposition of des,
                               to be used in
vel                          = velocity vector for future use.

The routine sends out some messages and plots the spectrum first and then the vector of singular values using the lin-log plot:  plot_io,ww

BFpro1.pro requires the following routines:
set_win.pro (a routine to switch back to a graph Windows window; it can be commented out if you always have another
graphic window
interpol.pro (IDL distribution)
svdc.pro (IDL distribution)

Concerning the selection of "m", the size of the BF: If you have no idea about the width of the BF, try the cross correlation function (CCF). The CCF is a poor substitute of the BF, but gives some idea about the necessary width. Read about that in the "Cookbook" itself (which has been used to call this text), Sec.3.

The smaller "m" is, the better. But it cannot be too small as you will not have enough baseline around the stellar feature. Use your judgment and be prepared to re-run your analysis several times with various values of "m".


Note, that the original DDO spectra are not in the heliocentric wavelengths: The centroid velocities will come in the geocentric system relative to the template. The velocities must be corrected for the heliocentric correction (for each program spectrum and for the template spectrum). Then, the velocity of the template must be added to reduce to the heliocentric absolute system.

If the original IRAF spectra are already in the heliocentric wavelengths, then the BF's will be relative to the standard in the heliocentric system.

C.1. Determination of time, phases and RV helio corrections.

To determine the heliocentric JD's and heliocentric phases (with assumed: T0,period), use: hjd_phase.pro

It is likely that T0 will be changed later in the stage of the orbit determination. Thus, the first set of the phases will be preliminary. When you find the new T0, return to this step. It is best to adjust T0 in such a way that it corresponds to the primary (deeper) eclipse.

For standards:


prg_lst                        = list of FITS files (same as in A4 above)
t0, period                     = T0 (initial epoch) and period,
                                 MUST be in double precision,
                                 eg: 2452607.5828d0, 0.3207186d0

hjd                            = helio JD of observation (calculated for centre of
                                 each exposure)
phase                          = phase of observation (as above)

hjd_phase.pro requires the following routines:

C.2. Determination of the heliocentric velocity corrections

Calculation of the heliocentric and rotational velocity correction using: hjd_vel.pro. Now, with the Earth rotation correction included, it should be good to a few m/sec.
For standards:


prg_lst                  = list of FITS files (same as in A4 above)
lon                      = longitude in degrees, double prec (DDO given above)
lat                      = latitude in degrees, double prec (DDO given above)

hjd                      = helio JD of observation, re-calculated here,
                           should come out
the same as from hjd_phase.pro
hvc                      = radial velocity heliocentric correction

hjd_vel.pro uses the same routines as hjd_phase.pro plus the IDL-astro routine baryvel.pro

To see which phases are well covered us the plot:
Such a plot is very useful when you will do measurements of the RV's.

You can use the above routines for standards as well. Then use any numbers for T0, P, as the phases will not have any meaning. The important quantities to evaluate will be hjd, hvc.

You may have to return to the above steps if you reject any observations in step D (BFpro2.pro), as below, with the modified list of spectra. This is to have all data correctly identified by same indices in vectors and arrays. The two routines run fast to run, just remember to do that...

Also, the final binary orbital solutions will probably result in adjustment of T0. This is another reason to run hjd_phase.pro again.

C.3. Geocentric to heliocentric corrections of velocities.

The radial velocities are normally measured by Gaussian fits but better results may come from the use of the rotational profile fits (this will be described in Sec.F). In a typical case this will be in the geocentric system. The results will be in the form of two vectors with measured RV's: rvpm and rvps.

If the measured RV's are rvpm (primary) and rvsm (secondary), they must be corrected in the following way:
rvp = rvpm - hvc[0] + hvc + RVstd
rvs = rvsm - hvc[0] + hvc + RVstd

That this assumes that the zero-th element of the vector hvc contains the heliocentric RV correction for the template. RVstd must be taken from the literature or catalogues.

For standards, later on, after the measurements, create a vector
rvexp=[.., .., ]
with the catalogue RV's and plot:
This will give you information about consistency of the data.


This routine is the workhorse of the approach. It permits rejection of poor spectra and an abort in the middle when things do not work out and you are confused or distracted. Each spectrum can be de-spiked for cosmic rays. If you have very good spectra, just keep on pressing the Return. The current spectrum with BF being determined is displayed.

If you this routine has been used to de-spike and read-in spectra, one can use an auxiliary routine BFpro2b.pro can be used to obtain a series of BF's with a different template or a different range of velocities. See the end of this section.

Typical input:
if no spectra expected to be removed, then you can name a new file as junk for later removal:
if a telluric feature is to be separately found:
if standards:

prg_lst = 'V2357Oph.lst' = list of FITS files, prepared as described above
                           in A.4.
w1                       = the log-wave vector, created with BFpro1.pro
ww                       = singular values determined by BFpro1.pro,
size of ww sets the length of the BF's
u,v                      = auxiliary arrays from the SVD step BFpro1.pro

prg_lst1 = 'V2357Oph1.lst' = a new list of FITS files which are acceptable
further processing; if you are sure that
spectra are very good, call it 'junk.lst'
discard later.
images                   = string array with names of files
spec                     = spectra, not used, but read in, just in case;
[Example: to see a 10th spectrum: plot,w1,spec[10,*]                  ]
bf                       = array of broadening functions, each bf[i,*]
                           can be plotted versus vel,
[eg.: plot,vel,bf[10,*],   but will - normally - look horrible without

Routines used, in addition to those used by BFpro1.pro:

If this routine gives many messages that the wavelength vector w1 is too long, it may be necessary to re-run BFpro1.pro with smaller number of pixels or change the velocity step. The best is to run BFpro4.pro before anything to avoid this problem. BFpro4 produces no variables, just the screen output. It can be run anytime.

After this step, the broadening functions will be ready, but will look rather terrible. One must run the next step to be sure that the shapes are well defined. If the BF's do not have well-defined baselines or, conversely, if they take very little of the velocity space (when K1+K2 are small and the BF definition would be better with a smaller m), it may be necessary to re-do the whole process setting m lower or higher. This would mean going back to step B.

Re-processing of spectra with a different template (BFpro2b.pro).

One can re-process the already de-spiked spectra stored in variable spec by using this auxiliary routine.
Typical input:
Note, that two variables, [w1,spec] are kept the same, but you can change the rest. The variables [des,ww,u,v,vel,bf_n] are new, so you may want to call them differently. Consult the header of this routine for detailed instructions.


The broadening functions derived as in step D. usually look pretty terrible as their individual points are treated as independent variables. They do not  "know" that there was a slit whose image usually corresponds to some 2 - 5 pixels (the optimally designed spectrograph provides 2 pixels per resolution element, but 2.5-3 may be a bit safer to allow for optical imperfections). Thus, we must re-introduce the smoothing. The details for our two systems are rationalized at the end of this section.

For sigma = 2.0.
Note that you may wish to run this routine a few times for a few values of sigma.

bf                        = the 2-D array of BF's from BFpro2.pro
bfss                      = smooth version of bf

You may call them by the various values of sigma:

Typically, we use 1.5 for the old CCD system and 2.0 for the new CCD system.

The old CCD system:
The Cassegrain spectrograph on the 1.88m used a Thomson CCD with 19 micron pixels. With the internal "minification" in the spectrograph of 4-times, the 300 micron slit is projected into a feature of FWHM~75 microns or about 4 pixels. If these were a Gaussian profile (FHWM=2.35*sigma), then the appropriate
sigma ~ 1.5 pix (bf15). For faint stars with poor S/N, bf20 may have to be used.

The new CCD system:
The pixels are smaller, 13.5 microns. The 250um slit gives images 62.5um wide, or 4.62 pix, hence the appropriate sigma=2.0 pix. Therefore use bf20 or bf25.

Note: If you re-binned the spectra to a lower resolution, and the re-binning was done correctly (preferably with some Gaussian smoothing first), then you may have already introduced some smoothing.



The velocities are determined interactively through a somewhat time consuming "manual" process. This could be done automatic, but this way one fully monitors the fitting process.

One can fit Gaussian or rotational profiles to the double peaked BF's. As pointed by Theo Pribulla in the Spring 2005, the latter is advantageous and gives better results.


Plot a few broadening functions to see how they look like. Typically:
i=1 & plot,vel,bf15[i,*] (i=0, std itself; i=1, 2, ... binary)
for poorer spectra, use:
i=1 & plot,vel,bf20[i,*]

Fit a 7-parameter function to each broadening function: 2x3 parameters of two Gaussians plus an arbitrary vertical offset. The parameters are:
a[0] & a[3]      = strengths of both components
                   (the BF will be displayed
normalized to max,
                   so numbers like 1.0 to 0.3.

a[1] & a[4]      = positions in velocity scale,
                   these are the numbers we are after!

a[2] & a[5]      = widths of the Gaussians; use two_gs2.pro
                   which is a version of the more general
                   two_gs.pro, but
with the widths fixed -
                   otherwise the routine can give
wrong results.
a[6]             = a vertical offset, not interesting,
                   but keep floating.

When two_gs2.pro is used, the fit is 5-parameter, and is more stable. All the routines two_gs* utilize simple LSQ iteration with derivatives of the Gaussian computed analytically.

The measurement instructions:
rvpm=fltarr(...) ; as many as there are phases
You may wish to declare weights:

A typical set:
a=[1, -50, 120, 0.5, 250, 80, 0]
; str1 pos1 wid1 str2 pos2 wid2 baseline

; when the primary is at left, i.e. more negative RV:
a=[1,-50,120,0.5,250,80,0] & da=0.
; when the primary is at right:
a=[1,50,120,0.5,-250,80,0] & da=0.

Now follow the cycle for the Gaussian fits ->
i=0   is the template itself, not interesting
bbb=reform(bf15[i,*]) & bbb=bbb/max(bbb) ; once per measurement
(use the appropriate bf* in place of bf15.)

Recall with up/down arrows & repeat the line below until all da corrections converge to zero (make sure that the line below does NOT wrap)
; once the differences stabilize, save the results in
rvpm[i]=a[1] & rvsm[i]=a[4]
wtp[i]=1. & wts[i]=1.
Set any of these to =0.5 for poor observations or =0.25 for very poor observations
<- end of the manual cycle

The weigths will be used in the RV orbit solutions, see section H.

It is important to start with a good set of the parameters a. Starting from a wrong set may give you entirely wrong iterated solution.

For weak companions, the Gaussian fits may go wrong as the 2-Gauss fitting may not "feel" the second star well. Use the cursor then to measure the RV in a simplified (but secure) way:

It is useful to identify a nice, double BF for a quadrature and measure it first. Then, for spectra observed close in time, the first result will be a good approximation for the neighbouring phases.

Rotational profiles:

This is a very similar approach generally to that of Gaussians. Also 7 parameters, but the widths seem to be better constrained than for the Gaussian fitting.

Sometimes, when the secondary is faint, the solutions are not satisfactory and the fitted profile is systematically shifted from what you see.. Then try to use the version with the fixed widths (Rot_two1.pro instead of Rot_two.pro). If you still get questionable results for the secondary, use the cursor measurements for less biased results, as described a few lines above here.

Generally, be very careful with widths. Do not expect them to converge to the right values automatically. It may be necessary to use the fixed width version for all of your BF's. The rotational profile fitting is a bit less forgiving than the Gaussian fitting, but when it works, it gives better results.

Normal fits:
Declare rvpm, rvsm as above and set values a the same way:
for as many as there are phases.
You may wish to declare weights:

A typical set:
a=[1, -50, 120, 0.5, 250, 80, 0]
; str1 pos1 wid1 str2 pos2 wid2 baseline

; primary at left:
a=[1,-50,120,0.5,250,80,0] & da=0.
; primary at right:
a=[1,50,120,0.5,-250,80,0] & da=0.

Cycle Rotational fitting->
i=1 ; any i>1, since i=0 is the template itself
bbb=reform(bf15[i,*]) & bbb=bbb/max(bbb) ; once per measurement
; use the appropriate bf*

Recall & repeat the line below until all corrections converge to zero (make sure that the line does NOT wrap)
              [   rot_two1  (for fixed widths)   ]
; once the differences stabilize, save the results in
rvpm[i]=a[1] & rvsm[i]=a[4]

wtp[i]=1. & wts[i]=1.
; set =0.5 for poor observations
; =0.25 for very poor observations

<- end of rotational fitting cycle

The end-product of the above measurements will be the measured velocities of the primary and secondary components, in the geocentric system:
rvpm( )
rvsm( )
for as many points as there are phases.

(Here repeat of the text from Sec. C above: The radial velocities are in the geocentric system. They must be
corrected in the following way:
rvp = rvpm - hvc[0] + hvc + RVstd
rvs = rvsm - hvc[0] + hvc + RVstd
That this assumes that the zero-th element of the vector hvc contains the heliocentric RV correction for the template. RVstd must be taken from the literature or catalogues.)


It may happen that the BF's show only one component. For measurements, use: one_gs.pro. Usually there is no problem with fitting all 4 parameters (strength, position, width, offset or zero-point) of the Gaussian. For especially stubborn cases, fix the width by using one_gs1.pro in the cycle below.

It may happen that after reductions of the visible component, a second look on the BF's at phases 0.25 and 0.75 will show weak signatures of a secondary component, partly buried in the BF noise. Then re-do the analysis as for two stars.

a=[1, 0, 100, 0] & da=0.
str pos width base
a=[1,0,30,0] & da=0.
for a standard

Declare for as many as phases or std star observations:
rvm=fltarr( )
rvm_s=fltarr( )

Measurement cycle ->
; change as needed
or for a standard

; recall & repeat the line below until all da=0
; no wrap of the above line
<- end cycle

In place of the Gaussian fits, one can use the rotational profiles:

If the above is a standard, give the new vectors some special names like rvm_std.

Later on:
rv_s = rvm_s - hvcs[0] + hvcs + RVstd

To check consistency of the RV std's:
where RV_std are the values from catalogues


This assumes the most typical case: the 3rd component shows a sharp spectral feature in the BF, due to its slow rotation (compared with the huge rotational velocities of components in a close binary). So far, this is what we have always seen.

You may want to establish the presence of the third component by using a contour plot to see if there is something stable in the BF's. The heliocentric system of velocities is the proper one here, but usually one sees the 3rd component also in geocentric system (if the dates of observations are close). Use:

w=sort(phase) ; sorts indices in phase

Possible changes, eg: ..., levels=[0.05,0.5,0.9,0.95]

Basically, the approach for 3-component BF's is similar to that for close binaries. However, because we want to determine the best velocities of the close pair, we subtract the 3rd body signature and measure the BF for the close pair itself.

Thus, three stages are needed:
1. Fitting of three Gaussian components, two (for the close binary) with fixed widths, and the 3rd with all parameters determined. This is going to be the only information of velocities of the 3rd body. The individual velocities of the close pair components will be discarded as this will be re-done in step 3.
2. Removal (subtraction) of the 3rd (sharp) component from the combined broadening function;
3. Measurement of the close binary velocities, as for a system of only two stars.

Typical set for three_gs2.pro:
a=[1,50,120, 0.5,-250,80, 1,-100,30, 0]

Define an array of measured data for the 3rd component:
rv3m=fltarr(3, )
Define an array of BF with subtracted 3rd component:
(can be bf20, and here is just a simple substitute declaration)

Work in the same way as for binaries, but then store the data for the 3rd component in rv3m.

Cycle ->
i=1 (change as needed)

a=a+da&y=three_gs2(vel,bbb,a,da)&plot,vel,bbb&oplot,vel,y,line=2&print,form='(10f8.2/10f8.4)',a,da ; no wrap of theses lines

For phase index i:
rv3m[*,i] = a[6:8]
; stores parameters for the 3rd comp

Subtract the 3rd comp:
<- end of cycle.

You can do the subtraction later:
Cycle ->
<- end of cycle

Heliocentric velocities of the 3rd comp:
rv3 = reform(rv3m[1,*]) - hvc[0] + hvc + RVstd

Analyze bf15_2 as for binary systems (see above).

Arguably, velocities of the 3rd body will be more affected by the blending with the close binary than in reverse, but this may not be true when the close system has a very small mass ratio so that its primary moves little on its orbit. Then the BF signature of the primary will almost always merge with that of the 3rd component.

It may be useful to plot the velocities of the 3rd component versus the phase of the close binary. This way any "cross-talk" could be detected. It will be present at some level. Analyse for systematic temporal trends.

The light contribution: L3/(L1+L2)

To find the light contribution of the 3rd component in the system, use the interactive routines: BFanal3.pro or BFanal3x.pro. They give systematically different results, with BFanal3x.pro probably more trustworthy. Take the difference as the level of uncertainty.


Repeat for as many phases as you clearly see the 3rd component, usually at quadratures, when both binary components are well visible. The scatter in L3/(L1+L2) is typically 0.02 - 0.05.

Note that the ratio L3/(L1+L2)  may be affected by the degree of broadening of the lines. If one is confronted with a (rather typical) situation of a  WUMa plus a sharp-line star, the  pseudo-continua are different for both stars, so that the ratio of luminosity will be also affected.


Make sure that the data do not require any "night corrections". Frequently we use just one observation of the template star. It is essential to check how this particular observation relates to other observations of standard stars and if it is a representative one, without a shift which would affect the whole series of observations. It may have happened that you selected a template spectrum with the star being guided on one side of the spectrograph slit.

Create an .lst file for the standard star spectra and run through all steps of the RV derivation. You may see deviations in radial velocities on individual nights. Or you may have a global shift due to the template observation selected in such a way that there is a constant shift for all other std observations; in the latter case all radial velocities should be shifted by some value, eg.

- create "std.lst" the same way as for the program star, with the template observation in front (zero-index position), it may be convenient to do that in the sub-dir "std" in your star's dir,
- run BFpro2.pro and BFpro3.pro programs,
- run:
(Note that you must calculate the phone phases in order to obtain the HelioJD for standards (hjd))
- measure the standards the same way as single-component binaries (Sec.F.2. above), create a vector rvm (substitute elements: rvm[i]=a[1])
- insert the corrections
rv = rvm - hvc[0] + hvc + RVstd
where RVstd is the catalogue radial velocity for the template

To find out of the standards really were in the common system:
create a vector with all catalogue data: rvexp=[.., .., ]
with the catalogue RV's and plot:

The vector rv should give values close to the catalogue ones. If not, find the corrections and appropriately correct your program star observations:
- create a vector ncor of the same length as rv (or rvp & rvs)
ncor = fltarr( )
- insert nightly corrections (the difference between the expected
velocity and observed), like say: ncorr[20:*]=2.5
- rvp = rvp + ncorr
- rvs = rvs + ncorr

H. Binary star RV orbits

Determination of V0,K1,K2,T0, where V0 is determined from both stars and Ki are semi-amplitudes which correspond to sinusoidal variations. T0 is the moment of conjunction of the primary star behind (i.e. the primary eclipse).

Use three routines as described below:

H1. Preparations.

You must form (this is IMPORTANT!) a double precision, 5 column array: xx(5,n). n is the number of phases for the binary itself (skip the 0th element containing the template RV), so this is by one element shorter than vectors discussed above.The columns are:
hjd, rvp, weightp, rvs, weights

xx should be selected as a characteristic, abbreviated name of the star; use anything like vw or v356, or anything which you will easily associate with a given star.

xx=dblarr(5,n) ; double to keep correctly the many-digit JD number
; n=number of BF's minus one
These columns 1 and 3 must be corrected for "night corrections", if need be.

Fill in the weights. If all the same:
You can edit some later by hand, eg. xx[2,30]=0.5
Or you can use the weights determined during the RV measurement process:

Prepare the 2-el. vector with the ephemeride. IMPORTANT: Enter both numbers as double precision, i.e. add d0 to the value entered or multiply by 1.d0.

H2. The first solution

This determines the first approximate V0,K1,K2


This will be an 8-el vector, first 4 are V0,K1,K2,place-holder; the second 4 are their errors

H3. The corrections

This determines corrections to V0,K1,K2,T0

The same format as above so that the first 4 can be added to give final elements:
Then disregard the last four results of the sum.

Make a note of the sigmas (errors per observation) printed on the screen.

The T0 correction must be a small number (say <0.02 day) for the routine to work well. It is better to adjust the phases by hand, by changing T0 or you may like to find T0 through successive corrections derived by two_star2.pro (see below). In any case, correct:
T0 = T0 + xx2[3],
and substitite into
as above.
Repeat the first solution with two_star1.pro and then find the corrections using two_star2.pro.

Usually, one such iteration is sufficient. The sigmas given by two_star2.pro will be smaller when this is done - so, for other corrections to be accurate, you must have a very small value of the T0 correction.

If you like, you can use the plotting routine plot_2star.pro in two ways:
plot_2star,xx,xx0,xx1,'XX Pho'
plot_2star,xx,xx0,xx1+xx2,'XX Pho corrected'

Remember to run again:
to update the phases (Section C).

H4. Errors

The errors from the two routines two_star*.pro above are usually unrealistic, because they come from the linearized LSQ. To obtain better estimates, use:
It runs a minute or so.

The 3 rows give the two 1-sigma ranges of elements and the median values estimated from the bootstrap sampling experiment. Strictly speaking, one should use xx1+xx2 in place of xx1, but this makes almost no difference when xx1 are the best possible. Note that the median values are usually different from the best xx1+xx2 which indicates skewness in the error distributions. However, we use two_star3.pro to obtain error estimates only and nothing more.

The best final elements: V0,K1,K2,delT, are the first 4 elements of xx1+xx2. The best final 1-sigma errors are from a simple symmetric approximation:

The routines for single stars (one_star*.pro) are similar, but not identical. One must create a 3-column array of data, say:
star, (star=[hjd,rv,weight]; star[3,*])
and process it with
where as before, star=[t0,period], and star1, star2, star3 are results.


This is a murky area. We interpreted some problems as flexure effects, but we are no longer sure if this should be our main concern. In any case, to evaluate the flexure corrections, one should know the hour angle, HA. The hour angle and airmass can be computed with:

For DDO:

INPUT: prg_lst,lon,lat
prg_lst  = the list of FITS files;
lon, lat = coordinates of the observatory (note the longitude sign convention as for DDO).

OUTPUT: hjd, ha, am.
HJD is re-calculated again, should be the same from other routines;
HA = hour angle, in degrees;
AM = the airmass.


A set of routines BFimage*.pro can be used to show BF's in two dimensions, with the phase as the X-axis and the heliocentric velocity as the Y-axis. The BF's form "intensities" in these two dimensions. Weak stellar components can be detected that way as such images re-introduce the temporal correlation between individual BF's.

The routines BFimage*.pro have their own comments, so please analyze these.

This is all, folks. Have a good time...