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

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行。

Windows运行参数分析

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

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

阅读全文 »

正则表达式(一)

元字符

元字符（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}	限定符
^, $,	位置和顺序
\|	“或”操作

最大间隙问题

题目

对于输入的一组数，个数为 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_1,a_2,...,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

立体角

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

球面坐标系及积分

球坐标系

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

\[ x=\rho\sin\phi\cos\theta\\ y=\rho\sin\phi\sin\theta\\ z=\cos\phi \]

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

球面积分

从上图可以看出，面积的微元，长为$\rho d\phi$，因为这个方向上的长度为$\rho\phi$，所以微元为$\rho d\phi$；宽为$\rho\sin\phi d\theta$，这是因为首先要投影到$xy$平面，在这个平面上的长度为$\rho \sin\phi$，从图上可以看出，这个长度也就是宽，所以微元的宽为$\rho \sin\phi d\phi$。所以球面积分的微元面积$\rho^2 \sin\phi d\phi d\theta$。

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

\[ S=\int_0^{2\pi}d\phi\int_0^{\theta}r^2\sin\theta d\theta \] $\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\pi r^2$，跟之前的球表面积也对上了。

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

角度

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

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

立体角

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

\[ \Omega=\frac{S}{r^2} \]

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

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

球的立体角球的表面积为$S=4\pi r^2$，球的半径为$r$，立体角为$\Omega=\frac{S}{r^2}=4\pi$，这也是最大的立体角。
球帽的立体角球帽的表面积为$2\pi 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=\Delta day=(year_{now}-2060)\times365.2425 \] 参数解释:

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

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

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

\[ \omega=0.017202786day^{-1} \]

后角度计算部分

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

第一种算法的误差

算法的计算步骤

算法的实现代码

% 主函数
% 第一种算法的计算误差
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倍，说明这个差距还是很大的，是有一定问题的。

第四种算法的误差

算法的计算步骤

.png)

算法的实现代码

% 主函数
% 第三种算法的计算误差
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计算的时候需要给表格中的参数。

.png)

算法的实现代码

% 主函数
% 第三种算法的计算误差
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

三种大气质量和一个公式的计算方法

发表于 2019-10-04

三种大气质量和一个公式的计算方法

简介

整理的文献为 Revised optical air mass tables and approximation formula 。

文章一开始介绍了一个由 Karsten 在1965年发表的并且广泛被世界所采用的关于大气质量的近似公式，并且讨论了一些由于各个学科对于不同的物理量符号和术语的不同使得读者经常由此而困惑。

其后介绍了一个计算大气光学质量的近似公式，然后说明了在公式中存在的一种不定情况，之后又对这个近似公式用非线性最小二乘法进行修正得到了一组新的系数。后面又根据索引文献[1]文中也多次提到这篇文献，很多都是从这篇文献里面来的。

还有从一篇《基于遥感与地面监测数据的城市气溶胶定量反演研究》，作者是王耀庭，南京师范大学博士论文。

一个通用的计算公式

\[ m(\gamma)=\frac{m_{abs}(\gamma)}{m_{abs}(90^{\circ})} \]

\[ m_{abs}(\gamma)=\rho_{0}\int^{\infty}_{0} \frac{\rho}{\rho_{0}}([1-[1+2\delta_{0}(1-\frac{\rho}{\rho_{0}})]]\times [\frac{\cos \gamma}{1+\frac{h}{R}}]^{2})^{-\frac{1}{2}}dh \] $h$是相对于海平面的平均高度；

$\rho=\rho(h)$，是在高度$h$处的大气质量；

根据以上的式子 (1),(2) 和已知的参数表1。要计算这个定积分，那就还需要知道 $ρ$, 也就是 $ρ(h)$ 在高度 h 处的大气密度，但是我在文献中找不到，这是个问题，不知道是不是需要再去别的地方找这个 $ρ$，看完了这篇文章之后，知道了这个 $ρ$ 还是没有找到，但是文章已经给出了计算得到的结果的表格。

近似计算公式和不同的系数

