DHCP协议分析

发表于 2024-01-14

协议流程

直接DHCP

DHCP允许系统管理员在服务器上集中定义设备配置，当这些设备启动时，可以请求这些配置参数。

当客户端发送广播DISCOVER包时，DHCP 进程开始，本质上是说“外面有任何 DHCP 服务器吗?” 这就是我正在寻找的信息。” 然后 DHCP 服务器回复一个包含所有信息的OFFER报文。客户机用一个REQUEST包进行响应，这个包的名称似乎很奇怪——实际上，客户机只是通过确认的方式，发送它刚刚从服务器返回的信息。然后，服务器发送最后的ACKNOWLEDGEMENT数据包，同样带有相同的信息，再次确认它。

这是通常称为DORA序列(发现、提供、请求、确认)，通常是这样描述的:

因为这些都是 UDP 数据包，请记住 UDP 协议中没有任何会话信息，那么是什么将这四个数据包绑定到一个“会话”中呢? 因此，最初的 Discover 报文有一个事务 ID，在随后的三个报文中匹配- Wireshark 跟踪如下所示:

实际上，客户端直到第四个包才有一个地址，所以 Discover 和 Request 包是从 IP 为0.0.0.0的客户端的 MAC 地址到255.255.255.255的广播地址(即到整个局域网)。

中继DHCP

在许多公司网络中，服务器在它们自己的子网中——分离服务器和工作站是相当普遍的做法。这种情况下 DHCP 顺序如何工作? DORA 序列的前三个数据包被发送到广播地址，因此它们只能到达同一 VLAN 上的其他主机。

我们通过在客户端子网中的主机上放置一个 DHCP“转发器”或“中继”进程来完成这项工作。该进程接收本地广播，然后将其以单播形式转发给 DHCP 服务器。当服务器应答时(以单播方式向转发器主机)，转发器将包“转换”为客户机所期望的广播应答。几乎总是，这个转发器功能是在客户端子网上的路由器或交换机 IP 地址上完成的——换句话说，接口最终将成为客户端的默认网关。这个函数在技术上不需要在那个接口上，但它是一个我们知道会在那里的接口，而且这个函数几乎总是可供我们使用。另外，如果我们将其作为一种不成文的惯例，当我们以后需要更改它时，它将更容易找到该命令! 在 Cisco 路由器或交换机上，这个命令看起来像这样:

1	interface VLAN <x> ip helper-address 10.10.10.10

这里，10.10.10.10是我们 DHCP 服务器的 IP 地址。

在操作中，这改变了我们在大多数家庭网络上拥有的简单的广播操作：这如何修改我们的 DORA 序列? 简单的回答是，它实际上不会修改任何数据包的 DHCP 内容。它所做的是修改数据包中的上层“IP 地址”字段-路由器和服务器之间修改的数据包有“真实的”源和目的 IP 地址。但是，客户端看到的包内容保持不变。如果你深入研究 DHCP 报文，你会发现不管是否有中继，DHCP 客户端的 MAC 地址和 DHCP 服务器的 IP 地址实际上包含在 7 层 DHCP 协议的数据字段中。

现在我们开始为基本配置 DHCP 服务器工作站操作,但在我们到达之前,我们要考虑我们需要的专用设备,如 iphone,无线访问点(WAP), 或者甚至是预执行环境(PXE)设备，可以从 DHCP 信息加载其整个操作系统。

仿真

组网

路由设置

<Huawei>sys
[Huawei]sysname AR1
[AR1]int g0/0/0
[AR1-GigabitEthernet0/0/0]ip adderss 1.1.1.200 24
[AR1]dhcp enable
[AR1]ip pool ip_pool1
[AR1-ip-pool-ip_pool1]network 1.1.1.0 24
[AR1-ip-pool-ip_pool1]gateway-list 1.1.1.200
[AR1-ip-pool-ip_pool1]lease day 3
[AR1]int g0/0/0
[AR1-GigabitEthernet0/0/0]dhcp select global

抓包

dhcp discover

dhcp offer

资源收集

发表于 2024-01-13

malloc: http://blog.codinglabs.org/articles/a-malloc-tutorial.html

duartes.org: http://duartes.org/gustavo/blog/archives/

github mannul: http://www.epubit.com.cn/article/844#what

liaoxuefeng.com: http://www.liaoxuefeng.com/

TCP/IP network: http://blog.packagecloud.io/eng/2016/10/11/monitoring-tuning-linux-networking-stack-receiving-data-illustrated/

netfilter.org: https://people.netfilter.org/pablo/netdev0.1/papers/

tuning-linux-sending: https://blog.packagecloud.io/eng/2017/02/06/monitoring-tuning-linux-networking-stack-sending-data/

tuning-linux-receiving: https://blog.packagecloud.io/eng/2016/06/22/monitoring-tuning-linux-networking-stack-receiving-data/

strace: https://blog.packagecloud.io/eng/2016/02/29/how-does-strace-work/

IBM-Tim: https://www.ibm.com/developerworks/cn/views/linux/libraryview.jsp?search_by=Linux+%E5%89%96%E6%9E%90

yeolar: http://www.yeolar.com/

http://www.yeolar.com/note/2012/03/29/virtual-memory/

vxlan: https://blogs.vmware.com/vsphere/2013/07/vxlan-series-how-vmotion-impacts-the-forwarding-table-part-6.html

http://www.cisco.com/c/en/us/products/collateral/switches/nexus-5000-series-switches/white-paper-c11-733618.html#_Toc439799767

http://www.cisco.com/c/en/us/products/collateral/switches/nexus-9000-series-switches/white-paper-c11-729383.html

ali-kernel: http://kernel.taobao.org/index.php?title=%E5%86%85%E6%A0%B8%E6%9C%88%E6%8A%A52017-02

linux-performance: http://www.brendangregg.com/linuxperf.html

systemTap: https://sourceware.org/systemtap/SystemTap_Beginners_Guide/

calico: http://docs.projectcalico.org/v2.0/introduction/

docker: http://www.infoq.com/cn/articles/docker-network-and-pipework-open-source-explanation-practice

