*! arc v1.1 (11 Oct 2024)
*! Asjad Naqvi (asjadnaqvi@gmail.com)

* v1.1 (11 Oct 2024): Minor bug fixes.
* v1.0 (08 Oct 2024): first release.


cap program drop arc

program define arc, rclass

version 11

	syntax, ///
		x1(numlist max=1) y1(numlist max=1) ///  //  from
		x2(numlist max=1) y2(numlist max=1) ///  //  to
		[ n(real 40) RADius(numlist max=1 >0) major swap ]	///
		[ genx(string) geny(string)	replace ]


    if `n' < 5 {
        display in red "The number of points n() must be greater than 5."
        exit
    }


	// mid points
    local xm = (`x1' + `x2') / 2
    local ym = (`y1' + `y2') / 2

    // slope
    local dx = `x2' - `x1'
    local dy = `y2' - `y1'
	local L = sqrt(`dx'^2 + `dy'^2)


    // radius
	if "`radius'" == "" {
		local radius = `L'   // radius = chord length
		display in yellow "Radius of `radius' assumed."
	}

    // validate the radius
    if `radius' <= `L'/2 {
		local mychord = `L' / 2
        di as error "Radius must be greater than half the chord length. Half chord length for given coordinates = `mychord'"
        exit
    }

	return local chord  `L'
	return local radius `radius'

    // distance from midpoint to center
    local h_dist = sqrt(`radius'^2 - (`L'/2)^2)


    // direction of the perpendicular bisector
    local ux_perp = -`dy' / `L'
    local uy_perp =  `dx' / `L'

	local cross_prod = `dx' * `uy_perp' - `dy' * `ux_perp'

	// swap left to right orientation
    if "`swap'" == "" {
		local h = `xm' + `h_dist' * `ux_perp'
		local k = `ym' + `h_dist' * `uy_perp'
	}
	else {
		local h = `xm' - `h_dist' * `ux_perp'
		local k = `ym' - `h_dist' * `uy_perp'
	}


   // get the angles in order
   local angle_start = atan2(`y1' - `k', `x1' - `h')
   local angle_end   = atan2(`y2' - `k', `x2' - `h')

    // bound angles between 0, 2pi
    if `angle_start' < 0	local angle_start = `angle_start' + 2 * _pi
    if `angle_end'   < 0	local angle_end   = `angle_end'   + 2 * _pi


	// correct for minor vs major arc

	if "`swap'"!= ""  {
		if "`major'" != "" {
			if (`angle_end' - `angle_start') < _pi local angle_end = `angle_end' + 2 * _pi
		}
		else {
			if (`angle_end' - `angle_start') > _pi local angle_end = `angle_end' - 2 * _pi
		}
	}
	else {
		if "`major'" != "" {
			if (`angle_end' - `angle_start') > _pi 		local angle_end = `angle_end' - 2 * _pi
		}
		else {
			if (`angle_end' - `angle_start') < _pi 		local angle_end = `angle_end' + 2 * _pi
		}
		if (`angle_end' - `angle_start') > 2 * _pi 	local angle_end = `angle_end' - 2 * _pi
	}

    local delta_theta = (`angle_end' - `angle_start') / (`n' - 1)

	if _N < `n' set obs `n'

	tempvar theta
	gen double `theta' = `angle_start' + (_n - 1) * `delta_theta'

	local xvar _x
	local yvar _y

	if "`genx'" != "" local xvar `genx'
	if "`geny'" != "" local yvar `geny'

	if "`replace'" != "" {
		capture drop `xvar'
		capture drop `yvar'
	}

	gen double `xvar' = `h' + `radius' * cos(`theta')
	gen double `yvar' = `k' + `radius' * sin(`theta')

	return local xcirc = `h'
	return local ycirc = `k'

end