% @(#)user.tex	19.1 (ESO-IPG) 02/25/03  13:28:26
%++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
%.COPYRIGHT (c)  1992  European Southern Observatory
%.IDENT     user.tex
%.LANGUAGE  LaTeX
%.AUTHOR    Andrew T. Young, ESO/IPG
%.PURPOSE   Documentation for PEPSYS package
%.VERSION   5.4  1993-JULY-03 : Modified, ATY
%---------------------------------------------------------------------
\documentstyle[11pt,twoside]{report}	% remove ^% to print
\input ../../../doc/user/manual
\hyphenation{your}
\begin{document}			% remove ^% to print

\day=31  \month=1  \year=1993          % release date for chapter

\setcounter{page}{1}
\chapter{PEPSYS general photometry package}
This chapter describes the PEPSYS package for general photometric reductions.
The function of photometric reductions is to measure and remove, so
far as possible, the instrumental signature from the data.
If we regard the Earth's atmosphere as a part of the
instrument not under the control of the observer,
we see that measuring and correcting for extinction is an important part of
this process.
PEPSYS performs
extinction correction and transformation to a standard system, if possible.

The package contains two main parts: a planning program and a reduction
program.
The planning program interacts with the user to find out the goals of the
observing program, and then produces a schedule of the observations needed
to meet those goals efficiently.
The reduction program is flexible enough to model most types of photometric
observations accurately, and produce reliable estimates of the uncertainties of
the parameter values it obtains from the raw data.

In addition, some auxiliary tools are provided to simplify the work of
making the required tables.

Because the data reduction must model the instrument and observational
procedure as accurately as possible, we must discuss photometric techniques.
The recent book by Sterken and Manfroid~\cite{SMbook} and
an earlier review~\cite{Y74} are good general references.
We also discuss the modelling process in some detail, to provide users with
some assurance that good models are used.

{\bf Please read through the documentation carefully before trying to use the
programs!}
Many problems can be avoided if you are thoroughly familiar with the contents
of this chapter.


\section{Introduction}
Good photometry requires careful attention to calibrations.
If there are no calibration observations of standard and extinction stars,
the program-star observations cannot be transformed to a standard system.
At the opposite extreme, if there are {\it only} observations of standard and
extinction stars, there is no time for program-star observations.
In either case, the data produce no useful information.
Obviously, somewhere in between these extremes lies an optimum distribution
of program and standard stars.

The choices of which standards to observe, and when to observe them, involve
not only balancing the gain in information by adding a calibration observation
against the loss of a program-star observation, but also assumptions about the
stability and linearity of the equipment used.
A good distribution of standards in each color index is also needed to
determine accurate transformations.
When all these details had to be attended to by hand, it was easy for observers
to neglect some essential item.
Much experience was required to distribute observations effectively.

Fortunately, the computer can be told to keep track of all these details,
and will not forget them.
The planning program helps you approach the optimal distribution of standards,
which gives the
most accurate and precise results possible for a given amount of observing
time.

\subsection{What is needed}

A basic principle of photometry goes back to Steinheil, who built the world's
first stellar photometer in the early 1830's.
Steinheil's principle is that only {\it similar} things should be compared.
That is, good results are obtained when as many variables as possible are
kept fixed.
This minimizes the number of instrumental parameters that have to
be calibrated,
and allows the necessary calibrations to be a small part of
the total data gathered.
All observations must be made with the same instrument, used
under the same conditions.

However, we cannot make all our observations at the same time, or at the same
place in the sky.
So we must be able to separate different effects that can vary with time, such
as extinction coefficients and instrumental zero points.
To do this, we must distribute the observations so that the main independent
variables (such as airmass, star color, and instrumental temperature)
are uncorrelated with each other.
In particular, they should all be uncorrelated with time, so far as possible.
As some (like temperature) are inherently likely to be correlated with time, it
is doubly important that others (like airmass) be uncorrelated.
Achieving this condition requires a certain amount of advance planning.

\subsection{How to get it}

To help you get the necessary calibration data in the least
observing time, the planning program generates a schedule
of required observations.
A considerable amount of photometric expertise is built into this program, so
that observers without much photometric experience can confidently adopt its
recommended defaults, and obtain a plan that should be satisfactory for
most photometric observing programs.
At the same time, the program is flexible enough to let experienced observers
design special-purpose programs.

This flexibility requires considerable interaction between the user and the
planning program.
The program needs information about the program stars, the calibration stars,
the peculiarities of the instrument, and the requirements of the observer.
This input information is described in detail below, and in the appendix on
data-file formats.

\subsection{What to do with it}

Once the observations are made, we must remove the instrumental signature,
including atmospheric effects, as fully as possible.
To use the observational data effectively,
the reduction program must make correct assumptions about the way the data were
gathered, and the way the instrument itself (including the atmosphere!)
really works.

The atmospheric part can be removed quite accurately, but some remaining
instrumental effects, due to a mismatch between the instrumental and standard
passbands, cannot be completely removed, because some information is missing
(especially in the conventional photometric systems).
Nevertheless, if the filters are quite close to the standard ones,
the missing information is fairly small, and good results are possible.

Again, interaction between the observer and the program is necessary, both to
obtain the necessary information and to decide how to proceed when choices are
not clear-cut.
Safe defaults are offered to the beginner; more experienced observers can
try various assumptions to see what works best for their purposes.
In both programs, the goal is to allow a wide variety of choices, but to
warn users of potential problems.


\section{Getting started}

The first step is to gather the information that will be needed to plan a
program and to reduce the data.
Information that is more or less permanent is stored in several MIDAS table
files:

\begin{itemize}

\item Star tables:  These give names and positions of stars to be used for
calibrations, as well as program stars.
Common supplementary information, such as spectral types, may be included.
Standard-star tables are provided for some popular systems, but you will
have to compile your own table of program stars.

\item Observatory table: This table contains the positions of various
telescopes, together with essential information such as apertures, and
subsidiary information.

\item Horizon table: This table describes the apparent horizon as seen from a
particular telescope.

\end{itemize}

All of the above information, except for program-star tables,
should already be available to you.
Only
if you are the first user of PEPSYS at an observatory will you need
to construct the Observatory and Horizon tables yourself.
They are described in the Appendix; see also section~\ref{installtbl} below.
The Observatory file is very simple, and data for
the Horizon file for a telescope can
be collected in a few hours
(the command {\tt MAKE/HORFORM} will help you compile the horizon data).
Program-star tables are described below, along with the {\tt MAKE/STARTABLE}
command provided to produce them.

In addition, you will need information about the instrument.
As instruments tend to evolve with time, you will probably need to put this
information together for each observing run separately.
The {\tt MAKE/PHOTOMETER} command will ask you questions, and make a
new file from your answers; it can
also show you the contents of an existing file.
If an instrument is fairly stable, and a previous instrument-configuration file
is available, you may be able to just edit the old file, using the
{\tt EDIT/TABLE} command.

Some of the instrumental information may not be available at the time you plan
your observing run.
You can give ``don't know'' responses to questions about such things,
or leave items blank if you have no information.
However, for the sake of accuracy, you should try to obtain as much
of the missing information as you can before reducing the observations.

\subsection{Star tables}
\label{startables}
The essential data needed for program stars are names, positions, and
the equinox to which the positions are referred.
You may also include proper motions and the epoch of the position; spectral
types; rough magnitudes, such as might be obtained from the HD; and comments.
The command {\tt MAKE/STARTABLE} helps you create a table of program stars.

Some telescopes record the position toward which the control system thinks the
telescope is pointing as part of the data stream, so it is tempting to
extract these data to make a positional catalog for the program stars.
In general, this is very unwise.
The telescopes usually used for photometry generally have neither
accurate nor precise pointing.

The {\it precision} (repeatability) of the recorded positions does not indicate
their {\it accuracy} --- which could be far worse, if the observer adjusted the
zero-points of the telescope's
coordinates to agree with mean places for some equinox
removed a few years in time from the date of observations.
A 2$'$ altitude error corresponds to an airmass error of 0.002 at 2 airmasses,
more than 0.005 at 3 airmasses, and 0.01 airmass at 4 airmasses.
Thus, such errors can cause appreciable systematic errors in extinction
determination and correction.
Although such recorded positions are available (by means of the
CONVERT/PHOT command) as a last resort, they are {\it not} recommended.
Good catalog positions should be used whenever they are available.

Standard-star tables also require columns for the standard values, which are
often not known for program stars.
Because of this difference in content, you should keep standard and program
stars in separate table files.
You can have several files of each type.
The programs will ask for the standard-star files first, and then the
program-star files.

You should also plan to keep extinction stars (and other constant objects, such
as comparison stars in variable-star programs) in a {\it separate} file from
variable stars.
The reduction program treats standard, constant, and variable stars
differently, so each group should be kept in one or more separate files
(see subsection~\ref{reductionproc}, ``Reduction procedure'').


\subsubsection{ASCII source files}
If you are lucky, your program stars are available in machine-readable form.
If so, simply copy out the data for your program stars
as fixed-format ASCII records, one star per record.
It doesn't matter if the ASCII file contains extraneous data; they can
be ignored in converting to MIDAS format.
The command {\tt MAKE/STARTABLE} helps you turn this flat ASCII file into a MIDAS
table of program stars.
It will first ask for the name you want to give your new program-star table.
Then it checks to make sure you have the necessary information, such as the
ASCII table of stellar data, and a format file.

It is important to make sure the star names in your source files
match those in your data files.
In building star files, keep in mind the need to have the {\it same}
designations appear in the observational data.
If your data-logging system makes you enter star names when they are observed,
try to keep names short to avoid typing errors.
That means using short, unique names in your program-star files, to match those
that will appear in the data.
If your data-logging system takes star names from files prepared in advance,
try to adhere to the standard IAU name format (see the
guidelines published in A\&A ``Indexes 1990 and
Thesaurus, Supplementary Issue, May, 1991'', pp.~A11-A13; PASP {\bf 102},
1231-1233, 1990; and elsewhere).

If possible, use at least two names for each star, as recommended by the IAU.
You can use two names in the star files, and just use one in data files; the
reduction program is smart enough to match them up properly, or will ask for
help if similar but not identical names occur.
It is a good idea to separate aliases with an ``='' sign; just leave a space
on either side of it, so the program doesn't take it as part of a name string.
Ordinary spacing is allowed in names to make things readable:
{\tt HR 8832}, {\tt Chi Cygni}, {\tt BD +4 4048}, etc.
Thus a name field might contain ``{\tt HD 24587 = CD -24 1945 = HR 1213}''.

Although you can use any naming system you like for program stars,
so long as the {\it same}
name appears in star files and data files, a system
of priorities is
suggested for making up catalogs of standard stars.
The basic principle is that small catalogs generally yield shorter names that
are easier to use than do bigger catalogs.
Generally, catalogs for brighter stars contain fewer entries; so the lists with
the brightest limiting magnitude are preferred.
Thus,
common designations like Bayer letters and Flamsteed numbers are
usually included for the brightest stars.
The bright stars are also listed by HR number,
and fainter ones by HD or DM number.
Don't forget that there is considerable overlap among the BD, CD, and CPD;
specify the catalog, not just a zone and number.
The HST Guide Star Catalog is recommended for still fainter objects.

If the reduction program finds stars with different names but nearly identical
positions, it will ask you if they are the same star.
Be prepared to answer such questions, if you enter names inconsistently.
Many users find that repeatedly answering such questions becomes tedious
and irritating; you can avoid this problem by using identical name strings in
both star files and data files.

Don't use names like {\tt STAR A} and {\tt STAR B},
as the matching algorithm will
spot the common word STAR and ask if they are the same; instead, just use
the letter.
If it is necessary to intermix several similar names, try to make unique
strings out of them.
For example, if you
are working on a group of clusters, and have local standards designated by
the same
letters in each, attach the letter to the field designation:
{\tt M67\_A}, {\tt M67\_B}, etc.

Although it is not recommended, the CONVERT/PHOT command can extract apparent
star positions from raw data files in ESO and Danish formats (see
section~\ref{startables}).
As intermediate output, this command will produce an ascii file that can be
edited.
Manual editing may be necessary if you do not use consistent naming conventions
while observing.

Finally, don't try to create star files with incomplete data.
Missing values cause problems when you try to reduce data.
Be sure every entry is filled in correctly.


\subsubsection{Format files}
If you have no format file ready,
{\tt MAKE/STARTABLE} copies a standard table-format template named
{\tt starfile.fmt}
into your working directory.
The ASCII format file
describes the input file by giving the column numbers (character positions)
of the input file in which each field (table column)
begins and ends; it describes the final display formats with
Fortran-like format specifications.
Note that the widths of the input and display formats do {\it not} have to
match.

You can edit this ASCII format file with any editor.
All you need to do is fill in the correct first and last column numbers of the
various fields, and delete references to fields you do not need.
You may also want to modify some formats; for example, you may have
Right Ascensions given only
to minutes and tenths instead of minutes and seconds.
The appropriate formats are explained in comment lines (beginning with
``!'') in the dummy format file.
After editing the format file, re-invoke {\tt MAKE/STARTABLE}.
This time, tell it you have a format file ready; it will use the formatting
information to convert your ASCII file to a MIDAS table file.

\subsubsection{Manual entry of data}
If you have no ASCII table ready, you can
enter stars manually from printed tables,
using the {\tt EDIT/TABLE} command to add them to the list.
This editor is briefly described under ``Interactive Editing of Tables'' in the
MIDAS Users Guide.
{\tt EDIT/TABLE} will automatically be invoked by {\tt MAKE/STARTABLE}
if you have no ASCII file to convert.
If you don't like the EDT-like editor provided by {\tt EDIT/TABLE},
create an ASCII file with whatever editor you like
and set up a *.fmt file to convert it, as described above.
Be {\it very} careful to double-check the results, as transposed digits and
other errors are almost certain to creep in when entering data by hand!

\subsubsection{Multiple star files}
Probably you will have several sources of program stars.
Each one can be turned into a separate MIDAS table file, as just described;
then, if you wish, the
tables can be combined into one big table, using the MERGE/TABLE command.
Note that the column
display formats will be taken from the {\it first} of the merged
tables; this may influence the order in which you do the merging.

If the formats are different, you may want to keep the program-star tables
separate.
In any case, you should keep standard, constant,
and variable stars in separate tables.
Then you just enter the names of the different tables
one by one, as you run the planning and reduction programs.

A more convenient way to handle multiple files is to make a catalog, and then
just supply the catalog name.
Only stars of a single type (standard, constant, or variable) should be kept in
the same catalog.
At present, only the reduction program can use star catalogs.

\subsubsection{Standard-star files}
Some standard-star files are supplied for the UBVRI and uvby-H-Beta systems.
These should be installed in the {\tt \$MIDASHOME/\$MIDVERS/calib/data/pepsys}
directory; see section~\ref{installtbl} if they are not.
The files are named {\tt UBVSTD.tbl}, {\tt UVBYST.tbl}, {\tt saaoUBVRI.tbl},
and {\tt saaouvbyHB.tbl}.
The latter two are the E- and F-region standards, established by
A.W.J.Cousins.

You can also make new tables of standard stars.
The easiest way to make standard-star tables is to begin as though they were
program-star tables, but to include columns of standard values in the ASCII
file, and field descriptions for them in the format file, before running
{\tt MAKE/STARTABLE}.
The standard column names are tabulated in the ``Standard values'' subsection
of ``Star tables'' in the Appendix.
Don't forget to add the SYSTEM descriptor when you are done!

In making new tables of standard stars, remember that they will be used
as extinction stars by the planning program.
Therefore, you should look for standards that pass within about 20 degrees
of your zenith; this picks a zone of declination centered on your latitude.
In right ascension, you want a fairly uniform distribution, so that a
standard will be available at large air mass frequently during the night.

The planning program will not select standards that are so near the Sun's
right ascension that they can only be observed at large air masses.
The instrumental magnitudes of stars that cannot be observed at
small air masses must be determined by transforming the standard values to the
instrumental system, and these transformed values have relatively large errors;
therefore, they contribute relatively little extinction information
(see pp.\ 162-163 of Young~\cite{Y74} for a discussion of this problem, and
Beckert and Newberry\cite{BN89} for examples of transformation errors).

If you must go far from the zenith to pick up standards of rare types,
put them in a special file, and treat them as program stars when running the
planning program.
Then observe them only near the meridian,
and don't try to use them as extinction stars.
They can still be treated as standards when you run the reduction program.


\subsection{Observatory table}
The standard ESO observatory table, called
{\tt esotel.tbl},
is stored in the MIDAS directory system at
{\tt \$MID\_PEPSYS/esotel.tbl} .
This directory is {\it not} distributed with the regular MIDAS updates,
but is part of the calib directory that can be reached by anonymous ftp
from mc3.hq.eso.org
(see the Dec.\ 1991 issue of the {\it Courier} for details).
Just cd to the midaspub directory, and see the README file there.
These table files must be installed on your local system before they can be
used; there are procedures for creating them in the calib directory.
These and other files in the calib directory subtree
can also be obtained on tape by special request.

This is a short table, so if you need to make up a new one, just copy the
ESO table and edit the copy with the {\tt EDIT/TABLE} command.
(If you do this, don't forget to change the descriptors as well as the table
entries!)
Or, using the dummy format file
{\tt \$MIDASHOME/\$MIDVERS/contrib/pepsys/lib/esotel.fmt}
to convert a *.dat file to MIDAS table form, you could create a whole new table.
The absolutely essential columns are the designation of the telescope focus,
accurate coordinates (including height above sea level), and the diameter of the
telescope.
Please consult the Appendix to clarify the meaning of the more obscure columns!
They are not yet required by PEPSYS, but may be in the future.

\subsection{Horizon tables}
Every telescope needs a horizon table, which tells where the telescope can look
in the sky, and also from which parts of the sky the Moon can shine on the
objective.
Without a horizon table, the programs will assume that the telescope has an
unobstructed view of the horizon --- an assumption that is probably
wrong.
The name of the table is keyed to the name of the telescope; see the Appendix.

If you need to make up a new table, the command {\tt MAKE/HORFORM} will provide you
with a form to fill in while gathering the data,
which can then be entered into the empty table provided.


\subsection{Instrument file}

Instruments are more complicated than star catalogs, but a MIDAS table file is
still useful for holding the description of an instrument.
The file is partitioned into two sub-tables.
The main sub-table contains information about the passbands used, and the
second sub-table describes the detector(s).
General information about the instrument is stored in descriptors.

The {\tt MAKE/PHOTOMETER} command will solicit the necessary information (if you
need to make up a new table), or display an existing table (so you can check its
contents).
Because the structure of instrument files is rather complicated, you should use
{\tt MAKE/PHOTOMETER} to build a new file when an instrument changes, rather than
trying to edit the file.
This will ensure that all the necessary information is included, and that the
file contains an internally consistent description of your instrument.

Basic data for an instrument file include the passbands used (and any
coding that is used to represent them in the data);
the detectors used, and their properties; and general information about
temperature control and optical cleanliness.
If your data include neutral-density filters or indicate which measuring
aperture was used, be sure to include the corresponding information in this
file.
See the ``File Formats'' Appendix for details.

One property that needs special attention in pulse-counting systems is the
{\it type} of dead time.
There are (in the textbook approximation) two kinds of counters; these are
known as paralyzable and non-paralyzable, and the corresponding dead
times are described as extending or non-extending.
Different analytical expressions relate observed and true counting rates for
the two types, so it is important to know which kind you have.
The matter is very clearly discussed by Evans~\cite{Evans}, which is a
standard reference on this subject.

If an instrument contains filters for more than one system, it may be useful to
maintain separate files for the different systems used (e.g., UBV and uvby).
This should certainly be done if different filter wheels are used for
different systems.


\subsection{General advice about table files}

The MIDAS table file system contains several commands for editing and modifying
table files.
Columns can be re-named, added, or deleted; various operations can be performed
on column contents, and the results stored in a new column; and so on.
You should read the on-line help for all the commands with the TABLE qualifier,
just to see what operations are supported.

One common problem arises with columns that contain character strings.
When you create a table from an ASCII file, the width of the table column is
(by default) the width of the field in the original *.dat file.
Later, you may want to add information to such a column, but find that the
existing table column is too narrow to hold the whole string.
For example, you might want to add aliases to the OBJECT names, or expand a
COMMENT column.

The column cannot be made wider by the {\tt EDIT/TABLE} command;
however, there is a way to widen it.
First, use the ``concat'' operation of {\tt COMPUTE/TABLE}
to add a string of blanks
to the existing column, and store the widened column under a temporary name.
Then delete the original column, and rename the new column with the old name.

You can avoid this problem by specifying a width wider than the original ASCII
field by using (say) C*32 instead of just C for the column in the format file,
when you create the table.
Specifying the width in the format file overrides the default width.
However, you may find that this wastes a lot of disk space if the column is
mostly blanks.


\section{Planning your observing run}
\subsection{Introduction}

Now suppose you have all the files set up to make MIDAS happy, and you are
ready to plan an observing run.
The proper distribution of extinction and standard stars was discussed by
Young~\cite{Y74}.
The program is based on this discussion, but includes other considerations as
well (such as the need to measure time-dependent extinction, and to measure
instrumental zero-point drifts).
It should help you get the right amount of good calibration data
with a minimum of observing effort.

The command {\tt MAKE/PLAN} will ask questions about what you are trying
to do, and produces a schedule of observations that will meet your goals, if
possible.
If you are trying to do something impossible, it will tell you so, and offer
suggestions.
In general, if you follow the suggestions of the planning program, you will get
the data you need --- assuming that the weather cooperates!

The planning program needs to know which photometric system you want to work
in, and will try to supply appropriate standard stars.
You can provide other table files of standard stars, if you wish; but be sure
they really are on the same system as the built-in standards!


\subsection{Preparing to use the planning program}
The planning program will ask you for a lot of information, so be
ready with the answers to its questions.
Obviously, it needs to know the location and size of the telescope you will
use, and the stars you want to observe;
so you must tell it the names of the files in which this
information is stored.
It also needs to know the accuracy you want to achieve, when you want to
observe, and so on.
Some questions ask you for choices between plausible alternatives, such as the
kind of time you want to use (UT, local zone, or sidereal).

You may need to look up or find out some of the information ahead of time.
If you are not using a common photometric system, you will need to supply the
central wavelengths and widths of the passbands you are using.
Even if you are using a common system, information on your actual
instrumental passbands will be helpful, if it is available.
If you get halfway through and discover you lack some vital piece of
information, just enter Q to quit.
Get what you need, and begin again.

Star coordinates are needed to a minute of arc or better for accurate data
reduction, so you may as well supply accurate coordinates for planning, too.
A few files of standard stars are already available.

If you are doing pulse-counting photometry, the uncertainty in the
pulse-overlap correction sets a limit to the brightest stars that can be used
as standards.
Therefore, you will need to have both the effective dead-time of your system,
and a realistic estimate of its uncertainty.

\subsection{Using the planning program}
The planning program asks a series of questions; you provide the answers.
Many are simple yes/no choices that can be answered as simply as
{\tt y} or {\tt n}.
In general, when a choice among several alternatives is required, you can
truncate your answer (usually, to just the initial letter)
so long as it uniquely specifies the choice.
The less typing you do, the less chance you have to make a mistake.

In many cases, the program will offer some guidance, or suggest reasonable
default values.  You will not go far wrong if you adopt its suggestions.
However, experienced observers may have good reasons to ``bend the rules''
a little at times.
The program will allow this, though it may complain if you are doing something
it thinks is unwise.

In general, input that is comprehensible to an astronomer is also understood by
the program.
For example, the first and last dates of your observing run
can be entered in any reasonable format, so long as the month is specified
either in full or by the usual 3-letter abbreviations.
Because European and North American conventions for dates differ in the
order of month, day, and year, the program will not accept ambiguous
forms like 3/8/95.
Because of the possibility of misinterpreting dates, the program will
tell you how it has interpreted what you typed in, and give you a chance to
correct any misunderstanding.

If you make a mistake when entering the name of a file, the program will give
you a chance to recover.
The program checks to make sure files exist and appear to be of the proper
type.
The suffix ``.tbl'' will be supplied automatically if you do not type it in.
If you mistakenly say there is another star file, and the program finds your
last filename does not exist, you can enter ``{\tt NONE}''
when it asks for another name, and it will go on to the next question.

\subsection{Selection criteria}
An overall criterion for standard-star selection is the accuracy needed in the
final results:
if you need more accuracy, you will need more standard stars.
The program asks what accuracy you are trying to reach, and tries to select
enough stars to meet your request without being excessive.

To determine the extinction accurately, you must observe some ``extinction
stars'' at both high and low altitudes.
While it is possible to get a rough estimate of the extinction by observing
{\it different}
standard stars at high and low airmasses, the relatively large conformity
errors in most photometric systems, together with the low accuracy of some
standards, make this method very inefficient (see
Young~\cite{Y74}, p.178).
The practical problem is that the extinction coefficient can change
considerably in the few hours needed for an extinction star to move from low to
high airmass (or vice versa).
To minimize this problem, the program selects stars that traverse a large range
in airmass in the shortest possible time; these are stars that pass near your
zenith.
Furthermore, it asks you
to observe extinction stars that are both rising and setting, to
avoid a correlation of airmass with time.

However, the airmass changes slowly when stars are near the zenith.
But to separate extinction drift from instrumental drift, you must
observe a wide range in air masses in a short period of time.
Therefore, the program selects times when the extinction stars cross an
almucantar a little removed from the zenith, rather than when they are on the
meridian.
This almucantar typically corresponds to about 1.1 airmasses.
The times of these crossings are denoted by {\tt 'EAST'} and {\tt 'WEST'}
in the output of the planner.
To optimize the precision of the extinction determination, the low-altitude
observations are placed at about $25^\circ$ altitude, near 2.36 airmasses.
Those scheduled observations are denoted as {\tt 'RISING'} and {\tt 'SETTING'}.

Furthermore, to track changes in the extinction accurately, you need an
extinction measurement (i.e., an observation of an extinction star at large
airmass) two or three times per hour.
This means that the stars used must be in the right places in the sky to
be at large airmasses when you need them.
In particular, although the Cousins E-region standards are excellent
secondary standards for {\it transformation} purposes, Southern-Hemisphere
observers should augment them with {\it extinction} stars more evenly
distributed on the sky.

Obviously, standard stars used for the transformation from instrumental to
standard system can also be used for the transformation from inside to outside
the atmosphere (traditionally called ``extinction correction'').
To minimize the number of calibration observations, the planning program makes
standard stars do double duty as extinction stars.
These stars should be bright enough that their photon noise is negligible;
the proper magnitude range depends on telescope size and the bandwidth of
the filter system used.
On large telescopes, bright stars are {\it too} bright, especially if you are
doing pulse counting.

Finally, the standard stars must have a good distribution in each of the
color indices of the system you are using.
Because transformations are generally non-linear (\cite{BN89}, \cite{Y92}),
a wide range of each color
should be covered rather uniformly; it is {\it not}
enough to observe a few very red
and a few very blue stars.
All these requirements impose constraints on the selection of standard
stars.

Notice that, in using standard stars to measure extinction, we need
not use the standard values transformed to the instrumental system
(though this is possible).
Instead, we use the actual observed instrumental values for these stars, which
are considerably more accurate than standard values transformed to the
instrumental system (cf.\ p.184 of \cite{Y74}),
because of conformity errors~\cite{MS92}.
We observe standards at both large and small airmasses, and determine the
extinction directly from these observations.
This matter is discussed more fully in connection with the reduction
program (see below).

\subsection{Influencing the plan}
The program will propose reasonable default choices, and will provide some
approximate estimates of the magnitudes where photon noise becomes excessive,
and where photon and scintillation noise are comparable.
Because the crossover between photon and scintillation noise depends on both
zenith distance and azimuth, a set of values will be presented for your
inspection.
Table~\ref{tbl:crossover} shows what a typical crossover table looks like.


\begin{table}[ht]
\begin{verbatim}
            SCINTILLATION = PHOTON NOISE   Photon Noise   Present
          at  secZ = 2.36   secZ = 1.10   of  5.sec.int.   FAINT
                between       between     is0.005 mag.at   limit
           U   9.9 &  8.9    7.1 &  7.0        10.5         5.5
           B  10.6 &  9.7    7.6 &  7.5        11.2         6.0
           V  10.7 &  9.8    7.5 &  7.4        11.3         5.9
\end{verbatim}
\caption[]{A sample crossover table}
\label{tbl:crossover}
\end{table}


There's a lot of useful information in this table, so let's go over it
carefully.
First, notice the similar pairs of columns at the left, under the title
{\tt SCINTILLATION = PHOTON NOISE}.
The left-hand pair of columns gives crossover magnitudes for an airmass near
2.36; the right-hand pair, for an airmass of 1.10.
These are the maximum and minimum airmasses at which the planning program
expects to schedule observations of extinction stars.
The smaller airmass will be adjusted by the program to provide more or
fewer extinction stars, as needed; you can push it back and forth a little
if you want.

For each airmass, there is one pair of columns.
These contain the magnitudes at which scintillation noise and photon noise are
expected to be equal, for two different lines of sight: one looking along the
projected wind vector in the upper atmosphere, and the other looking at right
angles to the wind.
When we look at right angles to the wind, the scintillation shadow pattern
moves with the full wind speed across the telescope aperture, and we have the
maximum possible averaging of scintillation noise in a given integration time,
the least scintillation noise, and thus the brightest possible crossover
magnitude.
But when we look directly along the wind azimuth, the motion of the shadow
pattern is foreshortened by a factor of sec $z$; then the scintillation noise
is maximized (for a given zenith distance), and the crossover does not occur
until some fainter magnitude, where the photon noise is big enough to match the
increased scintillation.

As we do not know the wind azimuth in advance, we can only say that the
scintillation and photon noises will be equal somewhere in the interval between
these two extremes.
We would generally prefer to have the extinction measurements limited only by
scintillation noise, so the initial faint limit for extinction stars is set
1.5 magnitudes brighter than the brightest of the crossover values.
This makes the photon noise half as big as scintillation, near the zenith.
(Our sample table shows these initial values.)
These are conservative values.

However, we may need to set the actual faint limit used in selecting standard
and extinction stars (shown in the rightmost column of
the table) somewhat fainter.
For example, we may be using such a large telescope that we cannot observe such
bright stars.
(In particular, if you are doing pulse counting, the program will impose a
bright limit as well as a faint one.)
In this case, both photon and scintillation noise will be quite small, and we
can safely use considerably fainter stars without compromising our requested
precision.

In the example above, the user has requested an accuracy of 0.01 magnitude.
The planning program divides this error budget into four parts:
scintillation, photon noise, transformation errors, and instrumental
instabilities.
If these are uncorrelated, each can have half the size of the allowed total; in
our example, that's 0.005 mag.
So the table gives the magnitude at which the photon noise reaches its allowed
limit (in the next-to-last column), for the adopted integration time (5
seconds, in our example).
This magnitude
should be regarded as an absolute limit for extinction and standard stars.

Obviously, we could actually push the extinction stars close to this
photon-noise limit, without exceeding the requested error tolerance.
However, between the photon-noise limit in the right half of the table and the
crossover values to the left, there is a substantial contribution of
photon noise to the total, and hence a substantial advantage to using brighter
stars.
If we use this advantage, we provide some ``insurance'' --- a little slack in
the error budget.

Whenever the crossover table appears, you will be given an opportunity to
change the actual planning limits, whose current values are given in the last
column.
The columns to its left provide you with the information you need to make a
good choice: the crossover magnitudes, and the pure photon limit, for each
band.

Although these values
are given to 0.1 mag precision, you should be aware that the
scintillation noise can fluctuate by a factor of 2 or more
within a few minutes, so that only a rough estimate is really possible.
Furthermore, the photon-noise estimates are only as good as the estimates
available to the program for the transmission of the instrument and the
detective quantum efficiency of your detector.
So all these numbers are a little ``soft''; you should not take that last
digit literally.
Just bear in mind that the photon noise varies with magnitude, and that the
scintillation varies with airmass and azimuth, by the amounts
shown in the table.

Now, let's consider adjusting the circumzenithal airmass.
If the high-altitude, low-airmass almucantar is too close to the zenith,
only a few standard stars will be available in the small zone it intercepts
as the diurnal motion carries the stars past it.
Then it may be necessary to choose a larger minimum airmass to expand the
region of sky available for choosing extinction stars.
Conversely, if too many extinction stars are selected, it makes sense to
reduce the zone width a little, thereby getting not only a more reasonable
number of stars, but also a little bigger range in airmass for each star.

The planning program will make coarse adjustments in the minimum airmass (i.e.,
in the width of this zone)
to get about the right number of extinction and standard stars, but you can
also make fine adjustments yourself.
The program gives you this option after making a preliminary selection of
standards.
When it asks whether you want to adjust the sky area or magnitude limits, reply
``sky,'' and it will make a small change and try again.
You may need to make several adjustments to get what you want.

If you expect from past experience that you will have excellent photometric
conditions, you may be able to reduce the number of extinction stars a little.
However, this is risky: if the weather turns against you, you may need more
stars than you expected!
Conversely, if you know the observing run is at a place or time of year that
usually has mediocre conditions, you will surely want to play it safe and add
some extra extinction stars.

You can also alter the number of candidate stars by adjusting the
magnitude limits.
By adjusting magnitude limits separately in different bands,
you can manipulate the range of colors available.
For example, because signals tend to be low in ultraviolet passbands, the
photon noise is high there, so
it often happens that the default faint-magnitude limits do not
allow enough very red stars.
In this case, making (say) the faint U limit fainter and the V limit brighter
biases the selection toward redder standards.

Thus, you can manipulate the region of sky and the magnitude limits to select
a reasonable number of standards with a good range of colors.
You can make
such adjustments iteratively until you are satisfied with the selection.
At each stage, the program will show you the distribution of the selected stars
in a two-color diagram, as well as on the sky, and ask if its selection is
satisfactory.
If it is not, you reply ``no'' to the question ``Are these stars OK?'', and
then have another chance to manipulate the zone width and magnitude range.
When you finally reply ``yes'', the program will make up an observing schedule
for each night of the observing run, using the set of stars you approved.

Keep in mind the possibility that a star previously certified as a photometric
standard can still turn out to be variable!
A surprising number of bright eclipsing binaries continue to be discovered
every year.
Furthermore, stars of extremely early or late spectral type, and stars of high
luminosity, tend to be intrinsically variable by small amounts.
Therefore, it is important to use a few more standards than would otherwise be
absolutely necessary, just in case of trouble.
A little redundancy is cheap insurance against potential problems.
We also need some redundancy to find data that should be down-weighted
(see section~\ref{reductions} to see why.)

\section{Getting the observations}

When you go to the telescope, follow the plan as closely as you can.
If the telescope and instrumentation are unfamiliar to you, you may have
trouble keeping to the schedule for the first night or two, so
try to work with a previous observer on the equipment for a night or
two of training before your observing run.
If you fall behind, try to keep as many low-altitude observations of extinction
stars as you can, and drop some circumzenithal ones if necessary.

Don't change your mind in the middle of the night!
If you think you will need to adopt a different plan in the middle of the run,
make sure you generated that alternate plan ahead of time.
The time to be creative is during the planning stage, not while observing.

Remember Steinheil's principle: only {\it similar} data can be compared.
That means: DON'T switch to a different filter set or detector in the middle of
a run.
DON'T use different focal-plane apertures for different stars or different
nights; they change the instrument's spectral response as well as its
zero-point.
DON'T alter temperature, high-voltage, or discriminator settings.

There are very great advantages to getting homogeneous data, because then the
determination of the instrumental parameters can be spread over several nights;
for example, see
\cite{YI67},
\cite{Y74},
\cite{P82},
\cite{MH84}, and
\cite{ASC91}.
Likewise, the extinction can be determined much more accurately from
observations spread over full nights than from partial nights.
If you are only interested in a small region of sky that is not up all night,
either use the rest of the night to get extinction and standard-star data, or
try to combine your run with another program that uses the same equipment;
then all the data from both programs can be reduced together.

If there were instrumental parameters you weren't sure of when you made up the
plan, make sure you find out their values during the observing run.
The printed schedule will provide blank spaces for you to write them in.
Many useful instrumental tests can also be run during the daytime, or during
cloudy nights.

Remember that the extinction changes more on long time scales than on short
ones; so arrange your observations to chop as much extinction
variation as possible out of the most important data.
If you need very good $(u-v)$ colors, observe those filters in direct
sequence.
If you need to determine an eclipsing binary's light curve in several filter
bands, you may do better chopping between program and comparison stars before
you change filters (see~\cite{Y91} for further discussion).

A very common problem in observing records is an incorrect clock setting.
Even a one-second error is worth correcting.
An hour off is a common and much more serious mistake.
Sometimes tired observers enter the wrong {\it kind} of time
(say, UT instead of LST) into the observing records.
It's a good idea to note both times in your observing log --- you {\it do}
keep a written log of your observations, don't you?
(It's easy to annotate the printed output of the planning program for this
purpose.)
Even if the
computer is supposed to keep track of these things, what happens if there is a
power failure, or the computer crashes in the middle of the night?

It doesn't hurt to take note of moonrise and moonset; they can be used to
double-check the recorded times for gross errors.
Be sure to write down anything that seems suspicious in the data, or possible
problems with the equipment.
Don't imagine you will remember all the details of what went on when you get
around to reducing the data, which might be weeks or months later.
For example, which night was it the mirror got snowed on?
You can expect a big zero-point shift (or worse!) after that.

If temperature and humidity readings are not recorded automatically, be sure
you write them down {\it at least} once an hour.
(Space for this information is provided on the planning schedules.)
The temperature that is important is the temperature of the photometer
itself, not the observing room!
But sometimes you only have a thermometer available on the wall; then record
that, as any information is better than none at all.
Be sure to record the {\it time} when each reading was taken.
And don't forget to check the dark level now and then; it is a useful
check on the health of the instrument.

It is very easy to forget that ``nuisance'' parameters, like extinction and sky
brightness, have effects that propagate directly into errors of the
measurements you care about.
But this means you must be just as careful in measuring sky as in measuring
stars, and as careful with extinction and standard stars as with program stars.
Accurate and careful centering of stars is very important, particularly with
photometers using opaque (inclined) photocathodes.
Machine-like regularity is the ideal, attainable with some (but not all)
automatic centering systems.

Likewise, one must always measure the sky the {\it same distance} from each
star, regardless of brightness.
Otherwise, you include a different fraction of the star's light in the sky
measurement.
A common mistake is to offset the telescope until the glow from the star is no
longer visible in the measuring aperture; this guarantees that bright stars are
measured with a different system than faint ones.
A fixed angular offset --- with due regard to the telescope's diffraction spikes
-- should always be used.
If possible, always measure the sky in the same {\it place} for each star,
to make sure you include the same unseen background stars in every
measurement.

Finally, because sky brightness fluctuates on short time scales, it's best
to measure sky for each star.
You can use much shorter integrations for sky than for bright stars; but don't
let more than a few minutes go by without re-measuring sky.
During lunar twilight, or when the Moon is near the horizon, when the sky is
changing rapidly, measure sky both before and after each star.

And, if you are new to photometry, you won't believe how sensitive the results
are to the slightest {\it hint} of a cloud anywhere in the sky.
Even a tiny cloud anywhere in the sky means you probably have ``incipient
clouds'' (variable patches of haze) all over the sky, because of the layered
structure of the atmosphere.
Even Steinheil, a century and a half ago, found
that the slightest trace of cirrus made photometric
measurements impossible --- that's how gross the effects are.

It's a good idea to look at the sky around sunset.
If there are ``layers'' visible in the sky near the western horizon,
or ``notches'' in the sides of the setting sun, you are in for trouble,
and will need to put extra effort into extinction determinations.
During the night, any indications of clouds should be noted in your log.
Often, cirrus is not visible until the Moon rises --- or until the data are
reduced!


\section{Reducing the observations}

\subsection{Preliminaries}

Now you have your data, and are ready to reduce them.
Table~\ref{tbl:checklist} gives a list of the files you will need.


\begin{table}[ht]
\begin{center}
\begin{tabular}{|l|c|c|}           \hline
 File type & Contents				& Made by ... \\ \hline
Observatory & Telescope information		& editing esotel.tbl  \\
Built-in Standard-star& Std.-star positions \& values & copying {\tt calib} files \\
Additional Standard-star& Std.-star positions \& values & edit files, {\tt MAKE/STARTABLE} \\
Program-star & Program-star positions		& {\tt MAKE/STARTABLE}  \\
Instrument & Instrument information		& {\tt MAKE/PHOTOMETER}  \\
Data		& Stellar measurements		& {\tt CONVERT/PHOT}  \\
	\hline
\end{tabular}
\end{center}

\caption[]{Checklist of necessary table files}
\label{tbl:checklist}
\end{table}


If you need to make supplemental standard-star files, refer to the
``Standard-star files'' subsection at the end of section~\ref{startables}.
If your data files use different names for standard stars than are used in
some existing standard-star file, you can just supply a ``program-star'' file
that contains your names for these stars; use {\tt MAKE/STARFILE} to make this
table.
The reduction program will try to make the cross-identifications,
using positional coincidences as well as names
to decide which entries to match up.

If you were able to find additional information about the instrument while you
were observing, be sure to add any missing information to the instrument file.
The reduction program needs to have an accurate model of the instrument, as
well as your procedure.
You can check an existing instrument file by running {\tt MAKE/PHOTOMETER}.

Next, you must convert your observational data to the format MIDAS can digest.
If other people have already used the same equipment you did, there are
probably programs available to do this conversion.
If not, you will have to devise the conversion yourself.

In any case, you should look over the data files, line by line, to make sure
there are no garbled records.
Often, either equipment malfunctions, operator errors, or
problems on reading and writing tapes convert some part of the data
into nonsense.
You cannot reasonably expect any file-conversion program to cope with nonsense,
so edit it out first.

The sooner mistakes are removed from data, the less trouble they cause.
Each data-logging system is likely to make its own peculiar kinds of
mistakes; and each observer tends to make some particular kinds of errors at
the telescope.
Therefore, it is {\bf strongly recommended} that users have their own
preprocessing software, to detect mistakes and inconsistencies in raw data
files, {\it before} running {\tt CONVERT/PHOT}.

A common problem is an incorrect time or date setting; be particularly careful
to correct any errors that were discovered belatedly while observing.
Other common problems include misidentified stars; incorrect designations of
star and sky measurements; missing or incorrect data in star files; and
non-standard data formats.

In fixing errors in the data, {\it always} save a copy of the original in case
you make a mistake while editing.

\subsection{Format conversion}

If you have a program to convert data to MIDAS input format, fine.
Use it and go on to the next section.

If you do not,
first see if there is {\it any} existing program that reads your data
format.
(Very likely there is, unless the instrument is new.)
All you need to do is cannibalize the part of the existing program that reads
the data, and add a short piece of code to
re-format the data first as a flat ASCII file, one measurement per line; and
then convert this to *.tbl format.
Sound familiar?
It's just like the process of creating a program-star table.
Once again, a dummy *.fmt file is provided to help you.

That's usually the easiest way to go.
You may decide to write a short MIDAS script to apply your format-conversion
program and the edited *.fmt file to each night's data.
That can be put into a procedure, given a name, and invoked as needed.
This has already been done for data in the formats produced by the photometers
at the 1-meter ESO and the Danish
telescopes on La Silla; the command CONVERT/PHOT will convert
their data files to the proper MIDAS tables.

If you are {\it very} lucky, your data may already be in the form of a flat
ASCII file of consistent format; then all you need to do is edit the *.fmt file
and apply it to generate the *.tbl files MIDAS needs.
More likely, you have most of your data in this form, but interspersed with
comments and/or records of other types.
In this case, it may be possible to write a simple pre-processor (or use
some UNIX utility, like {\tt awk}) to strip out the data as a fixed-format
ASCII file.

If there are just a few distinct record formats in your data files, it may be
most convenient to strip out each type separately; convert each one first to a
homogeneous ASCII file, and then to a MIDAS *.tbl file; and then use
MERGE/TABLE to combine the separate tables into the final file for input to the
reduction program.
If each line of each file is correctly time-tagged, a simple sort on the
MJD\_OBS column will then put everything back together in the right order.

In any case, you should automate this procedure as much as possible.
That means writing a short MIDAS or shell script.
Whatever you do, {\it don't} resort to copying things out by hand!
At worst, you might need to insert a little timing information by hand to label
comments in the original data.

You should be careful to make sure enough information is in the final
*.tbl file to tell the reduction program everything it needs to know.
For example, if you measure dark current, and have more than one detector
channel, the table must show which detector's dark current is being measured.
That is usually indicated by the filter position; so you will probably need to
have a non-blank filter position for ``dark'' measurements.
Note, however, that you do {\it not} need to distinguish among standard,
extinction, program, and variable stars in data files, as this information is
either
supplied when you read in the star files, or determined dynamically during the
reductions.

The minimum of essential columns in a data file are
object identification, band identification, signal, time, and integration time.
The section ``Required observational data'' in the Appendix describes these
columns.
For standard band names in commonly-used systems, see the Appendix
(particularly Table 2 in the section ``Standard values'' under ``Star
tables'', and the description of the {\tt BAND} column in the ``Passbands''
subsection under Section 5).
Usually, you will also need the {\tt STARSKY} column, to distinguish between
star and sky measurements.
But many other data can be useful: temperatures, relative humidity, various
instrumental settings, and error estimates.
The section ``Additional information'' in the Appendix describes these in
detail.

Note that temperature and humidity data can be interpolated, if they are not
routinely recorded with each observation.
If your instrument table shows that neutral-density filters may be used, you
need to include this information as part of the band designation in the data
files.
Conversely, if your data indicate that ND filters have been used, their
presence and nominal attenuation factors should be included in the instrument
table.

It would be nice to have a universal data-conversion utility.
Unfortunately, the amazing variety of formats produced by dozens of home-grown
data-logging systems seems to prevent it.
You know your data better than anyone else; so you can put them in standard form
better than anyone else.

\subsection{Reductions --- at last!}
\label{reductions}

Now you really are ready to run the reduction program.
You have all the necessary information in the right file formats.
Usually, it's convenient to keep each night's observations in a separate file.

The command {\tt REDUCE/PHOT} will start up the reduction program.
It starts out very much like the planning program, as it also needs the
telescope, standard-star, and instrument information.
But instead of asking about output formats and the date of the observing run,
it requests the names of the data files.
Normally, each data file is one night; you can reduce up to 30 nights
together.

To save typing, it is convenient to make a catalog that refers to all the data
files; use the MIDAS command {\tt CREATE/TCAT} to do this.
For example, if your data files have names like night1.tbl, night2.tbl, and so
on, the MIDAS command line

\vspace{1.0mm}

{\tt CREATE/TCAT \ data \ night*.tbl}

\vspace{1.0mm}

\noindent
will make a catalog named data.cat that refers to the whole group.

As with star files, if you say by mistake that there is another data file,
you can enter ``{\tt NONE}''
when it asks for another file name, and it will end the list and go on.

If the data are pulse counts, they are corrected for the nominal dead time
(using a formula appropriate to the type of counter used) as they are read in.
If they are raw magnitudes, they are converted to intensities at this stage, so
the reduction can proceed in the same way for all types of data.
At this point, observations are in intensity units
on some arbitrary instrumental scale.
The intensities are then corrected for the nominal values of any neutral
attenuators that may have been used.

Once the data are read, the reduction proceeds in two main stages: first,
filling in missing meteorological data, subtracting dark readings, and
subtracting sky; and
second, fitting the remaining stellar data to appropriate models.

If they are present,
the temperatures and relative humidities for each night are displayed
graphically.
After looking at each graph, you can choose one of three treatments: polygon
interpolation (i.e., just a straight line segment between adjacent data);
linear smoothing (a straight line fitted to the whole set); or adopting a
constant value for the whole night.
Usually, polygon interpolation is adequate.
However, it is sensitive to aberrant or bad points.
If you believe there are outliers in the data, and the rest are reasonably
linear, use the linear fit, which resists the effects of bad points.
If you think a bad point can be fixed, or should be removed before proceeding,
just enter ``q'' to quit, and deal with the bad datum.

\subsubsection{Robust fits and bad points}

Already, we must decide how to deal with discordant data.
As this is a subject of considerable importance, a short digression is in order.

Poincar\'e attributed to Lippmann the remark that ``Everyone believes in the
normal law of errors: the mathematicians, because they think it is an
experimental fact; and the experimenters, because they suppose
it is a theorem of mathematics.''
However, it is neither an experimental fact nor a theorem of mathematics.

Experimentally, numerous investigations have shown that real errors are rarely
if ever normally distributed.
Nearly always, large errors are much more frequent than would be expected for a
normal distribution (see \cite{MT77}, pp.~10 -- 12,
and \cite{H86}, pp.~20 -- 31).
Menzies and Laing \cite{ML85} show clear examples in photometric data.

Mathematically, the reason for this behavior
is well understood: although the Central Limit Theorem promises
a Gaussian distribution
{\it in the limit} as the number of comparable error sources
approaches infinity, the actual approach to this limit is agonizingly slow ---
especially in the tails, where a small number of large individual contributors
dominate.
In fact, if there are $n$ independent and identically distributed contributors,
the rate of convergence is no faster than $n^{-1/2}$ \cite{H80}.
If we wanted to be sure that our distribution was Gaussian to an accuracy of
1\%, we would need some $10^4$ elemental contributions --- clearly, an
unrealistic requirement.
In practice, a few large error sources dominate the sum.

Furthermore, the proportionality constant
in the convergence formula changes rapidly with distance from
the center of the distribution, so that
convergence is very slow in the tails.
This guarantees that the tails of real error distributions are always far from
Gaussian.

In the last 30 years, the implications of these deviations from ``normality''
for practical data analysis have become widely appreciated by
statisticians.
Traditionally, the excess of large errors was handled by applying the method of
least squares, after rejecting some subset of the data that appeared
suspiciously discordant.
There are several problems with this approach.

First, the decision whether to keep or reject a datum has an arbitrary
character.
A great deal of experience is needed to obtain reliable results.
But manual rejection may be impractical for large data sets; and
automated rejection rules are known to have inferior performance.
Second, rejection criteria based on some fixed number of standard deviations
result in no rejections at all when the number of degrees of freedom is small,
because a single aberrant point greatly inflates the estimated standard
deviation (\cite{H86}, pp.~64 -- 69).
The common ``3-$\sigma$'' rejection rule rejects nothing in samples smaller
than 11, no matter how large the biggest residual is; the inflation of the
estimated standard deviation by just one wild point outruns the largest
residual in smaller data sets.
There is no hope of rejecting a bad point this way in samples of 10 or smaller;
but one rarely measures the same star 10 times.
For the more typical sample sizes of 3 and 4, the largest possible residuals are
only 1.15 and 1.5 times the estimated standard deviation.
Third, including or rejecting a single point typically introduces
discontinuous changes in the estimated parameters that are comparable to their
estimated errors, so that the estimated values undergo relatively large jumps
in response to small changes in the data.
We would have more trust in estimators that are continuous functions of the
data.

Finally, the nature of most observed error distributions is not that
data are clearly either ``good'' or ``bad'', but that the few obviously wrong
points are accompanied by a much larger number of ``marginal'' cases.
Thus the problem of rejection is usually not clear-cut, and the data analyst is
left with doubts, no matter where the rejection threshold is set.
The reason for this situation is also well understood: most data are affected
by error sources that vary, so that the ``marginal'' cases represent data
gathered when the dominant error source was larger than average.
Such observations are not ``wrong'',
though they clearly deserve smaller weights than
those with smaller residuals.

In particular, we know that {\it photometric} data are afflicted with variable
errors.
For example, scintillation noise can vary by a factor of 2 on time scales of a
few minutes; and by an additional factor of $sec \, Z$ at a given air mass,
depending on whether one observes along or at right angles to the
upper-atmospheric wind vector.
Menzies and Laing~\cite{ML85} discuss other possible sources of error.
Therefore, we know we must deal with an error distribution that is
longer-tailed than a Gaussian.
Furthermore, both scintillation and photon noise are decidedly asymmetrical.
As these are the main sources of random error in photometric observations, we
can be sure that we {\it never} deal with truly Gaussian errors in photometry.

Unfortunately, the method of least squares, which is optimal for the Gaussian
distribution, loses a great deal of its statistical efficiency for even
slightly non-Gaussian errors.
(Statistical efficiency simply refers to the number of observations you need to
get a desired level of reliability.  If one estimator is twice as efficient as
another, it will give you the same information with half as many
observations.)
The classical example is Tukey's contaminated distribution.
Suppose all but some fraction $\epsilon$ of the data are drawn from a normal
distribution, and the remainder are drawn from another Gaussian that is three
times as broad.
Tukey~\cite{T60}
asked for the level of contamination $\epsilon$ that would make the
mean of the absolute values of the residuals (the so-called {\it average
deviation}, or A.D.) a more efficient estimator of the population width
than the standard deviation, which is the least-squares estimator of width.

Although the mean absolute deviation has only 88\% of the
efficiency of the standard deviation for a pure Gaussian, Tukey found that less
than 0.2\% contamination was enough to make the A.D. more efficient.
The reason is simply that least squares weights large errors according to the
squares of their magnitudes, which gives them an unreasonably large influence
on the results.

Similar, though less spectacular, results exist for position estimators.
For example,
about 10\% contamination is enough to make the median as efficient as
the mean (the least-squares estimator); while several ``robust''
estimators are some 40\% more efficient than the mean at this level of
contamination.
Real data seem to be somewhat longer tailed than this, so the mean (i.e., least
squares) is typically even worse than this simple example suggests.

Because convergence of the central limit theorem is much faster near the center
of the error distribution than in the tails, we can expect real error
distributions to be nearly Gaussian in the middle, and this is in fact
observed to be true.
A practical approach to data analysis is then to treat the bulk of the data
(in the middle of the distribution) as in least squares; but to reduce the
contribution of the points with large residuals, which would be rare in
a genuinely Gaussian distribution, in a smooth and continuous
fashion.

There is now a large literature on ``robust'' estimation --- that is, on methods
that are less critically dependent on detailed assumptions about the actual
error distribution than is least squares.
They can be regarded as re-weighted least squares, in which the weights of
data with moderate to large residuals are decreased smoothly to zero.
There are many ways to do this; all produce rather similar results.
The really ``wild points'' are completely rejected; the marginal cases
are allowed to participate in the solution, but with reduced weight.
The result is only a few per cent less efficient than least squares for
exactly Gaussian errors, and much better than least squares
--- typically, by a factor of the order of two ---for realistic error
distributions.
These methods
are also typically 10\% or so more efficient than results obtained by
experienced data analysts using careful rejection methods (\cite{H86}, pp.~67
-- 69).

The particular method used here for reducing photometric data is known to the
statisticians as ``Tukey's biweight''; it is easy to calculate, and produces
results of uniformly high efficiency for a range of realistic distributions.
To prevent iteration problems, it is always started with values obtained from
even more robust (but less efficient) estimators, such as the median and its
offspring, Tukey's robust line~\cite{HMT}.
The usual method starts with a very robust but inefficient estimator
such as the median or Tukey's robust line; switches to Huber's M-estimator for
initial refinement until scale is well established; and then iterates to the
final values using the biweight.
If you are unaware of the need to precede the biweight with an
estimator having an non-redescending influence function, don't worry.
This is known to be a numerically stable procedure.

As robust methods depend on ``majority logic'' to decide which data to
down-weight, they obviously require a certain amount of redundancy.
One cannot find even a single bad point unless there are at least three to
choose from (corresponding to the old rule about never going to sea with two
chronometers).
Therefore, it is better to obtain a large number of short integrations
than a smaller number of longer ones, provided that
the repetitions are separated in time enough to be independent.
The planning program will help you get the necessary data.

In summary, photometric data are known to have decidedly non-Gaussian error
distributions; so we use methods designed to be nearly optimal for
these distributions, rather than the less reliable method of least squares.
These methods are closely related to least squares, but are much less
sensitive to the bigger-than-Gaussian tails of real error distributions.
From the point of view of the average user, the methods employed here are
simply a more effective refinement of the old method of rejecting outliers.

The advantage of using these well-established, modern
methods is a gain in efficiency of some tens of
per cent --- exactly equivalent to increasing the amount of observing time by
such an amount.
It's about like getting an extra night per week of observing time.
This advantage is well worth having.

\subsubsection{Subtraction of dark and sky measurements}

After interpolating any meteorological values that exist,
you may have to subtract dark current and/or sky values.
If there are no data for dark current, the reduction program will complain but
continue.
Sky values may have already been subtracted (e.g., in CCD measurements or
other data marked as RAWMAG instead of SIGNAL.)

DARK CURRENT:
Dark current is usually a strong function of detector temperature.
If you have regulated the detector temperature, then only a weak time
dependence might be expected --- perhaps only a small difference from one night
to the next, or a weak linear drift during each night.
You will have the choice of the model to be used in interpolating dark values.

If the detector temperature is measured, you should look for a temperature
dependence that is the same for all nights.
The program will show you all the dark data, with a separate symbol for each
night, as a function of temperature.
If all nights seem to show the same behavior, you can then fit a smoothing
function of temperature to the dark values.
You can choose a constant, a simple exponential (i.e., a straight line in the
plot of log$\,($DARK) vs. temperature), or the sum of two exponential terms.

Although in principle these ought to be of the form
$D = a \, exp \, (-b/T)$, the range of temperature available is usually insufficient
to distinguish between this and a simple $a \, exp \, ( c T )$ term.
Furthermore, though the temperatures {\it ought} to be absolute temperatures in
Kelvins, you may have only some arbitrary numbers available, which might even
assume negative values.
In this case, an attempt to fit the correct physical form would blow up, but
the simple exponential term might still give reasonable results.
So the simpler form is actually used.

If the plot of log$\,$(DARK) vs. temperature bends up at the ends, or at least at
the right end, you should be able to get a good two-term fit.
If it looks linear, you can just fit a single line.
You also have the option of adopting a single mean value.

On the other hand, if the data are not consistent from night to night, or
show a temperature dependence that is different from the expected form, or if
you have no temperature information at all, you may have to interpolate the
dark data simply as a function of time.
As with the weather data, you have a choice of polygon, linear, or constant
fits.
Remember that a polygon fit uses every datum, right or wrong, and so is not
robust.

After removing a temperature-dependent fit, you will see the remaining
residuals plotted as a function of time.
This provides a double check on the adequacy of dark subtraction.

SKY SUBTRACTION:
Sky data must be treated separately for each night and passband.
Here, your options are more numerous.
You can choose the usual linear or constant fits; but those are likely to be a
poor representation of sky brightness.

More conventional choices are to use either the preceding or following sky for
each star observation, or the ``nearest'' sky (in which both time and position
separations are used to decide what ``nearest'' means).
Linear interpolation between successive sky measurements (i.e., a polygon fit)
is also an option.
These choices, while conventional, are not robust.
They are sensitive to gross errors in sky data, such as
star observations that have been marked as sky by mistake.

One might argue that bad sky data will stand out in the plots discussed below,
and that careful users will remove them and re-reduce their data.
One might also argue that really bad sky values will cause the stellar data to
be discarded or down-weighted in the later reductions, so that a robust fit at
this stage is not absolutely necessary.
However, such arguments are not completely convincing.
Therefore a more elaborate sky subtraction option is available, which tries to
model the sky brightness, discriminating
against outlying points in a robust regression.

To help you choose the best method, the program displays three plots of sky
brightness against different independent variables: time, airmass, and lunar
elongation.
In the time plot, the times of moonrise and moonset are marked,
and twilight data are marked t;
Figure~\ref{skyplot} shows an example.
In the other two plots, points with the Moon above the horizon are marked with
the letter M, points with the Moon below the horizon are marked by a minus
sign, and twilight data are marked t.
In these and other plots, the characters {\tt \^\ } on the top line or {\tt v}
on the lower edge indicate points
outside the plotting area; and \$ indicates multiple overlapping points.
You can re-display the plots if you want to look at them again before deciding
which sky-subtraction method to use.

\begin{figure}[htp]
\begin{center}
\begin{verbatim}

u  SKY on  FEB 14 1993    R = moonrise    t = twilight sky

  X  I  ^                                                ^      ^^
  10.I
     I  t                      R
   3.+                         |                                t
     I                         |
     I
     I
   2.+                                                     *   tt
     I          *                                          $
     I           $*                       *      * **$**   *
     I  t         $*                     **$$**$*$* $$ * *     t
   1.+   tt       *      *       $ $$*$*$$$**$*$* **  * $      t
     I   t *             *   $$$* $*$$**$
     I   t  ****$*  *          **
     I       *$    * $$$$ $
   0.+               $
     I
     I
     I
  -1.+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+--+---
     0  2  4  6  8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 4E-02
    TIME (MJD = 49032 + decimal of day) -->
\end{verbatim}
\caption[]{Plot of sky brightness as a function of time}
\label{skyplot}
\end{center}
\end{figure}

Note that no one method is best for all circumstances.
While modelling the sky should work well under good conditions, there
are certainly cases in which it will fail.

For example, when using DC or charge-integration equipment,
an observer commonly
uses the same gain setting for both (star+sky) and sky alone.
This is perfectly appropriate, as it makes any electrical zero-point error
cancel out in taking the difference.
But often the limited precision available --- for example, a 3-digit digital
voltmeter --- means that the sky brightness is measured with a precision of
barely one significant figure when bright stars are observed.
If a bright star reading is 782 and the sky alone is 3, one does not have much
information to use in modelling the sky.

Another case where one does better to subtract individual sky readings is
observations made during auroral activity.
While one would prefer not to use such data, because of the rapidity of sky
variations, they must sometimes be used.
Here again, subtraction of the nearest sky reading is better than using a
model, because the rapid fluctuations are not modelled.
Likewise, when terrestrial light sources around the horizon make the sky
brightness change rapidly with azimuth
and/or time, no simple sky model would be adequate.

If it is necessary to make measurements of some objects through two or more
different focal-plane diaphragms, these measurements cannot be combined
directly.
Ordinarily, all observations to be reduced together should be measured through
the same aperture, because the instrumental system changes in an unpredictable
way with aperture size.
Even the sky measurements are not exactly proportional to the diaphragm area.
However, it may be possible to reduce program objects observed with a
non-standard aperture {\it as if} they were measured through the standard
one, and then apply a suitable transformation after the fact.
This means that a sufficient number of calibration measurements of stars having
a considerable range in color {\it must} be taken, using both aperture sizes,
to determine the transformation between the two instrumental systems.
In such cases, individual sky readings taken through the same apertures must be
used in the reductions.
The reduction program will complain if you try to intermix data taken through
different diaphragms, and data taken with a peculiar aperture will be rejected
if there are no corresponding sky measurements.

Finally, when very faint stars are observed (as in setting up secondary
standards for use with a CCD), so that the sky is a large fraction of the star
measurement, it may be necessary to subtract individual sky readings simply
because the model used is not sufficiently accurate.
The model is reasonably good, but is not good enough to produce
estimates free of systematic error.

In any case, the plots of sky vs. time, airmass, and lunar elongation should
prove useful in assessing the quality of the sky data, and in choosing the
best subtraction strategy.
Furthermore, the residuals from the sky model may be useful in identifying bad
sky measures that should be removed; so it is a good idea to run the sky
model, even if you decide not to subtract its calculated values from the star
data.

\subsubsection{Sky models}
For the user interested in the details of the methods used, here are the
actual algorithms used in selecting the ``nearest'' sky, and in the
sky-brightness model:

In the ``nearest'' method, a distance estimator is computed that includes
separations both in time and on the sky.
The estimator is
$$
S = 20 | t_1 - t_2 |  +
    | AM_1 - AM_2 | + | ( AM_1 + AM_2 ) (AZ_1 - AZ_2 ) |
$$
where the $t$'s are times in decimal days,
the $AM$'s are air masses, and the $AZ$'s are
azimuths (in radians) of the two observations being compared.
(The azimuth difference is always taken to be less than $\pi$ radians.)
$S$ crudely takes account of the greater variation in sky brightness with
position near the horizon.
At moderate air masses, it makes a separation of a minute in time about
equivalent to a degree on the sky.

This estimator is computed for the two sky samples closest in time to each star
observation, one before the star and the other after it.
Only observations made with the same filter (and, if diaphragm size is
available, with the same diaphragm) are used.
The sky observation that gives the smaller value of $S$ is the one used in the
``nearest'' sample method.
Obviously, if the star is the first or last of the night, there may be no sky
sample on one side; then the one on the other side (in time) is used.

The angular part of $S$ is necessary to prevent problems when groups of
variable and comparison stars are observed.
It can happen that the sky is observed only after two stars in a group have
been observed, if only a single sky position is used for the whole group.
If sky was measured after the last star before the group, that
previous sky may be closer
to the first star of the group, in time alone, than the appropriate sky within
the group.
Then a purely time-based criterion would assign the (distant) previous star's
sky to the first star of the group, instead of the correct (following) sky.
Similar effects can occur at the end of a group, of course.

While this problem can be prevented by careful observers, not all observers
are sufficiently careful to avoid it.
The crude separation used here is adequate to resolve the problem without
going into lengthy calculations.
Note that both the airmass dependence of sky brightness and the possible
presence of local sky-brightness sources around the horizon (e.g., nearby
cities) make the horizon system preferable to equatorial coordinates for this
purpose.

This brings us to the sky model.
The general approach is to represent the sky brightness as the sum of two
terms, a general one due to scattered starlight, zodiacal light, and airglow;
and an additional moonlight term that applies only
when the Moon is above the horizon.

The airglow and scattered starlight are proportional to the airmass, to a first
approximation.
Actually, extinction reduces the brightness of the sky light near the horizon.
However, the full stellar extinction appears only in the airglow and
zodiacal components, not in the scattered light.
All three components are of comparable magnitude in the visible part of the
spectrum.

While Garstang~\cite{G1,G2,G3} has made models for the night-sky brightness,
these require much detailed information that is not usually available to the
photometric observer.
Garstang's models also were intended to produce absolute sky brightnesses, while
data at this preliminary stage of reduction do not have absolute calibrations.
Finally, they do not include moonlight.
Therefore, a simpler, parametric model is used.

Some guidance regarding the scattered starlight
can be obtained from the multiple-scattering results given
by van de Hulst~\cite{vdH}.
In photometric conditions, the aerosol scattering is at most a few per cent,
and the total optical depth of the atmosphere is less than unity.
In the visible, the extinction is dominated by Rayleigh scattering, which is
not far from isotropic, and nearly conservative.
Therefore, we are interested in cases with moderate to small optical
depth, and conservative, nearly isotropic scattering.
Because the light sources (airglow, stars, and zodiacal light) are widely
distributed over the sky, we expect small variations with azimuth, and
can use the values in van de Hulst's Table 12 to see that
azumuthally-averaged scattering has the following properties:

1.
For air masses such that the total optical depth along the line of sight is
less than 1, the brightness is very nearly proportional to air mass, regardless
of the altitude of the illuminating source.

2.
For vertical optical depths $\tau$
less than about 1.5, the sky brightness reaches a
maximum at an air mass on the order of $1/\tau$, and then declines to a fixed
limiting value as $ M \rightarrow \infty$ (remember that for the plane-parallel
model, the airmass does go to infinity at the horizon).

The decrease in the scattered light at the horizon is also to be expected in
the direct airglow and zodiacal components attenuated by extinction; so the
same general behavior is expected for all components.
The simplest function that has these properties is
$$
B_1 = (aM + bM^2)/(1 + dM^2) ,
$$
where M is the airmass; the limiting brightness at the horizon is just $b/c$.
Actually, a substantially better fit can be obtained to the values in van de
Hulst's Table 12 by including a linear term in the denominator; so
the approximation
$$
B_1 = (aM + bM^2)/(1 + cM + dM^2) ,
$$
is used to represent the airglow and scattered
light.

There is a problem in fitting such a function to sky data that cover a limited
range of airmass.
Except for optical depths approaching unity (i.e., near-UV bands),
the maximum in the function lies well beyond the range of airmasses usually
covered by photometric observations.
That means that the available data usually do not sample the large values of M
at which the squared terms become important.
Thus, one can usually choose rather arbitrary values of these terms,
and just fit the well-determined linear terms.
It turns out that choosing $b = 0$ is often satisfactory. If the data bend
over enough, $d$ can be determined; otherwise, it defaults to 0.01 times the
square of the largest airmass in the data.

An example of data that extend to large enough airmass to determine all four
parameters is the dark-sky brightness table published by Walker~\cite{W87},
which were used by Garstang to check his models.
These data extend to $Z = 85^{\circ}$.
The model above fits them about as well as do Garstang's models; typical errors
are a few per cent.
Various subsets, with the largest-$Z$ data omitted, give similarly good fits.
This indicates that the model is adequate for our purposes here.

In principle, one could add separate terms for zodiacal and diffuse galactic
light that depend on the appropriate latitudes; but this seems an excessive
complication, as these components vary with wavelength and longitude as well.
We have also neglected the ground albedo.  Unless the ground is covered with
snow, this is a minor effect except near the horizon.
Furthermore, the airglow can vary by a factor of 2 during the night; so we
cannot expect a very tight fit with any simple formula.

The moonlight term is more complicated.
In principle, it consists of single scattering, which in turn depends on
the size and height distributions of the aerosols, as well as Rayleigh
scattering from the molecular atmosphere; and additional terms due to
multiple scattering.
The radiative-transfer problem is complicated further by the large
polarization of the Rayleigh-scattered component, which can approach 100\%.
Rather than try to model all these effects in detail, we adopt a simple
parametric form that offers fair agreement with observation, but does not have
too many free parameters to handle effectively.

First of all, van de Hulst \cite{vdH} points out that the brightness of the
solar aureole varies nearly inversely with elongation from the Sun.
We assume the lunar aureole has the same property.
And, for
the small optical depths we usually encounter, and the moderate to small
airmasses at which we observe, we can simply assume the brightness of the
lunar aureole is nearly proportional to airmass.

Second, interchanging the directions of illumination and observation would
give geometries related by the reciprocity theorem if the ground were black.
For typical ground albedoes, we can still assume approximate reciprocity.
We can also assume $cos \, Z = 1/M$, where $M$ is the airmass, accurate to a few
per cent for actual photometric data, in calculating elongations from the Moon.

The adopted form is
$$
B_2 = M ( a/E + b + cE) \cdot [ \exp ( -dS) + e/P ] ,
$$
where $M$ is the airmass in the direction of observation, $E$ is the
angular elongation from the Moon, $S$ is the sum of observed and lunar
airmasses, and $P$ is their product.
(Note that the lower-case parameters here are different from the ones in the
dark-sky model.)

The factor in parentheses mainly represents the single-scattering phase
function, and should be nearly constant in good photometric conditions.
It is plotted as a ``normalized'' sky brightness.
Its parameter $a$ is a measure of the lunar aureole strength; if $a/b$  
is large, you probably have non-photometric conditions.
The factor in brackets handles the reciprocity effects.
The $e/P$ term produces the correct asymptotic behavior for a
homogeneous atmosphere; however, it cannot represent ``Umkehr'' effects at
wavelengths where ozone absorbs strongly.

Both factors have the symmetry required by the reciprocity theorem.
This condition may be violated by ground-albedo effects, and by photometers
that have a large instrumental polarization.

Unfortunately, when the telescope is pointed so that the Moon can shine on the
objective, the apparent aureole
is nearly always dominated by light scattered from
the telescope optics, not from the atmosphere.
Even if the mirror has been freshly aluminized, the scattered light may not
be negligible, because of surface scattering.
This scattering has different angular and wavelength dependences from those 
of the true sky brightness, and so is not represented by the sky model.
This means we should avoid observing so near the Moon that it can
shine directly on the primary mirror.

However, we cannot avoid having the {\it star} shine on the mirror; so we must
include a term in the ``sky'' model that is
proportional to the brightness of the measured star, to allow for the wings of
the star image that we include in our ``sky'' measurements.
Even if this fraction is so small (say, below $10^{-3}$) as to have no effect
on the sky-subtracted stellar data, it can be important in the sky {\it model},
if fairly bright stars are observed.
The fraction of starlight measured depends on the distance from the star to the
sky area chosen; as mentioned before,
one must keep this distance fixed for all stars when observing.

\subsubsection{Using the sky model}
If you elect to model the sky, you will see some graphs showing the 
progress of the fitting procedure.
First, if some bright stars were observed, you see a plot of sky data (for a
given band and night) as a function of the adjacent star's brightness.
A crude linear fit, indicated by + signs, shows the initial estimate of the
star contribution to ``sky'' data.
Second, if there are enough dark-sky data to model, you will see a plot of the
dark-sky data (corrected for star light) as a function of airmass, with the
fitted $B_1$ term drawn in as + signs.
Both these plots show instrumental intensity units on the vertical scale.

Next, if there are enough data with the Moon above the horizon to fit the $B_2$
term, you will see a plot of the moonlit sky, corrected for both the stellar
and the dark-sky terms, as a function of elongation from the Moon.
To show the importance of the lunar aureole (whose presence is an indicator of
considerable aerosol scattering, and hence probably non-photometric
conditions), this plot is normalized by multiplication by the factor in square
brackets (cf. the equation for $B_2$ in the previous subsection) and division
by the airmass; that is, it is simply a plot of the $E$-dependent factor in
parentheses.
Again, the fit is shown with + signs.  If the subtraction of the stellar and
dark-sky terms produced some apparently negative intensities, the calculated
zero level for the $B_2$ term will be drawn as a horizontal line.

If the sky is good, this plot will be nearly horizontal.
Usually, it bends up near 0 and 3 radians, with a minimum near 1.7; your data
may not cover this whole range, so pay attention to the numbers on the
horizontal scale.
The vertical scale is chosen to make the range of most measurements visible, so
the zero level may be suppressed; look at the numbers on the vertical scale.

It often happens that the range of the independent variables in these fits is
inadequate to allow a full fit of all the parameters.
Reasonable starting values will be used for the indeterminate parameters.
The fitting strategy is to adjust the best-determined parameters first, and
release parameters in turn until nothing more can be done.
At each stage, the results of a fit are examined to see whether the values
determined are reasonable.
For example, most of the parameters must be positive; and the dimensionless
ones should not be much larger than unity.

You will be informed of the progress of these fitting attempts.
Do not be alarmed by ``error messages'' like {\tt DLSQ FAILED} or
{\tt SINGULAR MATRIX} or {\tt Solution abandoned}.
It is quite common for solutions to fail when several parameters are being
adjusted, if the data are not well distributed over the sky.
Sometimes, when a solution fails, the fitting subroutine will check to make
sure the partial derivatives were computed correctly.
This is a safety feature built into this subroutine, and you should always find
that the error in the derivatives is $10^{-6}$ or less.
The program should then comment that the model and data are incompatible,
because we are usually trying to fit more parameters than the data can
support.
Remember that the program will adopt the best solution it can get; so watch for
messages like
{\tt 3-parameter fit accepted}, and don't worry about the failures.

Pay more attention to the graphs that compare the fits to the data.
Are there regions where they become widely separated?
If so, the fit is poor, and you can forget about the model.
If the fit looks good, the model is useful for detecting bad sky data, and may
even be useful for interpolating missing sky measurements.

When the fitting is complete, the program will print the terms used in the fit.
Then, three summary graphs display the quality of the fit.
The moonlit data are marked M in each of these three plots.
The first shows the observed sky brightness as a function of the model value.
This plot should be a straight line.

The second diagnostic
plot shows the {\it residuals} of the sky fit as a function of the
adjacent star brightness.
These points should be clustered about the axis indicated by dashes.
If you tend to measure sky farther from bright stars than from faint stars, as
some beginners tend to do, the points will show a downward trend toward the
right.
That's a sign that you need to be more careful in choosing sky positions.
Likewise, a large scatter in this plot probably means you have been careless in
choosing sky positions, sometimes measuring closer to the star and sometimes
farther away.
(Here, ``large'' means large compared to the typical sky values on previous
plots.)

The final plot in this group shows
the {\it ratio} of the observed to the computed (modelled) sky brightness,
as a function of time.
If the airglow changes with time, you will see waves and wiggles in the
dark-sky portion.
Likewise, if the aerosol is varying with time, you will see coherent variations
in the moonlit portion.
The upper and lower limits of this plot are fixed in absolute value, so the
scatter visible here is a direct indication of the quality of the overall fit.

Finally, if there are aberrant points that do not fit the model, they will be
tabulated.
If the data are not well represented by the model, there will be many entries
in this table.
Pay particular attention to the last column, which gives the ratio plotted in
the last diagnostic graph.
If the fits were generally satisfactory,
the few sky measurements tabulated here may
be in error, and may indicate a instrumental or procedural problem.
They should be examined carefully to determine the cause of the problem.

After all this information has been presented, you will be asked whether you
want to subtract the modelled sky from the stellar data.
You can reply Yes or No, or enter R to repeat the whole process, or H to ask
for help.
If you reply Yes, the model values will be subtracted from all the stellar data
that were not taken during twilight.
However, as twilight is not modelled, the nearest-neighbor method will be used
to correct stars observed during twilight.
If you reply No, you will be given the option of subtracting the nearest
neighbor in all cases.

The sky models do not work well if only a few observations are made with
the Moon either above or below the horizon.
They do not handle solar or lunar twilight.
They can have difficulties if the observations are not well distributed over
the sky.
In these, as in some other cases discussed above, one should choose some other
method instead of using a model for sky subtraction.


\subsection{Choosing a gradient estimator}
\label{2-color}

Now the program knows which stars have been observed, and can show you a
two-color diagram of the available standards.
It displays a plot of the two nearest color indices for each band.
For example, in the $uvby$ system,
plots of $(u-v)$ vs. $(v-b)$ are displayed for
the $u$ and $v$ passbands, and plots of $(v-b)$ against $(b-y)$
for both $b$ and $y$.
The degree of correlation between the two displayed color indices is also
printed.
The goal is to select an appropriate estimator of the stellar spectral gradient
across each passband.

You are now asked to choose whether to use the
single default color index (which would
be $(u-v)$ for $u$, and $(u-b)$ for $v$), or {\it both} colors displayed on the
plot, as the basis for transformations.
In general, if the plot shows a high degree of correlation, adding the second
index adds little information, and may weaken the solution by allowing an
additional degree of freedom to be adjusted.
If the points cover an appreciable two-dimensional area, or follow a curved
line rather than a straight one, you may gain substantial accuracy in
transforming by using both indices as independent variables.
When in doubt, ask the program for help.

Note that you make this choice separately for each passband.
(It might make sense to use two colors in reducing $b$, but only one in reducing
$y$, even though the same 2-color diagram is presented in each case.)
Once you have set this choice, it will be followed throughout the remainder of
the run.

Note that you will be given this choice only if there are enough standard stars
to make it a reasonable course of action.
With only a few standards, and too many adjustable parameters,
there is the danger that the program may simply fit a function to
a few points and ignore the rest.

In principle, these plots should employ the instrumental
rather than the standard colors of the standard stars.
Unfortunately, the choice of a gradient estimator
must be made before extinction-corrected instrumental colors are available,
so the standard colors must stand in for the instrumental ones.
In practice, this should not be a problem, unless the instrumental system is a
very poor match to the standard one.
For example, if the B band of a UBV photometer includes too much ultraviolet,
the instrumental B-V index may contain appreciable contamination from U that
would not be seen in the standard 2-color plot.


\subsection{Extinction and transformation models}

To determine the parameters we want (stellar magnitudes and colors) from the
observations we have, we fit the data to a model.
The model should represent the circumstances of the observations as
realistically as possible; otherwise, effects present in the data but not in
the model will be distributed among the available adjustable parameters, such
as stellar magnitudes, extinction coefficients, and transformation parameters.
The result may well be a good fit to the data, with small residuals; but the
parameters will be biased --- that is, they will have significant systematic
errors.
Such errors are called {\it reduction errors} by Manfroid and
Sterken~\cite{MS92}


\subsubsection{The bias problem}

Bias can also be introduced by improper fitting procedures.
For example, the reduction program known as SNOPY, used for many years at La
Silla and also at La Palma, uses a simple iterative technique in which some
subsets of the parameters are solved for, while others are held fixed; then the
ones last solved for are held fixed, while the ones previously fixed are
redetermined.
The earliest version of PEPSYS used this technique \cite{YI67}.
However, it required {\it hundreds} of iterations to approach the true minimum
in
the sum of squares, and consequently was abandoned as soon as computers became
large enough to allow simultaneous solutions for all the parameters.

A particularly treacherous aspect of such subspace alternation is that the sum
of squares decreases rapidly for the first few iterations and then levels off.
If the iterations are performed by hand, as they are with SNOPY, the user is
likely to think that the system has converged when both the residuals and the
parameter changes have become small.
Nevertheless, the current values of the parameters can be far from the desired
solution.

What happens is that the solution finds its way to the floor of a long, narrow
valley in the sum-of-squares hypersurface in the $n$-dimensional parameter
space.
Even with optimal scaling, such valleys occur whenever the independent
variables are partially correlated, so that the parameters themselves are
correlated.
This nearly always happens in multiparameter least-squares problems.
At each iteration,
the descent to the valley floor (i.e., the minimum in the mean-square
error) occurs within the subspace of the parameters that are adjusted at each
iteration, but is necessarily at the starting values of the parameters that
were held fixed.
At the next iteration, the solution moves a short distance in this
orthogonal subspace, and again finds a point on the valley floor; but, again,
it only finds a {\it nearby} point on the axis of the valley, and is unable to
progress {\it along} the axis, because of the parameters that are held fixed.
Succeeding iterations take very small steps, back and forth across the valley
axis, while making only very slow progress toward the true minimum at the
lowest point on that axis.
This behavior is notorious in the numerical-analysis community, where it is
known as ``hemstitching''.

The fatal slowness of convergence along coordinate
directions is explicitly described
and illustrated on pp.\ 294-295 of {\it Numerical Recipes}~\cite{NR}, although
the schematic diagram shown there does not show how slow the process really is.
Actual computations show that hundreds or
even many thousands of iterations may not be
enough to get acceptably close to the true minimum: see, for example, Figure 4j
of Gill et al.~\cite{GMW} or Figure 9.7 of Kahaner et al.~\cite{KMN}.
But in every iteration after the first few, the sum of the squared residuals
is small and changes very slowly.
One must not forget that obtaining small residuals is not an end in itself, but
is merely the means of finding what we really want, namely, the best estimates
of the parameters.

The most effective way to prevent hemstitching is to adjust all the parameters
simultaneously.
In photometric problems, this means solving a rather large system of
simultaneous equations.
For example, if we have 100 stars observed in 4 colors on 5 nights, there are
400 different stellar parameters, 25 extinction coefficients and nightly
zero-points, and a few additional transformation coefficients to be determined,
if they are estimated simultaneously.
In broadband photometry, there are also parameters proportional to the second
central moments of the passbands to be estimated, which introduce cross-terms
and hence nonlinearity.
Fortunately, the nonlinearity is weak, and can be very well handled by standard
techniques (cf.~\cite{DS}) that converge quickly.

The problem can be speeded up further by taking advantage of its special
structure: the matrix of normal equations is moderately sparse, with a
banded-bordered structure; the largest block (containing the stellar parameters)
is itself of block-diagonal form.
Thus matrix partitioning, with a block-by-block inversion of the
stellar submatrix, provides the solution much more quickly than simple
brute-force elimination of the whole matrix.
A single iteration takes only a few seconds, in typical problems.


\subsubsection{Equations of condition}

Assured of a reliable method of solving the normal equations, we now consider
the equations of condition.
The observed magnitude $m$ of a star at $M$ air masses
with true instrumental magnitude $m_0$ is modelled as (\cite{Y74}, \cite{Y91})
$$
m = m_0 + A_{eff} M + Z_n ,
$$
where $Z_n$ is a nightly zero-point.
The effective extinction coefficient
$A_{eff}$ depends on the effective spectral gradient \cite{Y92} of the star,
for which we use some color index, affected by half the atmospheric reddening
(\cite{YI67}, \cite{Y74}):
$$
A_{eff} = A_0 - W ( C_0 + RM/2 ) .
$$
Here $A_0$ is the extinction coefficient for the effective wavelength of the
passband;
$W$ is proportional to the second central moment of the instrumental passband;
$C_0$ is the extra-atmospheric color of the star (obtained by
differencing the $m_0$ values in the appropriate pair of passbands);
and $R$ is the atmospheric reddening per unit airmass, obtained by differencing
the corresponding pair of $A_0$ values.

On general numerical-analysis grounds, we initially assume the best pair of
passbands to use for $C_0$ and $R$ to be the pair that flank the passband in
which $m$ was observed.
Thus, for the B band of UBV, we use the (U-V) color.
This is not conventional practice; however, it provides a general rule that
works reasonably well for all photometric systems.
Such a rule is required in a general-purpose reduction program.
For a well-sampled system, this produces a strikingly accurate representation
of the extinction correction \cite{Y92b}.
For undersampled systems,
we can use a somewhat more accurate gradient estimator, using a linear
combination of the two adjacent color indices, as described above
(see subsection~\ref{2-color}, ``Choosing a gradient estimator'').
For bands at the extreme wavelengths of our system, we adopt just the
neighboring color index (e.g., (U-B) for the U magnitude and (B-V) for the V
magnitude of UBV).

Regarding the use of (U-V), we may note Bessell's remark~\cite{B90}
(in connection with
the problems created by the wide and variable B band) that ``in retrospect, it
would have been much better had (U-V) rather than (U-B) been used by Johnson in
his establishment of the UBV system.''
In any case, there are special problems in UBV, both due to the inadequate
definition of the original system, and the neglect of the transformation from
inside to outside the atmosphere entirely in the (U-B) index, which made the
``B'' of (B-V) and the ``B'' of (U-B) have different effective wavelengths; in
principle, it is incorrect to add the two indices to a V magnitude and come up
with a ``U magnitude'' as a result.
Consequently, one must be very careful in doing UBV photometry, no matter how
it is reduced; the only safe course is to observe at least the 20 standard
stars recommended by Johnson (more would be preferable), and to look very
carefully for systematic trends in the transformation residuals.

Although, for astrophysical reasons, there may be partial correlations
of the true stellar gradient at a given band with color indices that are
remote in wavelength
from the band in question, such correlations will depend on metallicity,
reddening, and other peculiarities of individual stars.
If correlations obtained from one group of stars are applied to another group,
the results may easily be worse than if the additional correlation had been
ignored in the first place.
Thus, it is exceedingly dangerous to employ distant
color indices unless the calibration and
program stars are very similar in all these respects.
We cannot rely on such good matching in general, so these
partial correlations are not used in the present package.
However, one can expect that very good results will be obtained if the program
and extinction/standard stars cover the same region of parameter space.
That is, they should have the same ranges of spectral types, reddening,
chemical composition, etc.

Because filter photometers observe different bands at different times, we have
to do the reduction in terms of magnitudes (which are measured),
rather than colors (which are not ---  see \cite{Y74},
pp.~152 -- 154).
This also allows isolated observations of extinction stars in a single filter,
or a partial subset of filters, to be used.
Furthermore, the best estimate of the extra-atmospheric color $C_0$ is used to
reduce every observation, so that errors in individual measurements have the
least effect on the reduction of each particular observation.

Of course, it may be necessary to allow the extinction $A_0$ to be a function
of time; we always solve for individual values on different nights, as many
careful investigations have demonstrated the inadequacy of ``mean extinction''.
The instrumental zero-point terms $Z$ may likewise be functions of time,
temperature, or relative humidity.
Finally, the zero-points may be kept the same in each passband for all nights
in a run if the instrument is sufficiently stable.

\subsubsection{Strategy}

Bear in mind that you should have {\it enough} standard stars, but not too
many.
It is better to have 5 observations each of 30 standard/extinction stars,
than to have 3 observations each of 50 stars; in the latter case, we are
spreading the same observational weight over nearly twice as many parameters.
Because one observation is ``used up'' in determining the stellar values, only
$(n-1)$ of $n$ observations per star are available for measuring extinction.
Thus a star observed 5 times contributes twice as much to the extinction
determination as a star observed only 3 times, which in turn contributes twice
as much as a star observed only twice.
The planning program will help you choose the right number of stars.

Furthermore, the running time of the reduction
program is roughly proportional to 
the number of observations, but to the {\it cube} of the number of parameters
to be adjusted.
Thus, the program will take almost 5 times as long to run if there are 50 stars
to be adjusted than if there are only 30.
The storage required is proportional to the square of the number of parameters,
and this may limit the number that can be handled on machines with small
memories.
At some point, any machine becomes overloaded; there is a practical limit to
the number of stars that can be adjusted in a single solution.
The reduction program automatically excludes non-standard stars that have been
observed only once from the extinction solution.

In sum, though one can hardly have too many {\it observations} of standard
stars, one should not go beyond the optimum {\it number of stars} recommended
by the planning program.
Some observers who follow Hardie's~\cite{H62} method observe a large number of
standard stars just once each; this is extremely inefficient use of telescope
time.

Although these questions ought to be considered before one goes to the
telescope, some observers only consider them when trying
to reduce data that have been badly distributed.
Even a sophisticated reduction program cannot generate information that was not
gathered at the telescope.

\subsubsection{Standard stars}

Standard stars can be included in the extinction solution, if one is careful.
The problem is that there are usually appreciable ``conformity errors'' due to
mismatch between the instrumental and the standard passbands \cite{MS92}.
These are due to the integrated influence of features in the spectra of
individual stars, which differ in age, metallicity, reddening, etc.\ at a given
color index value.
Thus, it is wrong to assume that the transformation between the instrumental
system and the standard system is exact, even if there were no measurement
error in either system \cite{BN89}.

Usually, one finds that the instrumental system can be reproduced
internally with much more
precision than one can transform between it and the standard system.
This means we should regard the standard values as imprecise estimates of the
extra-atmospheric instrumental values,
after transforming from standard to instrumental system.
In effect, we regard the standard values as noisy pseudo-observations of the
instrumental, extra-atmospheric magnitudes.

Therefore, a reasonable equation of condition to use for the standard values
is:
$$
m_{std} = m_0 + a C_0
$$
where we effectively regress the noisy (because of conformity errors) standard
values on the internally precise instrumental system.
No zero-point term appears, as we already took it out in conjunction with
the extinction.
As before, $m_0$ is the extra-atmospheric instrumental magnitude,
and $C_0$ is the estimator of stellar gradients across the passband
being transformed; the same gradient estimator should be used in both
extinction and transformation equations.
The transformation coefficient $a$ is then estimated in the general solution.
In general, we expect a higher-order polynomial in $C_0$ to be
needed~\cite{Y92}; it can be determined if enough standard stars are available.
This equation then gives an explicit transformation from the instrumental to
the standard system.

One problem is to know what weights to assign such pseudo-observations.
The conformity errors in the transformation are usually much larger than the
internal errors of measurement in the standard system, so that quoted
uncertainties in well-observed standard stars are usually
irrelevant to our problem.
For convenience, unit weight in the extinction solution is intended to
correspond to errors of 0.01 magnitude.
This is a typical order of magnitude of conformity errors as well, so we can
start with unit weight for these values, and adjust the weights after a trial
solution, followed by examination of the residuals.
This process assigns self-consistent weights to the standard values.

Alternatively, we can omit the standard stars from the extinction solution,
determine the extinction entirely from the observations, and then determine the
transformation in a separate step.
In this case, a separate set of zero-points must be determined in the
transformation solution, as the nightly zero-points in the extinction solution
are no longer coupled to the transformation.
However, this does not increase the number of unknowns, as we must then
fix one night's zero points to prevent indeterminacy.
(In principle, one should do this using Lagrange multipliers; in practice, it
hardly matters, because the nightly zero points are always determined much more
precisely than other parameters in the solution.)
This process is preferable if there are enough extinction observations, as it
does not propagate systematic errors from the transformation into the
extinction (see pp.~178 -- 179 of \cite{Y74}, and \cite{P82}).

However, observers sometimes do not get enough extinction data to permit solving
for extinction directly.
Those who follow Hardie's~\cite{H62} bad advice sometimes observe each
standard star just once; then, if there are no real extinction stars,
the standard values have to be used to obtain any extinction solution at all.
But usually,
in combined solutions, the extinction coefficients will be slightly
more precise, but less accurate, than in separate solutions.
The reduction program will tell you if one method or the other seems
preferable; in any case, it lets you decide which method to use.

There has been much splitting of hairs in the astronomical literature over the
direction in which the standard-star regression should be performed.
Unfortunately, the arguments have all assumed that the regression model is
functional rather than structural; that is, they assume there is an exact
relation connecting the two systems, in the absence of measurement noise.
In practice, conformity errors are usually larger than measurement errors, so
the functional-regression model is incorrect.
In any case, the problem we have here is a calibration problem: given
measurements of some stars in the standard system A and the instrumental system
B, we want to predict the values that {\it would} have been observed in A
from the values we {\it have} observed in B, for the program stars.
In this case, the regression of A on B provides the desired relationship
\cite{FB92}.

While the conformity errors make it reasonable to do the regression
in the sense described above, it is also clear that they involve significant
effects that are not accounted for in the usual treatments of transformation.
In particular, substantial {\it non-linear} terms are to be expected
in transformations, and the neglect of higher-order cross-product terms in
general makes transformations appear to be {\it multivalued}~\cite{Y92}.
The correct treatment of these problems requires a well-sampled system.
Trying to correct for missing information after the fact is a lost cause;
effort should instead go into minimizing the conformity errors in the first
place.

In general, one cannot emphasize too strongly the need to measure the
instrumental spectral response functions, and to choose filters --- by
actual transmission, not by thickness! --- that reproduce the standard
system as accurately as possible.
The difference between instrumental and standard response functions represents
a serious disparity that cannot be made up by any reduction algorithm, no
matter how cleverly contrived.


\subsection{Reduction procedure}
\label{reductionproc}

After subtracting dark current and sky brightness from your stellar data,
and asking you to select a gradient estimator for each band,
the reduction program converts intensities to magnitudes.
At this stage, it warns you of any observations that seem to have zero or
negative intensities (which obviously cannot be converted to magnitudes).
These are often an indication that you have confused sky and star readings.
The program then estimates starting values for the full solution.

The program treats four different categories of stars differently:
standard stars, extinction stars, program stars, and variable stars.
It asks for the category of each file of stars when the positions are
read in.
{\it Standard} stars are those having known standard values, which are used
only to determine transformation coefficients.
They may also be used to determine extinction.
Be cautious about using published catalog values as standards; these often
have systematic errors that will propagate errors into everything else.
{\it Extinction} stars are constant stars, observed over an appreciable range of
airmass, whose standard magnitudes and colors are unknown, or too poorly
known to serve as standards.
Ordinary stars taken from catalogs of photometric data are best used as
extinction stars; later, you can compare their derived standard values with the
published ones as a rough check.
{\it Program} stars may be re-classified as either extinction or variable stars
during the course of the extinction solution.
If they are not used as extinction stars, they are not included in the
solutions, but are treated as variable stars at the end.
{\it Variable} stars are excluded from the extinction and transformation
solutions.
Their individual observations will be corrected for extinction and
transformation after the necessary parameters have been evaluated.

You will do best to maintain separate star files for each category; however,
you can intermix extinction and variable stars in a ``program star'' file, and
PEPSYS will do the best it can to sort them out, if you ask it to use program
stars for extinction (see below).
Remember that only star files need to be separated this way; a data file
normally has all the observations for a given night, regardless of the type of
object.

You can also group related files in a MIDAS ``catalog'' file.
Use the MIDAS command {\tt CREATE/TCAT}
to refer to several similar *.tbl files as a catalog file.
Note that (a) you must enter the ``.cat'' suffix explicitly when giving PEPSYS a
catalog; and (b) a catalog {\it must} contain only tables of the same kind ---
e.g., only standard stars, or only program stars.
Catalogs can be used for both star files and data files.

When all the star files have been read, the program asks for the name of the
table file that describes the instrument.
As usual, you can omit the ``.tbl'' suffix and it will be supplied.

Then the program asks for the data files.
Remember to add the ``.cat'' suffix if you use a catalog.
While reading data, it may ask you for help in supplying cross-identifications
if a star name in the data did not occur in a star table.
If all your data files have been read, but it is still asking for the name of a
data file, reply {\tt NONE} and it will go on.
(This same trick works for star files too.)

The program
will display a plot of airmass vs.~time for the standard and extinction
stars on each night, so you can judge whether it makes sense to try to solve for
time-dependent extinction later on.
If the airmass range is small, it will warn you that you may not be able to
determine extinction.

If the data are well distributed, and there are numerous standard stars,
it will obtain starting values of extinction coefficients from the standard
values; otherwise, it assumes reasonable extinction coefficients and estimates
starting values of the magnitudes from them.
(These starting values are simple linear fits, somewhat in the style of SNOPY,
except that robust lines are used instead of least squares.)
From the preliminary values of magnitudes, it tries to determine transformation
coefficients, if standards are available.

This whole process is iterated a few times to obtain a self-consistent set of
starting values for all the parameters, except for bandwidth factors.
Each time the program loops back to refine the starting values, it adds a line
like ``{\tt BEGIN CYCLE \ n}'' to the log.
Don't confuse these iterations, which are just to get good starting values,
with the iterations performed later in the full solution.

One problem in this preliminary estimation is that faint stars with a lot of
photon noise might just add noise to the extinction coefficients.
The program tries to determine where stars become too faint and noisy to be
useful for estimating extinction.
The rough extinction and transformation
coefficients will be displayed as they are determined; if any fall outside a
reasonable range of values, the program will tell you and give you a chance to
try more reasonable values.

When reasonable starting values have been found for the parameters that will be
estimated in the full solution, the program still needs to decide how to do
the solution.
Should the standard values of the standard stars be used in determining
extinction (which requires estimating transformation coefficients
simultaneously), or should extinction be determined from the observations
alone, and the transformations found afterward?
The program will ask your advice.

If you have designated no extinction stars,
the reduction program will ask you whether you want to treat program stars as
extinction stars.
If you believe most of them are constant in light, you can try using all of
them as extinction stars.
If some turn out to be variable, you will have a chance to label them as such
later on.
If many of the program stars are faint, they may be too noisy to contribute
anything useful to the extinction solution; then you could leave them out and
speed up the solution.
On the other hand, values obtained from multiple observations in the general
solution will be a little more accurate than values obtained by averaging
individual data at the end.

When the program is ready to begin iterating, it will ask how much output you
want.
When you first use PEPSYS, you may find it useful to choose option 2 (display
of iteration number and variance).
After you get used to how long a given amount of data is likely to take to
reduce, you can just choose option 1 (no information about iterations).
The detailed output from the other options is quite long, and should
only be requested if the iterations are not converging properly and you want
to look at the full details.

If you ask for iteration output, you will always see the value of the weighted
sum of squares of the residuals (called WVAR in the output), and the value of
the typical residual for an observation of unit weight (called SCALE), which
should be near 0.01 magnitude.
The weights are chosen so that the error corresponding to unit weight is
intended to be 0.01 mag.
When SCALE changes, there can be considerable changes in WVAR --- in
particular, don't be alarmed if it occasionally increases.
A flag called MODE is also displayed, which is reset to zero when SCALE is
adjusted.
During normal iterations, MODE = 4; ordinarily, most of the iterations are
in mode 4, with occasional reversions to mode 1 when SCALE is adjusted.
More detailed output contains the values of all the parameters at each
iteration, and other internal values used in
the solution; these are mainly useful for debugging, and can be ignored by the
average user.

Typically 20 or 30 iterations are needed for convergence; if the data are
poorly distributed, so that some parameters are highly correlated,
more iterations will be needed.
If convergence is not reached in 99 iterations, the program will give up and
check the values of the parameters (see next paragraph).
This usually indicates that you are trying to determine parameters that are not
well constrained by the data.
It will also check the values of the partial derivatives used in the problem;
if there is an error here larger than a few parts in $10^6$, you have found a
bug in the program.

At the end of a solution, the program prints the typical error, and then
examines the parameters obtained.
If extinction or transformation coefficients or bandwidth parameters seem
{\it very} unreasonable, the program will simply fix them at reasonable values
and try again.
If they seem marginally unreasonable, the program will ask for your advice.

When reasonable values are obtained for all the parameters, the program will
check the errors in the transformation equations, if standard values have been
used in the extinction solution.
Then, if necessary, it will readjust the weights assigned to the standard-star
data, and repeat the solution.
Usually 3 or 4 such cycles are required to achieve complete convergence.
Thus, including the transformation parameters
in the extinction solution means the program
may take longer to reach full convergence.

Having reached a solution in which the weights are consistent with the
residuals, the program examines the stellar data again.
If you have ``program'' stars that might be used to strengthen the extinction
solution, it will ask if you want to use all, some, or none of them for
extinction.
If you reply SOME, it will go through the list of stars one by one and offer
those that look promising; you can then accept or reject each individual star
as an extinction star.
Only stars whose observations cover a significant interval of time and/or
airmass will be offered as extinction candidates.

In examining the individual stars, the program may find some that show signs of
variability.
For those that have several observations, a ``light-curve'' of residuals will
be displayed.
Pay close attention to the numbers on the vertical scale of this plot!
Each star's residuals are scaled to fill the screen.
If there are only a few data for a star, only the calculated RMS variation is
shown.

A star may show more variation than expected, but if this is under 0.01
magnitude, it may still be a useful extinction star --- indeed, the star may
not really be variable, but may have had its apparent variation enhanced by one
or two anomalous data points.
Another problem that can occur, if you have only a few extinction stars, or
only one or two nights of data, is a variation in the extinction with time.
This can produce a drift in the residuals of a standard or extinction star
that may look like some kind of slow variation.
Watch out for repeated clumps of anomalously dim measurements that occur about
the same time for one star after another; this often indicates the passage of a
patch of cirrus.

The program will show you several doubtful cases for every star that really
turns out to be variable, so be cautious in deciding to treat them as variable
stars.
If you decide that a star really looks variable, you can change its category to
``variable'', and it will be excluded from all further extinction solutions.

If you change the category of any star, the program goes back and re-does the
whole extinction solution again from the very beginning.
When you get through a consistent solution without changing any star's
category, the program announces that it has done all it can do, displays
some residual plots, and prints the final results.

All reductions are done in the instrumental system, even if
standard-star values (and hence transformations) are included as part of the
extinction solution.
It may be helpful to look at a schematic picture of the reduction process
(see Figure~\ref{flowchart} on the next page).


\begin{figure}[htp]
\begin{center}
\begin{picture}(120,180)

\put(0,170){\makebox{Pre-process data:}}
	\put(60,172){\makebox(0,0){\framebox{
				\begin{minipage}{32mm}
				Subtract Dark, Sky, and Red leaks
				\end{minipage}}}}
	\put(60,165){\vector(0,-1){11}}
\put(0,148){\makebox{Initialize solution:}}
				\put(105,149){\vector(-1,0){35}}
	\put(60,149){\makebox(0,0){\framebox{START}}}
	\put(60,145){\vector(0,-1){10}}
	\put(60,130){\makebox(0,0){\framebox{Include std./program stars?}}}
\put(0,110){\makebox{Get starting values:}}
	\put(60,125){\vector(0,-1){10}}
	\put(60,107){\makebox(0,0){\framebox{
				\begin{minipage}{25mm}
				Use std. values or mean ext.?
				\end{minipage}}}}
\put(50,100){\vector(-1,-1){10}}
\put(40,98){\makebox(0,0){std.}}
\put(0,71){\makebox{Cycle:}}
\put(40,83){\makebox(0,0){\framebox{
			\begin{minipage}{17.8mm}
			find ext. from stars
			\end{minipage}}}}
	\put(52,85){\vector(1,0){8}}
\put(40,67){\vector(0,1){10}}
\put(40,61){\makebox(0,0){\framebox{
			\begin{minipage}{22mm}
			check photon noise level
			\end{minipage}}}}
		\put(70,100){\vector(1,-1){10}}
		\put(80,98){\makebox(0,0){mean}}
		\put(80,85){\makebox(0,0){\framebox{correct stars for ext.}}}
		\put(80,80){\vector(0,-1){10}}
		\put(80,65){\makebox(0,0){\framebox{find approx.~transform.}}}
		\put(80,60){\vector(0,-1){10}}
		\put(80,45){\makebox(0,0){\framebox{good starting values?}}}
	\put(55,47){\makebox(0,0){no}}
	\put(60,45){\vector(-1,0){10}}
	\put(40,50){\vector(0,1){5}}
	\put(50,50){\oval(20,10)[bl]}
		\put(80,40){\vector(0,-1){10}}
		\put(85,35){\makebox(0,0){yes}}
\put(0,25){\makebox{Full solution:}}
		\put(80,25){\makebox(0,0){\framebox{Do full iterations}}}
		\put(80,20){\vector(0,-1){10}}
		\put(80,5){\makebox(0,0){\framebox{Results OK?}}}
	\put(55,7){\makebox(0,0){yes}}
	\put(65,5){\vector(-1,0){20}}
\put(21,5){\makebox(0,0){\framebox{
			\begin{minipage}{42mm}
			Reduce individual data, tabulate results. DONE.
			\end{minipage}}}}
				\put(105,7){\makebox(0,0){no}}
				\put(95,5){\vector(1,0){20}}
					\put(105,77){\oval(34,144)[r]}

\end{picture}
\caption[]{Schematic flowchart for reductions}
\label{flowchart}
\end{center}
\end{figure}


\subsection{Final results}

\subsubsection{Residual plots}

First, there is a set of residual plots against airmass for each night.
The plots for individual passbands are compressed so that two will fit on the
screen at once; this allows you to notice correlations in the behavior of
residuals in adjacent bands.

Next, there are residuals as functions of time for each night (again compressed
to make comparisons easy).
Examine these plots carefully for evidence of trends with time.
After you have inspected the run of the residuals with time, the program will
ask if you want to solve for time-dependent parameters.
If you reply ``YES'', it will show you some plots intended to help you
distinguish between time-dependent extinction and time-dependent zero-point
drift.
These plots are accompanied by some numbers that may also be helpful (see
Fig.~\ref{driftplot}).


\begin{figure}[htb]
\begin{verbatim}
                          v drift plot     slope = 1.50   32 data
       differential offset
         X  I
        .01 I
            I                          *
            I
         -1.+
            I                                *
            I                   *   *  *                   *
            I
            I                 *    *  $               +++++++++++++
          0.+   - - -  - - ++++++++++*+++++*++++++++++  - - -
            I  ++++++++++++ *                * *
            I   *   *                 $*       $
            I      *                    *        *
            I           *
          1.+
            I                  *
            I
            I        *
            I
          2.+---+---+---+---+---+---+---+---+---+---+---+---+---+--
            14  12  10   8   6   4   2   0  -2  -4  -6  -8 -10 E-04
            mean offset  -->     ratios:  0.089  0.083
\end{verbatim}
\caption[]{Sample drift plot}
\label{driftplot}
\end{figure}


The principle behind these plots is that time-dependent extinction produces
errors (and hence, we hope, residuals) that are proportional to airmass;
but instrumental drift produces residuals uncorrelated with airmass.
Suppose we have two successive observations at airmasses $M_1$ and $M_2$,
which produce residuals $\epsilon_1$ and $\epsilon_2$.
If the extinction changes slowly with time, we expect the ratios
$ (\epsilon_1 - \epsilon_2 ) / (M_1 - M_2 )$ and
$ (\epsilon_1 + \epsilon_2 ) / (M_1 + M_2 )$
to be equal, as each is (on the average) the deviation of the current
extinction coefficient from the mean extinction for the night, multiplied by
the airmass difference of the points being compared.

Unfortunately, a plot of the first ratio against the second produces an
uninformative scatter diagram, because the divisor $ (M_1 - M_2 )$
makes the scatter of the points very unequal.
If, instead, we plot
$ (\epsilon_1 - \epsilon_2 )$ against
$ (\epsilon_1 + \epsilon_2 ) (M_1 - M_2) / (M_1 + M_2 )$,
the values should all have comparable scatter.
This is the plot shown in Fig.~\ref{driftplot}.
The line drawn with {\tt +} signs has the fitted slope and passes through the
origin.
Again, we expect unit slope if the extinction is changing with time,
and zero slope if the instrument is at fault.
The actual slope is printed above the plot.

We can obtain additional information from the scatter of the values.
The vertical spread of the plotted points
is just proportional to a typical residual.
The horizontal spread has the same value, but diminished by a typical value of
the ratio
$ (M_1 - M_2) / (M_1 + M_2 )$.
This average ratio, followed by the ratio of the spreads, is shown at the
bottom of the plot.
If the errors are dominated by changes in the extinction, we expect the
ratio of spreads to be just this average ratio.
But if instrumental drift dominates the errors, the spread in the ordinates
should be much smaller than the spread in the abscissas, because the
correlated part of the residuals cancels out in the differences.

If the drifts are small, so that the data are dominated by random noise, we
expect the mean slope to be zero, and the ratio of spreads to be the
average ratio of airmass difference to airmass sum,
$ (M_1 - M_2) / (M_1 + M_2 )$.
We can therefore take a mean slope greater than 0.5 as evidence in favor of
extinction drift, but a smaller slope is only weak evidence of zero drift, as
it could be due to noise.
Similarly, a spread ratio near zero is strong evidence for zero-point drift,
but a value near the mean airmass ratio could be due to noise.
On this basis, the sample data shown in Fig.~\ref{driftplot} 
show good evidence for extinction drift.

After showing the individual drift plots, the program displays a table
that summarizes the results for all bands and nights.
In this table, cases that favor extinction drift are marked ``E'' and
those that favor zero-point drift, ``Z''.
Lower-case letters indicate weaker evidence, and a dash indicates nights with
too few data to make a useful plot.

At this point, you will be asked to choose whether to solve for time-dependent
extinction, time-dependent zero points, neither, or both.
In general, you should choose the situation that seems most typical;
help is available if you want it.
It is usually not wise to try to solve for both kinds of drift; unless the data
are very well distributed, one usually finds a strong anti-correlation of the
derived drifts in extinction and zero-point.
At present, only linear drifts with time are modelled; but this should suffice
for most good nights.

If you choose to change the model (i.e., either add or remove a time-dependent
term, compared to what was previously used), the program loops back to the
beginning of the solution, and repeats its work, using the new model.
Initially, solutions are always done without time-dependent terms.
This strategy lets you choose an optimal set of extinction stars before trying
to solve for any time-dependent terms, which are usually rather weakly
determined.
If you make no changes in the time-dependence of the model, the extinction
solution is finished at this point.

If you included standard-star values in the extinction solution, the
transformation coefficients were included in the parameter adjustments.
Then full-screen plots of the transformation residuals are shown, followed by
residuals plotted as a linearity check.
If you did not include standard-star values, the linearity check is plotted,
and the program then determines the transformation separately.

\subsubsection{Estimated parameters}

After the last residual plot, the final parameter values
(and their errors) are printed, including
the magnitudes and colors of the individual stars in the extinction solution.
The errors for the color indices take account of the covariances between the
magnitudes.
This is followed by a table of the individual residuals, which includes the
symbols used to identify different stars in the residual plots.

\subsubsection{Individual data}

Finally, the individual observations themselves are reduced, using the
parameters from the extinction solution.
This is a little tricky, as one needs color values to include the color terms
in both extinction and transformation equations.
If all colors are observed simultaneously, as is done with multichannel
spectrometers, these colors can easily be extracted from the neighboring
data.
But a filter photometer or CCD must observe passbands sequentially.
Again, this may not be a problem, if the star is constant in light.

However, for variable stars, the brightness is changing in all bands, and it is
difficult to define the correct color to use in the reduction.
In principle, one should construct the full light curve of the star, and then
interpolate the colors to the time at which each passband was observed.
In practice, a simpler compromise is used: the colors are simply interpolated
linearly in time, if there are earlier and later observations of the star on
the same night.
For the first and last observations of a star on a given night, the nearest
observations (in time) in the adjoining bands are used.
In every case, the values of the colors used to reduce each magnitude are shown
on the output.

Every observation is given individually, including observations of standard and
extinction stars, together with the U.T.\ day fraction and heliocentric Julian
Date.
These data are of interest for standard and extinction stars if they later
turn out to be variable.
Note that the values of stellar parameters
adopted in the extinction solution are to be preferred to
values obtained by averaging these individual observations
that have been corrected for extinction and transformation.

NOTE:
Because colors are needed to transform the observations, first to outside the
atmosphere and then to a standard system, only stars that were observed in all
passbands on a given night can be reduced in this section of the program.
This is in contrast to the treatment of constant stars in the general solution,
where even single-passband measurements are useful, provided that (a) every
star was observed at least once during the run in every passband;
and (b) every night has some observations in all passbands.

That is, the extinction solution can in principle produce results if we observe
extinction star 1 only in passband A and star 2 in passband B on night 1, and
star 1 only in passband B and star 2 in passband A on night 2.
This fulfills the requirements for data in both bands for each night, and for
each star.
Even if the passbands are linked by color terms in the extinction, the solution
is possible, because we do obtain mean colors for each star.
However, none of these data could be reduced individually, as there are no
observations of either star in both bands on the same night.

\subsubsection{Logfile output}

All the on-line output goes into the MIDAS logfile.
Normally, you would print this out after a reduction run, to have a record of
what you have done.
This may run a hundred or more pages in length, so you may want to use a
reduced-size or ``two-up'' option in printing, if you have it available.
The final values are also stored in a table named results.tbl, which can be
incorporated in the Archive (using a procedure yet to be written!)
You can extract values from this table to make tables or plots for publication.

To avoid cluttering up your reduction log with other MIDAS output, you can run
the {\tt DELETE/LOGFILE} command {\it before} starting the reductions.
Note that the logfile (normally kept in your {\tt \$HOME/midwork} directory)
is a flat ascii file that can be edited to remove unwanted matter before
printing.

\subsubsection{Table file output}
The individual data, reduced as described above, are put into a simple output
table file named ``results.tbl''.
Its columns give both MJD and heliocentric Julian Date (minus 2400000); the
name of the object; the band; and the corresponding standard magnitude.
This makes the results of the reductions available for further processing with
MIDAS utilities.

For example, although the table entries are written in the same order they are
printed in the logfile (star by star, so that the necessary color data are
adjacent), they can be rearranged into chronological order, or by passband,
with the {\tt SORT/TABLE} command.
Columns can be extracted to a new table, using {\tt PROJECT/TABLE}.
You could interpolate the values of a comparison star to every time in the
table by using {\tt FIT/TABLE} or {\tt REGRESS/TABLE},
and then use {\tt COMPUTE/TABLE} to
subtract these from the values for a variable star to produce
differential photometry.

Additional columns could be transferred from the original data to the new
table, using {\tt PROJECT/TABLE} and {\tt JOIN/TABLE} or {\tt MERGE/TABLE},
so that correlations and effects not looked for by the reduction program can be
detected.
In general, because the PEPSYS package makes heavy use of table files, the
user should try to become familiar with the utilities MIDAS provides to
deal with them.


\subsection{Interpreting the output}

When you look over the printed output, pay particular attention to the
following:

1.  Rejected stars:  These are marked with an asterisk (*) in the right margin
when the residuals are listed.
An occasional reject is normal; but clumps of rejected observations are not.
If {\it every} observation of a standard star is rejected, it may be
misidentified.
If the standard values of a standard star are all rejected, there may be a
catalog (or copying) error; or the wrong star may have been observed.

Abnormally faint observations of program stars are often due to the dome being
in the way.
Be careful to check the dome slit while observing, if it is not controlled
automatically.
Another cause of abnormally faint star readings is confusion of star and sky
measurements (see item 4 below).

2.  Reading plots:  Remember the conventions used in low-resolution plotting:
the \$ symbol marks overlapping points; {\tt \^\ } and {\tt v} mark points
beyond the upper and lower edges of the plot (think of them as arrowheads).
Some plots have fitted lines indicated by a series of dashes.
The residual plots identify different stars with different symbols; these are
given in the tables of residuals.

3.  Trends in residual plots:  Similar trends in the run of residuals with time
in different bands usually indicate either instrumental instability (zero-point
drift) or varying extinction.
Instrumental drifts are usually a function of temperature and/or relative
humidity; but a bad high-voltage supply can produce irregular variations.
Extinction variations are usually larger at shorter wavelengths,
and may show short-lived dips as wisps of cirrus cross the sky.

4.  Observations with negative intensities:  These can be caused by star
observations misidentified as sky, or vice versa.
Check your data very carefully for errors.
Ordinarily, sky data should be much fainter than the star measurements, for
stars brighter than about magnitude 15.
If you have fainter program stars than this, use a CCD to obtain the
considerable benefits of simultaneous measurements of star and sky.
The sky is considerably brighter in the infrared, or during auroral displays.
Here, large fluctuations in sky brightness can occur, occasionally producing
negative star intensities after sky subtraction.
If your sky fluctuations are appreciable, compared to your faintest program
star, be sure to chop back and forth quickly between star and sky to subtract
the fluctuations accurately.


\subsection{Special problems}

\subsubsection{Missing bands}

Because a full set of colors is required to account for color terms, it is not
possible to reduce data for stars that have not been observed in all filters.
If you have some stars that were observed in only a subset of bands, you should
extract from the entire table of observations just the subset of
bands, and run the whole set of stars in this band
subset to obtain reduced values for the stars with incomplete observations.
Then make a second set of data from which the incomplete observations are
removed, and reduce these to get good values for the stars with complete
observations.
Be aware that the reduced values will differ systematically in the two subsets.

For example, some stars might be observed only in B and V in a run when most
stars were observed in the full UBV set of bands.
You would then do two separate reductions: one for all stars, using only
the B and V data, and one with only the stars having full UBV data.
Notice that the B and V values for the stars in common
will differ slightly in the two solutions.
The values from the full UBV solution should be more accurate; but
you should not intermix them with values from the BV solution.
If homogeneity is more important than accuracy, you could try adopting the B
and V values for all stars from the BV solution, and the U-B colors from the
full solution.
The B's in the B-V colors then differ from the B's in the U-B colors;
but this is basically what Johnson did in setting up the system.
To avoid the problem, make sure you observe every star in every band.

Although the reduction program assumes you will observe a contiguous subset of
the bands in any standard system, you could still force it to reduce a
{\it non}-contiguous subset by replying {\tt OTHER} when asked for the system
name.
You would then supply, in order of increasing wavelength, the bands you
actually used.
However, you would also have to make up a special set of standard-star files,
in which the indices skip any missing bands.

Note that standard-star data that are missing some passbands may still be
useful.
In general, at least one magnitude is required; except for H-Beta standards,
stars with only indices and no magnitudes are not useful.


\subsubsection{Nonlinearity}

Although it is possible to solve for nonlinearity (e.g., dead-time)
corrections in both
pulse-counting and DC photometry, {\it extreme} care should be used in doing so.
The problem is that the nonlinearity is strongly coupled to other parameters in
the solution.
A common error is to include standard-star values in the solution; this
aliases conformity errors for the 2 or 3 brightest stars into the nonlinearity
parameter.
Nonsensical values of both the nonlinearity and the transformation parameters
are the usual result, accompanied by a misleadingly ``good'' fit (small
residuals, and small standard errors on the coupled parameters).

To determine nonlinearity accurately, a neutral attenuator should be used to
observe a considerable number (say, 15 or 20) of the brightest stars.
The fainter stars serve to calibrate the attenuator; the brightest stars then
determine the nonlinearity, through comparisons between their attenuated and
unattenuated observations.
A less precise determination of nonlinearity is possible by using the
atmospheric extinction as the attenuator; unfortunately, this is not neutral,
so the coupling between nonlinearity, extinction, and bandwidth parameters can
produce systematic reduction errors.
In either case, the transformation solution should be done separately from the
extinction-and-nonlinearity fit.

It is particularly dangerous to have a single star much brighter than the rest,
as its high leverage on the nonlinearity parameter will guarantee systematic
errors.
Such a star will stand out as isolated points on the right-hand side of the
linearity plots.
If the brightest star is extremely red or blue in any color index, the errors
will affect the transformations more strongly.
Try to find 2 or 3 bright stars of intermediate color within
half a magnitude of one another, to determine the nonlinearity.

While the best linearity is obtained with a good DC system, used at small anode
currents, an overloaded DC system can be just as nonlinear as an overloaded
pulse counter.
The best policy is to know and understand your equipment thoroughly, and (if
possible) to avoid observing in the range
where nonlinearity is known to be a problem.


\subsubsection{Special systems}

If you want to either plan or reduce observations made in your own non-standard
system, both the planning and the reduction programs let you declare the system
as {\tt NONE} or {\tt OTHER}.
Use {\tt NONE} when you want to work entirely in the instrumental system, and
{\tt OTHER} when you want to use a special system of standard stars.
In either case, you will need to specify central wavelengths, bandwidths, and
the relation between the passband magnitudes and whatever color indices you
use.

You can maintain your own files of standard stars for your private system, or
for your instrumental system.
Sometimes it is convenient to maintain instrumental-system standards, to avoid
the information loss of conformity errors, which should be negligible in this
case.
Instrumental mean values can then be used to determine extinction very
accurately.


\subsubsection{Marginal nights}

Often, one finds that some nights of a run were of marginal photometric
quality.
There may have been clouds or cirrus visible; or the residuals may simply be
anomalously large.
How should these marginal nights be treated?

The safest thing is to reduce all the nights of a run together at first.
If one or two nights then have a large number of rejected observations,
try to decide whether the problem is instrumental or transparency problems.
Nights with instrumental problems should probably be rejected completely.
Nights with variable transparency may sometimes be salvaged, either by throwing
out the worst part of the night, and reducing the rest together with the good
data; or by removing a dubious night from the solution, and treating it
separately.

In treating a bad night separately, one can either force the extinction
coefficients to have particular values, or force the standard and extinction
stars to have particular values (by creating a special standard-star file).
Much of the time, the best thing to do with bad data is throw them away.


\subsubsection{Sky problems}

If the program complains that you have negative intensities, this is usually
due to misidentifying star observations as sky.
It can also be caused by marking sky observations as star data.
Check the {\tt STARSKY} column of your data tables.

Another problem with sky occurs when the program cannot find a suitable
sky observation to subtract from a star datum.
This can occur when different measuring diaphragms have been used.
Remember that you cannot correct an observation for sky unless there are sky
data taken in the same band through the same aperture.

Finally, as the sky-modelling algorithm does not try to model twilight,
it uses the ``nearest-neighbor'' method to correct star observations made
during twilight for sky.
You may find that such observations cannot be corrected for sky by this method,
if the star
was measured during twilight, but the corresponding sky measurement occurs just
after the end of evening twilight or just before the start of morning twilight.
You can prevent this problem by measuring sky both before and after each star
measured during twilight; but you should do this anyway.
Note that the planning program tells you when twilight begins and ends.


\subsubsection{CCD data}

If CCD data are to be reduced, it is {\bf essential} that they all be on the
same instrumental system.
First, all the data for each night must be reduced with a common
average flat field, for a given filter.
(It is possible to use a different flat for each night, which will introduce a
zero-point shift from one night to the next.)

Second, all the data for a given night must be comparable, to satisfy
Steinheil's principle.
One can have problems with some image-extraction routines that use PSF fitting.
Because of seeing variations during to night --- and especially because of the
dependence of seeing on airmass --- there may be systematic errors introduced
by using different point-spread functions on different frames.
If a very detailed PSF model is available, so that the whole energy in a star
image is well extracted, with very small residuals, one may expect PSF-fitting
to work adequately.
However, one must be sure that the extracted magnitudes refer to the total
{\it energy} in the image, and are not just scaled to the peak.

If you use a PSF-fitting routine that leaves obvious ``blemishes'' when the
fitted profile is subtracted from the original frame, it is likely that there
will be systematic errors that depend on seeing.
In turn, this means systematic errors that depend on airmass, which will spoil
the determination of extinction coefficients.

In general, the safest approach with CCD data is to simulate ``aperture''
photometry, as it is often called --- just integrate the total signal in a box
(round or square) of fixed size centered accurately on each star.
This may give larger random errors than PSF-fitting, but smaller systematic
errors.
This balance between accuracy and precision is a common dilemma in stellar
photometry.


\subsubsection{Problems with star names}

If the star names used in data files are not {\it exactly} the same character
strings as those used in star files, the reduction program will try to make
cross-identifications.
These are based on heuristic rules that are generally observed in naming stars
--- for example, many stars have a catalog abbreviation followed by a number,
or a Greek letter followed by a constellation abbreviation.
The program recognizes some common catalog abbreviations (HR, BS, BD, CD, CPD,
HD, NGC), but will guess that two or three letters followed by a number
represents such a name.
It also can cope with common suffixes like A, B, and AB for multiple stars.

The program applies these rules in attempting to parse name strings
containing multiple designations that are not separated by
the `` {\tt =} '' string (surrounded by blanks).
If you sometimes write a name with an embedded blank and sometimes without
(e.g., ``{\tt HR 123}'' vs. ``{\tt HR123}''), it may be able to identify the
two as equivalent, or it may ask for help.
It is difficult to write a simple set of foolproof rules for recognizing star
names; for example, the program cannot simply squeeze out embedded blanks, as
it would then confuse BD +1 2345 with BD +12 345.

Occasionally a typing error can confuse the program, and it will print
{\tt Cannot parse:} followed by a name string.
For example, spelling errors in constellation abbreviations, Greek letters,
or the letter O for a digit 0 can cause such a message.
In these cases, the program asks for help, and asks for a replacement name
string (which can contain embedded blanks).
Afterward, the replacement string will be used instead of the one that caused
problems.

This situation can be avoided by making sure that all star names are spelled
consistently, and that names in star files agree with names in data files.





\section{Installation}
\label{install}

The first user of PEPSYS at a new site must make sure that the necessary table
files have been installed.
Also, it may be useful to change some {\tt PARAMETER} statements to suit local
circumstances.
Most users will never have to worry about these things; but here they are, just
in case you need them.

\subsection{Table files}
\label{installtbl}
MIDAS table files are binary files, to save the considerable overhead of
converting between ascii and the machine's internal binary representation of
numbers.
The price to be paid for fast response while reading a table is that
tables are not transportable between different machines.
This means that the standard tables have to be created on each new machine.

Fortunately, MIDAS also understands FITS formats, which are portable.
Therefore, the tables of standard stars (and some other sample tables) are
available as portable ascii files, and a script is provided to create the
tables.
The FITS files are (on a UNIX system) in the directory
{\tt \$MIDASHOME/calib/raw/pepsys},
which contains several files with the suffix {\tt .mt}.
The directory {\tt \$MIDASHOME/calib/install}
contains scripts for installing such files;
the script you want is {\tt pepsys.prg}.
The local MIDAS guru should run this script to create the table files in
{\tt \$MIDASHOME/calib/data/pepsys}.
After running the script, please check the newly-created table files to make
sure users can read them but not re-write or remove them.

If the ``calib'' directories and files are not part of your system, they can be
obtained by ftp, or (if necessary) on magnetic tape.

\subsection{Maximum limits}
The planning and reduction programs are distributed with array sizes set large
enough to let most observers work without difficulty.
However, you might need to reduce a very large data set, or to deal with a very
large number of passbands (9 is the nominal limit).
The programs are all written in FORTRAN, and the dimensions of these arrays are
set by {\tt PARAMETER} statements in ``include'' files.
All these short files are in the {\tt \$MIDASHOME/\$MIDVERS/contrib/pepsys/incl}
subdirectory.

The programs check for array overflows as they read data.
If you encounter a limit, the program will tell you which {\tt PARAMETER}
statement to change.
Ask your MIDAS guru to edit the  {\tt PARAMETER} statement and recompile the
programs.
To make sure everything is internally consistent, compile the subroutines
{\it before} compiling the main programs: move to the {\tt pepsys/libsrc}
directory and make the subroutines, then go to the {\tt pepsys/src} directory
to make the main programs.

The main limitation on what is practical is machine memory.  If your machine
does not have enough memory to store all the necessary arrays without swapping
pages back and forth to disk, the reduction program may run very slowly.
(Bear in mind that response on multi-user machines also depends on what else is
being run by other users at the same time.)
The main problem is holding the matrix of normal equations during the
iterations; on the other hand,
star catalogs can be made quite large without paying much of a
price in terms of performance.
(However, if the number of stars exceeds 2048, the binary-search subroutine
{\tt NBIN} will need to be modified.)

You can get a rough idea of where you will run into problems by noting that the
matrix of normal equations is made of double-precision numbers, and is a square
matrix (well... triangular, actually) with as many elements each way
as there are parameters to be evaluated.
For example, to solve for 120 extinction and standard stars in 4 colors takes
480 parameters; there will be some extinction and transformation coefficients,
too, so this problem is about 500 x 500, using a quarter of a million matrix
elements.
Each element is 8 bytes long, on most systems; so we need roughly 2 MB of
memory --- no problem these days.
Of course, only half the matrix is stored, as it is symmetric.
On the other hand, you need space for the right-hand-side vector, and the
program as well, not to mention other arrays that are needed.
So this rough calculation is good enough for astrophysical accuracy.

On the other hand, if you wanted to do 1000 stars simultaneously in 4 colors,
you'd need some 16 million elements of 8 bytes each, or 128 MB of storage.
That's probably too big for most machines to handle gracefully.
Or, if you wanted to reduce 100 channels of spectrophotometry simultaneously,
even 20 stars would give you 2000 parameters (plus 100 extinction coefficients
per night!); that would be over 32 MB, and heading for trouble.
In this latter case, it would probably make sense to reduce smaller subsets of
bands together, to allow more stars and nights to be combined.
In general, you should try to reduce at least 4 or 5 nights' data together, to
improve the extinction determinations.

If your machine proves to be a bottleneck, you might want to maintain two
variants of the reduction program: one small enough to run quickly on ordinary
data sets, and a monster version for the occasional monster problem.
In this case, the small, fast version would be kept in the regular MIDAS
directory tree, and the gigantic version could be owned by the
rare user who needs it.

If you find it necessary to increase some of the array parameters, please let
the Image Processing Group at ESO know what changes you make.
This will allow revised parameters to be distributed in future releases.



\section{A brief history of PEPSYS}

This package is descended from a program
originally written for the IBM 704 in the late 1950's~\cite{YI67}.
It used all constant stars as extinction stars, combined multiple nights
with separate extinction coefficients, and
employed King's correct representation of color terms in the extinction.
That program was written in assembly language after an initial try in FORTRAN II
would not fit into the 8192 words of the machine's memory.
Data were packed, two to a word, on magnetic drums.
Computers at that time were so small that the program needed every cell of the
machine; even so, it had to be split into overlays.

When the 7090 came along,
the program was stored on a loadable tape that displaced the FMS
operating system from core memory, and pulled the system back in when it
finished.
This special system tape was labelled PEPSYS (for Photo-Electric Photometry
SYStem).
This original PEPSYS died when the IBM 709x systems became obsolete.

A second attempt was made on a Cyber 175 many years later.
The planning program was added in the 1980's.
This was ported to a VAX, and later to other machines (AT\&T 3B2 and Sun 4).
These ports shook a lot of portability bugs out of the common subroutines.
The Manfroid -- Heck program~\cite{MH84}, and a program used at that time by
Peter Stetson were investigated, but found to be too specialized to particular
photometric systems and data formats to be generally useful,
so a proposal was made to NSF to develop something along the lines
of the present system.
The proposal was rejected on the grounds that (a) it was too expensive, and (b)
such a thing is impossible anyway.

The heart of the present reduction program --- the fast matrix inversion using
partitioning -- was developed during a simulation study in the late 1980's.
However, it was used only on phony data, and lacked the extensive user
interface of the present version.
Some user interface, especially material now embedded in the
{\tt MAKE/PHOTOMETER}
command, was developed on an Intel 80386 system running SCO Xenix System V.

The current system was sponsored by ESO; I thank Chris Sterken and Preben
Grosb\o l for setting up my visit to Garching.
This version was developed within ESO's Image Processing Group.
Although an image-processing system offers the photometrist
as many inconveniences to work
around as useful tools, it has been possible to hide most of the problems from
the user.


\section{Acknowledgements}

Many people have contributed help and advice in developing the current version
of PEPSYS.
I have learned a lot about photometry over the years from A.~W.~J.~Cousins,
whose papers I recommend to every observer.
John Menzies communicated the Cousins E-region standards in both UBV and uvby
systems.
Erik Olsen provided the latest uvby -- H-Beta standards, and advice on their
use.
Harri Lindgren and
Petr Harmanec provided much useful discussion of reduction methods.
Stan Stefl was the first user of PEPSYS at ESO, and uncovered a number of
bugs as well as suggesting helpful features.
Chris Sterken tested the reduction program, helped find bugs, and
made useful suggestions.
Josef Hron carefully read the documentation and suggested numerous
improvements.
Fionn Murtagh pointed out a useful reference~\cite{FB92}.
The members of the IPG all helped me with numerous MIDAS problems.
Above all, I must thank Chris Sterken and Jean Manfroid for thinking hard about
photometric reductions, and publishing a stimulating series of papers on the
subject, as well as their recent textbook~\cite{SMbook}.


\section{Summary of PEPSYS commands}

\begin{itemize}

\item MAKE/HORFORM: Makes a blank form to fill in with horizon data

\item MAKE/STARTABLE: Helps you make a *.tbl file for program stars

\item MAKE/PHOTOMETER: Makes or displays an instrument-description file.

\item MAKE/PLAN: Makes an observing schedule for you

\item CONVERT/PHOT: Converts some ascii data files to MIDAS tables.

\item REDUCE/PHOT: Reduces the data.

\end{itemize}



%REFERENCES
%References here:

\begin{thebibliography}{99}
\bibitem{BN89}  Beckert, D.~C., and Newberry, M.~V. : 1989, {\it Pub.~A.~S.~P.},
				{\bf 101}, 849.
\bibitem{B90}  Bessell, M.~S. : 1990, {\it Pub.~A.~S.~P.},
				{\bf 102}, 1181.
\bibitem{DS}  Dennis, J.~E., and Schnabel, R.~B. : 1983,
				{\it Numerical Methods for Unconstrained
				Optimization and Nonlinear Equations},
				Prentice-Hall, pp.~227-228.
\bibitem{Evans}  Evans, R.~D. : 1955, {\it The Atomic Nucleus},
				McGraw-Hill, New York, pp.~785-790.
\bibitem{FB92}  Feigelson, E.~D., and Babu, G.~J. : 1992,
				{\it Ap.~J.}, {\bf 397}, 55.
\bibitem{G1}  Garstang, R.~H. : 1986, {\it Pub.~A.~S.~P.},
				{\bf 98}, 364.
\bibitem{G2}  Garstang, R.~H. : 1989, {\it Pub.~A.~S.~P.},
				{\bf 101}, 306.
\bibitem{G3}  Garstang, R.~H. : 1991, {\it Pub.~A.~S.~P.},
				{\bf 103}, 1109.
\bibitem{GMW}  Gill, P.~E., Murray, W., and Wright, M.~H. : 1981,
				{\it Practical Optimization},
				Academic Press, New York, p.~104.
\bibitem{H62}  Hardie, R.~H., : 1962, in
				{\it Astronomical Techniques}, W.~A.~Hiltner,
				ed., Univ. of Chicago Press, Chicago, p.~178.
\bibitem{H80}  Hall, P. : 1980, {\it Rates of convergence in the central
				limit theorem}, Pitman, Boston.
\bibitem{H86}  Hampel, F.~R., Ronchetti, E.~M., Rouseeuw, P.~J., and
				Staehl, W.~A. : 1986, {\it Robust Statistics},
				Wiley, New York.
\bibitem{HMT}  Hoaglin, D.~C., Mosteller, F., and Tukey, J.~W. : 1983,
				{\it Understanding Robust and Exploratory
				Data Analysis}, Wiley, New York, Chap.~5.
\bibitem{KMN}  Kahaner, D., Moler, C., and Nash, S. : 1989,
				{\it Numerical Methods and Software},
				Prentice Hall, Englewood Cliffs, NJ, p.~368.
\bibitem{MH84}  Manfroid, J., and Heck, A. : 1984, {\it Astron.~Astrophys.},
                                 {\bf 132}, 110.
\bibitem{MS92}  Manfroid, J., and Sterken, C. : 1992, {\it Astron.~Astrophys.},
                                 {\bf 258}, 600.
\bibitem{ML85}  Menzies, J.~W., and Laing, J.~D. : 1985, {M.~N.},
                                 {\bf 217}, 563.
\bibitem{MT77}  Mosteller, F., and Tukey, J.~W. : 1977,
				{\it Data Analysis and Regression},
				Addison-Wesley, Reading.
\bibitem{P82}  Popper, D.~M. : 1982, {\it Pub.~A.~S.~P.}
                                 {\bf 94}, 204.
\bibitem{NR}  Press, W.~H., Flannery, B.~P., Teukolsky, S.~A., and
				Vetterling, W.~T. : 1986, {\it Numerical
				Recipes}, Cambridge University Press.
\bibitem{ASC91}  Schwarzenberg-Czerny, A. : 1991, {\it Astron.~Astrophys.},
                                 {\bf 252}, 425.
\bibitem{SMbook}  Sterken, C., and Manfroid, J. : 1992,
				{\it Astronomical Photometry: a Guide},
				Kluwer, Dordrecht.
\bibitem{T60}  Tukey, J.~W. : 1960, in {\it Contributions to Probability and
				Statistics}, I.~Olkin, ed.,
				Stanford Univ.~Press, p.~448.
\bibitem{vdH}  van de Hulst, H.~C. : 1980,
				{\it Multiple Light Scattering},
				Academic, New York.
\bibitem{W87}  Walker, M.~F. : 1987, in {\it Identification, Optimization, and
				Protection of Optical Telescope Sites},
				R.~L.~Millis et al., eds.,
				Lowell Observatory, Flagstaff, p.~128.
\bibitem{YI67}  Young, A.~T., and Irvine, W.~M. : 1967,
			     {\it Astron. J.}, {\bf 72}, 945.
\bibitem{Y74}  Young, A.~T. : 1974, {\it Methods of Experimental Physics},
				   {\bf 12}, Part A, Astrophysics: Optical and
				   Infrared,
				   ed. Carleton, N.,
				   Academic, New York, Chap.~3.
\bibitem{Y91}  Young, A.~T., et al. : 1991, {\t Pub.~A.S.P.}, {\bf 103}, 221.
\bibitem{Y92}  Young, A.~T. : 1992, {\it Astron.~Astrophys.},
                                 {\bf 257}, 366.
\bibitem{Y92b}  Young, A.~T. : 1992, ``High-Precision Photometry'', in
			{\it  Automated Telescopes for Photometry and Imaging},
			(A. S. P. Conference Series, vol. 28)
			S.~J.~Adelman, R.~J.~Dukes, and C.~J.~Adelman, eds.
			(Astronomical Soc.~of the Pacific, San Francisco, 1992)
			pp.~73-89.
\end{thebibliography}



%\cleardoublepage
\end{document}				% remove ^% to print