Flannel: http://dockone.io/article/618

SDN: https://www.opennetworking.org/

Linux Kernel Networking: https://wiki.linuxfoundation.org/networking/start

MacVtap: https://blog.kghost.info/2013/03/27/linux-network-tun/

http://blog.csdn.net/dog250/article/details/45788279

Neutron: http://blog.csdn.net/quqi99/article/details/22853403

Linux Bridge: http://blog.csdn.net/yeasy/article/details/50728243

Docker Networking: http://edgedef.com/docker-networking.html

Tun/Tap interface: http://backreference.org/2010/03/26/tuntap-interface-tutorial/

TC ifb: http://blog.csdn.net/dog250/article/details/40680765?utm_source=tuicool&utm_medium=referral

汇编:

内核转储文件--内核崩溃

发表于 2024-01-13

Ubuntu

安装crash工具

1	apt-get install linux-crashdump

检查kdump配置

kdump-config show

DUMP_MODE:              kdump
USE_KDUMP:              1
KDUMP_COREDIR:          /var/crash
crashkernel addr: 0x73000000
   /var/lib/kdump/vmlinuz: symbolic link to /boot/vmlinuz-5.15.0-89-generic
kdump initrd:
   /var/lib/kdump/initrd.img: symbolic link to /var/lib/kdump/initrd.img-5.15.0-89-generic
current state:    ready to kdump

kexec command:
  /sbin/kexec -p --command-line="BOOT_IMAGE=/boot/vmlinuz-5.15.0-89-generic root=UUID=92414257-97c5-46a0-9154-66c415ee7358 ro net.ifnames=0 consoleblank=600 console=tty0 console=ttyS0,115200n8 noibrs reset_devices systemd.unit=kdump-tools-dump.service nr_cpus=1 irqpoll nousb" --initrd=/var/lib/kdump/initrd.img /var/lib/kdump/vmlinuz

USE_KDUMP=1 代表kdump打开

KDUMP_COREDIR 代表生成的core文件在/var/crash下面

-> # dmesg | grep -i crash
[    0.000000] Command line: BOOT_IMAGE=/boot/vmlinuz-5.15.0-89-generic root=UUID=92414257-97c5-46a0-9154-66c415ee7358 ro net.ifnames=0 consoleblank=600 console=tty0 console=ttyS0,115200n8 noibrs crashkernel=512M-:192M
[    0.005698] Reserving 192MB of memory at 1840MB for crashkernel (System RAM: 2047MB)
[    0.013125] Kernel command line: BOOT_IMAGE=/boot/vmlinuz-5.15.0-89-generic root=UUID=92414257-97c5-46a0-9154-66c415ee7358 ro net.ifnames=0 consoleblank=600 console=tty0 console=ttyS0,115200n8 noibrs crashkernel=512M-:192M

dmesg 显示了crashkernel保留了192M内存，kexec命令。

手动触发sysrq

1 2	cat /proc/sys/kernel/sysrq echo c > /proc/sysrq-trigger

sysrq值不为0代表工作正常
echo写主动触发kdump

CentOS

安装crash工具

1 2	yum install kexec-tools yum install crash

安装内核调试信息包

1 2	yum install kernel-debuginfo-$(uname -r) yum install kernel-debuginfo-common-$(uname -r)

判断服务正常开启

1	systemctl status kdump

主动触发分析

这一步跟Ubuntu的一致

分析调试

1	crash /usr/lib/debug/lib/modules/$(uname -r)/vmlinux /var/crash/XX/vmcore

内核转储文件--应用程序崩溃分析

发表于 2024-01-07

Ubuntu

设置如下的程序

#include <iostream>
#include <string.h>
using namespace std;

void crash_test()
{
        char *str =0;
        strcpy(str,"test"); // segment fault
}
int main()
{
        cout << "crash test" << endl;
        crash_test();
        return 0;
}

以下就是生成了dump文件

查看coredump是否打开

查看命令：ulimit -c

打开命令：ulimit -c unlimited

0：关闭

unlimited：打开 ### 配置规则主要看内核参数kernel.core_pattern，路径为/proc/sys/kernel/core_pattern。

所以修改有两种方法

1、sysctl -w kernel.core_pattern=/root/core/core-%e.%s.%p.%t

2、echo -e "/root/core/core-%e.%s.%p.%t" > /proc/sys/kernel/core_pattern

%e：进程名称

%s：崩溃时收到的信号

%p：进程ID

%t：时间戳

注意关闭apport服务

这个服务会读取core文件分析，不会保存源文件

1	systemctl stop apport.service

CentOS

编辑`/etc/security/limits.conf`文件

末尾增加soft core unlimited，表示应用程序生成的core文件不受大小限制 ### 配置规则和Ubuntu一样 ### reboot

分析崩溃

1	gdb test /root/corefile/core-test.11.6284.1705062175

发现最后的出错函数位于strcpy这一处，源文件第8行。

Obsidian和github和hexo博客流配置

发表于 2023-12-21

hexo配置

node.js配置

hexo主题配置

Obsidian配置

下载Obsidian Git插件，需要打开第三方插件按钮，下载插件。

在对应的.git目录下打开仓库，下图的打开本地仓库，生成本地目录。

配置Git设置地址等，这样就可以调用后台的git自动推送了

需要注意的是：git需要首先被安装，并且配置到环境变量

github配置

workflows配置

dependabot配置

性能分析-Windows（一）

发表于 2020-10-20

Windows运行参数分析

这篇文章较为简单，主要是使用 Windows Performance Record 和 Windows Performance Analysis 来分析在运行一段时间内的各个参数的变化，记录详细的信息。

一般大家都是从任务管理器上面来看机器当前的运行状态，但是，如果我想要了解的更详细一些，并且需要一段特定时间的运行参数，任务管理器仿佛就无能为力了，这时候，可以使用 Windows 测试框架里面的两个小工具来记录分析机器的运行状态。

阅读全文 »

正则表达式(一)

发表于 2020-03-17

正则表达式(一)

元字符

元字符（metacharacters）是按照一定规则限制表示的字符。