\[ f(\gamma)=[\sin \gamma+a(\gamma+b)^{-c}]^{-1} \] $\gamma$是高度角，单位是$^\circ$;$f(\gamma)$是用近似公式计算的$m(\gamma)$;$a,b,c$是式子的常数，$a=0.1500，b=3.885^{\circ}，c=1.253;$

a,b,c这三个常数决定于最小二乘法的相对误差，也就是用前面的计算公式计算数据之后，用最小二乘法进行拟合，使用(3)的形式来计算三个常数。

文献后面又介绍了两个不同的参数组合，一个是根据非线性最小二乘法计算的 $a=0.50572,b=6.07995°,$ $c=1.6364$；一个是根据 $Bemporad$ 的经典大气质量表确定的，$a=0.6556,b=6.379°,c=1.757[1]$，其中文献的表中的 $r(γ)$ 是根据公式(4)计算的相对误差，用来衡量计算大气质量的相对误差。 \[ r(\gamma)=\frac{f(\gamma)-m(\gamma)}{m(\gamma)} \]

积分问题和解决

对于公式(2)，积分会在 $γ=0$ 和 $h$ 接近于0的地方不定，在这种情况下，这个积分可以通过执行一个特殊的程序来进行计算，在参考文献[1]中有写这个程序。但是在计算的时候有个错误会混入，在地平线上的值 36.2648 会比实际的小 5% 左右。

举例而言，在 Link 和 Neuzil[3] 文章的表中所给出的地平线上的在1962年美国的标准大气的大气质量是38.16，这跟 1959 年 Karsen 用的 ARDC 模型十分接近。Snider 和 Goldman[4] 给出的关于 1962 年的模型的38.10也是高度相似。Treve[5] 使用1959年的 ARDC 模型，得到了在地平线上的相对大气质量分别是 $0.55μm $ 的38.11和在 $0.70μm$ 的38.08。

还有就是采用一种新的标准来却确定式子 (2) 中的参数会优于旧的模型，也就是最新的国际标准化组织的大气模型 (ISO Standrad Atmophere) 代替 ARDC 模型大气(由国际民航组织 ICAO 提出的)，这个仅有的变化也就是名义地球半径变为 $R=6.356766×10^6m$。

一个计算公式

\[ m=\frac{1}{cos \frac{\pi \theta_0}{180^\circ}+0.15\times(93.885-\theta_0)^{-1.253}} \]

其中m是需要计算的大气质量，$\theta_0$ 是天顶角。

我要做的工作

在这篇文章里面，我要做的就是编写一个程序，根据文献中的大气质量近似公式(3)，并且用不同的参数组带入，将表格中自变量太阳高度角γ作为自变量带入近似公式计算，再与表格中所给的大气质量数进行作差比较，即验证这个算法是否真的符合实际，如果误差较小，则可以用到我们的项目中去。

上面的这张表格也就说明了在文章计算的数据中天顶角 $γ$ 的取值变化，也就是计算的时候自变量所采用的值。过计算得到了一些结果。

结果

图1 标准大气质量和用其他公式计算的大气质量

图2 四种计算方法与标准大气质量

图3 三种拟合系数计算的大气质量

图4 误差曲线

分析

从上面的图中可以看到，用三种不同的系数计算的相对大气质量以及三组拟合系数的误差曲线，从图中可以看到，三者在天顶角大于 30° 之后都是差不多的经度，主要就是在30°之前的差异。而且可以看到在起始点的时候，第一组和第三组都有很大的误差，特别是第三组，误差都接近于4%，回想文章中提到的积分会在 $γ=0$ 和 h 接近于 0 的地方不定，需要查阅参考文献[1]来寻找解决方法。但是我看到这个计算的第二组拟合系数表现的很好，不知道是否可以用第二组数据来计算，或者是这三组数据都是在不同的情况下表现的经度水平会不一样。但是有个问题，我们没有找到拟合系数1的这条曲线，在下面会进行说明，实际上它是和其他公式计算的这条曲线重合了。

