#version 3.8;
global_settings {assumed_gamma 1.0 }
#default {finish {emission 1}}

// Given 3 points:
// Solve a system of linear equations in 3 unknows (a, b, c)
// to find the equation of a parabola
// Bill "Bald Eagle" Walker 	July/2024

// Make render fill the screen
#declare Width = image_width;
#declare Height = image_height;
#declare ScreenMin = min (Width, Height);

camera {
	location <0, 0, -30>
	right     x*image_width/image_height
	up y
	look_at <0, 0, 0>
	
}

light_source {< 0, 0, -250> rgb 1}
#declare Line = 0.01;
#declare E = 0.0001;

#declare Gridline = Line*4;
#declare fmod = function (Value, Modulo) {select (Value, 1 - mod (abs (Value), Modulo), mod (abs (Value), Modulo))}
#declare Grid = function {(fmod (x,1) > Gridline) * (fmod (y,1) > Gridline)}
#declare Grid_PM = pigment_map {
	[0 rgb y+z]
	[0.01 rgb y+z]
	[0.01 rgb 1]
	[1 rgb 1]
}


#declare Extent = 15;
box {
	<-Extent, -Extent, 0>, <Extent+Line*2, Extent+Line*2, 0.0001> 
	texture {pigment {function {Grid(x, y, z)} pigment_map {Grid_PM} scale 1} finish {diffuse 1} } 
	translate z*Line*10
}
sphere {0, Line*10 pigment {rgb 0}}

#macro gcd (a, b)

	#if (b = 0)
		#local Result = a;
	#else
		#local Result = gcd (b, mod (a, b));
	#end
	
	Result
#end


#macro gcdlcm (a, b)
	#local GCD = gcd (a, b);
	#local LCM = (a*b)/GCD;
	#local Result = array [2] {GCD, LCM};
	
	Result
#end

#macro Str (N) str (N, 0, 0) #end
#declare Sgn = array [3] {"-", "0", "+"}
#declare P1 = <1, 3>;
#declare P2 = <4, 9>;
#declare P3 = <2, -10>;

// a * pow(x,2) + b * x + c = y

#declare x2_1 = P1.x*P1.x;
#declare  x_1 = P1.x;
#declare  y_1 = P1.y;

#declare x2_2 = P2.x*P2.x;
#declare  x_2 = P2.x;
#declare  y_2 = P2.y;

#declare x2_3 = P3.x*P3.x;
#declare  x_3 = P3.x;
#declare  y_3 = P3.y;

#macro Sign (X, Y)
	#declare sx = select (X, 0, 1, 2);
	#declare sy = select (Y, 0, 1, 2);
#end

#debug "Substitute x and y into general equation of a parabola and simplify: \n"
Sign (x_1, y_1)
#debug concat ("1) ", Str(x2_1), "a ", Sgn[sx], Str(x_1), "b + c = ", Sgn[sy], Str(y_1), "\n")
Sign (x_2, y_2)
#debug concat ("2) ", Str(x2_2), "a ", Sgn[sx], Str(x_2), "b + c = ", Sgn[sy], Str(y_2), "\n")
Sign (x_3, y_3)
#debug concat ("3) ", Str(x2_3), "a ", Sgn[sx], Str(x_3), "b + c = ", Sgn[sy], Str(y_3), "\n")
#debug "\n"

#debug "Eliminate c by subtracting 2 equations with smaller a values from equation with largest a value: \n"
#declare EqN = array [3] {0, 1, 2}
#declare Eq123  = array [3][3]
{
	{x2_1, x_1, y_1},
	{x2_2, x_2, y_2},
	{x2_3, x_3, y_3}
}

#macro Match (Array, Value)
	#local Elements = dimension_size (Array, 1)-1;
	#local Result = val ("nan");
	#for (E, 0, Elements)
		#if (Value = Array [E])
			#local Result = E;
		#else
			// keep looping
		#end
	#end
	Result
#end

#macro Arraymax (Array)
	#local Elements = dimension_size (Array, 1)-1;
	#local Max = 0;
	#for (E, 0, Elements)
		#local Max = max (Array [E], Max);
	#end
	Max
#end

#macro Arraymin (Array)
	#local Elements = dimension_size (Array, 1)-1;
	#local Min = Elements;
	#for (E, 0, Elements)
		#local Min = min (Array [E], Min);
	#end
	Min
#end

#declare A = array [3] {x2_1, x2_2, x2_3}
#declare Largest  = Match (A, Arraymax(A));
#declare Smallest = Match (A, Arraymin(A));
#declare Middle   = (0+1+2) - (Largest + Smallest); // gives the remaining index

#declare Eq45 = array [2][3]
{
	{Eq123 [Largest][0] - Eq123 [Smallest][0], Eq123 [Largest][1] - Eq123 [Smallest][1], Eq123 [Largest][2] - Eq123 [Smallest][2]},
	{Eq123 [Largest][0] - Eq123   [Middle][0], Eq123 [Largest][1] - Eq123   [Middle][1], Eq123 [Largest][2] - Eq123   [Middle][2]}
}