.匹配所有字符
[]匹配方括号内的所有字符,这个字符按照原样,即使是具有特殊意义的字符,如’,$,>,+’等都是保持字符本意\-除外
[^c1c2]c3c4不匹配’c1,c2’,即不带有c1,c2的c3c4字符串
[c1-c2]匹配c1到c2之间的所有字符
\w匹配任意单个数字，字符串或者下划线，等价于[a-zA-Z_1-9]，\w*匹配多个连续字符\w{N}匹配N个连续字符

1
2
3

OriginStr = 'spain remain contain aint retain ';
AnsStr = regexp(OriginStr,'\w*','match');
PosStr = regexp(OriginStr,'\w*');

PosStr =

     1     7    14    22    27    33

>> AnsStr

AnsStr =

  1×5 cell 数组

    {'spain'}    {'remain'}    {'contain'}    {'aint'}    {'retain'}

可以看出，这个得到的是所有的匹配。而且对于语句regexp(OriginStr,MatchStr,'option')，如果不加option的话，就是得到匹配的一个字符的位置，如果加上option为match的话，就是得到匹配的值。

1 2	OriginStr = 'spain remain contain aint retain '; AnsStr = regexp(OriginStr,'\w* ain','match');

>> AnsStr

AnsStr =

  1×1 cell 数组

    {'contain ain'}

可以看到，AnsStr输出的contain ain，即匹配ain的前缀都包括的一串字符串。

\W是匹配\w的补集，\W*跟\w*含义一致\W{N}和\w{N}含义一致。
\s等价于[ \f\t\r\n\v]，匹配空白字符
\S等价于[^\f\t\r\n\v]匹配非空白字符 tips：小写字母和大写字母的含义往往是补集的关系
\d匹配单个十进制数\d*匹配多个 tips*匹配多个
D是\d的补集，即匹配除了数字以外的所有字符，等价于[^\d]
\oN 和\o{N}匹配八进制数N
\xN匹配\x{N}匹配十六进制数N

表达式

expr，表达式，这个也就是对于一串字符串，即xp=x&p，[xp]=x|p，即不加[]是且关系，加了[]是或关系，表达式是限定符前面构成的所有字符构成的匹配字符串。

转义字符

char(10)=\n，即回车。

1
2
3

str = ['some',char(10),'text']
regexp(str,'\n')
regexp(str,'\n','match')

输出为

str =
 'some
  text'
ans =
  5
ans =
  1×1 cell 数组
 {'?'}

转义字符	描述	示例
`\a`	警告	`char(7)`
`\b`	退格	`char(8)`
`\f`	换页符	`char(12)`
`\n`	换行符	`char(10)`
`\r`	回车	`char(13)`
`\t`	水平制表符	`char(9)`
`\v`	垂直制表符	`char(11)`
`\char`	任意特殊字符	-

限定符

限定符	含义
expr *	多次expr：0~∞
expr ?	一次expr：0~1
expr +	一次到多次expr：1~∞
expr{m,n}	m到n次expr：m~n
expr{m,}	大于等于m次expr：m~∞
expr{m}	m次expr：=m

expr*代表匹配0次或者多次，与任意多个字符匹配

>> regexp('12--asdf--45dfg','\w*','match')

ans =

  1×3 cell 数组

    {'12'}    {'asdf'}    {'45dfg'}

expr?代表至多匹配1次

>> regexp('12--asdf--45dfg','\w?','match')

ans =

  1×11 cell 数组

    {'1'}    {'2'}    {'a'}    {'s'}    {'d'}    {'f'}    {'4'}    {'5'}    {'d'}    {'f'}    {'g'}

expr+匹配一次或者连续多次

>> regexp('12--asdf--45dfg','\w+','match')

ans =

  1×3 cell 数组

    {'12'}    {'asdf'}    {'45dfg'}

expr{m,n}匹配m~n次，{0,1}='?'

>> regexp('12--asdf--45dfg','\w{1,3}','match')

ans =

  1×5 cell 数组

    {'12'}    {'asd'}    {'f'}    {'45d'}    {'fg'}

expr{m,}匹配≥m次，{0,}='*';{1,}='+'
expr{n}匹配n次，n={n,n}

模式

表达式+限定词=模式

表达式+限定词	描述	示例
expr q	积极表达式：与尽可能多的expr匹配	给定文本 `<tr><td><p>text</p></td>`，表达式 `</?t.*>`与介于 `<tr` 和 `/td>` 之间的所有字符匹配：
`<tr><td><p>text</p></td>`
expr q？	消极表达式：与所需尽可能少的expr匹配	给定文本`<tr><td><p>text</p></td>`，表达式 `</?t.*?>` 在第一次出现右尖括号 (>) 时结束每个匹配项：
`<tr>` `<td>` `</td>`
expr q+	主动表达式：最大程度的匹配	给定文本 `<tr><td><p>text</p></td>`，表达式 `</?t.+>` 不返回任何匹配项，这是因为右尖括号是使用 `.` 捕获的且不进行重新扫描。空`{}`

积极表达式：碰到符合条件的，就从开始处匹配到结尾，格式为expr q，q为任意限定符。

>> regexp('<txr><txd><p>text</p></tpd>', '</?t.*>' ,'match')

ans =

  1×1 cell 数组

    {'<txr><txd><p>text</p></tpd>'}

消极模式：只要匹配就立即停止搜索，进行下一次匹配，格式为expr q?如

>> regexp('<txr><txd><p>text</p></tpd>', '</?t.*?>' ,'match')

ans =

  1×3 cell 数组

    {'<txr>'}    {'<txd>'}    {'</tpd>'}

主动模式：最大程度的匹配，格式为expr q+，如

>> regexp('<txr><txd><p>text</p></tpd>', '</?t.*+>' ,'match')

ans =

  空的 0×0 cell 数组
  
>> regexp('<txr><txd><p>text</p></tpd>', '<t.*+' ,'match')

ans =

  1×1 cell 数组

    {'<txr><txd><p>text</p></tpd>'}

