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.
FIRST READ THESE (HOPEFULLY) USEFUL
NOTES:
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.
Content:
GENERALITIES
AND HISTORY
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.
Methodology
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:
rectif_fits,file_name
as in:
rectif_fits,'K0003456.fits'
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.
IMPORTANT:
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.
In UNIX:
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:
BFpro4,'V523Cas.lst',5079.7,1018,11.8
====
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:
BFpro4,'list.txt',5050.0,2000,8.19
(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::
BFPro4,'V523Cas.lst',6165.0,1946,6.5
====
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:
BFpro1,'K0004000.fits',5079.7,2000,151,6.7,[0.,0.],w1,des,ww,u,v,vel
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):
====
BFpro1,'K0004000.fits',6285.,328,21,6.5,[0.,0.],w1t,dest,wwt,ut,vt,velt
====
Note "t" added to names to signify "telluric".
For the stellar BF, use two spliced sections 6165-6275A (110A) and
6330-6430 (100A):
====
BFpro1,'K0004000.fits',6165.,1946,151,6.5,[6275.,6330.],w1,des,ww,u,v,vel
====
You may have an old version of this
routine. The old BFPro1.pro routine was used as:
BFpro1,std_fts,w00,n,m,stepV,w1,des,ww,u,v,vel
the new routine is:
BFpro1,std_fts,w00,n,m,stepV,blank,w1,des,ww,u,v,vel
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.
INPUT:
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 using:
3e5/((longest-wavelenght)/(step-in-wavelength));
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.
OUTPUT:
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 Bfpro2.pro
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:
readfits.pro
parse_fits_head.pro
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)
map4.pro
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.
===
Typically:
hjd_phase,'v401cyg.lst',t0,period,hjd,phase
===
For standards:
hjd_phase,'std.lst',2448500.d0,1.d0,hjds,phony
===
hjd_phase,prg_lst,T0,period,hjd,phase
INPUT:
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
OUTPUT:
hjd
= helio JD of observation (calculated for centre of
each exposure)
phase
= phase of observation (as above)
hjd_phase.pro requires the following routines:
readfits.pro
parse_fits_head.pro
helio_jd.pro
jdtime.pro
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.
===
Typically:
hjd_vel,'v401cyg.lst',-79.421667d0,43.863333d0,hjd,hvc
===
For standards:
hjd_vel,'std.lst',-79.421667d0,43.863333d0,hjds,hvcs
===
hjd_vel,prg_lst,lon,lat,hjd,hvc
INPUT:
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)
OUTPUT:
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:
plot,phase,psym=1,yran=[0,1]
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:
plot,rv-rvexp,psym=1
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:
BFpro2,'V2357Oph.lst','V2357Oph1.lst',w1,ww,u,v,images,spec,bf
===
if no spectra expected to be removed, then you can name a new file as
junk for later removal:
BFpro2,'V2357Oph.lst','junk.lst',w1,ww,u,v,images,spec,bf
===
if a telluric feature is to be separately found:
BFpro2,'V2357Oph.lst','junk.lst',w1t,wwt,ut,vt,imagest,spect,bft
===
if standards:
BFpro2,'std.txt','std1.txt',w1,ww,u,v,img_s,spec_s,bf_s
===
INPUT:
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,
the size of ww sets the
length of the BF's
u,v
= auxiliary arrays from the SVD step BFpro1.pro
OUTPUT:
prg_lst1 = 'V2357Oph1.lst' = a new list of FITS
files which are acceptable
for further processing;
if you are sure that
the spectra are very good,
call it 'junk.lst'
and 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
smoothing].
Routines used, in addition to those used by BFpro1.pro:
trunc.prot
svsol.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.
BFpro2b,std_fts,m,stepV,w1,spec,des,ww,u,v,vel,bf
============
Typical input:
BFpro2b,'K0003000.fits',111,6.5,w1,spec,des,ww,u,v,vel,bf_n
============
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.
===
Typically:
bf20=BFpro3(bf,2.0)
For sigma = 2.0.
===
Note that you may wish to run this routine a few times for a few values
of sigma.
INPUT:
bf
= the 2-D array of BF's from BFpro2.pro
OUTPUT:
bfss
= smooth version of bf
You may call them by the various values of sigma:
bf10,bf15,bf20,bf25,bf30
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.
F.1. BINARY STARS, DOUBLE PEAKS
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.
Gaussians:
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:
Declare:
rvpm=fltarr(...) ; as
many as there are phases
rvsm=fltarr(...)
You may wish to declare weights:
wtp=fltarr(...)
wts=fltarr(...)
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=1
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)
a=a+da&y=two_gs2(vel,bbb,a,da)&plot,vel,bbb&oplot,vel,y,line=2&print,form='(7f10.4)',a,da
; 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:
plot,vel,bf15[i,*]
cursor,x,y,/data
print,x
rvsm[i]=x
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:
rvpm=fltarr(...)
rvsm=fltarr(...)
for as many as there are phases.
You may wish to declare weights:
wtp=fltarr(...)
wts=fltarr(...)
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)
a=a+da&y=rot_two(vel,bbb,a,da)&plot,vel,bbb&oplot,vel,y,line=2&print,form='(7f10.4)',a,da
[ 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.)
F2. SINGLE STARS OR ONLY ONE-COMPONENT
VISIBLE (SB1) OR STANDARDS
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 ->
i=1
; change as needed
bbb=reform(bf15[i,*])&bbb=bbb/max(bbb)
or for a standard
bbb=reform(bf15_s[i,*])&bbb=bbb/max(bbb)
; recall & repeat the line below until all da=0
a=a+da&y=one_gs(vel,bbb,a,da)&plot,vel,bbb&oplot,vel,y,line=2&print,form='(4f10.4)',a,da
; no wrap of the above line
rvm[i]=a[1]
or
rvm_s[i]=a[1]
<- end cycle
In place of the Gaussian fits, one can use the rotational profiles:
a=a+da&y=rot_one(vel,bbb,a,da)&plot,vel,bbb&oplot,vel,y,line=2&print,form='(4f10.4)',a,da
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:
plot,rv_s-RV_std,psym=1
where RV_std are the values from catalogues
F3. TRIPLE SYSTEMS
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
contour,bf15[w,*],levels=findgen(10)/10+0.1,/fill
or
contour,bf15[w,*],phase[w],vel,lev=[....]
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:
bf15_2=bf15
(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)
bbb=reform(bf15[i,*])&bbb=bbb/max(bbb)
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:
bf15_2[i,*]=bbb-a[6]*exp(-(a[7]-vel)^2/a[8]^2)
<- end of cycle.
You can do the subtraction later:
Cycle ->
i=1
bbb=reform(bf15[i,*])&bbb=bbb/max(bbb)
bf15_2[i,*]=bbb-rv3m[0,i]*exp(-(rv3m[1,i]-vel)^2/rv3m[2,i]^2)
<- 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.
Typically:
BFanal3,vel,bf15[20,*]
BFanal3x,vel,bf15[20,*],bf15_2[20,*]
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.
rvp[1:*]=rvp[1:*]+2.5
rvs[1:*]=rvs[1:*]+2.5
Therefore:
- 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:
hjd_phase,'std.lst',2448500.do,1.d0,hjd,phony
hjd_vel,'std.lst',hjd,hvc
(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:
plot,rv-rvexp,psym=1
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
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:
two_star1.pro
two_star2.pro
two_star3.pro
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
xx[0,*]=hjd[1:*]
; n=number of BF's minus one
xx[1,*]=rvp[1:*]
xx[3,*]=rvs[1:*]
These columns 1 and 3 must be corrected for "night corrections", if
need be.
Fill in the weights. If all the same:
xx[2,*]=replicate(1.0,n)
xx[4,*]=replicate(1.0,n)
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:
xx[2,*]=wtp[1:*]
xx[4,*]=wts[1:*]
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.
xx0=[T0,period]
H2. The first solution
two_star1,xx,xx0,xx1
This determines the first approximate V0,K1,K2
print,xx1
This will be an 8-el vector, first 4 are V0,K1,K2,place-holder; the
second 4 are their errors
H3. The corrections
two_star2,xx,xx0,xx1,xx2
This determines corrections to V0,K1,K2,T0
print,xx2
The same format as above so that the first 4 can be added to give final
elements:
final=xx1+xx2
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
xx0=[T0,period]
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:
hjd_phase,'star.lst',t0,period,hjd,phase
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:
two_star3,xx,xx0,xx1,xx3
It runs a minute or so.
print,xx3
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:
print,(xx3[*,2]-xx3[*,0])/2
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
one_star1,star,star0,star1
one_star2,star,star0,star1,star2
one_star3,star,star0,star1,star3
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:
hjd_ha,prg_lst,lon,lat,jdtime,ha,am
For DDO:
hjd_ha,'list.txt',-79.421667d0,43.863333d0,hjd,ha,am
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...
Slavek