How to plot arc that passes through 3 points?

Question

I am trying to plot a circular arc that passes through given 3 control points (i.e. [P1,P2,P3]). I found the circle equation and the angles of the points onto it. However, I could not handle determining the direction (CCW, CW). I was getting erroneous results such as this:

The yellow arc must be on the other side obviously. I am open to the solution for handling this angular problem or direct program that generates arc through 3 points.

Please note that the green curve is unrelated to the question.

Answer 1

After considering all the possibilities, I finally wrote a code that generates arc passes through 3 points. Here is the calculation of thetaStart and thetaEnd angles;

    function  obj = fitArc(obj,points)
        
        x1 = points(1,1);
        x2 = points(2,1);
        x3 = points(3,1);
        y1 = points(1,2);
        y2 = points(2,2);
        y3 = points(3,2);
        
        A = x1*(y2-y3)-y1*(x2-x3)+x2*y3-x3*y2;
        B = (x1^2+y1^2)*(y3-y2)+(x2^2+y2^2)*(y1-y3)+(x3^2+y3^2)*(y2-y1);
        C = (x1^2+y1^2)*(x2-x3)+(x2^2+y2^2)*(x3-x1)+(x3^2+y3^2)*(x1-x2);
        D = (x1^2+y1^2)*(x3*y2-x2*y3)+(x2^2+y2^2)*(x1*y3-x3*y1)+(x3^2+y3^2)*(x2*y1-x1*y2);
        
        obj.x = -B/(2*A);
        obj.y = -C/(2*A);
        obj.r = sqrt((B^2+C^2-4*A*D)/(4*A^2));
        
        th1 = wrapTo2Pi(atan2((points(3,2)-obj.y),(points(3,1)-obj.x)));
        th3 = wrapTo2Pi(atan2((points(1,2)-obj.y),(points(1,1)-obj.x)));
        th2 = wrapTo2Pi(atan2((points(2,2)-obj.y),(points(2,1)-obj.x)));
        
        if (th1 < th2 && th2 < th3)
            obj.thStart = th1;
            obj.thEnd = th3;
            
        elseif (th1<th3 && th3<th2)
            obj.thStart = th3;
            obj.thEnd = th1+2*pi;
            
        elseif (th2<th1 && th1<th3)
            obj.thStart = th3;
            obj.thEnd = th1+2*pi;
            
        elseif (th2<th3 && th3<th1)
            obj.thStart = th1;
            obj.thEnd = th3+2*pi;
            
        elseif (th3<th1 && th1<th2)
            obj.thStart = th1;
            obj.thEnd = th3+2*pi;
            
        elseif (th3 < th2 && th2 < th1)
            obj.thStart = th3;
            obj.thEnd = th1;
            
        end
    end

After finding the start and end angles plotting function should be like:

th = linspace(thStart,thEnd,n);
xunit = r * cos(th) + x;
yunit = r * sin(th) + y;
plot(xunit, yunit,'linewidth',3);

Here are random tests;

I hope this solution would help so many people as I have seen lots of unanswered questions about this problem.