主动模式较难理解，我觉得可以这样看，他就就是直接把检测到开头的<字符就开始匹配，知道匹配到无可匹配也就是>为止，这时候才到达匹配模式里面的?部分，所以如果?后面有任何的字符，这个匹配出来的都是空字符串。

这个匹配模式跟匹配出来的结果就是两个，一个是最后一位是？,这时候等价于expr *，最后一位不是?，这时候结果是空字符串。

操作符的优先级

操作符	描述
\	转义符
(), (?:), (?=), []	圆括号和方括号
*, +, ?, {n}, {n,}, {n,m}	限定符
^, $,	位置和顺序
\|	“或”操作

最大间隙问题

发表于 2019-11-24

最大间隙问题

题目

对于输入的一组数，个数为 n ，求这 n 个数中大小相邻的两个数之间最大差。假设对于任何实数的下取整函数耗时都是 O(1) ，设计此问题的最大间隙解法。例如。

input.txt output.txt

5 0.8

1.1 1.3 2.5 0.7 2.1

其中这个几个数依次排列 0.7<1.1<1.3<2.1<2.5

可以看到，这几个相邻数字之间的差距最大在 2.1-1.3=0.8，所以求得这几个数字之间的最大间距是 0.8

思考

办法一

对于这个题目，首先我想到的也就是上面这个思路，首先要判断整个数组的两数相邻与否，那么我要做的就是将其排序，排序完，我再将其相邻数做差，那么也就是的都一个相邻两数的数组，得到数组之后，我就寻找数组中的最大元素，找到最大元素，再对应其索引，我就能找到那两个数字了。

但是，我立刻又想到了复杂度问题，因为排序算法的复杂度最低的是分治法，普通的搜索算法复杂度会比排序低的

复杂度计算如下 $$ O(n\log n)+O(n)+O(n)=O(n\log n) \\ O(n\log n) > O(n) $$ 这个复杂度有点尴尬，因为他不是线性复杂度，不符合题目的要求，所以这个思路是不可行的。

那么我只能换其他的办法来做。

办法二

设输入的是n个数，分别为 a₁, a₂, ..., a_n ，设在这些数中最大的数为a_max，最小的是a_min；
在a_max与a_min之间均匀插入n − 2个等分点，将其分为n − 1段，如下图；
把需要剩下的n − 2个数按照大小放入这n − 1个段中，根据鸽巢原理，那么其中必定有一个段是空的，如果有一个段是空的，那么我们就可以知道两个相邻数字之间的差值肯定是在两个段中的；
计算每个段中的最大值最小值low[i]和high[i]，再使用后一段的low[i+1]减去前一段的high[i]，得到这些差值中的最大值就是相邻两数的最大间隙了。

分析一下复杂度，对于找最大最小值，也只需要遍历一次即可，所以复杂度为O(n)，等分点只需要得到单位长度即可，单位长度为 $$ l=\frac{max-mix}{n-1} $$ 每个数要判断在哪个段里面也只需要知道自己跟起始处距离几个单位长度即可 $$ seg[i]=\frac{x[i]-minx}{l}+1 $$ 数组seg[i]是一个纽带氏作用，它的索引与x[i]的数的索引是一样的，数组里面存储的数对于的第几个段。这个过程也是用这个公式就可以计算完成的，所以这个过程的复杂度也是O(n)，最后需要做差，使用low[i+1] − high[i]得到数组的间隙的最大值即可，最后这个过程的复杂吨也只需要O(n)。

根据前面的分析，每个过程的复杂度均为O(n)，那么有限个O(n)的累加复杂度还是O(n)，说明这个方法是可行的，满足题目要求。

伪代码

通过前面的分析，确定了方法二是可行的，那么就可以采用这个方法来进行伪代码书写。

//首先要求出最大值最小值。在这里求出它的索引，这样通过 x 数组访问就可得到它的值了
function maxi
maxi = 1;
for i=1:1:n
	if(x(i) > x[maxi])	maxi = i;
end
return maxi;

function mini
mini = 1;
for i=1:1:n
	if(x(i) > x[mini])	mini = i;
end
return mini;
//计算最大最小值，计算刻度，计算索引映射数据段数组 seg
maxx = x[maxi];minx = x[mini];
l = (maxx - minx)/(n-1);
for i=1:1:n
	seg[i] = (x[i] - minx)/l;
end
//初始化 low 和 high 数组
for i=0:1:n-2

立体角

发表于 2019-10-14

立体角

在谈立体角之前，我们先来复习一下球坐标及其面积分。

球面坐标系及积分

球坐标系

如上图所示，就是一个典型的球坐标系统，在坐标系中的每一点都可以使用(ρ,ϕ,θ)来描述，与笛卡尔坐标系之间的转化如下

$$ x=\rho\sin\phi\cos\theta\\ y=\rho\sin\phi\sin\theta\\ z=\cos\phi $$

上式就是球坐标和笛卡尔坐标之间的转化，这个转化很容从图上看出来。

球面积分

从上图可以看出，面积的微元，长为ρdϕ，因为这个方向上的长度为ρϕ，所以微元为ρdϕ；宽为ρsin ϕdθ，这是因为首先要投影到xy平面，在这个平面上的长度为ρsin ϕ，从图上可以看出，这个长度也就是宽，所以微元的宽为ρsin ϕdϕ。所以球面积分的微元面积ρ²sin ϕdϕdθ。

来看几个例子
1. 球的表面积计算
  
  表面积的求法，就是对于微元而言，ϕ的取值为0到2π，θ的取值也是0到2π。ρ就是球的半径。所以微元dS = ρ²sin ϕdϕdθ。表面积的求法如下： S = ∫₀^2πdθ∫₀^2πρ²sin ϕdϕ → 2πρ²(−cosθ)|₀^2π = 4πρ²，这就是很熟悉的圆表面积计算公式了。S_圆 = 4πr²，这个公式应该在高中时候就经常使用了。
2. 球帽表面积的计算
蓝色部分就叫做球帽，计算这一部分的表面积。首先面积分的微元跟之前的一样，为dS = ρ²sin θϕdϕdθ，不过这个图上的角度标的不太一样，那么改变一下，把xy平面的角度设为ϕ，把与z轴之间的夹角设为α。很容易从图上看出，ϕ ∈ [0,2ϕ], α ∈ [0,θ]，这边要注意这个α，它的取值是从z轴为0，到xy平面为$\frac{\pi}{2}$。所以其表面积为

