使用两节点四自由度梁单元模型，结合瑞利阻尼模型，对两端固定支撑的圆截面梁的结构动力学问题进行求解，要求使用C语言编程

以下是使用两节点四自由度梁单元模型，结合瑞利阻尼模型，对两端固定支撑的圆截面梁的结构动力学问题进行求解的C语言程序示例：

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

#define N 2    //每个单元节点数
#define M 10   //梁的数量
#define DOF 4  //每个节点的自由度数

void assemble(double *K, double *M, double *C, double *f, double *u, double *v, double *a, double t);
void solve(double *K, double *M, double *C, double *f, double *u, double *v, double *a, double dt, double t, double tf);

int main()
{
    //定义节点坐标
    double x[N * M] = {0.0, 0.0, 0.0, 0.0, 0.1, 0.0, 0.2, 0.0, 0.3, 0.0, 0.4, 0.0, 0.5, 0.0, 0.6, 0.0, 0.7, 0.0, 0.8, 0.0, 0.9, 0.0, 1.0, 0.0};
    //定义单元节点编号
    int e[N * (M - 1)] = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10};
    //定义杨氏模量、密度、阻尼系数和损耗系数
    double E = 2.1e11, rho = 7800, alpha = 1.0e-3, beta = 1.0e-3;
    //定义节点载荷
    double f[N * M * DOF] = {0.0};
    //定义初始位移、速度和加速度
    double u[N * M * DOF] = {0.0};
    double v[N * M * DOF] = {0.0};
    double a[N * M * DOF] = {0.0};

    //定义时间参数
    double t = 0.0, tf = 1.0, dt = 0.001;

    //定义刚度矩阵、质量矩阵、阻尼矩阵和瑞利阻尼矩阵
    double K[N * M * DOF][N * M * DOF] = {0.0};
    double M[N * M * DOF][N * M * DOF] = {0.0};
    double C[N * M * DOF][N * M * DOF] = {0.0};
    double R[N * M * DOF][N * M * DOF] = {0.0};

    //组装刚度矩阵、质量矩阵、阻尼矩阵和瑞利阻尼矩阵
    assemble(&K[0][0], &M[0][0], &C[0][0], &f[0], &u[0], &v[0], &a[0], t);
    //求解结构动力学问题
    solve(&K[0][0], &M[0][0], &C[0][0], &f[0], &u[0], &v[0], &a[0], dt, t, tf);

    return 0;
}

void assemble(double *K, double *M, double *C, double *f, double *u, double *v, double *a, double t)
{
    //定义局部刚度矩阵、质量矩阵和阻尼矩阵
    double k[N * DOF][N * DOF] = {0.0};
    double m[N * DOF][N * DOF] = {0.0};
    double c[N * DOF][N * DOF] = {0.0};
    double R[N * DOF][N * DOF] = {0.0};
    double L, c1, c2, s1, s2, r;

    for (int i = 0; i < M; i++)
    {
        //计算单元长度
        L = sqrt(pow(x[N * i + 2] - x[N * i], 2.0) + pow(x[N * i + 3] - x[N * i + 1], 2.0));
        //计算单元方向余弦
        c1 = (x[N * i + 2] - x[N * i]) / L;
        s1 = (x[N * i + 3] - x[N * i + 1]) / L;
        c2 = (x[N * i + 2] - x[N * i]) / L;
        s2 = (x[N * i + 3] - x[N * i + 1]) / L;

        //计算局部刚度矩阵、质量矩阵和阻尼矩阵
        k[0][0] = k[1][1] = E * 3.1416 * pow(L, 3.0) / (4 * pow(L / 2, 4.0));
        k[0][1] = k[1][0] = -E * 3.1416 * pow(L, 3.0) / (4 * pow(L / 2, 4.0));
        k[2][2] = k[3][3] = E * 3.1416 * pow(L / 2, 2.0) / L;
        k[0][2] = k[2][0] = -E * 3.1416 * pow(L, 3.0) / (8 * pow(L / 2, 3.0));
        k[1][3] = k[3][1] = E * 3.1416 * pow(L, 3.0) / (8 * pow(L / 2, 3.0));
        k[2][3] = k[3][2] = -E * 3.1416 * pow(L, 2.0) / (4 * (L / 2));
        k[0][3] = k[3][0] = -E * 3.1416 * pow(L, 2.0) / (4 * (L / 2));
        k[1][2] = k[2][1] = E * 3.1416 * pow(L, 2.0) / (4 * (L / 2));
        k[1][3] = k[3][1] = E * 3.1416 * pow(L, 2.0) / (4 * (L / 2));

        m[0][0] = m[1][1] = rho * 3.1416 * pow(L / 2, 2.0) * L / 2;
        m[2][2] = m[3][3] = rho * 3.1416 * pow(L / 2, 4.0) / 4;
        m[0][2] = m[2][0] = rho * 3.1416 * pow(L / 2, 3.0) / 6;
        m[1][3] = m[3][1] = rho * 3.1416 * pow(L / 2, 3.0) / 6;
        m[2][3] = m[3][2] = rho * 3.1416 * pow(L / 2, 3.0) / 3;

        c[0][0] = c[1][1] = alpha * k[0][0] + beta * m[0][0];
        c[0][1] = c[1][0] = alpha * k[0][1] + beta * m[0][1];
        c[2][2] = c[3][3] = alpha * k[2][2] + beta * m[2][2];
        c[2][3] = c[3][2] = alpha * k[2][3] + beta * m[2][3];
        c[0][2] = c[2][0] = alpha * k[0][2] + beta * m[0][2];
        c[1][3] = c[3][1] = alpha * k[1][3] + beta * m[1][3];
        c[1][2] = c[2][1] = alpha * k[1][2] + beta * m[1][2];
        c[0][3] = c[3][0] = alpha * k[0][3] + beta * m[0][3];

        r = 0.1 * (k[0][0] + k[1][1] + k[2][2] + k[3][3]) / 4;

        R[0