偶的AS1 Matrix类

博主回顾去年五月为3D引擎编写类的经历,相关资料离校时遗失。在论坛看到以前发的版本也不确定。回家看图形学书,虽线代知识遗忘但后来适应。表示以后有时间会再搞3D,且已收集很多资料。

这是去年五月份为偶的3d引擎写的一个类.可惜那些东西离校的时候遗失在学校工作室了.在帝国论坛上看到以前发过的, 也不知道是哪一版本. 回家猛看了一下图形学的书, 发现线代已经扔给老师了(不应该说老师,偶从来没上过线代课),呵呵。不过后面还是适应了,毕竟线代学得还可以   以后有时间一定再搞一搞3d的,最近收集了好多资料

 

 

function Matrix() {
        //构造一个二维数组
        function constructor() {
                var i;
                t = new Array();
                for (i=0; i<4; i++) {
                        t[i] = new Array();
                }
                return t;
        }
        //公有数据成员
        this.mat = constructor();
        //私有成员函数, 计算矩阵的代数余子式
        function ValueDim(a) {
                var num1, num2;
                num1 = a[0][0]*a[1][1]*a[2][2]+a[0][1]*a[1][2]*a[2][0]+a[0][2]*a[1][0]*a[2][1];
                num2 = a[0][2]*a[1][1]*a[2][0]+a[0][1]*a[1][0]*a[2][2]+a[0][0]*a[1][2]*a[2][1];
                return (num1-num2);
        }
}
//清零
Matrix.prototype.ZeroMatrix = function() {
        var i, j;
        for (i=0; i<4; i++) {
                mat[i] = new Array();
                for (j=0; j<4; j++) {
                        mat[i][j] = 0;
                }
        }
};
//单位化
Matrix.prototype.LoadIndetity = function() {
        var i;
        ZeroMatrix();
        for (i=0; i<4; i++) {
                mat[i][i] = 1;
        }
};
//旋转 m=1时,绕x轴;m=2时,绕y轴;m=3时绕z轴
Matrix.prototype.Rotate3d = function(m, theta) {
        var m1, m2;
        var c, s;
        LoadIdentity();
        mat[m-1][m-1] = 1;
        mat[3][3] = 1;
        m1 = (m%3)+1;
        m2 = (m1%3)+1;
        m1 -= 1;
        c = Math.cos(theta);
        s = Math.sin(theta);
        mat[m1][m1] = c;
        mat[m1][m2] = -s;
        mat[m2][m2] = c;
        mat[m2][m1] = s;
};
//平移变换矩阵
//参数tx,ty,tz分别表示x,y,z位移量
//         1         0         0         0
//         0         1         0         0
//         0          0         1          0
//         tx         ty         tz         1
Matrix.prototype.Translate3d = function(tx, ty, tz) {
        LoadIdentity();
        mat[0][3] = tx;
        mat[1][3] = ty;
        mat[2][3] = tz;
};
//缩放变换矩阵
//sx,sy,sz分别表示沿x,y,z方向的缩放比例
//         sx         0         0         0
//         0         sy         0         0
//         0          0         sz         0
//         0         0         0   1
Matrix.prototype.Scaled3d = function(sx, sy, sz) {
        LoadIdentity();
        mat[0][0] = sx;
        mat[1][1] = sy;
        mat[2][2] = sz;
};
//去掉矩阵平移量,和透视变换量(以后还得定义一些透视函数得用到)
//                                    0
//                                    0
//                                     0
//         0         0         0
Matrix.prototype.RotComponet = function() {
        mat[0][3] = 0;
        mat[1][3] = 0;
        mat[2][3] = 0;
        mat[3][0] = 0;
        mat[3][1] = 0;
        mat[3][2] = 0;
};
//求逆矩阵,因为这里的是齐次坐标矩阵是四维的,
//所以求代数余子式ValueDim()比较好求
//假如是要求别的什么N维的话, 那就用别的方法。
Matrix.prototype.VertDim = function(b) {
        var i, j, lin, col, i1, j1;
        var d, deta1;
        var c = new Array();
        for (i=0; i<4; i++) {
                for (j=0; j<4; j++) {
                        lin = 0;
                        col = 0;
                        for (i1=0; i1<4; i1++) {
                                if (i1 != i) {
                                        c[lin] = new Array();
                                        for (j1=0; j1<4; j1++) {
                                                if (j1 != j) {
                                                        c[lin][col] = mat[i][i];
                                                        col += 1;
                                                }
                                        }
                                        lin += 1;
                                        col = 0;
                                }
                        }
                        deta1 = ValueDim(c);
                        if ((i+j)%2 == 0) {
                                b.mat[j][i] = deta1;
                        } else {
                                b.mat[j][i] = -deta1;
                        }
                }
        }
        d = 0;
        for (i=0; i<4; i++) {
                d += mat[0][i]*b.mat[i][0];
        }
        if (d == 0) {
                return;
        }
        for (i=0; i<4; i++) {
                for (j=0; j<4; j++) {
                        b.mat[i][j] /= d;
                }
        }
};
//求两个矩阵的乘积
Matrix.prototype.Matrix4x4 = function(v1, v2) {
        var i, j, k;
        for (i=0; i<4; i++) {
                for (j=0; j<4; j++) {
                        mat[i][j] = 0;
                        for (k=0; k<4; k++) {
                                mat[i][j] += v1.mat[i][k]*v2.mat[k][j];
                        }
                }
        }
};
//复制一个矩阵
Matrix.prototype.CopyMatrix = function(v1) {
        var i, j;
        for (i=0; i<4; i++) {
                for (j=0; j<4; j++) {
                        mat[i][j] = v1.mat[i][j];
                }
        }
};
//绕空间任意轴线放置变换矩阵
//这里用到了我还没有定义的一个类CVector,它是矢量类
//从CPointer类继承的,今天晚上再搞(反的工作顺序?)
//pbeg表示任意轴线的起点
//pend表示任意轴线的终点或者就是轴线的方向向量,这要看key的取值了
//key=0时pend表示终点
//key=1时pend就表示轴线方向向量(起点为默认为原点了)
Matrix.prototype.MakeRotateAxis = function(pbeg, pend, angle, key) {
        var r = 0, spsi, cpsi;
        var i;
        p = new Cvector();
        mat1 = new Matrix();
        ma = new Matrix();
        rx = new Matrix();
        ry = new Matrix();
        rz = new Matrix();
        rx1 = new Matrix();
        ry1 = new Matrix();
        mt1 = new Matrix();
        if (key != 1) {
                p.VectorPointMinus(pend, pbeg);
                p.Norvec();
        } else {
                p.Copy(pend);
                p.Norvec();
        }
        //平移矩阵
        //         1                 0                 0                 0
        //         0                 1                 0                 0
        //         0                  0                 1                  0
        //-pbeg[0] -pbeg[1] -pbeg[2] 1
        for (i=0; i<4; i++) {
                ma.mat[i][i] = 1;
        }
        ma.mat[0][3] = -pbeg[0];
        ma.mat[1][3] = -pbeg[1];
        ma.mat[2][3] = -pbeg[2];
        //逆平移矩阵
        //         1                 0                 0                 0
        //         0                 1                 0                 0
        //         0                  0                 1                  0
        // pbeg[0] pbeg[1]  pbeg[2]  1
        for (i=0; i<4; i++) {
                mt1.mat[i][i] = 1;
        }
        mt1.mat[0][3] = pbeg[0];
        mt1.mat[1][3] = pbeg[1];
        mt1.mat[2][3] = -pbeg[2];
        //绕x轴旋转矩阵
        //         1                 0                 0                 0
        //         0          cosθ        sinθ         0
        //         0                  -sinθ        cosθ         0
        //         0                 0                 0                 1
        spsi = 0;
        cpsi = 1;
        r = Math.sqrt(p.y*py+p.z*p.z);
        if (r>=1.e-5) {
                spsi = p.y/r;
                cpsi = p.z/r;
        } else {
                r = 0;
        }
        for (i=0; i<4; i++) {
                rx.mat[i][i] = 1;
        }
        rx.mat[1][1] = cpsi;
        rx.mat[1][2] = -spsi;
        rx.mat[2][1] = spsi;
        rx.mat[2][2] = cpsi;
        //绕x轴逆旋转矩阵
        //         1                 0                 0                 0
        //         0          cosθ        -sinθ         0
        //         0                  sinθ        cosθ         0
        //         0                 0                 0                 1
        for (i=0; i<4; i++) {
                rx1.mat[i][i] = 1;
        }
        rx1.mat[1][1] = cpsi;
        rx1.mat[1][2] = spsi;
        rx1.mat[2][1] = -spsi;
        rx1.mat[2][2] = cpsi;
        //绕y轴旋转矩阵
        //         cosθ         0                 -sinθ         0
        //         0          1           0                 0
        //         sinθ        0                cosθ         0
        //         0                 0                 0                 1
        for (i=0; i<4; i++) {
                ry.mat[i][i] = 1;
        }
        spsi = -p.x;
        cpsi = r;
        ry.mat[0][0] = cpsi;
        ry.mat[0][2] = -spsi;
        ry.mat[2][0] = spsi;
        ry.mat[2][2] = cpsi;
        //绕y轴逆旋转矩阵
        //         cosθ         0                 sinθ         0
        //         0          1           0                 0
        //         -sinθ        0                cosθ         0
        //         0                 0                 0                 1
        for (i=0; i<4; i++) {
                ry1.mat[i][i] = 1;
        }
        ry1.mat[0][0] = cpsi;
        ry1.mat[0][2] = spsi;
        ry1.mat[2][0] = -spsi;
        ry1.mat[2][2] = cpsi;
        //绕z轴旋转矩阵
        //         cosθ         sinθ 0                0
        //         -sinθ cosθ 0                 0
        //         0                 0     1                0
        //         0                 0          0                1
        for (i=0; i<4; i++) {
                rz.mat[i][i] = 1;
        }
        spsi = Math.sin(angle);
        cpsi = Math.cos(angle);
        rz.mat[0][0] = cpsi;
        rz.mat[0][1] = -spsi;
        rz.mat[1][0] = spsi;
        rz.mat[1][1] = cpsi;
        //mt1=ma·rx·ry·rz·rx1·ry1·mt1
        //最后得到的结果矩阵存在mat里
        mat1.Matrix4x4(mt1, rx1);
        mt1.Matrix4x4(mat1, ry1);
        mat1.Matrix4x4(mt1, rz);
        mt1.Matrix4x4(mat1, ry);
        mat1.Matrix4x4(mt1, rx);
        mt1.Matrix4x4(mat1, ma);
        CopyMatrix(mt1);
};
//以平面任意轴线对称变换矩阵
//这是一个二维矩阵, 主要是配合我要搞的教材,要不然就不会定义这个函数了
//pbeg,pend,key代表的意义和上面旋转函数一样
//结果放在mat里.不做什么解释了
Matrix.prototype.MakeReflectaxis = function(pbeg, pend, key) {
        var r = 0, spsi, cpsi;
        var i, j;
        p = new CVector();
        mat = new Matrix();
        ma = new Matrix();
        rx = new Matrix();
        rx1 = new Matrix();
        mt1 = new Matrix();
        if (key != 1) {
                p.VectorPointMinus(pend, pbeg);
                p.Norvec();
        } else {
                p.Copy(pend);
                p.Norvec();
        }
        for (i=0; i<4; i++) {
                ma.mat[i][i] = 1;
        }
        ma.mat[0][3] = -pbeg[0];
        ma.mat[1][3] = -pbeg[1];
        ma.mat[2][3] = -pbeg[2];
        for (i=0; i<4; i++) {
                mt1.mat[i][i] = 1;
        }
        mt1.mat[0][3] = pbeg[0];
        mt1.mat[1][3] = pbeg[1];
        mt1.mat[2][3] = pbeg[2];
        spsi = 0;
        cpsi = 1;
        r = Math.sqrt(p.x*p.x+p.y*p.y);
        if (r>1.e-5) {
                spsi = p.y/r;
                cpsi = p.x/r;
        } else {
                r = 0;
        }
        for (i=0; i<4; i++) {
                rx.mat[i][i] = 1;
        }
        rx.mat[0][0] = cpsi;
        rx.mat[0][1] = spsi;
        rx.mat[1][0] = -spsi;
        rx.mat[1][1] = cpsi;
        for (i=0; i<4; i++) {
                rx1.mat[i][i] = 1;
        }
        rx1.mat[0][0] = cpsi;
        rx1.mat[0][1] = -spsi;
        rx1.mat[1][0] = spsi;
        rx1.mat[1][1] = cpsi;
        for (i=0; i<4; i++) {
                rs.mat[i][i] = 1;
        }
        rs.mat[0][0] = 1;
        rs.mat[0][1] = 0;
        rs.mat[1][0] = 0;
        rs.mat[1][1] = -1;
        mat.Matrix4x4(mt1, rx1);
        mt1.Matrix4x4(mat, rs);
        mat.Matrix4x4(mt1, rx);
        mt1.Matrix4x4(mat, ma);
        CopyMatrix(mt1);
};
import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # 设置中文字体和负号显示 plt.rcParams['font.sans-serif'] = ['SimHei'] plt.rcParams['axes.unicode_minus'] = False # 物理常数 c = 3e8 # 光速 (m/s) def design_fda(N, M, f0, delta_f, r0, theta0, d=None): """ 设计频控阵并计算波束汇聚所需的相位补偿 """ # 1. 计算默认阵元间距 (半波长) lambda0 = c / f0 d = d if d is not None else lambda0 / 2 # 2. 生成阵元位置 [-Nd, ..., 0, ..., Nd] n_indices = np.arange(-N, N + 1) positions = n_indices * d # 3. 生成频率矩阵 (2N+1) x (M+1) m_range = np.arange(M + 1) n_abs = np.abs(n_indices) # Δf_n = Δf * log(|n| + 1) delta_f_n = delta_f * np.log(n_abs + 1 + 1e-12) # +1e-12 避免 log(0) # Δf_m = Δf * log(m+1) delta_f_m = delta_f * np.log(m_range + 1) # f_{n,m} = f0 + Δf_m + Δf_n freq_matrix = f0 + delta_f_n[:, np.newaxis] + delta_f_m[np.newaxis, :] # 4. 计算目标点直角坐标 target_x = r0 * np.cos(theta0) target_y = r0 * np.sin(theta0) # 5. 计算距离和传播时延 distances = np.sqrt((positions - target_x) ** 2 + target_y ** 2) time_delays = distances / c # 6. 计算相位补偿矩阵 (补偿传播相位) phase_comp = 2 * np.pi * freq_matrix * time_delays[:, np.newaxis] # 7. 添加额外补偿确保目标点同相 ref_phase = phase_comp[N, 0] # 中心阵元第一个载波的相位 phase_comp -= ref_phase # 使参考相位为零 return positions, freq_matrix, phase_comp def plot_beam_patterns(N, M, f0, delta_f, r0, theta0): """ 绘制频控阵的波束方向图(含二维俯视图) """ # 设计频控阵 positions, freq_matrix, phase_comp = design_fda(N, M, f0, delta_f, r0, theta0) # 创建计算网格 r_values = np.linspace(0.5 * r0, 2 * r0, 200) # 距离范围 theta_values = np.linspace(-np.pi / 2, np.pi / 2, 200) # 角度范围:-90°到90° R, Theta = np.meshgrid(r_values, theta_values) # 转换为直角坐标用于场强计算 X = R * np.cos(Theta) Y = R * np.sin(Theta) # 计算网格上每个点的场强 print("开始计算场强分布...") field_intensity = np.zeros_like(R, dtype=complex) num_elements = len(positions) num_freqs = freq_matrix.shape[1] wavelengths = c / freq_matrix # 预计算所有波长 # 优化计算:向量化操作 for i in range(num_elements): elem_x = positions[i] elem_y = 0 # 假设阵列沿x轴放置 # 预计算距离网格 dist_grid = np.sqrt((elem_x - X) ** 2 + (elem_y - Y) ** 2) for m in range(num_freqs): # 计算传播相位 prop_phase = 2 * np.pi * dist_grid / wavelengths[i, m] # 计算场强贡献 field_intensity += np.exp(1j * (phase_comp[i, m] - prop_phase)) # 取场强的模值 beam_pattern = np.abs(field_intensity) print("场强计算完成!") # 将角度转换为度 theta_deg = np.degrees(theta_values) r_km = r_values / 1000 # ================== 2D 距离和角度切面图 ================== plt.figure(figsize=(14, 6)) # 角度维切面 (固定距离r0) plt.subplot(1, 2, 1) r_idx = np.argmin(np.abs(r_values - r0)) plt.plot(theta_deg, beam_pattern[:, r_idx]) plt.axvline(np.degrees(theta0), color='r', linestyle='--', label='目标角度') plt.xlabel('角度 (度)') plt.ylabel('场强幅度') plt.title(f'固定距离 r={r0 / 1000:.1f}km 的角度维方向图') plt.legend() plt.grid(True) plt.xlim(-90, 90) # 距离维切面 (固定角度theta0) plt.subplot(1, 2, 2) angle_idx = np.argmin(np.abs(theta_values - theta0)) plt.plot(r_km, beam_pattern[angle_idx, :]) plt.axvline(r0 / 1000, color='r', linestyle='--', label='目标距离') plt.xlabel('距离 (km)') plt.ylabel('场强幅度') plt.title(f'固定角度 θ={np.degrees(theta0):.1f}° 的距离维方向图') plt.legend() plt.grid(True) plt.tight_layout() plt.show() # ================== 3D 场强距离角度图 ================== fig = plt.figure(figsize=(10, 8)) ax = fig.add_subplot(111, projection='3d') # 创建角度和距离网格 THETA, R_GRID = np.meshgrid(theta_deg, r_km) # 绘制3D曲面图 surf = ax.plot_surface(THETA, R_GRID, beam_pattern.T, cmap='viridis', edgecolor='none', alpha=0.8, antialiased=True) # 标记目标点 ax.scatter([np.degrees(theta0)], [r0 / 1000], [np.max(beam_pattern)], s=50, c='red', marker='o', label='目标点') ax.set_xlabel('角度 (度)') ax.set_ylabel('距离 (km)') ax.set_zlabel('场强幅度') ax.set_title('频控阵三维距离-角度场强分布') ax.legend() # 添加颜色条 fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10, label='场强幅度') # 设置视角 ax.view_init(elev=30, azim=45) plt.tight_layout() plt.show() # ================== 2D 俯视图 ================== plt.figure(figsize=(10, 8)) # 创建俯视图:角度-距离平面上的场强分布 plt.pcolormesh(theta_deg, r_km, beam_pattern.T, shading='gouraud', cmap='viridis') # 添加轮廓线 levels = np.linspace(0.2 * np.max(beam_pattern), np.max(beam_pattern), 10) plt.contour(theta_deg, r_km, beam_pattern.T, levels=levels, colors='white', linewidths=0.5, alpha=0.7) # 设置坐标轴和标签 plt.xlabel('角度 (度)') plt.ylabel('距离 (km)') plt.title('频控阵二维俯视场强分布') plt.colorbar(label='场强幅度') plt.legend() plt.grid(True, linestyle='--', alpha=0.5) # 设置坐标范围 plt.xlim(-90, 90) plt.ylim(0.5 * r0 / 1000, 2 * r0 / 1000) plt.tight_layout() plt.show() # 运行可视化 if __name__ == "__main__": # 参数设置优化 N = 5 M = 3 f0 = 2.4e9 # 10 GHz 中心频率 delta_f = 300e3 # 300 kHz 频率步进 # 目标点参数 r0 = 4000 # 1 km 距离 theta0 = np.pi / 4 # 30度方位角 print("开始计算波束方向图...") # 绘制方向图 plot_beam_patterns(N, M, f0, delta_f, r0, theta0) print("计算完成!")修改相位加权矩阵,使得只在目标位置形成峰值,不会在角度对称位置也形成一个峰值
最新发布
09-14
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值