S = ∫₀^2πdϕ∫₀^θr²sin θdθ $\rightarrow 2\pi r^2(-\cos\theta)|^{\theta}_0=2\pi r^2(1-\cos\theta)=2\pi r^2(1-\frac{r-h}{r})=2\pi rh$，如果h = r，那么求的是半球的面积，面积S = 2πr²，跟之前的球表面积也对上了。

球面坐标复习到这儿，下面就进入正题了。 ## 立体角

角度

在介绍立体角之前，也先做个铺垫，讲一下角度。

如上图所示，平面角，简称角度定义为圆的弧长与半径之间的比值，单位为弧度(rad)。 $$ \theta=\frac{l}{r} $$

立体角

参考平面角的定义，立体角的定义为表面积与半径平方的比值，即

$$ \Omega=\frac{S}{r^2} $$

反映的是从该点出发，向球面区域张成的视野大小，是平面角的三维扩展。

接上面计算的表面积的例子。

球的立体角球的表面积为S = 4πr²，球的半径为r，立体角为$\Omega=\frac{S}{r^2}=4\pi$，这也是最大的立体角。
球帽的立体角球帽的表面积为2πrh，半径为r，立体角为$\Omega=\frac{S}{r^2}=2\pi h$。

需要注意的是，立体角计算也可以理解为是所形成表面在以原点为圆心的球的球面上的投影除以半径的平方。因为要计算的立体角表面不一定是球面的一部分，所以需要先投影到球面再进行计算。

几种算法误差报告

发表于 2019-10-10

几种算法误差报告

本文主要探讨了几种计算天顶角的算法，并且根据已经测得的标准数据进行了误差对比，以此来判断几种算法的精度。

原算法的误差

算法的计算步骤

算法的详情我不是很了解，但是我得到了计算后的数据，将其与标准数据做对比。

算法的实现代码

MATLAB代码，new1是通过.csv文件导入得到的.mat数据。

%连接日期与时间
new1.Date.Format = 'dd.MM.uuuu HH:mm';
new1.Time.Format = 'dd.MM.uuuu HH:mm';
x = new1.Date + timeofday(new1.Time);
%计算误差
y = new1.Real-new1.Calculate;
%制图
plot(x,y);
xlabel('时间');ylabel('误差');
title('天顶角实际值与计算值的误差');
box off
grid on

算法的最终结果

得到的结果如图:

根据上图可以看出来，在所给的测试数据中，天顶角的计算误差在-0.006~0.01。

五种算法的公共部分

前时间处理部分

Δτ的计算是根据上图用插值法计算的一个线性表达式

$$ \Delta\tau=96.4+0.00158t\\ $$

t = Δday = (year_now−2060) × 365.2425 参数解释:

根据实际情况，确定了一些参数

参数名称	参数符号	单位	参数值
经度	θ	rad	1.647765346807846
纬度	φ	rad	0.699702497124527
压强	P	atm	0.85862324
温度	T	℃	25

在计算过程中，如果月份m<2，那么把月份加12，年份减1，这里的INT是表示向0取整。前处理过程如下：

ω = 0.017202786day⁻¹

后角度计算部分

计算时只需要知道经纬度，赤经赤纬，时角就可以了。

第一种算法的误差

算法的计算步骤

算法的实现代码

% 主函数
% 第一种算法的计算误差
clear;clc;
% 合并日期，提取信息
load('table.mat', 'new1');
new1.Date.Format = 'dd.MM.uuuu HH:mm';
new1.Time.Format = 'dd.MM.uuuu HH:mm';
x = new1.Date + timeofday(new1.Time);
[y,m,d,h,min,sec] = readDateTime(x);
%调用算法函数计算
data = [y,m,d,h,min,sec];
[GammaAngle,zAngle] = reportFun1(data);
Real = new1.Real;
% 计算差值
y = Real-zAngle';
% 画图
plot(x,y);
xlabel('时间');ylabel('误差');
title('天顶角实际值与计算值的误差');
box off
grid on

% reportFun1函数
% 主要是根据文献的步骤写的
function [GammaAngle,zAngle] = reportFun1(data)
[y,m,d,h,min,sec,lon,lat,pre,tem] = readFile(data);
%需要转化为的数值或者其他常量
[theta,phi,omega,t,te]=Preprocess(y,m,d,h,min,sec,lon,lat);
%算法第一步，计算s1=sin(omega*te),c1=cos(omega*te);
s1 = sin(omega*te);c1=cos(omega*te);
%算法第二步，计算s2=2*s1*c1,c2=(c1+s1)*(c1-s1)
s2=2*s1.*c1;c2=(c1+s1).*(c1-s1);
%算法第三步，计算赤经α
alpha = -1.38880+1.72027920*10^(-2)*te+3.199*10^(-2)*s1...
        -2.65*10^(-3)*c1+4.050*10^(-2)*s2+1.525*10^(-2)*c2;
%算法第四步，把α转换到方便的范围:α→mod(α,2π)
alpha = mod(alpha,2*pi);                                 %根据算法用方便的范围代替，此时还是弧度
alphaAngle = alpha*180/pi;                                %转换到角度
%算法第五步，计算赤纬δ
delta = 6.57*10^(-3)+7.347*10^(-2)*s1-3.9919*10^(-1)*c1...
        +7.3*10^(-4)*s2-6.60*10^(-3)*c2;
deltaAngle = delta*180/pi;
%算法第六步，注意都要用弧度                              
H = 1.75283+6.3003881*t+theta-alpha;
%算法第七步,转到方便的范围H→mod(H+π，2π)-π
H = mod(H+pi,2*pi)-pi;
HAngle = H*180/pi;                                        %转角度输出
% %Final Step,几个算法都一样的
[GammaAngle,zAngle]=finalStep(phi,delta,H,tem,pre);
end