对比三种方法和一个计算公式，发现计算公式的误差在几个计算方法折中的位置，在角度 >10° 之后，这个计算值的偏差与第二组拟合系数计算的误差一样，都是非常小的。

可以看看在高度角大于10°时候的表现。

图 5 在高度角大于10°时候计算大气质量的误差

这里没有找到拟合系数1这条曲线，是因为他的变化与公式计算的是一模一样的，两条线是重合的。

图 6 拟合系数1和公式计算的天顶角大于10°的误差曲线

这说明，其实拟合系数1也就是将天顶角计算的公式做了稍微的变化，就得到了太阳高度角的，本质上，这两个公式是一模一样的，只是取得系数不同罢了。

新给的数据的计算

在之后使用已经写好的这几种计算方法来计算新的数据值，数据可以在 ‘1.xlsx’ 表格中找到。

图 7 根据所给的数据计算的三天的大气质量数值

可以看到，其实一天内的天顶角并不是全是 0~90° 的，这几天维持在 40° 以下，这时候，查阅标准数据也是差不多在这样的数据范围。

参考文献

[1] Kasten,"A New Table and Approximation Formula for the Relative Optical Air Mass",Arch.Meteorol.Geophys.Bioklimatol.Ser.B 14,206-223(1965). [2] R.A Miner,K.S.W.Chamption,and H.L.Pond,The ARDC Model Atmosphere,1959,Air Force Surveys in Geophysics 11(AFCRL,1959) [3] F.Link and L.Neuzil,Table of Light Trajectories in the Terrestrial Atmosphere(Hermann,Paris,1969) [4] D.E Snier and A. Goldman,Refractive Effects in Remote Sensing of Atmosphere with Infrared Transmission Spectroscopy,(Ballistic Research Labratories,June 1975) [5] Y. M. Treve, New Values of the Optical Air Mass and the Refraction and Comparison with Previous Tables," in Proceed-ings, Second Tropospheric Refraction Effects Technical ReviewMeeting, Technical Documentary Rep. ESD-TDR 64-103, May1964 (National Technical Information Service Order AD-442626), pp.5-391. [6] International Organization for Standardization,Standard Atmosphere,International Standard ISO253(1972) [7] S.L.Valley,Handbook of Geophysics and Space Physics (AFCRL,1965), pp.23. [8] W.H.Jefferys,M.J.Fitzpatrick,B.E.McArthur,andJ.E. McCartney, GaussFit:A System for Least Squares and RobustEstimation (U. Texas at Austin, 1989). [9] A.T.Young,Observational Technique and Data Reduction," inle to Methods of Experimental Physics(Vol. 12, Astrophysics; Partrmly A:Optical and Infrared),N,Carleton,Ed.(Academic, New York, 1974),p.150.

代码

% 主函数
clear;clc;
%x的步长选取
x1=0:0.1:20;
x2=20.2:0.2:30;
x3=30.5:0.5:55;
x4=56:1:90;
x=[x1,x2,x3,x4];
%定义Latxe字符
gamma=texlabel('gamma');flambda=texlabel('f(lambda)');
txt = texlabel('f(lambda)=[sin gamma+a(gamma+b)^{-c}]^{-1}');
%输入标准数据，画标准数据图像
data=xlsread('datain.xlsx');
datax=data(:,1);px=datax;
datay=data(:,2);py=datay;
subplot(2,1,1)
plot(datax,datay);grid on;box off;
xlabel(['高度角',gamma]);ylabel({'标准的相对大气质量';flambda});
% print('Standard','-deps');
%调用三个计算函数
airMass1=massCal1(x);
airMass2=massCal2(x);
airMass3=massCal3(x);
airMass4=massCal4(x);
subplot(2,1,2)
plot(datax,airMass4);
xlabel(['高度角',gamma]);ylabel({'其他公式计算的相对大气质量';flambda});
suptitle('相对大气质量与高度角的关系');grid on;box off

