*  D A O _ S T A T   - REDUCTION PROGRAM *

					by Friedel Loinger

(I) M E T H O D

   The dao_stat program gets as input a set of images, containing
frame-dependent magnitudes. The purpose of the program is to decide 
which stars are standards and to compute light-curves of object stars.
   Absolute scaling can be done, if the user has a reference image,
containing photometric standard magnitudes, obtained by the 
the STDP reduction program.

   We start by checking each potential comparison star for either genuine
variability or for bad data points (due to cosmic rays, bad pixels, etc.).

   At the beginning (iteration 0), we take the first frame as
a reference image, and choose one of the stars as a reference
star. This star should appear in all the other frames.
For example,the observed magnitudes of the comparison star in
the first frame are: 
  13.0,14.0,15.0,16.0,17.0 and the second star is the reference star.
In the second frame the observed magnitudes of the same stars are
  13.5,14.5,15.5,16.5,17.5, then
the program decreases these magnitudes by 0.5, in order
to obtain initial values for "real" magnitudes. The value of linear shift of
each frame depends implicitly on conditions like clouds, airmass, etc.
The next step is to assign to each comparison star the mean, over all frames, 
of its magnitude relative to the chosen reference star.

   Then for each comparison star we calculate an improved mean by comparing 
with the mean magnitudes (as obtained in the previous stage) of all other 
comparison stars in the same frame. For each comparison star in turn,
except of the reference star, we construct a light curve and
compute a new mean value over the frames and its standard deviation. 
This allows us to identify discrepant points or real variability, and any
bad points are removed from the data set.
We repeat the procedure, until all the mean magnitudes converge.
Now we finished the first step, where the magnitude of the reference star was
kept fixed. Only after the process converges, also a light curve 
for the reference star is computed, and its mean value over the
frames is taken as a new reference magnitude. Afterwards the light curves
of the other comparison stars are computed again. We continue the procedure
until all mean magnitudes converge. 
The process usually converges in a few iterations. If it does'nt,
then you have chosen a bad set of comparison stars. The reference star
might be variable.

  We now have a set of comparison stars with their computed mean magnitudes
and standard deviation. The user can now remove 'bad' stars, 
according to their standard deviations, and compute new mean magnitudes,
based on the new set of comparison stars. We continue until we get a 
high quality set of comparison stars from which one can define an zero point 
for the magnitude scale of each image, including error estimations.
The larger the number of frames and number of comparison stars, the better
is the accuracy achieved. It is also important to have a large enough
range of magnitudes, and the objects star[s] being in that range.
   If the user has calibration data from the STDP reduction program, 
then dao_stat can now calculate real apparent magnitudes for the standard

   We then use the magnitudes to calculate light curve of either 
one object star or all variable stars. 
The normal procedure is to use a weighted linear regression of
instrumental versus apparent magnitudes of the local standards and to use
this relation to obtain magnitudes of the object star[s].
In some cases, where the number of local standards is very small (2 or 3), 
or when the magnitude of the object is outside of the standard stars scale,
it is preferable to force a fit to a slope 1.0 relationship. 

(II) I N P U T    F I L E S

(1) filename, containing a list of frame-names

recommended filename: xxxF.yyy (no obligation)
where :
xxx - name of object 
  F - name of filter (1 character)
yyy - extension 

example: s3c279b.dat

Important: Write the filenames without the 'nst' extension into the list-file.

Enter the *.nst - files in following order:

reference image (for arranging coordinates only)
field 1
field 2
field n
-----------   ( this line must contain at least 3 char ---)
reference image for standard magnitude (after STDP-reduction)

(a) the first image is only used for arranging the coordinates.
    If you want it to be considered as one of the other fields,
    then you should enter its name twice.
(b) The program will also run without the last 2 lines.
(c) field 1,...,field n, should be ordered by date
(d) maximal number of images:      8000          



(2)  names and contents of *.nst-files:

You should prepare the frames, so that they will match each other, i.e. each 
fixed star should have the same coordinate (+- eps) in all the frames.
In order to shift the images, you can use the procedure "trans" in
the vista-reduction program.

(a) If the observations were done in the WISE-Observatory.
    and the *.nst - files were created by DAOPHOT in our local VISTA-version,
    then they will contain 2 additional lines at the beginning. 
    The total number of header-lines is then 5, instead of 3 of the original

    The first one is a comment-line, containing the object-name,filtername 

    If you are interested in an accurate Julian-day, then you have 2 options:

    I) Insert the Julian-day at the second column in the input file of 
       framenames. The program will read it in format F11.5.

    II) The Julian day can be computed from a mean value of the ut-time,
        and the observation date. 

    If you do'nt want an accurate value, then 0.5 is added to the julian-day..

Example: s3c27902jul92b.nst

'3C 279          R  250s UT18.41 02/07/92 AM1.653 02:40W       '       #line 1
								       #line 2
 NL  NX  NY  LOWBAD HIGHBAD  THRESH     AIR   ITIME     HJD     AP1    #line 3
  1 320 512   115.3 12000.0    24.0    0.00    0.00    0.00     5.0    #line 4
								       #line 5
     1    33.18   481.19   14.421    0.004  139.289       3.     1.01    0.005
     2   299.31    98.65   15.380    0.007  138.769       3.     0.99   -0.059
     3   132.41   274.13   15.401    0.006  138.827       3.     0.85    0.000
     4   171.88   217.86   15.946    0.053  138.509       3.     1.02    0.117
     5   266.54   482.23   16.102    0.012  137.873       3.     1.08    0.055
     6   119.05   411.54   16.617    0.014  138.916       3.     0.83   -0.075
     7   192.72   281.58   17.630    0.003  137.825       5.     0.94    0.115
     8   208.85   337.49   17.640    0.003  138.822       4.     1.06    0.402

(b) If the frames come from another observatory, and do not contain our 'wise'
    comment-line, then give following names for the nst-files:


    xxx - any name
    DDMMMYY - date
    F   - name of filter

    If the user wants to get accurate Julian-Date, then add in the list-
    of filenames the UT-time after the name of the nst-file, for example:

    s3c27902jul92b 19:10

    or choose solution (2aI)

The dao_stat is able to handle both input-types (either containing 5 or 3
header-lines), since the program always searchs for the line 
starting with ' NL', and then skips the next 2 lines, 
in order to reach the input of the the magnitudes.


Type: dao_stat

(1) Do you want to weigh frames according to scatter (y/n)?
Usually answer 'y'
When averaging over frames, then the STD-error of the frame will be taken
as its weight. A lower limit of 1.0e-02 will be taken for the weight.

(2) Enter ID for unknown star ( 0 - all stars considered unknown ) 

If the user is interested in a light-curve for one particalur star, then
he should enter its ID-number here. The ID-number appears in the first 
column of the *.nst - file of the reference image.

If the user answers 0, then all variabled stars will be taken as unknown.

(3) Enter filename, containing frame-names
   (first frame is used only as coordinate reference for all others)

Enter here an input filename, as explained in (II)(1).

(4) Do you have already arranged stars according to ID numbers (y/n)?
If your answer is 'n', then the program will take the *.nst - files,
and create *.arr output-files. The arrange will be done, so that stars with
the same (+- eps) coordinates, will have the same ID-number in all the fields.
As reference the first *.nst - file will be taken. A star, which is missing
in a field will, be added to the file with a value 99.000 for its magnitude
and 9.999 for its error.

The *.arr will contain the original header of the *.nst - files, and another
header-line will be added. Thus the *.arr-files in case (2a) will contain
6 header-lines, and 4 in case (2b). After the title, there will
be 5 input columns for the stars, which will be read in free format by 
The program can handle up to 8000 stars.

Example: s3c27902jul92b.arr