% ReadFile函数
% 读取年月日时分秒
function [y,m,d,h,min,sec,lon,lat,pre,tem] = readFile(data)
y = data(:,1);
m = data(:,2);
d = data(:,3);
h = data(:,4);
min = data(:,5);
sec = data(:,6);
[~,col] = size(data);
if(col>6)
lon = data(:,7);
lat = data(:,8);
pre = data(:,9);
tem = data(:,10);
else
    [lon,lat,pre,tem] = definePara();
end
y=y';
m=m';
d=d';
h=h';
min=min';
sec=sec';
lon=lon';
pre=pre';
tem=tem';
end

% Preprocess函数
% 主要做一些前处理工作，如计算年积日等，也是根据文献来的
function [theta,phi,omega,t,te]=Preprocess(y,m,d,h,min,sec,lon,lat)
%调用格式[theta,phi,omega,t,te]=Proprocess(y,m,d,h,min,sec,lon,lat)
%输出参数[经度弧度θ，纬度弧度φ，ω常数，与2060相距天数t，t加上Δτ后的te]
%输入参数(年year，月month，日day，时hour，分minute，秒second，经度longitude(°)，纬度latitude(°))
h = h+min/60+sec/3600;              %精确时间
theta = lon*pi/180;                 %theta表示经度的弧度
phi = lat*pi/180;                   %phi表示纬度的弧度
omega = 0.017202786;                %omega是算法给的，单位是day^(-1)
tau = 96.4+0.00158*(y-2060)*365;    %tau是前面用插值曲线计算的，代表论文中的Δτ
if(m<=2)                            %根据前面的描述，如果月份小于2，那么月份+2，年份-1
    m=m+12;
    y=y-1;
end
%根据算法1中的公式计算
%根据论文，也通过验证，代码没问题，下式的t是距离2060的天数，小的为-，大的为+
t = fix(365.25*(y-2000))+fix(30.6001*(m+1))-fix(0.01*y)+d+h/24-21958;
%te就是基于TT独立于地球转动，下面都用这个te来建立全局太阳位置
te = t+1.574*10^(-5)*tau;
end

% finalStep函数
% 几个算法的最终的步骤都是一样的
function [GammaAngle,zAngle]=finalStep(phi,delta,H,tem,pre)
%函数调用格式:[GammaAngle,zAngle]=FinalStep(phi,delta,H,tem,pre)
%输出参数[方位角GammaAngle(°),天顶角zAngle(°)]
%输入参数(纬度phi,赤纬delta,时角H,温度tem,压强pre)
e0 = asin(sin(phi).*sin(delta)+cos(phi).*cos(delta).*cos(H));      %计算升角e0=arcsin(sinφsinδ+cosφcosδcosH),计算都要弧度
%计算Δpe应该是修正量
deltape = -4.26*10^(-5)*cos(e0);                               %以地心为中心
ep = e0+deltape;
%计算方位角Γ=atan2(...)见下
Gamma = atan2(sin(H),cos(H).*sin(phi)-tan(delta).*cos(phi));
GammaAngle = Gamma*180/pi;
deltare = 0.08422*pre./((273+tem).*tan(ep+(0.003138./(ep+0.08919))));
%计算天顶角z
z = pi/2-ep-deltare;
zAngle = z*180/pi;
end

算法的最终结果

根据文献的第一种算法计算的误差不容乐观，在-0.2_0.3之间，参考黄冬师兄算的误差，这个扩大了30倍。说明这个算法的精度不是很高，在文献中，这个算法最后也使用了计算数据进行了比较，他给出天顶角的误差范围在-0.190.19，实际计算的比这个范围稍大些，只能说是我们这个数据不够精确或者文献的数据凑得比较好。

第二种算法的误差

算法的计算步骤

算法的实现代码

% 主函数
% 第二种算法的计算误差
clear;clc;
% 合并日期，提取信息
load('table.mat', 'new1');
new1.Date.Format = 'dd.MM.uuuu HH:mm';
new1.Time.Format = 'dd.MM.uuuu HH:mm';
x = new1.Date + timeofday(new1.Time);
[y,m,d,h,min,sec] = readDateTime(x);
data = [y,m,d,h,min,sec];
[GammaAngle,zAngle] = reportFun2(data);
Real = new1.Real;
y = Real-zAngle';
plot(x,y);
xlabel('时间');ylabel('误差');
title('天顶角实际值与计算值的误差');
box off
grid on

% reportFun2函数
function [GammaAngle,zAngle] = reportFun2(data)
[y,m,d,h,min,sec,lon,lat,pre,tem]=readFile(data);
[theta,phi,omega,t,te]=Preprocess(y,m,d,h,min,sec,lon,lat);
[~,col]=size(te);
%Step1
s1=sin(omega*te);c1=cos(omega*te);
%Step2
s2=2*s1.*c1;c2=(c1+s1).*(c1-s1);
%Step3
s3=s2.*c1+c2.*s1;c3=c2.*c1-s2.*s1;
%Step4
s4=2*s2.*c2;c4=(c2+s2).*(c2-s2);
%Step5
alpha = -1.38880*ones(1,col)+1.72027920*10^(-2)*te+3.199*10^(-2)*s1...
        -2.65*10^(-3)*c1+4.050*10^(-2)*s2+1.525*10^(-2)*c2...
        +1.33*10^(-3)*s3+3.8*10^(-4)*c3+7.3*10^(-4)*s4+6.2*10^(-4)*c4;
%Step6
alpha = mod(alpha,2*pi);
%Step7
delta = 6.57*10^(-3)*ones(1,col)+7.347*10^(-2)*s1-3.9919*10^(-1)*c1...
        +7.3*10^(-4)*s2-6.60*10^(-3)*c2+1.50*10^(-3)*s3...
        -2.58*10^(-3)*c3+6*10^(-5)*s4-1.3*10^(-4)*c4;
%Step8
H = 1.75283*ones(1,col)+6.3003881*t+theta-alpha;
H = mod(H+pi*ones(1,col),2*pi)-pi*ones(1,col);
alphaAngle = alpha*180/pi;
deltaAngle = delta*180/pi;
HAngle = H*180/pi;
[GammaAngle,zAngle]=finalStep(phi,delta,H,tem,pre);
end

