This program takes three input pairs, A,B,C [in the form “x,y” when prompted], B being the point from which segments are joined to the other two forming tangents to the parabola constructed by forming the points P(t) = (1-t)Q(t) + t R(t), where Q(t)=(1-t)A+tB and R(T)=(1-t)B + tC for t in [0,1]. The proof that this forms a parabolic arc is given in oriolescience.blogspot.com entry 16SEP15. The program outputs the parabola which contains A and C and to which lines AB and CB are tangent. The program follows the logic of the proof.
To run the program just copy the program into a plain text file with extension .c and then compile it using "gcc name.c" and run "./a.out" which is the executable produced by gcc.
#include <stdio.h>
#include <math.h>
int main(void) {
char s[255],c;
long int i,j;
float x,y,z,a1,a2,b1,b2,c1,c2,d1,d2,mv,mu,mv1,mv2,mu1,mu2;
float A1,A2,B1,B2,C1,C2,D1,D2,ba1,ba2,ba,bc1,bc2,bc,xtemp,xp,aa,bb,cc,A,B,C,D,E,F,check;
float costheta,sintheta,theta;
int k,m,nr;
short n;
float xx,yy,zz;
/*Enter the three points*/
printf("Enter coordinates of Point A: x,y: ");
scanf("%f,%f", &a1,&a2);
printf("Enter coordinates of Point B x,y: ");
scanf("%f,%f", &b1,&b2);
printf("Enter coordinates of Point C: x,y: ");
scanf("%f,%f", &c1,&c2);
printf("A (%f,%f), B (%f,%f), C(%f,%f)\n",a1,a2,b1,b2,c1,c2);
/*Calculate midpoint of AC and determine unit vectors for positive u,v-axes*/
d1=(a1+c1)/2; d2=(a2+c2)/2; nr=-1;
if ((d1==b1) || (d2==b2))
{ /*Check for D with same 1st or 2nd coordinate as B*/
if (d1==b1) { /*D is directly above or below B*/
if(b2<d2) {nr=0;mu1=1;mu2=0;mv1=0;mv2=1;}
else {nr=180;mu1=-1;mu2=0; mv1=0; mv2=-1;};
}
else { /*D is directly to right or left of B*/
if (b1<d1) {nr=270; mu1=0;mu2=-1;mv1=1;mv2=0;}
else {nr=90; mu1=0;mu2=1;mv1=-1;mv2=0;};
}
}
else { /*There is a rotation which is not a multiple of 90*/
mv1=d1-b1; mv2=d2-b2; xtemp = sqrt( mv1*mv1 + mv2*mv2 );
mv1=mv1/xtemp; mv2=mv2/xtemp;
mu1=mv2; mu2=-mv1; /*make unit u = unit v cross unit k*/
}
/*get uv-coordinates for Point A and Point C as components along u, u dot BA, and v, v dot BA, etc*/
ba1=a1-b1;ba2=a2-b2;
A1= mu1*ba1 + mu2*ba2;
A2= mv1*ba1 + mv2*ba2;
bc1=c1-b1;bc2=c2-b2;
C1= mu1*bc1 + mu2*bc2;
C2= mv1*bc1 + mv2*bc2;
printf("unit u is (%f,%f)\n",mu1,mu2);
printf("unit v is (%f,%f)\n",mv1,mv2);
printf("Point A in {A,B,C} is (%f,%f)\n",A1,A2);
printf("Point C in {A,B,C} is (%f,%f)\n",C1,C2);
/*Now can write the coefficients of the parabola in standard form wrt uv-coordinates, i.e. in {A,B,C}
coordinate system as P(u) = au^2 + bu + c ,for u in A1 to -A1=C1 */
aa=(A2+C2)/(4*A1*A1);
bb=(A2-C2)/(2*A1);
cc=aa*A1*A1;
printf("P(u) = v = %6.3f u^2 + %6.3f u + %6.3f\n",aa,bb,cc);
/*Check that points are on this parabolic arc*/
xx=A1; yy=A2;
check= aa*A1*A1 + bb*A1 + cc - A2;
printf(" At (%5.2f, %5.2f) check = %5.4f \n",A1,A2,check);
xx=C1;yy=C2;
check= aa*C1*C1 + bb*C1 + cc - C2;
printf(" At (%5.2f, %5.2f) check = %5.4f \n",C1,C2,check);
/*Find the Vertex, Focus and Directrix with respect to the coefficients of P(u)*/
xtemp = cc - (bb*bb)/(4*aa); xp= 1/(4*aa);
printf("Vertex (%6.3f,%6.3f)\n", -bb/(2*aa), xtemp);
xtemp = xtemp + xp;
printf("Focus (%6.3f,%6.3f)\n", -bb/(2*aa), xtemp );
xtemp = xtemp - 2*xp;
printf("The directrix is the line v = %6.3f\n", xtemp);
/* Now get the equation in original XY-coordinates by applying appropriate rotation and translation*/
/*Rotation takes unit u-vector to unit i vector of XY-system*/
/* Theory: For simplicity of notation, identify the uv-coordinate system with the X"Y"-coordinate system.
Let Y"=aX"^2 + bX" + c be the equation in the {A,B,C} or X"Y"-coordinate system. This system is obtained
from the XY-coordinate system by first translating by B to the X'Y'-coordinate system, where X'=X-b1,
and Y'=Y-b2, which in turn is rotated by angle theta with costheta=mu1, sintheta=mu2, the coordinates
of the unit vector u of the X"Y"-coordinate system. That is, u=(mu1,mu2) is the positive unit vector in
the direction of the X"-axis. Then X" = X'costheta + Y'sintheta and Y" = -X'sintheta + Y'costheta. For
example the point with X'Y'-coordinates (1,0) will have X"Y"-coordinates (costheta, -sintheta). (Same point
has different "names" in the different coordinate systems.)
In our case we are going in reverse, from the X"Y"-system back through the X'Y'-system to the XY-system.
Starting with Y"=aX"^2 + bX" + c, we substitute (-X'sintheta + Y'costheta) for Y" and (X'costheta + Y'sintheta)
for X" to get an equation in X' and Y'. Now substitue X-b1 for X' and Y-b2 for Y' to get an equation in X and Y.
Expanding and collecting terms we get the expression A X^2 + B Y^2 + C XY + D X + E Y + F = 0, as below*/
costheta = mu1;
sintheta = mu2;
printf("Rotation theta has cosine %f and sine %f\n", costheta,sintheta);
printf("Translate by B and Rotate by theta (Y) = %6.3f (X)^2 + %6.3f (X) + %6.3f \nto get \n",aa,bb,cc);
A = aa*costheta*costheta; B = aa*sintheta*sintheta; C = aa*2*sintheta*costheta;
D = sintheta + (bb*costheta) - (2*aa*b2*sintheta*costheta) - (2*aa*b1*costheta*costheta);
E = -costheta + (bb*sintheta) - (2*aa*b1*sintheta*costheta) - (2*aa*b2*sintheta*sintheta);
F = (aa*b1*b1*costheta*costheta) +(aa*b2*b2*sintheta*sintheta) + (2*aa*b1*b2*sintheta*costheta);
F = F - b1*bb*costheta - b1*sintheta - b2*bb*sintheta + b2*costheta + cc;
printf(" %6.3f X^2 + %6.3f Y^2 + %6.3f XY + %6.3f X + %6.3f Y + %6.3f = 0\n", A,B,C,D,E,F);
/*CHECK original points in this parabolic equation*/
xx=a1;yy=a2;
check = (A*xx*xx) + (B*yy*yy) + (C*xx*yy) + (D*xx) + (E*yy) + F;
printf(" At (%5.2f, %5.2f) check = %5.4f \n",xx,yy,check);
xx=c1;yy=c2;
check = (A*xx*xx) + (B*yy*yy) + (C*xx*yy) + (D*xx) + (E*yy) + F;
printf(" At (%5.2f, %5.2f) check = %5.4f \n",xx,yy,check);
return 0;
}
/* Sample Program Prompts, Input and Output
Enter coordinates of Point A: x,y: 3,9
Enter coordinates of Point B x,y: 2.5,6
Enter coordinates of Point C: x,y: 1,4
A (3.000000,9.000000), B (2.500000,6.000000), C(1.000000,4.000000)
unit u is (0.707107,0.707107)
unit v is (-0.707107,0.707107)
Point A in {A,B,C} is (2.474874,1.767767)
Point C in {A,B,C} is (-2.474874,-0.353553)
P(u) = v = 0.058 u^2 + 0.429 u + 0.354
At ( 2.47, 1.77) check = -0.0000
At (-2.47, -0.35) check = 0.0000
Vertex (-3.712,-0.442)
Focus (-3.712, 3.889)
The directrix is the line v = -4.773
Rotation theta has cosine 0.707107 and sine 0.707107
Translate by B and Rotate by theta (Y) = 0.058 (X)^2 + 0.429 (X) + 0.354
to get
0.029 X^2 + 0.029 Y^2 + 0.058 XY + 0.520 X + -0.895 Y + 2.338 = 0
At ( 3.00, 9.00) check = 0.0000
At ( 1.00, 4.00) check = 0.0000
*/
