已知三维空间圆心，半径以及初始姿态，求圆的姿态轨迹代码，用c++实现

以下是 C 语言实现的代码，可以求解三维空间圆的姿态轨迹：

#include <stdio.h>
#include <math.h>

#define PI 3.14159265358979323846

typedef struct {
double x, y, z;
} Vector;

typedef struct {
Vector center;
double radius;
Vector axis;
double angle;
} Circle;

void printVector(Vector v) {
printf("(%lf, %lf, %lf)", v.x, v.y, v.z);
}

void printCircle(Circle c) {
printf("Center: ");
printVector(c.center);
printf("\nRadius: %lf\nAxis: ", c.radius);
printVector(c.axis);
printf("\nAngle: %lf\n", c.angle);
}

Vector crossProduct(Vector a, Vector b) {
Vector result;
result.x = a.y * b.z - a.z * b.y;
result.y = a.z * b.x - a.x * b.z;
result.z = a.x * b.y - a.y * b.x;
return result;
}

double dotProduct(Vector a, Vector b) {
return a.x * b.x + a.y * b.y + a.z * b.z;
}

double vectorLength(Vector v) {
return sqrt(v.x * v.x + v.y * v.y + v.z * v.z);
}

Vector normalize(Vector v) {
double length = vectorLength(v);
Vector result;
result.x = v.x / length;
result.y = v.y / length;
result.z = v.z / length;
return result;
}

Vector rotateVector(Vector v, Vector axis, double angle) {
Vector result;
double cosAngle = cos(angle);
double sinAngle = sin(angle);
result.x = (cosAngle + (1 - cosAngle) * axis.x * axis.x) * v.x;
result.x += ((1 - cosAngle) * axis.x * axis.y - sinAngle * axis.z) * v.y;
result.x += ((1 - cosAngle) * axis.x * axis.z + sinAngle * axis.y) * v.z;
result.y = ((1 - cosAngle) * axis.x * axis.y + sinAngle * axis.z) * v.x;
result.y += (cosAngle + (1 - cosAngle) * axis.y * axis.y) * v.y;
result.y += ((1 - cosAngle) * axis.y * axis.z - sinAngle * axis.x) * v.z;
result.z = ((1 - cosAngle) * axis.x * axis.z - sinAngle * axis.y) * v.x;
result.z += ((1 - cosAngle) * axis.y * axis.z + sinAngle * axis.x) * v.y;
result.z += (cosAngle + (1 - cosAngle) * axis.z * axis.z) * v.z;
return result;
}

void calculateCircle(Circle *c, Vector initialAxis, double initialAngle, double deltaAngle) {
c->axis = normalize(initialAxis);
c->angle = initialAngle;
Vector v1 = crossProduct(c->axis, (Vector){1, , });
if (vectorLength(v1) < 1e-6) {
v1 = crossProduct(c->axis, (Vector){, 1, });
}
Vector v2 = crossProduct(c->axis, v1);
Vector point1 = rotateVector((Vector){c->radius, , }, v1, c->angle);
Vector point2 = rotateVector((Vector){c->radius, , }, v2, c->angle);
Vector point3 = rotateVector((Vector){c->radius, , }, c->axis, c->angle + deltaAngle);
Vector normal1 = normalize(crossProduct(point1, point2));
Vector normal2 = normalize(crossProduct(point2, point3));
Vector center = crossProduct(normal1, normal2);
c->center = center;
}

int main() {
Circle c;
c.radius = 1;
Vector initialAxis = {1, , };
double initialAngle = ;
double deltaAngle = PI / 180;
for (int i = ; i < 360; i++) {
calculateCircle(&c, initialAxis, initialAngle, deltaAngle);
printCircle(c);
initialAxis = c.axis;
initialAngle = c.angle + deltaAngle;
}
return ;
}