高斯定理的证明（三重积分的C/C++实现）（C++）（大学物理）
2024-04-03 23:14:45
特别鸣谢：无私帮助我的ycc同学
#include <iostream>
#include <cstdlib>
#include <ctime>
#include <corecrt_math_defines.h>const int max_eletrons = 256;struct Eletrons {double x, y, z, q;
} eletrons[max_eletrons];
int num_eletrons = 0;
double rand_double() {return double(rand()) / 65536.0;
}
void init_eletrons(int range_min, int range_max) {srand((unsigned)time(0));num_eletrons = range_min + rand() % (range_max - range_min);for (int i = 0; i < num_eletrons; i++) {eletrons[i] = { rand_double(), rand_double(), rand_double(), rand_double() * 100.0 };}
}
struct Vec {double x, y, z;inline Vec operator + (Vec const& r) const { return { x + r.x, y + r.y, z + r.z }; }inline Vec operator - (Vec const& r) const { return { x - r.x, y - r.y, z - r.z }; }inline Vec operator * (double const& s) const { return { x * s, y * s, z * s }; }inline double dot(Vec const& r) const { return x * r.x + y * r.y + z * r.z; }inline double length2() const { return dot(*this); }inline double length() const { return std::sqrt(length2()); }inline double distance(Vec const& r) const { return (*this - r).length(); }inline Vec normalize() const { auto len = length(); return { x / len, y / len, z / len }; }
};const double pi = M_PI;
//const double eps0 = 8.854187817 * 1e-12f
const double eps0 = 1;
const double k = 1 / (4 * pi * eps0);double integral_dS(Vec const& pos, Vec const& normal) {double sum = 0.0;for (int i = 0; i < num_eletrons; ++i) {Eletrons e = eletrons[i];Vec dir = pos - Vec{ e.x, e.y, e.z };sum += dir.normalize().dot(normal) * k * e.q / dir.length2();}return sum;
}
struct Sphere {Vec p;double r;
};
bool sphere_inside(const Sphere& s, const Eletrons& e) {return s.p.distance(Vec{ e.x, e.y, e.z }) <= s.r;
}
double sphere_gauss(const Sphere& s) {double Q = 0;for (int i = 0; i < num_eletrons; ++i) {if (sphere_inside(s, eletrons[i])) Q += eletrons[i].q;}return Q / eps0;
}
double sphere_integral(const Sphere& s, int step) {double dphi = pi / step, dtheta = 2 * pi / step;double sum = 0.0;for (double phi = 0; phi < pi; phi += dphi)for (double theta = 0; theta < 2 * pi; theta += dtheta) {double dS = s.r * s.r * sin(phi) * dtheta * dphi;Vec norm = { sin(phi) * cos(theta), sin(phi) * sin(theta), cos(phi) };sum += integral_dS(s.p + norm * s.r, norm) * dS;}return sum;
}
struct Rectangle {Vec p, l;
};
bool rectangle_inside(const Rectangle& rec, const Eletrons& e) {return e.x >= rec.p.x - rec.l.x && e.x <= rec.p.x + rec.l.x &&e.y >= rec.p.y - rec.l.y && e.y <= rec.p.y + rec.l.y &&e.z >= rec.p.z - rec.l.z && e.z <= rec.p.z + rec.l.z;
}
double rectangle_gauss(const Rectangle& rec) {double Q = 0;for (int i = 0; i < num_eletrons; ++i) {if (rectangle_inside(rec, eletrons[i])) Q += eletrons[i].q;}return Q / eps0;
}
double rectangle_integral(const Rectangle& rec, int step) {double sum = 0.0, sum_xy = 0.0, sum_yz = 0.0, sum_xz = 0.0;double dx, dy, dz, dS;dx = 2 * rec.l.x / step, dy = 2 * rec.l.y / step, dS = dx * dy;sum_xy = 0.0;for (double x = rec.p.x - rec.l.x; x < rec.p.x + rec.l.x; x += dx)for (double y = rec.p.y - rec.l.y; y < rec.p.y + rec.l.y; y += dy) {sum_xy += integral_dS({ x + dx / 2, y + dy / 2, rec.p.z - rec.l.z }, { 0, 0, -1 });sum_xy += integral_dS({ x + dx / 2, y + dy / 2, rec.p.z + rec.l.z }, { 0, 0, 1 });}sum += sum_xy * dS;dy = 2 * rec.l.y / step, dz = 2 * rec.l.z / step, dS = dy * dz;sum_yz = 0.0;for (double y = rec.p.y - rec.l.y; y < rec.p.y + rec.l.y; y += dy)for (double z = rec.p.z - rec.l.z; z < rec.p.z + rec.l.z; z += dz) {sum_yz += integral_dS({ rec.p.x - rec.l.x, y + dy / 2, z + dz / 2 }, { -1, 0, 0 });sum_yz += integral_dS({ rec.p.x + rec.l.x, y + dy / 2, z + dz / 2 }, { 1, 0, 0 });}sum += sum_yz * dS;dx = 2 * rec.l.x / step, dz = 2 * rec.l.z / step, dS = dx * dz;sum_xz = 0.0;for (double x = rec.p.x - rec.l.x; x < rec.p.x + rec.l.x; x += dx)for (double z = rec.p.z - rec.l.z; z < rec.p.z + rec.l.z; z += dz) {sum_xz += integral_dS({ x + dx / 2, rec.p.y - rec.l.y, z + dz / 2 }, { 0, -1, 0 });sum_xz += integral_dS({ x + dx / 2, rec.p.y + rec.l.y, z + dz / 2 }, { 0, 1, 0 });}sum += sum_xz * dS;return sum;
}
struct Cylinder {Vec p;double r, h;
};
bool cylinder_inside(const Cylinder& clin, const Eletrons& e) {return e.z >= clin.p.z - clin.h / 2 && e.z <= clin.p.z + clin.h / 2 && clin.p.distance(Vec{ e.x,e.y,clin.p.z }) <= clin.r;//miaomiaomiao
}
double cylinder_gauss(const Cylinder& clin) {double Q = 0;for (int i = 0; i < num_eletrons; ++i) {if (cylinder_inside(clin, eletrons[i])) Q += eletrons[i].q;}return Q / eps0;
}
double cylinder_integral(const Cylinder& clin, int step) {double dz = clin.h / step, dtheta = 2 * pi / step, dr = clin.r / step;double sum = 0.0;for (double theta = 0; theta < 2 * pi; theta += dtheta){for (double r = 0; r < clin.r; r += dr){double ds = dtheta * dr * r;Vec norm1 = { 0, 0, 1 }, norm2 = { 0, 0, -1 };sum += integral_dS(clin.p + Vec{ r * cos(theta), r * sin(theta), clin.h / 2 }, norm1) * ds;sum += integral_dS(clin.p + Vec{ r * cos(theta),r * sin(theta), -clin.h / 2 }, norm2) * ds;}}for (double theta = 0; theta < 2 * pi; theta += dtheta) {for (double z = -clin.h / 2; z < clin.h / 2; z += dz) {double ds = clin.r * dtheta * dz;Vec norm = { cos(theta),sin(theta), 0 };sum += integral_dS(clin.p + Vec{ clin.r * cos(theta), clin.r * sin(theta), z }, norm) * ds;}}return sum;
}
struct Circular_cone {Vec p; //默认圆锥的顶点朝z轴正方向double r, h;
};
double circular_cone_inside(const Circular_cone& cir, const Eletrons& e) {double r0 = (cir.r / cir.h) * (cir.h - (e.z - cir.p.z));return e.z >= cir.p.z && e.z <= cir.p.z + cir.h && r0 >= 0 &&(e.x - cir.p.x) * (e.x - cir.p.x) + (e.y - cir.p.y) * (e.y - cir.p.y) <= r0 * r0;
}
double circular_cone_gauss(const Circular_cone& cir) {double Q = 0;for (int i = 0; i < num_eletrons; ++i) {if (circular_cone_inside(cir, eletrons[i])) Q += eletrons[i].q;}return Q / eps0;
}
double circular_cone_integral(const Circular_cone& cir, int step) {double dr = cir.r / step, dtheta = 2 * pi / step, dy = cir.h / step;double sum = 0.0;for (double r = 0; r < cir.r; r += dr) {for (double theta = 0; theta < 2 * pi; theta += dtheta) {double ds = r * dr * dtheta;Vec norm = { 0, 0, -1 };sum += integral_dS(cir.p + Vec{ r * cos(theta), r * sin(theta), 0 }, norm) * ds;}}double alpha = atan(cir.r / cir.h);for (double y = 0; y < cir.h; y += dy)for (double theta = 0; theta < 2 * pi; theta += dtheta) {double r0 = (cir.r / cir.h) * (cir.h - y - dy);double ds = dtheta * r0 * dy / cos(alpha);Vec norm = { cos(theta), sin(theta), cir.r / cir.h };sum += integral_dS(cir.p + Vec{ r0 * cos(theta), r0 * sin(theta), y }, norm.normalize()) * ds;}return sum;
}int main() {/*init_eletrons(5, 7);Sphere sphere = { {0, 0, 0}, 1 };printf("Sphere use Gauss: %lf\n", sphere_gauss(sphere));printf("Sphere use Integral: %lf\n", sphere_integral(sphere, 1000));Rectangle rectangle = { {0, 0, 0}, {0.9, 0.9, 0.9} };printf("Rectangle use Gauss: %lf\n", rectangle_gauss(rectangle));printf("Rectangle use Integral: %lf\n", rectangle_integral(rectangle, 1000));Cylinder cylinder = { {0,0,0},1000,2000 };printf("Cylinder use Gauss: %lf\n", cylinder_gauss(cylinder));printf("Cylinder use Integral: %lf\n", cylinder_integral(cylinder, 1000));Circular_cone cir = { {0, 0, -1}, 1, 2 };printf("Circular_cone use Gauss: %lf\n", circular_cone_gauss(cir));printf("Circular_cone use Integral: %lf\n", circular_cone_integral(cir, 2000));system("pause");return 0;*/printf("%d,%d,%d,",267.377556 - 267.030334, 237.123285 - 236.817932, 138.437950 - 138.261414);
}
学习经验（按照代码顺序）
- struct Vec{}“面向对象”
这个结构帮助实现多元向量值积分中的一些常用操作
1)inline函数&内联函数
关键字 inline 必须与函数定义体放在一起才能使函数成为内联，仅将 inline 放在函数声明前面不起任何作用。
tip:当函数只有10行或者更少时将其定义为内联函数
内联函数让代码更简单，好用get；）
2）std::
是标准库的前缀，C++要用标准库函数就要用std::
主要是起到了资源管理的作用，下面是一个例子：
有两个软件公司A公司和B公司，他们都是用C++语言开发他们的产品。那么，他们分别编写了a.h和b.h两个自己的头文件，这两个文件中都有一个叫func()的函数。他们各自使用也没什么问题。
假设你公司也是一个软件公司，你现在要开发一个软件，必须同时用到A公司和B公司的头文件，同时会调用他们的func（）函数。这个时候问题就来了，你调用的func()函数，编译器不知道应该选用A公司的还是B公司的。为解决这个问题，C++采用了命名空间，这样，你调用A公司的func()函数，就使用A::func()，B公司亦然。
原文链接：https://blog.csdn.net/u013152895/article/details/44490923
- init_eletrons()函数
高斯定理的证明（三重积分的C/C++实现）（C++）（大学物理）相关推荐
- 电磁学-高斯定理的证明和应用
## **高斯定理的理解** ## 电场线密度:该点的电场强度,即E的大小 ## 电场线性质: 1.起源于正电荷,止于负电荷,不会中断 正电荷是静电场的源头,电荷是静电场的尾闾 所以说 2.电场线能形 ...
- 有介质的高斯定理详细证明（电偶极子模型）以及例题讲解
目录 静电场中的电介质 电极化强度的引入 电偶极子模型的计算 电介质极化过程 极化电荷引入 推导 各向同性和线性的电介质 例题 静电场中的电介质 电介质与导体的区别:所有的粒子被束缚在原子核周围(限制 ...
- 【大学物理·恒定电流的磁场】恒定磁场的高斯定理与安培环路定理
一.恒定磁场的高斯定理 通过任一闭合曲面的总磁通量总等于零 激发静电场的场源(电荷)足电场线的源头或尾闾,所以静电场 是属于发散状的场,可称作有源场:而磁场的磁感应线无头无尾,恒是闭合的,所以磁场可称 ...
- 【大学物理· 恒定电流的磁场】有磁介质时的安培环路定理和高斯定理 磁场强度
一.磁化强度 磁化强度:表征磁介质磁化的程度 设为圆柱形磁介质表面上单位长度的磁化面电流,即磁化面电流的线密 度,S为磁介质的截面积,l为所选取的一段磁介质的长度 M为单位体积内的磁矩 磁介质表而某处 ...
- 【大学物理·静止电荷的电场】静电场的高斯定理
一.静电场的高斯定理 如果闭合曲面内没有包围电荷,电荷在闭合曲面外面,那么进人 闭合曲面的电场线等于穿出该闭合曲面的电场线,所以总的E通量为零 (q在闭合曲面的外面) 高斯定理 在静电场中,通过任一闭 ...
- 【大学物理·静止电荷的电场】有电介质时的高斯定理和环路定理 电位移
一.有电介质时的高斯定理 电位移 电位移 有电介质时的高斯定理: 在有电介质的 静电场中作电位移线,使线上每一点的切线方向和该点电位移D的方向相同, 并规定在垂直于电位移线的单位面积上通过的电位移线数 ...
- 【大学物理】磁场的高斯定理
- 大学物理----旋转矢量法证明同方向同频率简谐振动的合运动公式
- c语言高斯定理图片,大学物理上册所有公式-20210516063148.docx-原创力文档
第一章质点运动学和牛顿运动定律 1.1 平均速度v = △r t △r dr 1.2瞬时速度 v= lim△t= dt △t 0 △r lim ds 1. 3 速度 v= lim dt △t 0 △t ...
- 知识见解关于高斯定理
知识见解 电磁学中的高斯定理 高斯定理的定义 指真空中的任何静电场中,穿过任一闭合曲面的电通量,在数值上等于该闭合曲面内包围的电量的代数和乘以1/ε. 首先学习一个新知识,必须牢记并理解其概念,牢记的 ...
最新文章
- 双链表（删除节点操作）
- R语言ggplot2可视化：置信区间与分组具有相同色彩、自定义置信区间带的色彩、Make confidence intervals the same color as line by group
- Unity3d疑难问题解决
- R开发(part5)--导数计算
- 掌握 React 与 React Native
- 轻松使用终端开启macOS系统的隐藏功能，小白都能看得懂
- 亮屏变“黄”，暗屏变“绿”，iPhone 12用户太难了
- Forrester报告：人工智能将取代6%的工作岗位
- 修改windows默认远程管理端口
- Java使用BufferedImage修改图片内容
- 银行突发事件演练方案_湘阴星龙村镇银行开展防抢劫应急预案实战演练
- linux卸载java rpm_详解Linux中查看jdk安装目录、Linux卸载jdk、rpm命令、rm命令参数...
- Redis基础--使用treeNMS管理及监控Redis
- win7桌面我的计算机打不开怎么回事,win7系统双击我的电脑打不开的解决方法
- 让传感器数据在三维地图上显示，更直观，更震撼！
- Windows如何编辑hosts
- MT6797处理器怎么样?Helio X20处理器资料介绍
- How browsers work翻译
- 【真相】ChatGPT和OpenAI的API KEY
- 【python文件读取】加密数据的读取