; $Id: waverage.pro,v 1.2 2005/02/03 18:04:12 azuah Exp $
;###############################################################################
;
; NAME:
;  waverage
;
; PURPOSE:
;   Determines the weigted average (weigts are the inverse
;   square of the errors) and associated errors of two input values. 
;   The input can be scaler or vector. However, both inputs should 
;   be of the same dimensions and size.
;   *** calculations are done in double precision
;
; CATEGORY:
;   DAVE Framework
;
; AUTHOR:
;   Richard Tumanjong Azuah
;   NIST Center for Neutron Research
;   100 Bureau Drive, Gaithersburg, MD 20899
;   United States
;   azuah@nist.gov; (301) 9755604
;   May, 2002
;
; LICENSE:
;  The software in this file is written by an employee of 
;  National Institute of Standards and Technology 
;  as part of the DAVE software project.
;
;  The DAVE software package is not subject to copyright protection
;  and is in the public domain. It should be considered as an
;  experimental neutron scattering data reduction, visualization, and
;  analysis system. As such, the authors assume no responsibility
;  whatsoever for its use, and make no guarantees, expressed or
;  implied, about its quality, reliability, or any other
;  characteristic. The use of certain trade names or commercial
;  products does not imply any endorsement of a particular product,
;  nor does it imply that the named product is necessarily the best
;  product for the stated purpose. We would appreciate acknowledgment
;  if the DAVE software is used of if the code in this file is
;  included in another product.
;
;###############################################################################



;===============================================================================
;  vartype
;
; PURPOSE:
;  Deterimines the type name of an IDL variable.
;
; PARAMETERS:
;   x1 [in] - The first value to be averaged. Can be a scaler or
;             vector
;   
;   sig1 [in] - The error in x1. Same dimensions and size as x1
;   
;   x2 [in] - Second value. Should have same dimensions and size as x1
;   
;   sig2 [in] - The error in x2.
;   
;   result [out] - The weighted average of x1 and x2 taking into
;                  account the error in x1 (sig1) and x2 (sig2) 
;   
;   sig [out] - The error in result.
;   
; KEYWORDS:
;   err [out] - Set this keyword to obtain an error number on exit
;               from this routine. err will have one of the following
;               values on exit:
;               0  Successful completion of average - no errors
;               1  Invalid number of (or undefined) parameters specified.
;               2  The size and dimensions of the inputs do not match.
;               3  Either sig1 or sig2 contains values that are <= 0.
;
;   emsg [out] - Set this keyword to obtain a descriptive text of the
;                error encountered within this subprogram. The
;                posisble values are as described above for errnos.
;
pro waverage, x1,sig1,x2,sig2,result,sig,err=err,emsg=emsg

compile_opt idl2

; check number of non-keyword parameters is 6
if (n_params() ne 6) then begin
    err = 1
    emsg = 'Invalid number of procedure parameters specified'
    return
endif

; check input parameters are all defined
arg1 = n_elements(x1)
arg2 = n_elements(sig1)
arg3 = n_elements(x2)
arg4 = n_elements(sig2)
if ((arg1 eq 0) or (arg2 eq 0) or (arg3 eq 0) or (arg4 eq 0)) then begin
    err = 1
    emsg = 'At least one input parameter is undefined/uninitialised'
    return
endif

; check input parameters have same dimensions/sizes
if ((arg1 ne arg2) or (arg2 ne arg3) or (arg3 ne arg4)) then begin
    err = 2
    emsg = 'The input parameters do not have the same dimension and/or size'
    return
endif

; check sigma's have no zero or -ve values
void = where(sig1 le 0,count1)
void = where(sig2 le 0,count2)
if ((count1 gt 0) or (count2 gt 0)) then begin
    err = 3
    emsg = 'Input error parameters contain values that are zero or -ve!'
    return
endif

sig1Sq = double(sig1)^2
sig2Sq = double(sig2)^2
invSigSq = 1.0D/sig1Sq + 1.0D/sig2Sq
result = (double(x1)/sig1Sq + double(x2)/sig2Sq)/invSigSq
sig = sqrt(1/invSigSq)
err = 0
emsg = ''

end 
