randomfox (randomfox) wrote,
randomfox
randomfox

PARI/GP script to find the center of a circle, given 3 points on its circumference

PARI is a high-precision math library and computer algebra system. GP is a scripting interface to that library.

Coordinates should be converted to UTM for use with the script.

To run a script from GP, enter \r filename on the command line. If the script fails to run because some of the symbols already exist, enter kill(circle), for example, before running the script, or start a fresh GP session.

dist(a, b) = sqrt(sqr(a[1]-b[1]) + sqr(a[2] - b[2]));

circle(b, c, d) = 
{
    local(temp, bc, cd, det, circ);

    temp = sqr(c[1]) + sqr(c[2]);
    bc = (sqr(b[1]) + sqr(b[2]) - temp) / 2.;
    cd = (temp - sqr(d[1]) - sqr(d[2])) / 2.;
    det = (b[1]-c[1]) * (c[2]-d[2]) - (c[1]-d[1]) * (b[2]-c[2]);

    circ = [
	(bc * (c[2]-d[2]) - cd * (b[2]-c[2])) / det,
	((b[1]-c[1]) * cd - (c[1]-d[1]) * bc) / det
    ];

    [ circ, [ dist(circ, b), dist(circ, c), dist(circ, d) ] ];
}

b = [515939, 4413282];
c = [519054, 4405649];
d = [516437, 4399896];

circle(b, c, d)

Subscribe
  • Post a new comment

    Error

    Anonymous comments are disabled in this journal

    default userpic

    Your reply will be screened

    Your IP address will be recorded 

  • 0 comments