*! arc v1.2 (20 Nov 2024) *! Asjad Naqvi (asjadnaqvi@gmail.com) * v1.2 (20 Nov 2024): append/stack added * 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) genid(string) genorder(string) ] /// [ replace append stack ] if `n' < 5 { display in red "The number of points n() must be greater or equal to 5." exit } quietly { // prepare the variables local xvar _x local yvar _y local idvar _id local ordervar _order if "`genx'" != "" local xvar `genx' if "`geny'" != "" local yvar `geny' if "`genid'" != "" local idvar `genid' if "`genorder'" != "" local ordervar `genorder' if "`replace'" != "" { cap drop `xvar' cap drop `yvar' cap drop `idvar' cap drop `ordervar' } cap generate double `xvar' = . cap generate double `yvar' = . cap generate `idvar' = . cap generate `ordervar' = . if "`stack'" != "" | "`append'" !="" { tempvar temp gen `temp' = _n if !missing(`ordervar') summ `temp', meanonly if r(N)!= 0 { local start = `=`r(max)' + 1' local end = `=`r(max)' + `n' + 1' if _N < `=`r(max)'+`n'+1' set obs `end' } else { local start = 1 local end = `=`n' + 1' if _N < `n' set obs `end' } } else { local start = 1 local end = `=`n' + 1' if _N < `n' set obs `end' } summ `idvar', meanonly if `r(N)'==0 { local k = 1 } else { local k = `r(max)' + 1 } // id summ `ordervar', meanonly replace `idvar' = `k' in `start'/`end' // order tempvar _temp _temp2 _seq _angle gen `_temp' = 1 replace `ordervar' = sum(`_temp') if `idvar'==`k' **** core routines below **** // mid points local xm = (`x1' + `x2') / 2 local ym = (`y1' + `y2') / 2 // slope local dx = `x2' - `x1' local dy = `y2' - `y1' local L = sqrt(abs(`dx')^2 + abs(`dy')^2) // check why this square is not working without abs() *noisily display "`dx', `dy', `L'" // 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 } // 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' + (`ordervar' - 1) * `delta_theta' in `start'/`=`end'-1' replace `xvar' = `h' + `radius' * cos(`theta') in `start'/`=`end'-1' replace `yvar' = `k' + `radius' * sin(`theta') in `start'/`=`end'-1' } return local chord = `L' return local radius = `radius' return local xcirc = `h' return local ycirc = `k' end