'3C 279          R  250s UT18.41 02/07/92 AM1.653 02:40W       ' #line 1
							         #line 2	
 NL  NX  NY  LOWBAD HIGHBAD  THRESH     AIR   ITIME     HJD      #line 3   
  1 320 512   115.3 12000.0    24.0    0.00    0.00    0.00      #line 4   
								 #line 5
    ID     X        Y     DAOPHOT-MAG.   ERR                     #line 6
     1    33.18   481.19    14.421     0.004  
     2   299.31    98.65    15.380     0.007  
     3   132.41   274.13    15.401     0.006  
     4   171.88   217.86    15.946     0.053  
     5   266.54   482.23    16.102     0.012  
     6   119.05   411.54    16.617     0.014  
     7   192.72   281.58    17.630     0.003  
     8   208.85   337.49    17.640     0.003 
     9   154.11   150.45    99.000     9.999

If your answer is 'n', then the program assumes that the *.arr - files already

If the user did not create the *.nst-files using DAOPHOT, and has only
a table of magnitudes of the stars, then he should prepare the *.arr - files,
by adding header-lines. The dao_stat program will always search for a line
starting with ' NL', and then skip over the next 3 lines, in order to
reach the input containing the magnitudes. Therefore it is enough to
create a ' NL' line, and then add 3 another dummy header-lines
(might be empty). If you do'nt have coordinates (X,Y) for the stars,
then add numbers in columns 2,3 so that for each star ID,X,Y will be
identical in all the frames. 
The 4.-th column should contain the observed magnitudes: 25 - 2.5*log(counts)
If you do'nt have an error-estimation for the stars, then the 5.-th column
should contain the same error for all stars.

Example: s3c27902jul92b.arr

 NL                                              #line 1
                                                 #line 2
                                                 #line 3
    ID     X        Y        MAG         ERR     #line 4
     1    10       10       14.123        1
     2    20       20       14.551        1
     3    30       30       15.045        1
     4    40       40       15.915        1
     5    50       50       16.003        1
     6    60       60       17.671        1

(5) Enter min. no. of stars in a frame required for linear regression
   ( 0 for no regression, ie. mean ) 
If you answer 0, then mean is computed instead of linear regression,
when comparing a star with the others stars in the same frame.
If you are interested in linear regression, then you should enter
the min. of stars required for linear regression. Enter at least 4, otherwise
you might loose accuracy. 
If there is a frame, which contains less than the chosen min., then the mean 
will be taken instead for this frame.

(6) Enter ID for reference star
The reference star should be a non-variable star, which appears 
in all the fields. It should not be too bright, to avoid saturation problems.
The standard magnitudes of all the others stars will be computed relative
to this reference star.

The program will tell you now, how many frames are in the input - list
( the first image, which is only used for arrange the coordinates is not
counted here).

(7) Enter minimum presence for standards ( 0 for default = 0.8 of all frames ) 
    A candidate for standard star, should appear in most of the fields.

(8) If you have chosen all variable stars to be unknown 
    ( ID = 0 in paragraph 2), then the next question will be:

    Enter minimum presence for unknowns ( 0 for default = 0.5 of all frames )
(9) Enter ID for stars not to be included except of unknowns
    (-1 for next step)
    If you know ID's of variable stars, insert them here, one star
    each time, ending the list with -1 .
Iterations are now done, in order to compute the "real" magnitudes of the
candidates of the standards stars (the mean magnitudes as explained in (I)).
At the beginning, the real mag. of the reference star is kept fixed, and only
after the process converges, also iterations are done in order to compute the
real mag. of the reference star.

The number of standards stars is limited to 900.

The program types out the iterations, and shows following RESULTS
after the process converges:

   ID   X     Y      presence   AV +- STD    AV_D +-STD_D  CHI2/N


ID       - identification no. of the star
X,Y      - coordinates     "  "   "
presence - relative occurence of the star in the frames. If it appears
           in all of them, this value will be 1.
AV       - real magnitude of the star, averaged over all those frames,
           where the star appears. This will be weighted average, if you
           answered 'y' to qestion (1), where the weight is the error
           of the frame.