其实最大的不同只是reportFun()函数是不一样的，因为这个是算法的核心，主函数也是类似的，写入数据，调用函数，得到输出，最后画一下误差图。

其他的前处理函数，后处理函数都是一模一样的。

算法的最终结果

根据文献的第二种算法计算的误差也是不太OK，而且有个严重的问题。从图上看出来，这个集散的误差范围在-0.2_0.2之间，参考黄冬师兄算的误差，这个误差其实相对于第一种算法减小的不是很多。而且对比文献，他给出天顶角的误差范围在-0.0340.034，实际计算的比这个范围大了6倍，这个误差还是很大的。

第三种算法的误差

算法的计算步骤

算法的实现代码

% 主函数
% 第三种算法的计算误差
clear;clc;
% 合并日期，提取信息
load('table.mat', 'new1');
new1.Date.Format = 'dd.MM.uuuu HH:mm';
new1.Time.Format = 'dd.MM.uuuu HH:mm';
x = new1.Date + timeofday(new1.Time);
[y,m,d,h,min,sec] = readDateTime(x);
data = [y,m,d,h,min,sec];
[GammaAngle,zAngle] = reportFun3(data);
Real = new1.Real;
y = Real-zAngle';
plot(x,y);
xlabel('时间');ylabel('误差');
title('天顶角实际值与计算值的误差');
box off
grid on

function [GammaAngle,zAngle] = reportFun3(data)
[y,m,d,h,min,sec,lon,lat,pre,tem]=readFile(data);
%需要转化为的数值或者其他常量
[theta,phi,~,t,te]=Preprocess(y,m,d,h,min,sec,lon,lat);
[~,col]=size(te);
omegaa = 0.0172019715;
%Step1
lambda = -1.388803*ones(1,col)+1.720279216*10^(-2)*te+3.3366*10^(-2)...
         *sin(omegaa*te-0.06172)+3.53*10^(-4)*sin(2*omegaa*te-0.1163);
%Step2
epsilon = 4.089567*10^(-1)*ones(1,col)-6.19*10^(-9)*te;
%Step3
slambda=sin(lambda);clambda=cos(lambda);
%Step4
sepsilon=sin(epsilon);cepsilon=sqrt(1-sepsilon.^(2));
%Step5
alpha = atan2(slambda.*cepsilon,clambda);
%Step6
if(alpha<0)
    alpha = alpha+2*pi;
end
%Step7
delta = asin(slambda.*sepsilon);
%Step8
H = 1.75283*ones(1,col)+6.3003881*t+theta-alpha;
H = mod(H+pi*ones(1,col),2*pi)-pi*ones(1,col);
alphaAngle = alpha*180/pi;
deltaAngle = delta*180/pi;
HAngle = H*180/pi;
[GammaAngle,zAngle]=finalStep(phi,delta,H,tem,pre);
end

算法的最终结果

从图上看出来，这个计算的误差范围在-0.17_0.17之间。对比文献，他给出天顶角的误差范围在-0.00930.0093，实际计算的比这个范围大了18倍，说明这个差距还是很大的，是有一定问题的。

第四种算法的误差

算法的计算步骤

算法的实现代码

% 主函数
% 第三种算法的计算误差
clear;clc;
% 合并日期，提取信息
load('table.mat', 'new1');
new1.Date.Format = 'dd.MM.uuuu HH:mm';
new1.Time.Format = 'dd.MM.uuuu HH:mm';
x = new1.Date + timeofday(new1.Time);
[y,m,d,h,min,sec] = readDateTime(x);
data = [y,m,d,h,min,sec];
[GammaAngle,zAngle] = reportFun4(data);
Real = new1.Real;
y = Real-zAngle';
plot(x,y);
xlabel('时间');ylabel('误差');
title('天顶角实际值与计算值的误差');
box off
grid on

function [GammaAngle,zAngle] = reportFun4(data)
[y,m,d,h,min,sec,lon,lat,pre,tem]=readFile(data);
[theta,phi,~,t,te]=Preprocess(y,m,d,h,min,sec,lon,lat);
[~,col]=size(te);
omegaa = 0.0172019715;
omegan = 9.282*10^(-4);
%Step1
L = 1.752790*ones(1,col)+1.720279216*10^(-2)*te+3.3366*10^(-2)...
    *sin(omegaa*te-0.06172)+3.53*10^(-4)*sin(2*omegaa*te-0.1163);
%Step2
nu = omegan*te-0.8;
%Step3
deltalambda = 8.34*10^(-5)*sin(nu);
%Step4
lambda = L.*ones(1,col)+pi*ones(1,col)+deltalambda;
%Step5
epsilon = 4.089567*10^(-1)*ones(1,col)-6.19*10^(-9)*te+4.46*10^(-5)*cos(nu);
%Step6
slambda=sin(lambda);clambda=cos(lambda);
%Step7
sepsilon=sin(epsilon);cepsilon=cos(epsilon);
%Step8
alpha = atan2(slambda.*cepsilon,clambda);
%Step9
if(alpha<0)
    alpha = alpha+2*pi;
end
%Step10
delta = asin(slambda.*sepsilon);
%Step11
H = 1.7528311*ones(1,col)+6.300388099*t+theta-alpha+0.92*deltalambda;
H = mod(H+pi*ones(1,col),2*pi)-pi*ones(1,col);
alphaAngle = alpha*180/pi;
deltaAngle = delta*180/pi;
HAngle = H*180/pi;
[GammaAngle,zAngle]=finalStep(phi,delta,H,tem,pre);
end

算法的最终结果

算法计算的误差在-0.17_0.17之间。与算法三计算的差不多，文献中写的误差范围是-0.00910.0091，差了17.5倍左右。

第五种算法的误差

算法的计算步骤

使用算法5计算的时候需要给表格中的参数。

算法的实现代码

