*! NJC 2.0.2 10 May 2004 
* NJC 2.0.1 3 May 2004 
* NJC 2.0.0 28 April 2004 
* NJC 1.0.0 1 August 2001 
* von Mises F(theta | kappa, mu)                                      
* ignore user type
program _gvm
	version 8  
	syntax newvarname =/exp [if] [in] ///
	[ , Kappa(numlist max=1 >=0) Mu(str) radians ]  
	marksample touse, novarlist
	
	tempname z p  
	tempvar theta u y v s c sn cn r zz 

	if "`kappa'" == "" {
		tempname kappa 
		scalar `kappa' = r(kappa) 
		if `kappa' == . { 
			di as err "kappa needed" 
			exit 198 
		} 	
	} 
	
	if "`mu'" == "" {
		local mu = r(vecmean) 
		if `mu' == . { 
			di as txt "mu assumed " as res 0 
			local mu 0 
		} 
		else if "`radians'" != "" { 
			local mu = `mu' * _pi / 180 
		} 	
	} 

	local a1 = 12 
	local a2 = 0.8 
	local a3 = 8 
	local a4 = 1 
	local ck = 10.5 
	local c1 = 56

	quietly { 
		scalar `z' = `kappa'
		
		if "`radians'" != "" { 
			gen double `theta' = ///
			mod((`exp') - `mu' + 2 * _pi, 2 * _pi) 
		} 	
		else { 
			gen double `theta' = mod((`exp') - `mu' + 360, 360) 
			replace `theta' = `theta' - 360 if `theta' > 180  
			replace `theta' = _pi * `theta' / 180 
		}	
		
		markout `touse' `theta' 
		count if `touse' 
		if r(N) == 0 error 2000 
				
		gen double `u' = mod(`theta' + _pi, 2 * _pi) 
		replace `u' = `u' + 2 * _pi if `u' < 0 
		gen double `y' = `u' - _pi 

		if `z' <= `ck' { 
			gen double `v' = 0 
			if `z' > 0 { 
				scalar `p' = /// 
					`z' * `a2' - `a3' / (`z' + `a4') + `a1'
				local ip = int(`p')
				scalar `p' = `ip' 
				gen double `s' = sin(`y') 
				gen double `c' = cos(`y') 
				replace `y' = `p' * `y' 
				gen double `sn' = sin(`y') 
				gen double `cn' = cos(`y') 
				gen double `r' = 0 
				scalar `z' = 2 / `z'
				forval i = 2 / `ip' { 
					scalar `p' = `p' - 1 
					replace `y' = `sn' 
					replace `sn' = `sn' * `c' - `cn' * `s' 
					replace `cn' = `cn' * `c' + `y' * `s' 
					replace `r' = 1 / (`p' * `z' + `r') 
					replace `v' = (`sn' / `p' + `v') * `r'
				} 
			}
			gen double `varlist' = /// 
				(`u' * 0.5 + `v') / _pi if `touse' 
		} 
		else { 
			scalar `c' = 24.0 * `z' 
			gen double `v' = `c' - `c1' 
			gen double `r' = ///
			sqrt((54 / (347 / `v' + 26 - `c') - 6 + `c') / 6)
			gen double `zz' = sin(`y' * 0.5) * `r' 
			gen double `s' = `zz' * `zz' 
			replace `v' = `v' - `s' + 3 
			replace `y' = (`c' - `s' - `s' - 16) / 3
			replace `y' = ((`s' + 1.75) * `s' + 83.5) / `v' - `y' 
			gen double `varlist' = ///
				norm(`zz' - `s' / (`y' * `y') * `z') if `touse' 
		} 

		replace `varlist' = 0 if `varlist' < 0 
		replace `varlist' = 1 if `varlist' > 1  & `varlist' < . 
	} 	
end