Sign (Eq45[0][1], Eq45[0][2])
#debug concat ("4) ", Str(Eq45[0][0]), "a ", Sgn[sx], Str(Eq45[0][1]), "b = ", Sgn[sy], Str(Eq45[0][2]), "\n")
Sign (Eq45[1][1], Eq45[1][2])
#debug concat ("5) ", Str(Eq45[1][0]), "a ", Sgn[sx], Str(Eq45[1][1]), "b = ", Sgn[sy], Str(Eq45[1][2]), "\n")
#debug "\n"

#declare GCDLCM = gcdlcm (Eq45[0][1], Eq45[1][1]);
#declare LCMb = GCDLCM[1];
#declare Eq67 = array [2][3]
{
	{Eq45[0][0]*LCMb/Eq45[0][1], LCMb, Eq45[0][2]*LCMb/Eq45[0][1]},
	{Eq45[1][0]*LCMb/Eq45[1][1], LCMb, Eq45[1][2]*LCMb/Eq45[1][1]}
}

#debug "Eliminate b by multiplying by Lowest Common Multiple and subtracting: \n"
Sign (Eq67[0][1], Eq67[0][2])
#debug concat ("6) ", Str(Eq67[0][0]), "a ", Sgn[sx], Str(Eq67[0][1]), "b = ", Sgn[sy], Str(Eq67[0][2]), "\n")
Sign (Eq67[1][1], Eq67[1][2])
#debug concat ("7) ", Str(Eq67[1][0]), "a ", Sgn[sx], Str(Eq67[1][1]), "b = ", Sgn[sy], Str(Eq67[1][2]), "\n")
#debug "\n"

#declare A2 = array [2] {Eq67[0][0], Eq67[1][0]}
#declare Largest  = Match (A2, Arraymax(A2));
#declare Smallest = 1-Largest;

#declare AM = Eq67[Largest][0] - Eq67[Smallest][0];
#declare YM = Eq67[Largest][2] - Eq67[Smallest][2];


Sign (AM, YM)
#debug concat (Str(AM), "a = ", Sgn[sy], Str(YM), "\n")
#declare A = YM/AM;
#debug concat ("a = ", Str(A), "\n")
#debug "\n"

#declare Eq45a = array [2][3]
{
	{Eq45[0][0]*A, Eq45[0][1], Eq45[0][2]},
	{Eq45[1][0]*A, Eq45[1][1], Eq45[1][2]}
}

#debug "Substitute a into Equations 4 & 5 to get b: \n"
Sign (Eq45a[0][1], Eq45a[0][2])
#debug concat ("4a) ", Str(Eq45a[0][0]), "a ", Sgn[sx], Str(Eq45a[0][1]), "b = ", Sgn[sy], Str(Eq45a[0][2]), "\n")
Sign (Eq45a[1][1], Eq45a[1][2])
#debug concat ("5a) ", Str(Eq45a[1][0]), "a ", Sgn[sx], Str(Eq45a[1][1]), "b = ", Sgn[sy], Str(Eq45a[1][2]), "\n")
#debug "\n"

#declare B1 = (Eq45a[0][2] - Eq45a[0][0]) / Eq45a[0][1];
#declare B2 = (Eq45a[1][2] - Eq45a[1][0]) / Eq45a[1][1];
#debug concat ("b = ", Str(B1), "\n")
#debug concat ("b = ", Str(B2), "\n")
#if (B1=B2) #debug "(OK: B-values from both 4 & 5 match) \n" #else "B-value mismatch! \n" #end
#debug "\n"

#debug "Substitute a and b into Equation 1 to get c: \n"
#declare C = y_1 - (x2_1*A + x_1*B1);
#debug concat ("c = ", Str(C), "\n")
#debug "\n"

Sign (B1, C)
#debug concat ("The equation of the parabola is: y = ", Str(A), "x^2 ", Sgn[sx], Str(B1), "x ", Sgn[sy], Str(C), "\n")

//#error "No objects in scene"

#declare Parabola = function (_X, _A, _B, _C) {_A*_X*_X + _B*_X + _C}

union {
	sphere {P1, Line*10}
	sphere {P2, Line*10}
	sphere {P3, Line*10}
	pigment {rgb x}
}

#macro DrawParabola (Ext, L)
	union {
		#local Counter = 0;
		#for (X, -Ext, Ext, 0.1)
			#local Y = Parabola (X, A, B1, C);
			#local Current = <X, Y, 0>;
			sphere {Current, L}
			#if (Counter > 0)
				cylinder {Last, Current, L}
			#end
			#local Last = Current;
			#local Counter = Counter + 1;
		#end
		pigment {rgb z}
	}
#end

DrawParabola (Extent, Line*5)