% 主函数
% 第三种算法的计算误差
clear;clc;
% 合并日期，提取信息
load('table.mat', 'new1');
new1.Date.Format = 'dd.MM.uuuu HH:mm';
new1.Time.Format = 'dd.MM.uuuu HH:mm';
x = new1.Date + timeofday(new1.Time);
[y,m,d,h,min,sec] = readDateTime(x);
data = [y,m,d,h,min,sec];
[GammaAngle,zAngle] = reportFun4(data);
Real = new1.Real;
y = Real-zAngle';
plot(x,y);
xlabel('时间');ylabel('误差');
title('天顶角实际值与计算值的误差');
box off
grid on

% reportFun5函数
function [GammaAngle,zAngle] = reportFun5(data)
[y,m,d,h,min,sec,lon,lat,pre,tem]=readFile(data);
%需要转化为的数值或者其他常量
[theta,phi,~,t,te]=Preprocess(y,m,d,h,min,sec,lon,lat);
%确定常数
L0 = 1.7527901;
L1 = 1.7202792159*10^(-2);
omegaa = 0.0172019715;
beta = 2.92*10^(-5);
omegan = 9.282*10^(-4);
omega = [1.49*10^(-3),4.31*10^(-3),1.076*10^(-2),1.575*10^(-2),...
         2.152*10^(-2),3.152*10^(-2),2.1277*10^(-1)];
a = [3.33024*10^(-2),3.512*10^(-4),5.2*10^(-6)];
dbeta = -8.23*10^(-5);
d = [1.27,1.21,2.33,3.49,2.67,1.28,3.14]*10^(-5);
b = [-2.0582*10^(-3),-4.07*10^(-5),-9*10^(-7)];
Phi = [-2.337,3.065,-1.533,-2.358,0.074,1.547,-0.488];
     
[~,volume]=size(te);
s1=sin(omegaa*te);c1=cos(omegaa*te);
s2=2*s1.*c1;c2=(c1+s1).*(c1-s1);
s3=s2.*c1+c2.*s1;c3=c2.*c1-s2.*s1;
Sigma7=0;
for i=1:7
    Sigma7=Sigma7+d(i)*sin(omega(i)*te+Phi(i)*ones(1,volume));
end
L=L0*ones(1,volume)+L1*te+(a(1)*s1+a(2)*s2+a(3)*s3+b(1)*c1+b(2)*c2+b(3)*c3)...
    +dbeta*s1.*sin(beta*te)+Sigma7;
nu=omegan*te-0.8;
deltalambda=8.34*10^(-5)*sin(nu);
lambda=L+pi*ones(1,volume)+deltalambda;
slambda=sin(lambda);clambda=cos(lambda);
epsilon = 4.089567*10^(-1)*ones(1,volume)-6.19*10^(-9)*te+4.46*10^(-5)*cos(nu);
sepsilon=sin(epsilon);cepsilon=sqrt(1-sepsilon.^(2));
alpha=atan2(slambda.*cepsilon,clambda);
for i=1:volume
if alpha(i)<0
    alpha(i)=alpha(i)+2*pi;
end
end
delta=asin(slambda.*sepsilon);
H = 1.7528311*ones(1,volume,1)+6.300388099*t+theta-alpha+0.92*deltalambda;
H = mod(H+pi*ones(1,volume),2*pi)-pi*ones(1,volume);
deltaAngle=delta*180/pi;
HAngle=H*180/pi;
alphaAngle=alpha*180/pi;
[GammaAngle,zAngle]=finalStep(phi,delta,H,tem,pre);
end

算法的最终结果

图上看误差大概也是-0.17_0.16，文献的误差范围给的是-0.00250.0027，差距是64倍左右，这个差距是很大的。

分析一下8月7日

看了上面两张图，第一张其实看不出来什么，因为误差相对于测量角度的实际值差了100倍左右，能看出来的只是说一天中太阳天顶角是这么分布的。

第二张图和第一张图结合看其实是能看出来一点东西的，看到其实到了6点左右也是天顶角最小，这也是说明这时候高度角最大，太阳在比较高的地方，这时候计算的误差是小的，在00：00和12：00处误差也较大。

看了上面的三幅图，看到这些计算的结果，可以发现误差的分布也有点稀奇古怪，看不出来是什么类型的分布，因为我觉得数据的总数也不是很大，来计算的点数也不够多。可能这样子反应的规律也不是跟明确。总体上来看，还是误差小的占的比例大，所以可能这也是可以稍微有点欣慰的一件事吧。

总结

几种计算方法的结果

说明

算法	文献给的误差范围	计算的误差范围	相差倍数
1	[-0.19,0.19]	[-0.2,0.3]	1.05~1.5
2	[-0,034,0.034]	[-0.2,0.2]	5.8
3	[-0.0093,0.0092]	[-0.17,0.17]	18.3
4	[-0.0091,0.0093]	[-0.17,0.17]	18.7
5	[-0.0025,0.0027]	[-0.16,0.16]	59.3

从这个表格中看出，几种计算出来的结果都在只能最高保证在0.16。

参考文献

[1] Roberto Grena Five new algorithms for the computation of sun position from 2010 to 2110 Solar Energy

https://www.sciencedirect.com/science/article/pii/S0038092X12000400

协议流程

直接DHCP

中继DHCP

仿真

组网

路由设置

抓包

dhcp discover

dhcp offer

Ubuntu

安装crash工具

检查kdump配置

手动触发sysrq

CentOS

安装crash工具

安装内核调试信息包

判断服务正常开启

主动触发分析

分析调试

Ubuntu

查看coredump是否打开

注意关闭apport服务

CentOS

编辑/etc/security/limits.conf文件

分析崩溃

hexo配置

node.js配置

hexo主题配置

Obsidian配置

github配置

workflows配置

dependabot配置

Windows运行参数分析

正则表达式(一)

元字符

表达式

转义字符

限定符

模式

操作符的优先级

最大间隙问题

题目

思考

办法一

办法二

伪代码

立体角

球面坐标系及积分

球坐标系

球面积分

角度

立体角

几种算法误差报告

原算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

五种算法的公共部分

前时间处理部分

后角度计算部分

第一种算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

第二种算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

第三种算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

第四种算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

第五种算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

编辑`/etc/security/limits.conf`文件