STD      - standard deviation the star's magnitude.
AV_D	-  weighted average, where the weight takes into account the
           error for the frame and the dao_phot - error of the star 
           (column 5 of the *.arr - file)
STD_D	-  standard deviation in the case above.
CHI2/N	-  weighted CHI2, where the weight includes the error of the
           frame and dao_phot - error. N is the number of occurences
           of the star.

The light curves of the standards are displayed together with errorbars.
These errorbars are the std-errors of the frames.
Usually the first display is bad (not all the data is displayed),
due to a bug in the display software.
If this happens, choose DS in the menu, to display again.
The options of the menu are:

DS - DISPLAY, B - BLOWUP , XYB -  X,Y - BOUNDS, PSL - psland,
PSP - psport, S - STOP 

DS -   display again with original bounds
       If you choose this option, move the mouse the display window.
       A crosshair cursor will appear. Now choose 2 opposite corners of
       a rectangle, by pressing any character (you can also press the mouse).
       The data is then displayed again with the new bounds.
       You will be asked for XMIN,XMAX,YMIN,YMAX.
       This options enables the user either to decrease or increase
       the bounds.
PSL -  landscape hardcopy
PSP -  portrait hardcopy
S   -  stop the display menu

(10) Do you want to remove bad stars (y/n)? 
According to the results in the previous paragraph, you can decide if there
are variable stars. If you answer 'y' then you will be asked following
questions :

   (10.1) Enter max. of STD for standards  ( 0 for default = 1 )
   If you answer for example 0.03, then all stars with STD > 0.03, will
   be removed from the list of standards stars.

   (10.2)  Enter max. of CHI2/N for standars ( 0 for default = 1,000,000 )
   If you answer for example 2, then all stars with CHI2/N > 2 will
   be removed from the list of standards stars.

   The program returns now to paragraph (7), and starts the iterations
   again for the new list of standards stars. You can continue
   to remove bad stars, until you get a good set of standard stars.
   Then answer 'n' to question (10).


If your list of frame-names contains a reference image of 
standard magnitudes (achieved by STDP-reduction), then the
program continues as follows:

For each standard star from the previous step, if it appears also in the
STDP-file, a comparison of their magnitudes will be written out, with
following header:

   ID    MAG-stdp    MAG-dao_stat  difference   daophot-err

Then you are asked:

   (11.1) Do you want weighted average of standards (y/n)? 
          If you answer 'y', then the weight will be based on the
          daophot error of the standard magnitude.

   (11.2) Enter ID no. not to be included in average (-1 for next step)

An average value of the magnitude differences is now computed and 
added to the dao_stat - magnitudes, and we finally obtain the
apparent magnitudes of the standard stars. 



Do you want accurate Julian-Date in output-files (y/n)?

If your answer is 'y' then either the Julian day is taken from
the second column of the framenames input-file, or the UT-time from our 
wise-header will be taken into account to compute the Julian-day.
otherwise 0 is taken as UT-time.
(see 2a)

The last step in dao_stat is obtaining the light curve of the
unknown object[s].

If you have chosen a positive ID-number for the unknown (see paragraph 2), 
then its light curve together with error-bars will be displayed twice.
The first light curve is computed by taking weighted mean in each frame,
and the second light curve by taking weighted regression in each frame
The progam allows you to inspect an individual frame in both cases.

dao_stat will prepare following OUTPUT files:

filename.var (only in case of ID=0 in paragraph 2)

where filename is the original input filename without its extension.

list of common stars in the *.nst-files, and also those which are
missing only in one frame.

contains all results before and after STDP-correction

contains for each standard star the computed real magnitude 
and its observed magnitudes and daophot-errors.

output-file in xvgr-format:
comment-lines will start with #
The output columns are the result of the unknown(s).


See also the article:
Optical Monitoring of Luminous AGN: I. Radio-loud Quasars, MNRAS 1994,
published by Prof. H. Netzer, Friedel Loinger, Ana Heller, Tal Alexander.