*! 1.0.0 NJC 9 April 2004
program circovmean, sort rclass
        version 8.0
        syntax varname(numeric) [if] [in]  [, noLIST * ]

	quietly {
		marksample touse
		count if `touse'
		if r(N) == 0 error 2000
		else local N = r(N)

		tempvar xsum ysum vmean tag
		tempname XSUM YSUM vecmean veclng vecstr

		bysort `touse' (`varlist') : ///
			gen `xsum' = sum(sin((`varlist' * _pi) / 180))
		by `touse' : replace `xsum' = `xsum'[_N]
		by `touse' : gen `ysum' = sum(cos((`varlist' * _pi) / 180))
		by `touse' : replace `ysum' = `ysum'[_N]

		scalar `XSUM' = `xsum'[_N]
		scalar `YSUM' = `ysum'[_N]

		// Stata atan routine takes a single argument
		// and gives the wrong answer in three out of four quadrants
		Atan2 `XSUM' `YSUM'
		scalar `vecmean' = `r(angle)'
		local vecshow : di %2.1f `r(angle)'
		scalar `veclng' = sqrt((`XSUM')^2 + (`YSUM')^2)
		scalar `vecstr' = `veclng' / `N'

		replace `xsum' = `xsum' - sin((`varlist' * _pi) / 180)
		replace `ysum' = `ysum' - cos((`varlist' * _pi) / 180)
		egen `vmean' = atan2(`xsum' `ysum') if `touse'
		bysort `touse' `varlist' : gen byte `tag' = _n == 1
		char `vmean'[varname] "mean without"
		label var `vmean' "vector mean omitting value"
		format `vmean' %2.1f
	}

	circscatter `vmean' `varlist' if `tag' & `touse' ///
		,  yli(`=`vecmean'', lp(dash))           ///
		subtitle("vector mean for all `vecshow'{c 176}", place(w)) ///
		`options'

	if "`list'" == "" {
		list `varlist' `vmean' if `tag' & `touse' ///
		, subvarname abb(12) noobs
	}

	return scalar N = `N'
	return scalar vecmean = `vecmean'
	return scalar veclng = `veclng'
	return scalar vecstr = `vecstr'
end

program Atan2, rclass
* 1.3.0 NJC 18 Dec 2003
* 1.2.0 NJC 14 July 1998
version 8.0
	tempname at

	local sign1 = sign(`1')
	local sign2 = sign(`2')

	if (`sign1' == 1 & `sign2' == 1) | ((`sign1' == 0) & `sign2' == 1) {
		scalar `at' = atan(`1'/`2')
	}
	else if `sign1' == 1 & `sign2' == 0 {
		scalar `at' = _pi / 2
	}
	else if `sign1' == -1 & `sign2' == 0 {
		scalar `at' = 3 * _pi / 2
	}
	else if `sign2' ==  -1 {
		scalar `at' = _pi + atan(`1'/`2')
	}
	else if `sign1' == -1 & `sign2' == 1 {
		scalar `at' = 2 * _pi + atan(`1'/`2')
	}
  	else if `sign1' == 0 & `sign2' == 0 {
		scalar `at' = .
	}

	return scalar angle = (180 / _pi) * `at'
end