figure(2);
set(figure(2),'PaperSize',[42,60]);
suptitle(['\fontsize{14}',txt]);
subplot(3,1,1)
plot(x,airMass1);grid on;box off;
text(60,30,'a=0.1500,b=3.885,c=1.253');
xlabel(['高度角',gamma]);ylabel({'计算的相对大气质量';flambda});
subplot(3,1,2)
plot(x,airMass2);grid on;box off;
text(60,30,'a=0.50572,b=6.07995,c=1.6364');
xlabel(['高度角',gamma]);ylabel({'计算的相对大气质量';flambda});
subplot(3,1,3)
plot(x,airMass3);grid on;box off;
text(60,30,'a=0.6556,b=6.379,c=1.757');
xlabel(['高度角',gamma]);ylabel({'计算的相对大气质量';flambda});
set(gcf, 'position', [1920/4 1080/4 1920/2 1080/1.5]);
% print('Calculate','-deps');
%画三个函数的误差曲线
datay=datay';
delta1 = (airMass1-datay)./datay*100;
delta2 = (airMass2-datay)./datay*100;
delta3 = (airMass3-datay)./datay*100;
delta4 = (airMass4-datay)./datay*100;
figure(3)
plot(x,delta1);grid on;box off;hold on
plot(x,delta2,'--');grid on;box off;hold on
plot(x,delta3,':');grid on;box off;hold on
plot(x,delta4);grid on;box off;
xlabel(['高度角',gamma]);ylabel('三种大气质量的误差(%)');
legend('第1组拟合系数','第2组拟合系数','第3组拟合系数','其他公式计算');
title('几种计算方法的误差');
% print('Error','-deps');
figure(4)
plot(datax,datay);grid on;box off;hold on
plot(x,airMass1);grid on;box off;hold on
plot(x,airMass2);grid on;box off;hold on
plot(x,airMass3);grid on;box off;hold on
plot(x,airMass4);grid on;box off;
xlabel(['高度角',gamma]);ylabel('相对大气质量');
legend('实际的','方法1','方法2','方法三','其他方法');
title('几种方法和实际的相对大气质量');

[~,Date,~] = xlsread('1.xlsx','A2:A105');
Date = datetime(Date,'InputFormat','dd/MM/yyyy');
Date.Format = 'yyyy-MM-dd';
Time = days(xlsread('1.xlsx','B2:B105'));
Time.Format = 'hh:mm:ss';
datetime = Date+Time;
datetime.Format = 'yyyy-MM-dd hh:mm:ss';
xxx=xlsread('1.xlsx','C2:C105');
xx=90-xxx;
airMassx1=massCal1(xx);
airMassx2=massCal2(xx);
airMassx3=massCal3(xx);
airMassx4=massCal4(xx);
output=[airMassx1,airMassx2,airMassx3,airMassx4];
xlswrite('output.xlsx',output);

figure
datetime = datenum(datetime);
plot(datetime,airMassx1);grid on;box off;hold on
plot(datetime,airMassx2);grid on;box off;hold on
plot(datetime,airMassx3);grid on;box off;hold on
plot(datetime,airMassx4);grid on;box off;hold on
dateFormat = 'yy-mm-dd HH:MM:SS';
datetick('x',dateFormat)
%plot(px,py);grid on;box off;
legend('计算方法1','计算方法2','计算方法3','其他方法计算的');
xlabel('时间');ylabel('大气质量');
title('新给的数据的高度角计算的值')

p = (x>10);
x = x(p);
delta1 = delta1(p);
delta2 = delta2(p);
delta3 = delta3(p);
delta4 = delta4(p);
figure
plot(x,delta1);hold on
plot(x,delta2);hold on
plot(x,delta3);hold on
plot(x,delta4);grid on;box off;
xlabel('高度角');ylabel('大气质量计算偏差');
title('高度角 >10° 时大气质量计算偏差');
legend('拟合系数1','拟合系数2','拟合系数3','公式计算');

figure
subplot(2,1,1)
plot(x,delta1);box off;grid on
xlabel('高度角');ylabel('大气质量计算偏差');
title('拟合系数1计算');
subplot(2,1,2)
plot(x,delta4);box off;grid on
xlabel('高度角');ylabel('大气质量计算偏差');
title('公式计算');
suptitle('高度角 >10° 时大气质量计算偏差')

