Найти координаты центра окружности

#include "stdafx.h"
#include <iostream>
 
const int N=6;
struct TPoint {
    double x, y;
};
 
struct TPoint PS[N]={{0,0},{2,0},{4,1},{4,3},{2,5},{1,2}};
double Dist(TPoint a, TPoint b, TPoint c);
double MinDist(void);
 
int main() {
    double RMin;
    RMin=MinDist();
    printf("R=%6.2f\n", RMin );
    system("pause");
}
 
double Dist(TPoint a, TPoint b, TPoint c){
    double R, x0, y0;
    R=(sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)))*     /*длина стороны 1 треугольника*/
    (sqrt((a.x-c.x)*(a.x-c.x)+(a.y-c.y)*(a.y-c.y)))*       /*длина стороны 2 треугольника*/
    (sqrt((c.x-b.x)*(c.x-b.x)+(c.y-b.y)*(c.y-b.y)))/       /*длина стороны 3 треугольника*/
 
    (sqrt((8*((sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)))+ /*восемь*/
    (sqrt((a.x-c.x)*(a.x-c.x)+(a.y-c.y)*(a.y-c.y)))+           /*умноженное на*/
    (sqrt((c.x-b.x)*(c.x-b.x)+(c.y-b.y)*(c.y-b.y))))*          /*периметр треугольника*/
    
    (((sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)))+    /*полупериметр*/
    (sqrt((a.x-c.x)*(a.x-c.x)+(a.y-c.y)*(a.y-c.y)))+      /*минус*/
    (sqrt((c.x-b.x)*(c.x-b.x)+(c.y-b.y)*(c.y-b.y))))/2-   /*сторона*/ 
    (sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y))))*     /*1*/
 
    (((sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)))+    /*полупериметр*/
    (sqrt((a.x-c.x)*(a.x-c.x)+(a.y-c.y)*(a.y-c.y)))+      /*минус*/
    (sqrt((c.x-b.x)*(c.x-b.x)+(c.y-b.y)*(c.y-b.y))))/2-   /*сторона*/ 
    (sqrt((a.x-c.x)*(a.x-c.x)+(a.y-c.y)*(a.y-c.y))))*     /*2*/
 
    (((sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y)))+    /*полупериметр*/
    (sqrt((a.x-c.x)*(a.x-c.x)+(a.y-c.y)*(a.y-c.y)))+      /*минус*/
    (sqrt((c.x-b.x)*(c.x-b.x)+(c.y-b.y)*(c.y-b.y))))/2-   /*сторона*/
    (sqrt((c.x-b.x)*(c.x-b.x)+(c.y-b.y)*(c.y-b.y)))))));  /*3*/
    x0=(a.x+b.x+c.x)/3;
    y0=(a.y+b.y+c.y)/3;
    printf("x0=%6.2f y0=%6.2f\n", x0, y0);
    return R;
}
double MinDist(void) {
    int i, j, k;
    double r, RMin;
    RMin=Dist(PS[0],PS[1],PS[2]);
    for (i=0;i<N-2;i++)
        for (j=i+1;j<N-1;j++)
            for (k=i+2; k<N; k++){
            r=Dist(PS[i],PS[j],PS[k]);
            if (r<RMin) {
                RMin=r;
            }
        }
    return RMin;
}