I am solving this problem-
U and K are vectors is R3. U rotated around K by an angle θ in radian produces a new vector V. Find the components of V as functions of U , K and θ.
I solved the problem using array and trigonometry -
First I converted maths functions like sin, cos, sqrt etc into C functions by using if and for. I use Taylor Series for sine and cosine and use the Newton-Raphson Method for sqrt. I rationalize the numbers by truncating the series after 32 runs.
I define vector products like scalar product, dot product, cross product. I also define sums of vectors.
Then I perform rotation by using the method prescribed in here. I do this by calling functions in void rotate ().
#include <stdio.h>
#include <math.h>
double norm (double *vec)
{
return (vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
}
double dot (double *vec1, double *vec2)
{
return (vec1[0]*vec2[0] + vec1[1]*vec2[1] + vec1[2]*vec2[2]);
}
void cross (double *vec1, double *vec2, double *res)
{
res[0] = vec1[1]*vec2[2] - vec1[2]*vec2[1];
res[1] = vec1[2]*vec2[0] - vec1[0]*vec2[2];
res[2] = vec1[0]*vec2[1] - vec1[1]*vec2[0];
}
void sca (double scaler, double *vec, double *res)
{
for (int i =0 ; i < 3; i++)
{
res[i] = scaler * vec[i];
}
}
void sum (double *vec1, double *vec2, double *vec3, double *V)
{
for (int i =0 ; i < 3; i++)
{
V[i] = vec1[i] + vec2[i] + vec3[i];
}
}
void rotate (double *K, double *U, double angle, double *V)
{
double proj[3], rej_para[3], rej_perp[3], temp1[3], temp2[3];
sca ((dot (K, U) / norm(K)), K, proj);
cross (K, U, temp1);
sca ((sin(angle) / sqrt(norm(K))), temp1, rej_perp);
cross (temp1, K, temp2);
sca ((cos(angle) / norm(K)), temp2, rej_para);
sum (proj, rej_para, rej_perp, V);
}
int main (void)
{
for (;;)
{
double U[3], K[3], V[3], angle;
char axes[] = "xyz", vectors[] = "UKV";
for (int i = 0; i < 2; i++)
{
for (int j = 0; j < 3; j++)
{
printf ("Give me the %c component of vector %c : ", axes[j], vectors[i]);
if (i == 0)
{
if (scanf ("%lf", U+j) != 1)
{
printf ("ERROR!!");
break;
}
}
else
{
if (scanf ("%lf", K+j) != 1)
{
printf ("ERROR!!");
break;
}
}
}
}
printf ("Give me an angle of rotation in radian : ");
if (scanf ("%lf", &angle) != 1)
{
printf ("ERROR!!");
break;
}
rotate (K, U, angle, V);
printf ("\n\n");
for (int i = 0; i < 3; i++)
{
printf ("The %c component of the rotated vector V = %.8lf\n", axes[i], V[i]);
}
printf ("\n\n");
}
return 0;
}
There is however another way, a more elegant and computationally more efficient way of getting this job done, which is to use struct and quaternion. So my question is this - how do I rotate vector using struct and quaternions just like it can be done using arrays and trig?
You can test this - no trigonometry
#include <stdio.h>
#include <math.h>
#include <float.h>
typedef struct
{
double x, y, z;
} Vec3;
/* ===== Vector helpers (ASCII only) ===== */
static inline Vec3 v_add(Vec3 a, Vec3 b)
{
return (Vec3){ a.x + b.x, a.y + b.y, a.z + b.z };
}
static inline Vec3 v_sub(Vec3 a, Vec3 b)
{
return (Vec3){ a.x - b.x, a.y - b.y, a.z - b.z };
}
static inline Vec3 v_scale(Vec3 a, double s)
{
return (Vec3){ a.x * s, a.y * s, a.z * s };
}
static inline double v_dot(Vec3 a, Vec3 b)
{
return a.x*b.x + a.y*b.y + a.z*b.z;
}
static inline Vec3 v_cross(Vec3 a, Vec3 b)
{
return (Vec3)
{
a.y*b.z - a.z*b.y,
a.z*b.x - a.x*b.z,
a.x*b.y - a.y*b.x
};
}
static inline double v_norm(Vec3 a)
{
return sqrt(v_dot(a, a));
}
static inline Vec3 v_project(Vec3 a, Vec3 khat)
{
return v_scale(khat, v_dot(khat, a));
}
static inline Vec3 v_perp(Vec3 a, Vec3 khat)
{
return v_sub(a, v_project(a, khat));
}
/* ===== Quaternion rotation (expects UNIT quaternion q = (w, qv)) ===== */
static inline Vec3 rotate_by_unit_quaternion(Vec3 V, double w, Vec3 qv)
{
/* t = 2 * (qv x V); V' = V + w*t + qv x t */
Vec3 t = v_scale(v_cross(qv, V), 2.0);
Vec3 Vp = v_add(V, v_add(v_scale(t, w), v_cross(qv, t)));
return Vp;
}
/* Axis-angle via tau = tan(theta/2). No trig functions used.
If K is degenerate (|K| ~ 0), returns U unchanged.
*/
static inline Vec3 rotate_axis_tau_quat(Vec3 U, Vec3 K, double tau)
{
const double eps = 1e-30;
Vec3 V = U;
double kn2 = v_dot(K, K);
if(kn2 > eps)
{
double inv_kn = 1.0 / sqrt(kn2);
Vec3 khat = v_scale(K, inv_kn);
/* Unit quaternion from tau:
w = 1 / hypot(1, tau)
s = tau / hypot(1, tau)
qv = khat * s
*/
double h = hypot(1.0, tau); /* robust sqrt(1 + tau^2) */
double w = 1.0 / h;
double s = tau / h;
Vec3 qv = v_scale(khat, s);
V = rotate_by_unit_quaternion(U, w, qv);
}
return V;
}
/* ===== Oracle: Rodrigues with tau (no trig)
V = U + s*(k x U) + t2*(k x (k x U)),
where k is unit axis, s = 2*tau/(1+tau^2), t2 = 2*tau^2/(1+tau^2).
*/
static inline Vec3 rodrigues_tau(Vec3 U, Vec3 K, double tau)
{
const double eps = 1e-30;
Vec3 V = U;
double kn2 = v_dot(K, K);
if(kn2 > eps)
{
Vec3 khat = v_scale(K, 1.0 / sqrt(kn2));
double denom = 1.0 + tau*tau;
double s = (2.0 * tau) / denom;
double t2 = (2.0 * tau * tau) / denom;
Vec3 kxU = v_cross(khat, U);
Vec3 kx_kxU = v_cross(khat, kxU);
V = v_add(U, v_add(v_scale(kxU, s), v_scale(kx_kxU, t2)));
}
return V;
}
/* ===== Close check ===== */
static int vec3_close(Vec3 a, Vec3 b, double atol, double rtol)
{
Vec3 d = v_sub(a, b);
double na = v_norm(a), nb = v_norm(b);
double tol = atol + rtol * fmax(na, nb);
return (fabs(d.x) <= tol) && (fabs(d.y) <= tol) && (fabs(d.z) <= tol);
}
/* ===== Tests ===== */
typedef struct
{
Vec3 U;
Vec3 K;
double tau; /* tau = tan(theta/2), no trig used */
const char *label;
} TestCase;
int main(void)
{
const double ATOL = 1e-12;
const double RTOL = 1e-12;
/* Precomputed tau values (all ASCII comments) */
const double TAU_15 = 0.26794919243112270647; /* tan(15 deg) */
const double TAU_22_5 = 0.19891236737965800691; /* tan(22.5 deg) */
const double TAU_30 = 0.57735026918962576451; /* tan(30 deg) = 1/sqrt(3) */
const double TAU_45 = 1.0; /* tan(45 deg) */
const double TAU_60 = 1.7320508075688772935; /* tan(60 deg) = sqrt(3) */
const double TAU_HUGE = 1e12; /* approx 180 deg (limit) */
TestCase cases[] =
{
{{ 1, 0, 0}, {0, 0, 1}, TAU_45, "x about +Z ~90deg (tau=1)"},
{{ 0, 1, 0}, {0, 0, 1}, TAU_45, "y about +Z ~90deg"},
{{ 1, 0, 0}, {0, 1, 0}, TAU_45, "x about +Y ~90deg"},
{{ 0, 0, 1}, {1, 0, 0}, TAU_45, "z about +X ~90deg"},
{{ 1, 0, 0}, {0, 0, 1}, TAU_HUGE, "x about +Z ~180deg (tau>>1)"},
{{ 1, 0, 0}, {0, 0, 1}, 0.0, "identity (tau=0)"},
{{ 1, 2, 3}, {0, 0, 1}, TAU_22_5, "U(1,2,3) about Z ~45deg"},
{{ 1, 2, 3}, {0, 1, 0}, TAU_22_5, "U(1,2,3) about Y ~45deg"},
{{ 1, 2, 3}, {1, 0, 0}, TAU_22_5, "U(1,2,3) about X ~45deg"},
{{ 1, 2, 3}, {1, 1, 1}, TAU_30, "diag axis (1,1,1) ~60deg"},
{{ 1, 2, 3}, {2, 0, 0}, -TAU_22_5, "non-unit X, ~-45deg"},
{{ 1, 2, 3}, {0, 3, 0}, -TAU_45, "non-unit Y, ~-90deg"},
{{ 1, 2, 3}, {0, 0, 5}, TAU_60, "non-unit Z, ~120deg"},
{{ 0, 0, 0}, {0, 0, 1}, TAU_30, "zero vector U"},
{{ 1, 0, 0}, {0, 0, 0}, TAU_30, "degenerate axis K=0"},
{{ 1e-9,-2e-9,3e-9},{1,0,0},TAU_15, "tiny U, about X ~30deg"},
{{ 1, 0, 0}, {1, 0, 0}, TAU_45, "U parallel K"},
{{ 0, 1, 0}, {1, 0, 0}, TAU_15, "U perpendicular K (+) ~30deg"},
{{ 0, 1, 0}, {1, 0, 0}, -TAU_15, "U perpendicular K (-) ~-30deg"},
{{ 0, 0, 1}, {0, 1, 0}, TAU_15, "z about Y ~30deg"},
{{-1, 0, 0}, {0, 0, 1}, TAU_30, "-x about Z ~60deg"},
{{ 0,-1, 0}, {0, 0, 1}, TAU_30, "-y about Z ~60deg"},
{{ 0, 0,-1}, {0, 1, 0}, TAU_30, "-z about Y ~60deg"},
{{ 3, 4, 0}, {0, 0, 1}, TAU_22_5, "U(3,4,0) about Z ~45deg"},
{{ 5,-2, 7}, {-1,2,-3}, 0.7, "general case #1 (tau=0.7)"},
{{ 1,-1, 1}, { 2,-1,0.5},2.5, "general case #2 (tau=2.5)"},
{{-3,0.5, 2}, {0.1,0.2,0.3},0.001, "tiny angle (tau=1e-3)"},
{{ 8,-6, 2}, {0, 0, 1}, -TAU_45, "~270deg about Z (tau=-1)"},
{{ 8,-6, 2}, {0, 0, 1}, TAU_45, "~-270deg about Z (tau=+1)"},
{{ 1, 2, 3}, {0, 0, 1}, 0.0, "full turn ~360deg (tau=0)"},
{{ 1, 2, 3}, {1, 1, 0}, TAU_HUGE, "~pi about (1,1,0)"},
{{ 1, 2, 3}, {1,-1, 0}, 0.23, "approx ~26deg (tau=0.23)"},
{{-1,-2,-3}, {0, 1, 1}, -0.23, "approx ~-26deg (tau=-0.23)"},
{{0.3,0.4,0.5},{0.6,-0.7,0.8},0.9, "general case #3 (tau=0.9)"},
{{ 1, 0, 0}, {0, 1, 0}, 0.0, "~720deg about Y (tau=0)"},
{{ 2, 3, 4}, {0, 1, 0}, 3.0, "large tau about Y (3.0)"},
{{ 2, 3, 4}, {1, 0, 0}, -3.0, "large tau about X (-3.0)"},
{{ 2, 3, 4}, {1, 2, 3}, -1.0, "general case #4 (tau=-1)"},
{{-2,-3,-4}, {1, 2, 3}, 1.0, "general case #5 (tau=1)"},
{{ 0, 0, 1}, {1, 1, 1}, 2.2, "diag axis vs z (tau=2.2)"},
{{ 0, 1, 0}, {1, 1, 1}, -2.2, "diag axis vs y (tau=-2.2)"},
{{ 1, 0, 0}, {1, 1, 1}, 2.2, "diag axis vs x (tau=2.2)"},
{{ 1, 1, 0}, {0, 1, 0}, TAU_45, "xy about Y ~90deg"},
{{ 1, 1, 1}, {0, 0, 1}, TAU_15, "sum vec about Z ~30deg"},
{{ 1, 1, 1}, {0, 1, 0}, TAU_15, "sum vec about Y ~30deg"},
{{ 1, 1, 1}, {1, 0, 0}, TAU_15, "sum vec about X ~30deg"},
{{ 0, 1, 1}, {1, 0, 1}, 0.333, "skew axis #1 (tau=0.333)"},
{{ 1, 0, 1}, {0, 1, 1}, -0.333, "skew axis #2 (tau=-0.333)"},
{{ 1e6,-1e6,5e5},{0,0,1},TAU_30, "large magnitude U, ~60deg"},
{{ 1, 2, 3}, {1e-16,0,0},TAU_22_5, "nearly-zero axis (degenerate)"}
};
const int N = (int)(sizeof(cases) / sizeof(cases[0]));
int passed = 0, failed = 0;
double worst_L2 = 0.0, worst_Linf = 0.0;
int worst_i_L2 = -1, worst_i_Linf = -1;
for(int i = 0; i < N; ++i)
{
Vec3 U = cases[i].U;
Vec3 K = cases[i].K;
double tau = cases[i].tau;
double kn = v_norm(K);
int degenerate = (kn <= 1e-30);
Vec3 khat = degenerate ? (Vec3){0,0,0} : v_scale(K, 1.0/kn);
/* Build quaternion (for reporting) from tau */
double h = hypot(1.0, tau);
double w = degenerate ? 1.0 : (1.0 / h);
double s = degenerate ? 0.0 : (tau / h);
Vec3 qv = degenerate ? (Vec3){0,0,0} : v_scale(khat, s);
double qnorm = sqrt(w*w + v_dot(qv, qv)); /* should be 1 */
Vec3 vq = rotate_axis_tau_quat(U, K, tau);
Vec3 vr = rodrigues_tau(U, K, tau);
Vec3 d = v_sub(vq, vr);
double L2 = v_norm(d);
double Linf = fmax(fabs(d.x), fmax(fabs(d.y), fabs(d.z)));
if(L2 > worst_L2) { worst_L2 = L2; worst_i_L2 = i; }
if(Linf > worst_Linf) { worst_Linf = Linf; worst_i_Linf = i; }
/* Invariants */
double normU = v_norm(U);
double normV = v_norm(vq);
double norm_err = fabs(normV - normU);
double parU = degenerate ? 0.0 : v_dot(khat, U);
double parV = degenerate ? 0.0 : v_dot(khat, vq);
double par_delta = degenerate ? 0.0 : (parV - parU);
double perpU = degenerate ? 0.0 : v_norm(v_perp(U, khat));
double perpV = degenerate ? 0.0 : v_norm(v_perp(vq, khat));
double perp_delta = degenerate ? 0.0 : (perpV - perpU);
int ok_core = vec3_close(vq, vr, ATOL, RTOL);
int ok_norm = (norm_err <= ATOL + RTOL * fmax(normU, normV));
int ok_par = degenerate || (fabs(par_delta) <= ATOL + RTOL * fmax(fabs(parU), fabs(parV)));
int ok_perp = degenerate || (fabs(perp_delta) <= ATOL + RTOL * fmax(perpU, perpV));
int ok = ok_core && ok_norm && ok_par && ok_perp;
if(ok) ++passed; else ++failed;
/* Detailed printout per test (ASCII only) */
printf("----- Test %02d: %s -----\n", i+1, cases[i].label ? cases[i].label : "(unnamed)");
printf("U = (%.16g, %.16g, %.16g) |U|=%.16g\n", U.x, U.y, U.z, normU);
printf("K = (%.16g, %.16g, %.16g) |K|=%.16g\n", K.x, K.y, K.z, kn);
printf("tau = %.16g (tau = tan(theta/2))\n", tau);
if(!degenerate)
{
printf("khat = (%.16g, %.16g, %.16g)\n", khat.x, khat.y, khat.z);
printf("quat = [w=%.16g, qv=(%.16g, %.16g, %.16g)] |q|=%.16g\n",
w, qv.x, qv.y, qv.z, qnorm);
}
else
{
printf("khat = (n/a, degenerate axis)\n");
printf("quat = identity (degenerate axis)\n");
}
printf("V_quat= (%.16g, %.16g, %.16g) |V_quat|=%.16g\n", vq.x, vq.y, vq.z, v_norm(vq));
printf("V_rod = (%.16g, %.16g, %.16g) |V_rod |=%.16g\n", vr.x, vr.y, vr.z, v_norm(vr));
printf("diff = (%.3e, %.3e, %.3e) L2=%.3e Linf=%.3e\n",
d.x, d.y, d.z, L2, Linf);
if(!degenerate)
{
printf("axis projection: (k.U)=%.16g, (k.V)=%.16g, delta=%.3e\n", parU, parV, par_delta);
printf("perp magnitudes: |U_perp|=%.16g, |V_perp|=%.16g, delta=%.3e\n", perpU, perpV, perp_delta);
}
printf("norm preservation: |V|-|U| = %.3e %s\n",
norm_err, ok_norm ? "[OK]" : "[WARN]");
printf("Result: %s (core=%s, axis-par=%s, perp=%s)\n\n",
ok ? "PASS" : "FAIL",
ok_core ? "OK" : "FAIL",
ok_par ? "OK" : "FAIL",
ok_perp ? "OK" : "FAIL");
}
printf("===== Summary =====\n");
printf("Total: %d Passed: %d Failed: %d\n", N, passed, failed);
if(worst_i_L2 >= 0) printf("Worst L2 : %.3e at test %02d\n", worst_L2, worst_i_L2+1);
if(worst_i_Linf >= 0) printf("Worst Linf : %.3e at test %02d\n", worst_Linf, worst_i_Linf+1);
return (failed == 0) ? 0 : 1;
}