% 计算方法1
function y=massCal1(gamma)
%函数调用格式：airMass=massCal1(gamma)
%输入参数说明：gamma是天顶角，单位是°
%输出参数说明：airMass是大气质量
a=0.1500;b=3.885;c=1.253;
% gamma=gamma*pi/180;%b=b*pi/180;
y=(sin(gamma*pi/180)+a*(gamma+b).^(-c)).^(-1);
end

% 计算方法2
function y=massCal2(gamma)
%函数调用格式：airMass=massCal2(gamma)
%输入参数说明：gamma是天顶角，单位是°
%输出参数说明：airMass是大气质量
a=0.50572;b=6.07995;c=1.6364;
y=(sin(gamma*pi/180)+a*(gamma+b).^(-c)).^(-1);
end

% 计算方法3
function y=massCal3(gamma)
%函数调用格式：airMass=massCal(gamma)
%输入参数说明：gamma是天顶角，单位是°
%输出参数说明：airMass是大气质量
a=0.6556;b=6.379;c=1.757;
y=(sin(gamma*pi/180)+a*(gamma+b).^(-c)).^(-1);
end

考拉兹猜想

发表于 2019-09-25

考拉兹猜想

考拉兹猜想定义

考拉兹函数定义如下

\[ f(x)=\left\{ \begin{array}{**lr**} 3n+1&x为奇数且x\neq1\\ n/2&x为偶数\\ 1&x=1 \end{array} \right. \] 通过对 $x$ 取不同的值，发现最后都会收敛到 1。求该函数构成算法的上下界。

当然，下界是很容易求出来的，如果输入 $n$ ，下降最快的也就是每次下降 $\frac{1}{2}$，这个下降速度对于的计算时间是 $\log n$ 。对于上界，用 $MATLAB$ 带入一些数值计算。得到的结果如下：

部分 $n_0$ 的步长与 $n$ 值变化

$n_0$ 与计算次数分布曲线

论文笔记-场地自动化定标方法研究及应用

发表于 2019-03-01

论文笔记-场地自动化定标方法研究及应用

辐射定标的意义

卫星获取地物图像时候，由于存在各种各样的干扰，造成观测值与实际值的偏差，定标就是来（尽可能）消除这些偏差，使得卫星得到的图像与真实的物体之间（尽可能）没有这些偏差，即反映客观事物。

辐射定标

辐射校正

遥感辐射定标方法

因为卫星所在环境为外太空，会受到强烈的电磁干扰和强辐射，所以器件的老化程度是很快的，这也就是卫星存在使用寿命的一个原因，受限于载荷的使用寿命。那么要延长载荷的工作年限，就需要实时测出其老化程度，对于传感器的灵敏度做出评估，改变具体参数使得其能够继续工作。比如之前对于输入为 10 就能输出 10 ，现在只能输出 5，就乘以系数 2 就可以完成跟之前一样的功能。这样器件就还可以接着使用。

阅读全文 »

Ubuntu

安装crash工具

检查kdump配置

手动触发sysrq

CentOS

安装crash工具

安装内核调试信息包

判断服务正常开启

主动触发分析

分析调试

Ubuntu

查看coredump是否打开

注意关闭apport服务

CentOS

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

分析崩溃

Windows运行参数分析

正则表达式(一)

元字符

表达式

转义字符

限定符

模式

操作符的优先级

最大间隙问题

题目

思考

办法一

办法二

伪代码

立体角

球面坐标系及积分

球坐标系

球面积分

角度

立体角

几种算法误差报告

原算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

五种算法的公共部分

前时间处理部分

后角度计算部分

第一种算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

第二种算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

第三种算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

第四种算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

第五种算法的误差

算法的计算步骤

算法的实现代码

算法的最终结果

分析一下8月7日

总结

几种计算方法的结果

说明

参考文献

三种大气质量和一个公式的计算方法

简介

一个通用的计算公式

近似计算公式和不同的系数

积分问题和解决

一个计算公式

我要做的工作

结果

分析

新给的数据的计算

参考文献

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