冲毕业结束！

转移一下blog 回到我小破电脑

---

还有部署的bug 同步gitee的bug

在根目录_config.yml

# Deployment

## Docs: https://hexo.io/docs/one-command-deployment

deploy:

type: git

# repo: https://github.com/tsuiwade/tsuiwade.github.io

repo:

# gitee: https://gitee.com/tsuiwade/tsuiwade

# github: https://github.com/tsuiwade/tsuiwade.github.io

gitee: https://tsuiwade:私人令牌@gitee.com/tsuiwade/tsuiwade.git

github: https://tsuiwade:私人令牌（token）@github.com/tsuiwade/tsuiwade.github.io.git

branch: master

设置aactions -》settings -》 Actions secrets and variables -》 GITEE_RSA_PRIVATE_KEY -》为私钥.ssh id_rsa

搞了很久！

---

还是非常感谢 https://www.zhihu.com/question/373117145

---

接下来需要考虑二维的问题了。

突然发现一个这个算法

https://stackoverflow.com/questions/18227014/selecting-evenly-distributed-points-algorithm

算法是遍历前面的最小cost

double tempcost = solutions[j].cost + std::abs(dis[i] - dis[j] - 0.1);

0.1为1/10 但问题是凑不到10个点

在下面的回答中看到了

提到了ANMS算法自适应非极大值抑制算法

有点像那么回事

---

应该会有一些对比实验： 非极大值抑制算法ANMS、bucketing算法 https://blog.csdn.net/frozenspring/article/details/117174032

https://blog.csdn.net/u013453604/article/details/45350653 ANMS

对比实验

EFFICIENTLY SELECTING SPATIALLY DISTRIBUTED KEYPOINTS FOR VISUAL TRACKING 2011

2005年 Matthew Brown在论文Multi-Image Matching using Multi-Scale Oriented Patches提出了自适应非极大值抑制算法（Adaptive Non-Maximum Suppression，ANMS

最近有一篇论文讨论了图像上均匀关键点分布的问题。这个资源库提供了 C + + 、 Python 和 Matlab 接口

There is a recent paper that tackles the problem of homogeneous keypoint distribution on the image. C++, Python, and Matlab interfaces are provided in this repository

---

L. F. Kozachenko and N. N. Leonenko. Sample estimate of entropy of a random vector. Problems of Information Transmission, 23:95–101, 1987.

https://kaba.hilvi.org/tim-1.3.0/tim/core/differential_entropy_kl_details.htm

---

学习b站kd树

https://www.bilibili.com/video/BV1d5411w7f5/?spm_id_from=333.337.search-card.all.click&vd_source=6c92aa3e5d0f2e0347ec135013a906d8

后来在思考我的矩形是首尾相连的，该怎么办呢。

学习kd树的搜索原理，自己画了个图就能理解了，当点的最小半径 < 该点的x，说明会和x=0相交，那么将该点的x轴加360，再进行搜索，最终比较两个半径，取最小的半径作为distance。

from scipy.special import gamma, digamma
import numpy as np
from entropy_estimators import continuous
X = np.loadtxt(open("J2000_5_hor2.csv", "rb"), delimiter=",", skiprows=0)
print(X.shape)
kozachenko = continuous.get_h(X, k=1, norm='euclidean')

print(f"K-L estimator: {kozachenko:.5f}")

所有点的 K-L estimator: 11.33121（魔改distance前）

魔改效果为

改善第51个点，魔改算法为

n, d = x.shape

if norm == 'max':  # max norm:
    p = np.inf
    log_c_d = 0  # volume of the d-dimensional unit ball
elif norm == 'euclidean':  # euclidean norm
    p = 2
    log_c_d = (d/2.) * log(np.pi) - log(gamma(d/2. + 1))
else:
    raise NotImplementedError(
        "Variable 'norm' either 'max' or 'euclidean'")
kdtree = cKDTree(x)
# query all points -- k+1 as query point also in initial set
# distances, _ = kdtree.query(x, k + 1, eps=0, p=norm)
distances, _ = kdtree.query(x, k + 1, eps=0, p=p)
distances = distances[:, -1]
print("distances11 ", len(distances))
print("x11 ", x[0, 1], " distances ", distances[0])
for cur in range(len(x)):
    if (x[cur, 0] < distances[cur]):
        value = x[cur, 0] + 360.0
        distances1, _ = kdtree.query([value, x[cur, 1]], k, eps=0, p=p)
        print(" cur ", cur, " va ", value, " distances1 ",
              distances1, " qian ", distances[cur])
        # print(" distances1 ", distances[cur])
        distances[cur] = min(distances[cur], distances1)
    elif (x[cur, 0] > 360.0-distances[cur]):
        value = x[cur, 0] - 360.0
        print(" va ", value)
        distances1, _ = kdtree.query([value,  x[cur, 1]], k, eps=0, p=p)
        distances[cur] = min(distances[cur], distances1)

# print("distances ", distances)
# enforce non-zero distances
distances[distances < min_dist] = min_dist

# where did the 2 come from? radius -> diameter
sum_log_dist = np.sum(log(2*distances))
h = -digamma(k) + digamma(n) + log_c_d + (d / float(n)) * sum_log_dist

return h

h值为11.33087（魔改后，虽然变小 ，考虑到左右宽度）

---

接下来改为matlab

网上找的代码

https://blog.csdn.net/john_xia/article/details/107563005 不好用 直接用内置函数

data = [2 3;
    5 4;
    9 6;
    4 7;
    8 1;
    7 2];
Mdl = KDTreeSearcher(data)
[n,d] = knnsearch(Mdl,[6,3.1],'k',10);

n =

 2     6     5     1     3     4

d =

1.3454    1.4866    2.9000    4.0012    4.1725    4.3829

clear;close all;clc;
tic;
M = csvread('J2000_5_hor.csv');
points = M(:, 2:3);
k = 1;
Mdl = KDTreeSearcher(points);
[n,d] = knnsearch(Mdl,points,'k',k+1);
distances = d(:,end);
for i = 1 : length(points)
    if ( points(i,1) < distances(i) ) 
        value = points(i,1) + 360.0;
        [~, dd] = knnsearch(Mdl,[value, points(i, 2)],'k',1);
        distances(i) = min(distances(i), dd);
    elseif (points(i, 1) > 360.0-distances(i) ) 
        value = points(i,1) - 360.0;
        [~, dd] = knnsearch(Mdl,[value, points(i, 2)],'k',1);
        distances(i) = min(distances(i), dd);        
    end
end
Copy_2_of_totalH(distances,length(points))
toc;

function  h =  Copy_2_of_totalH(distances, len)
    d = 2;
    log_c_d = (d/2.) * log(pi) - log(gamma(d/2 + 1));
    sum_log_dist = sum(log(2*distances));
    h = - psi(1) + psi(len) + log_c_d + (d/len) * sum_log_dist;
end

验证是对的

接着

%% 模拟退火求解选点问题 直接计算h 更慢  % 466
clear all; clc; close all;
tic
% load points;
M = csvread('J2000_5_hor.csv');
points = M(:, 2:3);
% points = sort(points);
rng('shuffle')  % 控制随机数的生成，否则每次打开matlab得到的结果都一样
len = 50;
map = [1:len,1:len];

%% 参数初始化
T0 = 1000;   % 初始温度
T = T0; % 迭代中温度会发生改变，第一次迭代时温度就是T0
maxgen = 800;  % 最大迭代次数
Lk = 200;  % 每个温度下的迭代次数
alfa = 0.95;  % 温度衰减系数

%%  随机生成一个初始解
A = int32(1:length(points));
random_num = A(randperm(numel(A),len));
index0 = sort(random_num);
h0 = Copy_2_of_totalH(points, index0);

disp('初始方案是：'); disp(mat2str(index0))
disp('此时最优值是：'); disp(h0)
min_money = h0;     % 初始化找到的最佳的解对应的花费为money0
MONEY = zeros(maxgen,1); % 记录每一次外层循环结束后找到的min_money (方便画图）

%% 模拟退火过程
for iter = 1 : maxgen  % 外循环, 我这里采用的是指定最大迭代次数
    for i = 1 : Lk  %  内循环，在每个温度下开始迭代
        index = index0;
        x_num = index(randperm(numel(index),1));
        x = find(index == x_num);

        
        shengxia = setxor(A, index);
        y_num = shengxia(randperm(numel(shengxia),1));
        index(x) = y_num;
        index = sort(index);

        index1 = index;  % 调用我们自己写的gen_new_way函数生成新的方案

        h1 = Copy_2_of_totalH(points, index1);

        if h1 > h0    % 如果新方案的花费小于当前方案的花费
            index0 = index1; % 更新当前方案为新方案
            h0 = h1;
        else
            p = exp((h1 - h0)/T); % 根据Metropolis准则计算一个概率
            if rand(1) < p   % 生成一个随机数和这个概率比较，如果该随机数小于这个概率
                index0 = index1; % 更新当前方案为新方案
                h0 = h1;
            end
        end
        % 判断是否要更新找到的最佳的解
        if h0 > min_money  % 如果当前解更好，则对其进行更新
            min_money = h0;  % 更新最大的花费
            best_way = index1;  % 更新找到的最佳方案
        end
    end
    MONEY(iter) = min_money; % 保存本轮外循环结束后找到的最小花费
    T = alfa*T;   % 温度下降
end

disp('最佳的方案是：'); disp(mat2str(best_way))
disp('此时最优值是：'); disp(min_money)

%% 画出每次迭代后找到的最佳方案的图形
figure
plot(1:maxgen,MONEY,'b-');
xlabel('迭代次数');
ylabel('最小花费');
toc

function  h =  Copy_2_of_totalH(points,index)
    len = length(index);
    point = points(index, :);
    k = 1;
    Mdl = KDTreeSearcher(point);
    [n,d] = knnsearch(Mdl,point,'k',k+1);
    distances = d(:,end);
    for i = 1 : length(point)
        if ( point(i,1) < distances(i) ) 
            value = point(i,1) + 360.0;
            [~, dd] = knnsearch(Mdl,[value, point(i, 2)],'k',1);
            distances(i) = min(distances(i), dd);
        elseif (point(i, 1) > 360.0-distances(i) ) 
            value = point(i,1) - 360.0;
            [~, dd] = knnsearch(Mdl,[value, point(i, 2)],'k',1);
            distances(i) = min(distances(i), dd);        
        end
    end
    d = 2;
    log_c_d = (d/2.) * log(pi) - log(gamma(d/2 + 1));
    sum_log_dist = sum(log(2*distances));
    h = - psi(1) + psi(len) + log_c_d + (d/len) * sum_log_dist;
end

初始方案是：
[34 38 53 80 81 111 121 134 140 143 194 203 217 230 233 242 246 253 260 271 272 273 287 288 364 370 371 373 411 413 420 472 479 488 500 508 527 571 593 603 615 624 628 632 641 642 647 651 673 684]
此时最优值是：
11.6156

最佳的方案是：
[16 20 22 31 36 40 41 50 52 53 66 68 75 77 85 93 131 149 150 156 164 194 213 244 274 323 345 351 355 388 394 403 415 426 447 490 493 504 510 557 565 566 616 637 645 658 660 664 665 703]
此时最优值是：
13.4517

历时 57.334977 秒。

vscode python 验证

选择了这些点

---

2023/2/3 续

但是耗时过长

由于每次只替换一个点，因此可以快速计算，函数(一堆点集，要替换的点，新点)

现在还需要增加一个功能，就是记录那50个点中最近的点的序号

找最近的点，如果是要被删的，那么就重新计算这些点，但是要增加的点呢

还是需要全部进行查询最近点

所以该改进点只得作罢

---

开始写文章吧

---

画圆

https://ww2.mathworks.cn/help/matlab/ref/cylinder.html

https://blog.csdn.net/weixin_40525909/article/details/106636093

clc;clear all;
t = linspace(0,2*pi,1);
a = size(t); 
r = 1 * ones(a);
figure(1);
% polar(t,r,'w');
polar([0],[0],'w');
hold on;
% cylinder()
[X,Y,Z] = cylinder(1,400);
Z = Z+0.3;
surf(X,Y,Z);
view(90,30)
shading flat 
colorbar
opengl('save','hardwarebasic')

画好红 插入灯光，修改颜色，去除坐标点，加上°，加上colorbar,右击打开颜色图编辑器，指定颜色，删除颜色 设置两头黑，P图。

---

clf;
% box on;
axis equal;
xticks([0 60 120 180 240 300 360])
yticks([0 30 60 90])
xlim([-10 370]);
ylim([-10 100]);
% grid on;
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;
hold on;
scatter([10,20,350,370],[45,20,50,45],12);
xt=get(gca,'xtick');
for k=1:numel(xt);
	xt1{k}=sprintf('%d°',xt(k));
end
set(gca,'xticklabel',xt1);
yt=get(gca,'ytick');
for k=1:numel(yt);
	yt1{k}=sprintf('%d°',yt(k));
end
set(gca,'yticklabel',yt1);
line([0,360],[90,90],'color','k','linewidth',1);
line([0,0],[0,90],'color','k','linewidth',1);
line([0,360],[0,0],'color','k','linewidth',1);
line([360,360],[0,90],'color','k','linewidth',1);
% xlabel('\bf 方位轴');
% ylabel('\bf 高度轴');
r = pdist([[10,45];[20,20]],'euclidean');
rectangle('position',[10-r,45-r,2*r,2*r],'curvature',[1,1],'LineStyle',':');
rectangle('position',[370-r,45-r,2*r,2*r],'curvature',[1,1],'LineStyle',':');
theta=0:0.001:2*pi;
Circle1=10+r*cos(theta);
Circle2=45+r*sin(theta);
Circle1(Circle1 < 0) = 0;
plot(Circle1,Circle2,'k','linewidth',0.5);

Circle1=370+r*cos(theta);
Circle2=45+r*sin(theta);
Circle1(Circle1 > 360) = 360;
plot(Circle1,Circle2,'k','linewidth',0.5);
txt = ' M';
text(10,45,txt);
txt = ' N';
text(20,20,txt)
txt = ' M^{''}';
text(360,45,txt);
txt = ' P';
text(350,50,txt)

https://blog.csdn.net/luochao5862426/article/details/89219311

Matlab图形高清插入word文档的几种方法。

---

在选取恒星点的k近邻时，需要考虑的问题是如何对星点进行快速k近邻搜索。k近邻最简单的实现是枚举搜索，即对选取星点与其他星点计算距离，时间复杂度为 。当恒星点集很大时，计算非常耗时。可以使用特殊的数据结构存储星点位置，减小计算次数，以提高k近邻搜索效率。

​ KD树是一种对K维空间中的数据进行组织和存储以便进行快速搜索查找的二叉树形数据结构，用垂直于坐标轴的超平面将K维空间划分成一系列的K维超矩形区域。在本文中，KD树的每个节点对应于一个由方位轴和高度轴组成的二维矩形区域，构造得到的KD树是平衡二叉树，搜索查找的时间复杂度为 。

​

对于输入的恒星坐标点集 ，其中 ， ， ，构造**KD**树的方法如下：

BuildKDTree(star pointsT, tree_depth)
1) // 根据树深度的奇偶性来选择AE维度
2) axis = tree_depth & 1
3) // 选择第axis维度的中间值为根节点进行KD树的构造 4) Select median by axis from points 5) // 构造节点和子树 6) node = median 7) node.left = KDTree(points in subList(0, median), axis, tree_depth + 1) 8) node.right = KDTree(points in subList(median+1, points.end()), axis, tree_depth + 1) 9) return node

（1）构造根节点：根节点的深度为0，对应于望远镜两轴指向上半天区所组成的二维空间，其中方位轴范围为0°至360°，高度轴范围为0°至90°。接着，重复通过下面（2）的递归方法，不断地对二维空间进行切分，生成子节点。

（2）重复：对深度为奇数/偶数j的节点重复E/A轴切分：选择E/A轴为坐标轴，以该节点的区域中所有恒星的E/A坐标中位数为切分点，将该节点对应的二维区域切分为两个子区域。由该节点生成深度为j+1的左右子节点，左子节点对应坐标E/A小于切分点的子区域，右子节点对应坐标E/A大于切分点的子区域。并将落在该二维平面上的恒星点保存在该节点。

（3）直到两个子区域内没有恒星时，构造KD树停止，完成划分。

---

​ 利用KD树可以减少搜索的计算量。在望远镜两轴指向形成的圆柱面上，为一个恒星点搜索其k近邻与传统的KD树搜索策略有所差异，本文提出改进的KD树搜索策略IKDS，首先找到包含选取点的叶子节点，然后由该节点出发，依次回退到父节点；不断查找所选星点的最近邻，当不存在更近的节点时中止。具体方法如下：

---

利用KD树可以省去对大部分数据点的搜索，从而减少搜索的计算量。给定一个恒星点，搜索其最近邻，首先找到包含选取点的叶子节点，然后由该节点出发，依次回退到父节点；不断查找所选星点的最近邻，当不存在更近的节点时中止。具体方法如下

---

适应于圆柱面上的**KD**树搜索策略

​ 利用KD树可以省去对大部分数据点的搜索，从而减少搜索的计算量。给定一个恒星点，搜索其最近邻，首先找到包含选取点的叶子节点，然后由该节点出发，依次回退到父节点；不断查找所选星点的最近邻，当不存在更近的节点时中止。由于恒星位于高度轴和方位轴两个垂直轴所组成的圆柱面上，KD树的搜索策略有所不同，具体方法如下：

​ （1）在KD树种找到包含选取点的叶子节点，并以此叶子节点为“当前最近点”。递归地向上回退，在每个节点进行以下（2）（3）操作：

​ （2）如果该节点保存的恒星点比当前最近点更接近选取点，则该节点为当前最近点。

（3）当前最近点一定存在于该节点的一个子节点所对应的区域。检查该子节点的父节点的另一个子节点所对应的区域是否有更近的点。具体的，检查另一个子节点对应的区域是否与目标点为圆心，与目标点与当前最近点间的距离为半径的圆相交。如果相交，可能在另一个子节点对应的区域内存在距离目标更近的点，移动到另一个子节点。接着，递归地进行最近邻搜索。如果不相交，往上回退。

（4）当回退到根节点时，搜索结束。搜索得到的“当前最近点”即为选取点的最近邻。当选取点的方位轴 坐标减去最近邻距离小于0°时，以选取点为中心，最近邻距离为半径的最近邻圆与0°方位轴相切，可能存在方位轴为 （ 为一个较小的数）的星点距离选取星点更近。如上图所示，使用原始KD树搜索找到目标星点M的最近邻N，但方位轴尾端可能存在星点P，距离选取星点M更近。因此需要重新计算选取虚拟星点 的最近邻距离，与 的最近邻距离相比取较小值；同样的，当方位轴 坐标加上最近邻距离大于360°时，还需要计算选取虚拟星点 的最近邻距离。

---

正交网格法

clear all;  close all;clc;
% tic
% load points;
M = csvread('J2000_5_hor.csv');
points = M(:, 2:3);
% points = sortrows(points,1);
rng('shuffle')  % 控制随机数的生成，否则每次打开matlab得到的结果都一样
len = 50;
%%
jiange = 5;
x = 0;
y = 0;
ret = [];

while x < 360 && y < 90
    tmp = find( (points(:,1)> x) & (points(:,2) >y) & (points(:,1)<= x+ jiange) & (points(:,2)<= y+ jiange) );
    if ( ~isempty(tmp) )
        ret = [ret; points(tmp(1),:)];
    end
    y = y + jiange;
    if ( y >= 90 ) 
        y = 0;
        x = x + jiange;
    end
end
% scatter(points(:,1), points(:,2));
scatter(ret(:,1), ret(:,2));

选自 Vedder, J. D., STAR TRACKERS, STAR CATALOGS, AND ATTITUDE DETERMINATION - PROBABILISTIC ASPECTS OF SYSTEM-DESIGN. J. Guid. Control Dyn. 1993, 16, (3), 498-504.

例如间隔为5°，选择了487个点

间隔26.6°可以选择50个星

选择索引 ret = [ret; tmp(1)];

计算Copy_2_of_totalH(points,ret)熵值，12.5235，所有星是11.3309

找寻最近点 确实很快 0.05秒但存在问题 就是找寻50个点需要的jiange

用二分法试试

clear all;  close all;clc;
tic
% load points;
M = csvread('J2000_5_hor.csv');
points = M(:, 2:3);
% points = sortrows(points,1);
rng('shuffle')  % 控制随机数的生成，否则每次打开matlab得到的结果都一样
len = 50;
%%
% jiange = 26.6;
left = 0;
right = 90;

while left < right
    ret = [];
    x = 0;
    y = 0;
    jiange = double( left + right ) / 2
    while x < 360 && y < 90
        tmp = find( (points(:,1)> x) & (points(:,2) >y) & (points(:,1)<= x+ jiange) & (points(:,2)<= y+ jiange) );
        if ( ~isempty(tmp) )
    %         ret = [ret; points(tmp(1),:)];
                ret = [ret; tmp(1)];
        end
        y = y + jiange;
        if ( y >= 90 ) 
            y = 0;
            x = x + jiange;
        end
    end
    length(ret)
    if (length(ret) == 50 ) 
        break;
    elseif ( length(ret) > 50 ) 
        left = jiange;
    else 
        right = jiange;
    end
end
toc;
% scatter(points(:,1), points(:,2));
% scatter(ret(:,1), ret(:,2));

jiange =

26.3672

ans =50

历时 0.044573 秒。 很快

---

画结果图

plot([30,40,50,60],[11.4426	11.5679	11.5762	11.6659],'-bo');hold on;
plot([30,40,50,60],[12.5960	12.5311	12.6415	12.6410],'-v','Color',[0 0.6 0.1]);hold on;
plot([30,40,50,60],[12.8165	12.7812	12.8322	12.8459],'-^','Color',[0.5 0 0.8]);hold on;
plot([30,40,50,60],[13.5425,13.4567,13.4374,13.3992],'-rsquare');hold on;
plot([30,60],[11.3309,11.3309],'--k');hold on;
xticks([30 40 50 60])
yticks([11 12 13 14])
xlim([27 63]);
ylim([11 14.5]);
legend( '1:Random','2:RGM', '3:SO', '4:MKLE','5:All');
xlabel('n');
ylabel('KL 熵');



plot([27,30,44,48, 52],[11.0260	11.5331	11.4096	11.0756	11.8410],'-bo');hold on;
plot([27,30,44,48, 52],[12.2689	12.5960	12.5225	12.5195	12.6379],'-v','Color',[0 0.6 0.1]);hold on;
plot([27,30,44,48, 52],[13.5609	13.5329	13.4409	13.4240	13.3890],'-^','Color',[0.5 0 0.8]);hold on;
plot([27,30,44,48, 52],[13.5537	13.3811	13.3753	13.4138	13.2940],'-rsquare');hold on;

xticks([27 30 44 48 52])
yticks([11 12 13 14])
xlim([25 55]);
ylim([11 14.5]);
legend( 'Random','RGM', 'MKLE', 'ETRP');
xlabel('n');
ylabel('KL 熵');

仿真平台为Windows11操作系统下的Matlab，计算机配置为Intel(R) Core i5-11400F@ 2.6GHz，16GB内存。

没敢写时间

---

初始方案是：
[30 43 86 126 134 146 152 183 214 216 232 263 265 278 283 299 307 309 354 363 376 377 379 396 406 424 425 429 435 439 478 485 490 499 522 526 535 556 561 576 602 611 615 625 638 639 692 698 702 705]
此时最优值是：
11.1946

最佳的方案是：
[8 20 22 36 40 41 42 43 50 53 54 68 74 75 93 97 104 133 150 158 162 189 206 248 274 325 331 343 378 379 380 394 399 413 424 432 460 490 510 514 517 546 576 612 624 647 660 664 699 702]
此时最优值是：
13.4316

历时 88.230554 秒。

---

画结果图。

figure;
subplot(2,1,1);
scatter(points(:,1), points(:,2), 8)
box on; 
axis equal;
xlim([0 360]);
ylim([0 90]);
xticks([0 60 120 180 240 300 360])
yticks([0 30 60 90])
xlabel('方位轴');
ylabel('高度轴');
% title('Distribution of all stars');
subplot(2,1,2);
scatter(points(best_way,1), points(best_way,2), 12)
box on; 
axis equal;
xlim([0 360]);
ylim([0 90]);
xticks([0 60 120 180 240 300 360])
yticks([0 30 60 90])
xlabel('\bf Azimuth axis');
ylabel('\bf Elevation axis');
title('Distribution of selected stars');

或者

figure;
scatter(points(:,1), points(:,2), 8)
box on; 
axis equal;
xlim([0 360]);
ylim([0 110]);
xticks([0 60 120 180 240 300 360])
yticks([0 30 60 90])
xlabel('Azimuth axis');
ylabel('Elevation axis');
hold on;
scatter(points(best_way,1), points(best_way,2), 20,'filled','MarkerFaceColor','#FF8000')
xtickformat('usd')
legend('All stars', 'Selected stars')
xtickformat('degrees')
ytickformat('degrees')

---

画环形图

ax = polaraxes;
ax.RAxis.Label.String = 'Elevation';
ax.ThetaAxis.Label.String = 'Azimuth';
ax.ThetaZeroLocation = 'bottom';
hold on;
polarscatter(double(points(:,1))./180.0 * pi,points(:,2),20)
thetalim([0 360])
rticks(0:22.5:90)
hold on;
polarscatter(double(points(best_way,1))./180.0 * pi,points(best_way,2), 30,'r','filled')
legend('All stars', 'Selected stars')



ax = polaraxes;
ax.RAxis.Label.String = 'Elevation axis';
ax.ThetaAxis.Label.String = 'Azimuth axis';
ax.ThetaZeroLocation = 'bottom';
hold on;
polarscatter(double(points(:,1))./180.0 * pi, 90-points(:,2),20)
thetalim([0 360])
rticks(0:22.5:90)
hold on;
polarscatter(double(points(best_way,1))./180.0 * pi, 90-points(best_way,2), 30,'filled','MarkerFaceColor','#FF8000')
legend('All stars', 'Selected stars')
ax.ThetaTickLabel = string(ax.ThetaTickLabel) + char(176)
rticks([0 22.5 45  67.5 90])                        % 在r = 50,100,200处显示刻度
rticklabels({'90','67.5','45','22.5','0'})            % 在刻度线处加标记
ax.RTickLabel = string(ax.RTickLabel) + char(176)
ax.RAxisLocation = 90;




ax = polaraxes;
ax.RAxis.Label.String = '高度轴';
ax.ThetaAxis.Label.String = '方位轴';
ax.ThetaZeroLocation = 'bottom';
hold on;
polarscatter(double(points(best_way,1))./180.0 * pi, 90-points(best_way,2), 20)
thetalim([0 360])
rticks(0:22.5:90)
hold on;

ax.ThetaTickLabel = string(ax.ThetaTickLabel) + char(176)
rticks([0 22.5 45  67.5 90])                        % 在r = 50,100,200处显示刻度
rticklabels({'90','67.5','45','22.5','0'})            % 在刻度线处加标记
ax.RTickLabel = string(ax.RTickLabel) + char(176)
ax.RAxisLocation = 90;

---

接着做实验

画一个图

clear all;  close all;clc;
tic
jiange = 22.5;
ret = [];
x = 0;
y = 0;
while x < 360 && y <= 90
    ret = [ret; x, y];
 
    y = y + jiange;
    if ( y > 90 ) 
        y = 0;
        x = x + jiange;
    end
end
toc;

3个图拼一起

figure(1)
jiange = 30;
ret = [];
x = 0;
y = 0;
while x < 360 && y <= 90
    ret = [ret; x, y];
 
    y = y + jiange;
    if ( y > 90 ) 
        y = 0;
        x = x + jiange;
    end
end
jiange = 22.5;
ret1 = [];
x = 0;
y = 0;
while x < 360 && y <= 90
    ret1 = [ret1; x, y];
 
    y = y + jiange;
    if ( y > 90 ) 
        y = 0;
        x = x + jiange;
    end
end
jiange = 18;
ret2 = [];
x = 0;
y = 0;
while x < 360 && y <= 90
    ret2 = [ret2; x, y];
 
    y = y + jiange;
    if ( y > 90 ) 
        y = 0;
        x = x + jiange;
    end
end
subplot(3,1,2)

clf;
% box on;
axis equal;
xticks([0 60 120 180 240 300 360])
yticks([0 30 60 90])
xlim([-10 370]);
ylim([-10 100]);
% grid on;
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;
hold on;

xt=get(gca,'xtick');
for k=1:numel(xt);
	xt1{k}=sprintf('%d°',xt(k));
end
set(gca,'xticklabel',xt1);
yt=get(gca,'ytick');
for k=1:numel(yt);
	yt1{k}=sprintf('%d°',yt(k));
end
set(gca,'yticklabel',yt1);


xlabel('方位轴');
ylabel('高度轴');
line([0,360],[90,90],'color','k','linewidth',1);
line([0,0],[0,90],'color','k','linewidth',1);
line([0,360],[0,0],'color','k','linewidth',1);
line([360,360],[0,90],'color','k','linewidth',1);

scatter(ret1(:,1), ret1(:,2))

n = 120, distances = 18; KL = 18.0668

去掉任意一个星 distance不变 kl = 18.05834851010958

n = 80, distances = 22.5; KL = 18.551770401636414

去掉任意一个星 distance不变 kl = 18.53911217378831

n = 48 dis = 30; KL = 19.187483254228805

19.16620665848412

---

试试三角形

clear all;  close all;clc;
jiange = 22.5;
ret = [];
x = 0;
y = 0;
while x < 360 && y <= 90
	if x >= 0
	    ret = [ret; x, y];
 	end
    x = x + jiange;
    if ( x >= 360 ) 
        x = x - 360 -jiange/2.0;
        y = y + jiange/2.0*sqrt(3);
    end
end

存在这种情况

---

试试六边形

clear all;  close all;clc;
jiange = 18;
ret = [];
x = 0;
y = 0;
while x < 360 && y <= 90
	if x >= 0
	    ret = [ret; x, y];
 	end
    x = x + jiange;
    if ( x >= 360 ) 
        x = x - 360 -jiange/2.0;
        y = y + jiange/2.0*sqrt(3);
    end
end
ret1 = [];
for i = 1 : 120
    if ret(i,2) > 31 && ret(i,2) < 32
        if mod(i, 3) ~= 1
            ret1 = [ret1;ret(i,:)];
        end
    elseif ret(i,2) < 31
        if mod(i, 3) ~= 0
            ret1 = [ret1;ret(i,:)];
        end
    elseif ret(i,2) > 60
        if mod(i, 3) ~= 2
            ret1 = [ret1;ret(i,:)];
        end
    elseif ret(i,2) > 32
        if mod(i, 3) ~= 1
            ret1 = [ret1;ret(i,:)];
        end
    end
end
% for i = 1 : 120
%     if mod(i, 3) ~= 0
%         ret1 = [ret1;ret(i,:)];
%     end
% 
% end

clear all;  close all;clc;
jiange = 22.5;
ret = [];
x = 0;
y = 0;
while x < 360 && y <= 90
	if x >= 0
	    ret = [ret; x, y];
 	end
    x = x + jiange;
    if ( x >= 360 ) 
        x = x - 360 -jiange/2.0;
        y = y + jiange/2.0*sqrt(3);
    end
end
ret1 = [];
for i = 1 : 80
    if ret(i,2) < 3
        if mod(i, 3) ~= 0
            ret1 = [ret1;ret(i,:)];
        end
    elseif ret(i,2) > 19 && ret(i,2) < 20
        if mod(i, 3) ~= 2
            ret1 = [ret1;ret(i,:)];
        end
    elseif ret(i,2) < 39
        if mod(i, 3) ~= 2
            ret1 = [ret1;ret(i,:)];
        end
    elseif ret(i,2) > 60
        if mod(i, 3) ~= 1
            ret1 = [ret1;ret(i,:)];
        end
    elseif ret(i,2) > 32
        if mod(i, 3) ~= 1
            ret1 = [ret1;ret(i,:)];
        end
    end
end









jiange = 30;
ret = [];
x = 0;
y = 0;
while x < 360 && y <= 90
	if x >= 0
	    ret = [ret; x, y];
 	end
    x = x + jiange;
    if ( x >= 360 ) 
        x = x - 360 -jiange/2.0;
        y = y + jiange/2.0*sqrt(3);
    end
end
ret2 = [];
for i = 1 : 48
    if ret(i,2) < 3
        if mod(i, 3) ~= 0
            ret2 = [ret2;ret(i,:)];
        end
    elseif ret(i,2) < 39
        if mod(i, 3) ~= 1
            ret2 = [ret2;ret(i,:)];
        end
    elseif ret(i,2) < 60
        if mod(i, 3) ~= 0
            ret2 = [ret2;ret(i,:)];
        end
    else 
        if mod(i, 3) ~= 1
            ret2 = [ret2;ret(i,:)];
        end
    end
end

81个点 17.6717

随机点

tic
ret = [];
sum = 0;
len = 100000;
for i = 1 : len
    x = 360.*rand(48, 1);
    y = 90.*rand(48, 1);
    sum = sum + Copy_3_of_totalH([x,y]);
end
av = sum / len
toc

随机加扰动

tic
% ret = [];
sum = 0;
len = 100000;
len1 = 48;
minv = 60;
maxv = 0;
for i = 1 : len
    ret1 = ret;
    ret1 = ret1 + 5*rand(len1,2)-2.5;
    tmp = Copy_3_of_totalH(ret1);
    minv = min(minv, tmp);
    maxv = max(maxv, tmp);
    sum = sum + tmp;
end
av = sum / len
toc

随机加点

tic
jiange = 22.5;
ret = [];
x = 0;
y = 0;
while x < 360 && y <= 90
    ret = [ret; x, y];
 
    y = y + jiange;
    if ( y > 90 ) 
        y = 0;
        x = x + jiange;
    end
end

sum = 0;
len = 100000;
minv = 60;
maxv = 0;
for i = 1 : len
    ret1 = ret;
    x = 360 * rand;
    y = 90 * rand;
    ret1 = [ret1;x,y];
    tmp = Copy_3_of_totalH(ret1);
    minv = min(minv, tmp);
    maxv = max(maxv, tmp);
    sum = sum + tmp;
end
av = sum / len
toc

---

画图 随机扰动

hold on;
plot([48,80,120],[13.7714,13.7110,13.6724],'-ksquare')
hold on
plot([48,80,120],[13.7436,13.6732,13.6244],'-o','Color',[0 0.6 0.1])
hold on
plot([48,80,120],[13.6305,13.5185,13.4272],'-^','Color',[0.5 0 0.8])
hold on
%%
xlim([43, 125])
line([48-1.000,48+1.000],[13.7311,13.7311],'Color',[0 0.6 0.1],'linewidth',1);
hold on;
line([48-1.000,48+1.000],[13.7549,13.7549],'color',[0 0.6 0.1],'linewidth',1);
hold on;
line([48,48],[13.7311,13.7549],'color',[0 0.6 0.1],'linewidth',1);

line([80-1.000,80+1.000],[13.6591,13.6591],'color',[0 0.6 0.1],'linewidth',1);
hold on;
line([80-1.000,80+1.000],[13.6844,13.6844],'color',[0 0.6 0.1],'linewidth',1);
hold on;
line([80,80],[13.6844,13.6591],'color',[0 0.6 0.1],'linewidth',1);

line([120-1.000,120+1.000],[13.6115,13.6115],'color',[0 0.6 0.1],'linewidth',1);
hold on;
line([120-1.000,120+1.000],[13.6364,13.6364],'color',[0 0.6 0.1],'linewidth',1);
hold on;
line([120,120],[13.6115,13.6364],'color',[0 0.6 0.1],'linewidth',1);

%%
line([48-1.000,48+1.000],[13.5590,13.5590],'color',[0.5 0 0.8],'linewidth',1);
hold on;
line([48-1.000,48+1.000],[13.6914,13.6914],'color',[0.5 0 0.8],'linewidth',1);
hold on;
line([48,48],[13.6914,13.5590],'color',[0.5 0 0.8],'linewidth',1);

line([80-1.000,80+1.000],[13.4549,13.4549],'color',[0.5 0 0.8],'linewidth',1);
hold on;
line([80-1.000,80+1.000],[13.5855,13.5855],'color',[0.5 0 0.8],'linewidth',1);
hold on;
line([80,80],[13.5855,13.4549],'color',[0.5 0 0.8],'linewidth',1);

line([120-1.000,120+1.000],[13.3566,13.3566],'color',[0.5 0 0.8],'linewidth',1);
hold on;
line([120-1.000,120+1.000],[13.4927,13.4927],'color',[0.5 0 0.8],'linewidth',1);
hold on;
line([120,120],[13.3566,13.4927],'color',[0.5 0 0.8],'linewidth',1);

box on;
legend('正方形网格点分布','网格点附近扰动±1°','网格点附近扰动±5°')
xlabel('n')
ylabel('KL熵');
xticks([48 80 120]);

hold on;
plot([48,80,120],[13.7714,13.7110,13.6724],'-ksquare')
hold on
plot([48,80,120],[13.7501,13.6984,13.6639],'-o','Color',[0 0.6 0.1])
hold on
plot([48,80,120],[13.6847,13.6586,13.6372],'-^','Color',[0.5 0 0.8])
hold on
%%
xlim([43, 125])
line([48-1.000,48+1.000],[13.1914,13.1914],'Color',[0 0.6 0.1],'linewidth',1);
hold on;
line([48-1.000,48+1.000],[13.7214,13.7214],'color',[0 0.6 0.1],'linewidth',1);
hold on;
line([48,48],[13.1914,13.7214],'color',[0 0.6 0.1],'linewidth',1);

line([80-1.000,80+1.000],[13.3497,13.3497],'color',[0 0.6 0.1],'linewidth',1);
hold on;
line([80-1.000,80+1.000],[13.6806,13.6806],'color',[0 0.6 0.1],'linewidth',1);
hold on;
line([80,80],[13.3497,13.6806],'color',[0 0.6 0.1],'linewidth',1);

line([120-1.000,120+1.000],[13.4646,13.4646],'color',[0 0.6 0.1],'linewidth',1);
hold on;
line([120-1.000,120+1.000],[13.6520,13.6520],'color',[0 0.6 0.1],'linewidth',1);
hold on;
line([120,120],[13.4646,13.6520],'color',[0 0.6 0.1],'linewidth',1);

%%

box on;
legend('正方形网格点分布','网格点附近扰动±1°','网格点附近扰动±5°')
xlabel('n')
ylabel('KL熵');
xticks([48 80 120]);

hold on;
plot([48,80,120],[13.7714,13.7110,13.6724],'-ksquare')
hold on
plot([48,80,120],[13.6630,13.6456,13.6286],'-v','Color',[0.6 0.1 0.1])
hold on

%%
line([48-1.000,48+1.000],[13.1646,13.1646],'Color',[0.6 0.1 0.1],'linewidth',1);
hold on;
line([48-1.000,48+1.000],[13.7697,13.7697],'color',[0.6 0.1 0.1],'linewidth',1);
hold on;
line([48,48],[13.1646,13.7697],'color',[0.6 0.1 0.1],'linewidth',1);

line([80-1.000,80+1.000],[13.3706,13.3706],'color',[0.6 0.1 0.1],'linewidth',1);
hold on;
line([80-1.000,80+1.000],[13.7095,13.7095],'color',[0.6 0.1 0.1],'linewidth',1);
hold on;
line([80,80],[13.3706,13.7095],'color',[0.6 0.1 0.1],'linewidth',1);

line([120-1.000,120+1.000],[13.4553,13.4553],'color',[0.6 0.1 0.1],'linewidth',1);
hold on;
line([120-1.000,120+1.000],[13.6715,13.6715],'color',[0.6 0.1 0.1],'linewidth',1);
hold on;
line([120,120],[13.4553,13.6715],'color',[0.6 0.1 0.1],'linewidth',1);



box on;
legend('正方形网格点分布','随机替换一个点')
xlabel('n')
ylabel('KL熵');
xticks([48 80 120]);

随机替换一个点

clear all;  close all;clc;tic;
jiange = 22.5;
ret = [];
x = 0;
y = 0;
while x < 360 && y <= 90
    ret = [ret; x, y];
 
    y = y + jiange;
    if ( y > 90 ) 
        y = 0;
        x = x + jiange;
    end
end
ret(end,:)=[];

 
 
 
sum = 0;
len = 100000;
minv = 60;
maxv = 0;
for i = 1 : len
    ret1 = ret;
    x = 360 * rand;
    y = 90 * rand;
    ret1 = [ret1;x,y];
    tmp = Copy_3_of_totalH(ret1);
    minv = min(minv, tmp);
    if ( tmp > maxv ) 
    	maxv = tmp;
    	x
    	y
    end

    sum = sum + tmp;
end
av = sum / len
toc

我本章的实验基于的是同n值的，因为n值大，KL熵小 ，不同n值不能说明均匀性。

---

正三角形网格点旋转

clear all;  close all;clc;
jiange = 22.5;
rec = [];
x = 0;
y = 0;
while x < 360 && y <= 90
	if x >= 0
	    rec = [rec; x, y];
 	end
    x = x + jiange;
    if ( x >= 360 ) 
        x = x - 360 -jiange/2.0;
        y = y + jiange/2.0*sqrt(3);
    end
end
d = 15;
rec1(:,1) = rec(:,1)*cosd(d) - rec(:,2)*sind(d);
rec1(:,2) = rec(:,1)*sind(d) + rec(:,2)*cosd(d);

subplot(2,1,1)
axis equal;
xticks([0 60 120 180 240 300 360])
yticks([0 30 60 90])
xlim([-30 400]);
ylim([-30 140]);
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;
hold on;
xt=get(gca,'xtick');
for k=1:numel(xt);
	xt1{k}=sprintf('%d°',xt(k));
end
set(gca,'xticklabel',xt1);
yt=get(gca,'ytick');
for k=1:numel(yt);
	yt1{k}=sprintf('%d°',yt(k));
end
set(gca,'yticklabel',yt1);
xlabel('方位轴');
ylabel('高度轴');
line([0,360],[90,90],'color','k','linewidth',1);
line([0,0],[0,90],'color','k','linewidth',1);
line([0,360],[0,0],'color','k','linewidth',1);
line([360,360],[0,90],'color','k','linewidth',1);
scatter(rec1(:,1), rec1(:,2))

jiange = 22.5;
rec2 = [];
x = 0;
y = -jiange *sqrt(3) ;
while x <= 360+jiange && y <= 90
	if x >= 0
	    rec2 = [rec2; x, y];
 	end
    x = x + jiange;
    if ( x >= 360+jiange ) 
        x = x-360 -jiange/2.0 - jiange;
        y = y + jiange/2.0*sqrt(3);
    end
end
d = 15;
rec3(:,1) = rec2(:,1)*cosd(d) - rec2(:,2)*sind(d);
rec3(:,2) = rec2(:,1)*sind(d) + rec2(:,2)*cosd(d);  


subplot(2,1,2)
axis equal;
xticks([0 60 120 180 240 300 360])
yticks([0 30 60 90])
xlim([-30 400]);
ylim([-30 140]);
% grid on;
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;
hold on;

xt=get(gca,'xtick');
for k=1:numel(xt);
	xt1{k}=sprintf('%d°',xt(k));
end
set(gca,'xticklabel',xt1);
yt=get(gca,'ytick');
for k=1:numel(yt);
	yt1{k}=sprintf('%d°',yt(k));
end
set(gca,'yticklabel',yt1);


xlabel('方位轴');
ylabel('高度轴');
line([0,360],[90,90],'color','k','linewidth',1);
line([0,0],[0,90],'color','k','linewidth',1);
line([0,360],[0,0],'color','k','linewidth',1);
line([360,360],[0,90],'color','k','linewidth',1);

scatter(rec3(:,1), rec3(:,2))

---

clear all;  close all;clc;
jiange = 22.5;
ret = [];
x = 0;
y = 0;
while x < 360 && y <= 90
	if x >= 0
	    ret = [ret; x, y];
 	end
    x = x + jiange;
    if ( x >= 360 ) 
        x = x - 360 -jiange/2.0;
        y = y + jiange/2.0*sqrt(3);
    end
end

clf;
% box on;
axis equal;
xticks([0 60 120 180 240 300 360])
yticks([0 30 60 90])
xlim([-10 370]);
ylim([-10 100]);
% grid on;
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;
hold on;

xt=get(gca,'xtick');
for k=1:numel(xt);
	xt1{k}=sprintf('%d°',xt(k));
end
set(gca,'xticklabel',xt1);
yt=get(gca,'ytick');
for k=1:numel(yt);
	yt1{k}=sprintf('%d°',yt(k));
end
set(gca,'yticklabel',yt1);


xlabel('方位轴');
ylabel('高度轴');
line([0,360],[90,90],'color','k','linewidth',1);
line([0,0],[0,90],'color','k','linewidth',1);
line([0,360],[0,0],'color','k','linewidth',1);
line([360,360],[0,90],'color','k','linewidth',1);

scatter(ret(:,1), ret(:,2))

选取点建KD树， 对600个星点找最近邻 ，算出每个星的权值，再遍历一遍每个星，为每个星的空间不断替换最优值。遍历一遍也行

从哪个点找哪个点更快

---

clear all;  close all;clc;
tic;
M = csvread('J2000_5_hor.csv');
points = M(:, 2:4);

jiange = 360.0/9;
ret = [];
x = 0;
y = 0;
while x < 360 && y <= 90
	if x >= -0.001
	    ret = [ret; x, y];
 	end
    x = x + jiange;
    if ( x >= 359 ) 
        x = x - 360 -jiange/2.0;
        y = y + jiange/2.0*sqrt(3);
    end
end

n = length(ret);
tmp = ones(n,1)*360;
rev = ones(n,2)*0;

Mdl = KDTreeSearcher(ret);
kk = 0;
for i = 1 : length(points)
    [nn,d] = knnsearch(Mdl,points(i,1:2),'k',kk+1);
    if d < tmp(nn)
        tmp(nn) = d;
        rev(nn,:) = points(i,1:2);
    end
end




clf;
% box on;
axis equal;
xticks([0 60 120 180 240 300 360])
yticks([0 30 60 90])
xlim([-10 370]);
ylim([-10 100]);
% grid on;
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;
hold on;

xt=get(gca,'xtick');
for k=1:numel(xt);
	xt1{k}=sprintf('%d°',xt(k));
end
set(gca,'xticklabel',xt1);
yt=get(gca,'ytick');
for k=1:numel(yt);
	yt1{k}=sprintf('%d°',yt(k));
end
set(gca,'yticklabel',yt1);


xlabel('方位轴');
ylabel('高度轴');
line([0,360],[90,90],'color','k','linewidth',1);
line([0,0],[0,90],'color','k','linewidth',1);
line([0,360],[0,0],'color','k','linewidth',1);
line([360,360],[0,90],'color','k','linewidth',1);

scatter(ret(:,1), ret(:,2))
scatter(points(:,1), points(:,2),5)
scatter(rev(:,1), rev(:,2),50,'g')
toc;	

纯距离

共27个

随机KL熵 11.7958

距离最近加权 13.4318 历时 0.078488 秒。

星等加权 随机 13.4006 ，总权值 -247.0324

优化后 13.5629 总权值 -195 最好的位移是 32.5583 6.7 历时25s

粒子50 迭代30 13.5537 -190.6025 位移30 7 30.2979600360942 7.02950188386670 25s

粒子15迭代10 13.5537 -190.6025 位移30.1858 6.9537 历时 2.705201 秒。

-196.325020765185 -196.325020765185 -196.325020765185 -196.325020765185 -191.733361560326 -190.655450867045 -190.655450867045 -190.655450867045 -190.655450867045 -190.655450867045

暴力点 13.5557 -190.6024 30.299 7.036 300粒子500迭代 历时 2638.095214 秒。

模拟退火最大化 13.6253 历时 74.406457 秒。

模拟点 13.7632

---

权值 = -距离 + -星等 ，星等越亮，越小，-星等越大，距离越近，越好，-距离越大，越好 权值越大

Mdl = KDTreeSearcher(ret);
kk = 0;
for i = 1 : length(points)
    [nn,d] = knnsearch(Mdl,points(i,1:2),'k',kk+1);
    value = -d - points(i,3);
    if value > tmp(nn)
        tmp(nn) = value;
        rev(nn,:) = points(i,1:2);
    end
end

toc;
Copy_3_of_totalH(rev)

clear all;  close all;clc;
tic;
M = csvread('J2000_5_hor.csv');
points = M(:, 2:4);
jiange = 360.0/10;
ret = [];
x = 0;
y = 0;
while x < 360 && y <= 90
	if x >= -0.0000001
	    ret = [ret; x, y];
 	end
    x = x + jiange;
    if ( x >= 359.9999999 ) 
        x = x - 360 -jiange/2.0;
        y = y + jiange/2.0*sqrt(3);
    end
end

N = 50;                           % 初始种群个数
d = 2;                             % 空间维数（参看上述的函数表达式）
ger = 30;                         % 最大迭代次数     
plimit = [0,jiange;0,90-ret(end)];              % 设置位置参数限制(矩阵的形式可以多维)，现在2X2矩阵
vlimit = [-1.5, 1.5;-1.5, 1.5];    % 设置速度限制
w = 0.8;                           % 惯性权重,个体历史成绩对现在的影响0.5~1之间
c1 = 0.5;                          % 自我学习因子
c2 = 0.5;                          % 群体学习因子 
for i = 1:d
    xx(:,i) = plimit(i, 1) + (plimit(i, 2) - plimit(i, 1)) * rand(N, 1);%初始种群的位置
end                              %rand(N,1)产生N行一列范围在1之内的随机数
                              %第一列,第二列：x=0+（20-0）*（1之内的随机数）                                    
v = rand(N, d);                   % 初始种群的速度,500行2列分别在两个维度上
xm = xx;                           % 每个个体的历史最佳位置
ym = zeros(1, d);                 % 种群的历史最佳位置，两个维度，设置为0
fxm = ones(N, 1) * (-inf);                % 每个个体的历史最佳适应度，设置为0
fym = -inf;                       % 种群历史最佳适应度,求最大值先设置成负无穷
kk = 0;
iter = 1;  
n = length(ret);
tmp = ones(n,1)*(-500);
rev = ones(n,2)*(0);

fx = ones(n,1);

while iter <= ger

%     fx = f(xx(:,1),xx(:,2));         % 代入x中的二维数据，算出个体当前适应度,为500行1列的数据   

    for i = 1:N                   %对每一个个体做判断
        ret1 = ret + xx(i,:);
        ret1(ret1(:,1) > 360) = ret1(ret1(:,1) > 360) - 360;
        
        Mdl = KDTreeSearcher(ret1);
        tmp = ones(n,1)*(-500);
        rev = ones(n,2)*(0);
        for ii = 1 : length(points)
          
            [nn,ddd] = knnsearch(Mdl,points(ii,1:2),'k',kk+1);
            value = -ddd - points(ii,3);
            if value > tmp(nn)
                tmp(nn) = value;
                rev(nn,:) = points(ii,1:2);
            end
        end
        fx(i) = sum(tmp);

        if fxm(i) < fx(i)           %如果每个个体的历史最佳适应度小于个体当前适应度
            fxm(i) = fx(i);         % 更新个体历史最佳适应度,第一轮就是把小于零的清除
            xm(i,:) = xx(i,:);       % 更新个体历史最佳位置
        end 
     end
 

    if fym < max(fxm)                  %种群历史最佳适应度小于个体里面最佳适应度的最大值
            [fym, nmax] = max(fxm);    % 更新群体历史最佳适应度,取出最大适应度的值和所在行数即位置
            ym = xm(nmax, :);          % 更新群体历史最佳位置
    end
    
     v = v * w + c1 * rand *(xm - xx) + c2 * rand *(repmat(ym, N, 1) - xx); % 速度更新公式,repmat函数把ym矩阵扩充成N行1列
    
     %%边界速度处理
     for i=1:d 
            for j=1:N
            if  v(j,i)>vlimit(i,2);      %如果速度大于边界速度，则把速度拉回边界
                v(j,i)=vlimit(i,2);
            end
            if  v(j,i) < vlimit(i,1)     %如果速度小于边界速度，则把速度拉回边界
                v(j,i)=vlimit(i,1);
            end
            end
     end      
       
     xx = xx + v;                          % 位置更新
     
     
      for i=1:d 
            for j=1:N
            if  xx(j,i)>plimit(i,2)
                xx(j,i)=plimit(i,2);
            end
            if  xx(j,i) < plimit(i,1)
                xx(j,i)=plimit(i,1);
            end
            end
      end
      
      record(iter) = fym;            %记录最大值
      
%    if times >= 10
%         cla;                    %清除轴线图形
%         mesh(x0_1, x0_2, y0);
%         plot3(x(:,1),x(:,2),f(x(:,1),x(:,2)), 'ro');title('状态位置变化');
%         pause(0.5);
%        times=0;
%     end
    iter = iter+1;
end

toc;

ret2 = ret + ym;
ret2(ret2(:,1) > 360) = ret2(ret2(:,1) > 360) - 360;
Mdl = KDTreeSearcher(ret2);
tmp = ones(n,1)*(-500);
rev = ones(n,2)*(0);
for ii = 1 : length(points)
    [nn,ddd] = knnsearch(Mdl,points(ii,1:2),'k',kk+1);
    value = -ddd - points(ii,3);
    if value > tmp(nn)
        tmp(nn) = value;
        rev(nn,:) = points(ii,1:2);
    end
end
aaan = sum(tmp);

Copy_3_of_totalH(rev)



clf;

p1 = scatter(ret2(:,1), ret2(:,2))
hold on
p2 = scatter(points(:,1), points(:,2),5)
hold on
p3 = scatter(rev(:,1), rev(:,2),50,'g')
hold on


% box on;
axis equal;
xticks([0 60 120 180 240 300 360])
yticks([0 30 60 90])
xlim([-10 370]);
ylim([-10 120]);
% grid on;
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;
hold on;

xt=get(gca,'xtick');
for k=1:numel(xt);
	xt1{k}=sprintf('%d°',xt(k));
end
set(gca,'xticklabel',xt1);
yt=get(gca,'ytick');
for k=1:numel(yt);
	yt1{k}=sprintf('%d°',yt(k));
end
set(gca,'yticklabel',yt1);


xlabel('方位轴');
ylabel('高度轴');
line([0,360],[90,90],'color','k','linewidth',1);
line([0,0],[0,90],'color','k','linewidth',1);
line([0,360],[0,0],'color','k','linewidth',1);
line([360,360],[0,90],'color','k','linewidth',1);

legend([p1, p2,p3],{'三角形网格点','上半天区所有点','距离最近点'})

距离-星等加权法

基于正三角形网格点的均匀选取方法，寻求变换

基于粒子群算法的最大化总体评价权值

总体评价权值V

综合评价权值

贡献度

---

计算实时最远距离

dddd = 0;
tmp =[];
iii = 0;
tic;
Mdl = KDTreeSearcher(rev);

for ii = 1 : length(points)
    [nn,ddd] = knnsearch(Mdl, points(ii,1:2),'k',kk+1);
	if ( ddd > points(ii,1) )
    	[~, dd] = knnsearch(Mdl,[points(ii,1) + 360.0, points(ii, 2)],'k',1);
    	ddd = min(ddd, dd);
    end
  	if (ddd > 360 - points(ii,1) )
    	[~, dd] = knnsearch(Mdl,[points(ii,1) - 360.0, points(ii, 2)],'k',1);
    	ddd = min(ddd, dd);
	end
    if ddd > dddd
        dddd = ddd;
        iii = ii;
        tmp = points(ii,1:2);
    end
end
toc;
hold on;
p4 = scatter(tmp(1), tmp(2),'k');

legend([p1, p2,p3, p4],{'三角形网格点','上半天区所有点','距离最近点','增补的一个星'})

pdist([137.3420  , 87.7710;117.085,73.8926],'euclidean')
pdist([137.3420  , 87.7710;158.055,72.5033],'euclidean')

---

15(cosd(weidu) cosd(A) / tand(90-h) + sind(weidu))

---

重新做一下星图软件

假设五等星 1471个星

纯计算方位角 需要10ms，显示成表格 12ms 显示出星图 229ms （加上地平380ms） 加个判断显示在地平上的 95ms

显示表格+视赤道 14ms 235ms

6等星 纯计算方位角 需要30ms，显示成表格 36ms 显示出星图

---

方位角误差

% function  requation()

clc;clear;
load data1.csv

azimuth=data1(:,1)/180*pi; %#ok<*NODEF>
altitude=data1(:,2)/180*pi;
deltaA=0-data1(:,3)/3600.0; % °
deltaH=0-data1(:,4)/3600.0; 

X=[0.*altitude+1 0.*altitude -cos(azimuth).*tan(altitude) -sin(azimuth).*tan(altitude) sec(altitude) -tan(altitude) sin(azimuth) cos(azimuth) 0.*altitude sin(2.*azimuth) cos(2.*azimuth) 0.*altitude 0.*altitude 0.*altitude 0.*altitude sin(2.*azimuth).*sec(altitude) cos(2.*azimuth).*sec(altitude) 0.*altitude 0.*altitude];
y=[0.*altitude 0.*altitude+1 sin(azimuth) -cos(azimuth) 0.*altitude 0.*altitude 0.*altitude 0.*altitude cot(altitude) 0.*altitude 0.*altitude sin(azimuth) cos(azimuth) altitude.*sin(azimuth) altitude.*cos(azimuth) 0.*altitude 0.*altitude sin(2.*azimuth) cos(2.*azimuth)];

A=[X;y];
Y=[deltaA;deltaH]; % °

%Res 为方程系数 Y° = A(弧度) * Res -> 
Res = inv(A'*A)*A'*Y;
%E 为残差
E = (A*Res - Y)*3600;
%MSE 为平均均方误差
RMS = sqrt(mse(E))
%输出 到TXT文件 科学计算法
%save ResultXISHU.txt Res -ascii; 
file1=fopen('ResultXISHU.txt','w');
for i=1:length(Res)
   fprintf(file1,'%5.5f\r\n',Res(i));
end
EA = E(1:64);
EE = E(65:end);
RMSA = sqrt(mse(EA));

RMSE = sqrt(mse(EE));
% aa = 

fclose(file1);
%输出 到TXT文件 科学计算法
%save ResultRMS.txt RMS -ascii; 

% end


for i = 1 :64
    if data1(i,1) < 0 
        data1(i,1) = data1(i,1) + 360;
    end
end


subplot(2,1,1,'Position',[0.1 0.5 0.8 0.5]);
title('模型校正前的方位角偏差')
grid on;
box on; 
axis equal;
xlim([0 360]);
ylim([-80 50]);
xticks([0 60 120 180 240 300 360])
% yticks([0 30 60 90])
xlabel('方位轴');
ylabel('误差（角秒）');


ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;


hold on;

scatter(data1(:,1), data1(:,3), 28)

% scatter(points(best_way,1), points(best_way,2), 20,'filled','MarkerFaceColor','#FF8000')
% xtickformat('usd')
% legend('All stars', 'Selected stars')
xtickformat('degrees')
% ytickformat('degrees')


subplot(2,1,2,'Position',[0.1 0.05 0.8 0.45]);
title('模型校正后的方位角残差')
box on; 
axis equal;
xlim([0 360]);
ylim([-50 50]);
xticks([0 60 120 180 240 300 360])
yticks([-50 -40 -30 -20 -10 0 10 20 30 40 50])
xlabel('方位轴');
ylabel('残差（角秒）');



ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;


hold on;

scatter(data1(:,1), EA, 28)

% scatter(points(best_way,1), points(best_way,2), 20,'filled','MarkerFaceColor','#FF8000')
% xtickformat('usd')
% legend('All stars', 'Selected stars')
xtickformat('degrees')
% ytickformat('degrees')

% function  requation()

clc;clear;
load data1.csv

azimuth=data1(:,1)/180*pi; %#ok<*NODEF>
altitude=data1(:,2)/180*pi;
deltaA=0-data1(:,3)/3600.0; % °
deltaH=0-data1(:,4)/3600.0; 

X=[0.*altitude+1 0.*altitude -cos(azimuth).*tan(altitude) -sin(azimuth).*tan(altitude) sec(altitude) -tan(altitude) sin(azimuth) cos(azimuth) 0.*altitude sin(2.*azimuth) cos(2.*azimuth) 0.*altitude 0.*altitude 0.*altitude 0.*altitude sin(2.*azimuth).*sec(altitude) cos(2.*azimuth).*sec(altitude) 0.*altitude 0.*altitude];
y=[0.*altitude 0.*altitude+1 sin(azimuth) -cos(azimuth) 0.*altitude 0.*altitude 0.*altitude 0.*altitude cot(altitude) 0.*altitude 0.*altitude sin(azimuth) cos(azimuth) altitude.*sin(azimuth) altitude.*cos(azimuth) 0.*altitude 0.*altitude sin(2.*azimuth) cos(2.*azimuth)];

A=[X;y];
Y=[deltaA;deltaH]; % °

%Res 为方程系数 Y° = A(弧度) * Res -> 
Res = inv(A'*A)*A'*Y;
%E 为残差
E = (A*Res - Y)*3600;
E = E / 2.0;
%MSE 为平均均方误差
RMS = sqrt(mse(E))
%输出 到TXT文件 科学计算法
%save ResultXISHU.txt Res -ascii; 
file1=fopen('ResultXISHU.txt','w');
for i=1:length(Res)
   fprintf(file1,'%5.5f\r\n',Res(i));
end
EA = E(1:64);
EE = E(65:end);

RMSA = sqrt(mse(EA));

RMSE = sqrt(mse(EE));
% aa = 

fclose(file1);
%输出 到TXT文件 科学计算法
%save ResultRMS.txt RMS -ascii; 

% end


for i = 1 :64
    if data1(i,1) < 0 
        data1(i,1) = data1(i,1) + 360;
    end
end


subplot(2,1,1,'Position',[0.1 0.5 0.8 0.5]);
title('模型校正前的方位角偏差','FontWeight','bold','FontSize',12)
grid on;
box on; 
axis equal;
xlim([0 360]);
ylim([-80 50]);
xticks([0 60 120 180 240 300 360])
% yticks([0 30 60 90])
xlabel('方位轴');
ylabel('误差（角秒）');


ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;


hold on;

scatter(data1(:,1), data1(:,3), 28)

% scatter(points(best_way,1), points(best_way,2), 20,'filled','MarkerFaceColor','#FF8000')
% xtickformat('usd')
% legend('All stars', 'Selected stars')
xtickformat('degrees')
% ytickformat('degrees')


subplot(2,1,2,'Position',[0.1 0.1 0.8 0.3]);
title('模型校正后的方位角残差','FontWeight','bold','FontSize',12)
box on; 
% axis equal;
xlim([0 360]);
ylim([-10 10]);
xticks([0 60 120 180 240 300 360])
yticks([-50 -40 -30 -20 -10 0 10 20 30 40 50])
xlabel('方位轴');
ylabel('残差（角秒）');



ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):2:ax.YAxis.Limits(2);
grid minor;


hold on;

scatter(data1(:,1), EA, 28)

% scatter(points(best_way,1), points(best_way,2), 20,'filled','MarkerFaceColor','#FF8000')
% xtickformat('usd')
% legend('All stars', 'Selected stars')
xtickformat('degrees')
% ytickformat('degrees')

---

% function  requation()

clc;clear;
load data1.csv

azimuth=data1(:,1)/180*pi; %#ok<*NODEF>
altitude=data1(:,2)/180*pi;
deltaA=0-data1(:,3)/3600.0; % °
deltaH=0-data1(:,4)/3600.0; 

X=[0.*altitude+1 0.*altitude -cos(azimuth).*tan(altitude) -sin(azimuth).*tan(altitude) sec(altitude) -tan(altitude) sin(azimuth) cos(azimuth) 0.*altitude sin(2.*azimuth) cos(2.*azimuth) 0.*altitude 0.*altitude 0.*altitude 0.*altitude sin(2.*azimuth).*sec(altitude) cos(2.*azimuth).*sec(altitude) 0.*altitude 0.*altitude];
y=[0.*altitude 0.*altitude+1 sin(azimuth) -cos(azimuth) 0.*altitude 0.*altitude 0.*altitude 0.*altitude cot(altitude) 0.*altitude 0.*altitude sin(azimuth) cos(azimuth) altitude.*sin(azimuth) altitude.*cos(azimuth) 0.*altitude 0.*altitude sin(2.*azimuth) cos(2.*azimuth)];

A=[X;y];
Y=[deltaA;deltaH]; % °

%Res 为方程系数 Y° = A(弧度) * Res -> 
Res = inv(A'*A)*A'*Y;
%E 为残差
E = (A*Res - Y)*3600;
%MSE 为平均均方误差
RMS = sqrt(mse(E))
%输出 到TXT文件 科学计算法
%save ResultXISHU.txt Res -ascii; 
file1=fopen('ResultXISHU.txt','w');
for i=1:length(Res)
   fprintf(file1,'%5.5f\r\n',Res(i));
end
EA = E(1:64);
EE = E(65:end);
RMSA = sqrt(mse(EA));

RMSE = sqrt(mse(EE));
% aa = 

fclose(file1);
%输出 到TXT文件 科学计算法
%save ResultRMS.txt RMS -ascii; 

% end


% for i = 1 :64
%     if data1(i,1) < 0 
%         data1(i,1) = data1(i,1) + 360;
%     end
% end


subplot(2,1,1,'Position',[0.1 0.55 0.8 0.4]);
title('模型校正前的高度角偏差')
grid on;
box on; 
% axis equal;
xlim([0 90]);
ylim([0 300]);
xticks([0 15 30 45 60 75 90]);
% yticks([0 30 60 90])
xlabel('高度角');
ylabel('误差（角秒）');


ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):5:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):5:ax.YAxis.Limits(2);
grid minor;


hold on;

scatter(data1(:,2), data1(:,4), 28)

% scatter(points(best_way,1), points(best_way,2), 20,'filled','MarkerFaceColor','#FF8000')
% xtickformat('usd')
% legend('All stars', 'Selected stars')
xtickformat('degrees')
% ytickformat('degrees')


subplot(2,1,2,'Position',[0.1 0.1 0.8 0.35]);
title('模型校正后的高度角残差')
box on; 
% axis equal;
xlim([0 90]);
ylim([-20 20]);
xticks([0 15 30 45 60 75 90]);
% yticks([0 30 60 90])
xlabel('高度角');
ylabel('残差（角秒）');



ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):5:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):5:ax.YAxis.Limits(2);
grid minor;


hold on;

scatter(data1(:,2), EE, 28)

% scatter(points(best_way,1), points(best_way,2), 20,'filled','MarkerFaceColor','#FF8000')
% xtickformat('usd')
% legend('All stars', 'Selected stars')
xtickformat('degrees')
% ytickformat('degrees')

% function  requation()

clc;clear;
load data2.csv

azimuth=data2(:,1)/180*pi; %#ok<*NODEF>
altitude=data2(:,2)/180*pi;
deltaA=0-data2(:,3)/3600.0; % °
deltaH=0-data2(:,5)/3600.0; 

X=[0.*altitude+1 0.*altitude -cos(azimuth).*tan(altitude) -sin(azimuth).*tan(altitude) sec(altitude) -tan(altitude) sin(azimuth) cos(azimuth) 0.*altitude sin(2.*azimuth) cos(2.*azimuth) 0.*altitude 0.*altitude 0.*altitude 0.*altitude sin(2.*azimuth).*sec(altitude) cos(2.*azimuth).*sec(altitude) 0.*altitude 0.*altitude];
y=[0.*altitude 0.*altitude+1 sin(azimuth) -cos(azimuth) 0.*altitude 0.*altitude 0.*altitude 0.*altitude cot(altitude) 0.*altitude 0.*altitude sin(azimuth) cos(azimuth) altitude.*sin(azimuth) altitude.*cos(azimuth) 0.*altitude 0.*altitude sin(2.*azimuth) cos(2.*azimuth)];

A=[X;y];
Y=[deltaA;deltaH]; % °

%Res 为方程系数 Y° = A(弧度) * Res -> 
Res = inv(A'*A)*A'*Y;
%E 为残差
E = (A*Res - Y)*3600;
E = E/2.0;
%MSE 为平均均方误差
RMS = sqrt(mse(E))
%输出 到TXT文件 科学计算法
%save ResultXISHU.txt Res -ascii; 
file1=fopen('ResultXISHU.txt','w');
for i=1:length(Res)
   fprintf(file1,'%5.5f\r\n',Res(i));
end
EA = E(1:64);
EE = E(65:end);
RMSA = sqrt(mse(EA));

RMSE = sqrt(mse(EE));
% aa = 

fclose(file1);
%输出 到TXT文件 科学计算法
%save ResultRMS.txt RMS -ascii; 

% end


% for i = 1 :64
%     if data2(i,1) < 0 
%         data2(i,1) = data2(i,1) + 360;
%     end
% end


subplot(2,1,1,'Position',[0.1 0.55 0.8 0.4]);
title('模型校正前的高度角偏差','FontWeight','bold','FontSize',12)
grid on;
box on; 
% axis equal;
xlim([0 90]);
ylim([-100 200]);
xticks([0 15 30 45 60 75 90]);
% yticks([0 30 60 90])
xlabel('高度角');
ylabel('误差（角秒）');


ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):5:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;


hold on;

scatter(data2(:,2), data2(:,5), 28)

% scatter(points(best_way,1), points(best_way,2), 20,'filled','MarkerFaceColor','#FF8000')
% xtickformat('usd')
% legend('All stars', 'Selected stars')
xtickformat('degrees')
% ytickformat('degrees')


subplot(2,1,2,'Position',[0.1 0.1 0.8 0.33]);
title('模型校正后的高度角残差','FontWeight','bold','FontSize',12)
box on; 
% axis equal;
xlim([0 90]);
ylim([-10 10]);
xticks([0 15 30 45 60 75 90]);
% yticks([0 30 60 90])
xlabel('高度角');
ylabel('残差（角秒）');



ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):5:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):5:ax.YAxis.Limits(2);
grid minor;


hold on;

scatter(data2(:,2), EE, 28)

% scatter(points(best_way,1), points(best_way,2), 20,'filled','MarkerFaceColor','#FF8000')
% xtickformat('usd')
% legend('All stars', 'Selected stars')
xtickformat('degrees')
% ytickformat('degrees')

% function  requation()

clc;clear;
load data2.csv

azimuth=data2(:,1)/180*pi; %#ok<*NODEF>
altitude=data2(:,2)/180*pi;
deltaA=0-data2(:,3)/3600.0; % °
deltaH=0-data2(:,5)/3600.0; 

X=[0.*altitude+1 0.*altitude -cos(azimuth).*tan(altitude) -sin(azimuth).*tan(altitude) sec(altitude) -tan(altitude) sin(azimuth) cos(azimuth) 0.*altitude sin(2.*azimuth) cos(2.*azimuth) 0.*altitude 0.*altitude 0.*altitude 0.*altitude sin(2.*azimuth).*sec(altitude) cos(2.*azimuth).*sec(altitude) 0.*altitude 0.*altitude];
y=[0.*altitude 0.*altitude+1 sin(azimuth) -cos(azimuth) 0.*altitude 0.*altitude 0.*altitude 0.*altitude cot(altitude) 0.*altitude 0.*altitude sin(azimuth) cos(azimuth) altitude.*sin(azimuth) altitude.*cos(azimuth) 0.*altitude 0.*altitude sin(2.*azimuth) cos(2.*azimuth)];

A=[X;y];
Y=[deltaA;deltaH]; % °

%Res 为方程系数 Y° = A(弧度) * Res -> 
Res = inv(A'*A)*A'*Y;
%E 为残差
E = (A*Res - Y)*3600;
%MSE 为平均均方误差
RMS = sqrt(mse(E))
%输出 到TXT文件 科学计算法
%save ResultXISHU.txt Res -ascii; 
file1=fopen('ResultXISHU.txt','w');
for i=1:length(Res)
   fprintf(file1,'%5.5f\r\n',Res(i));
end
EA = E(1:64);
EE = E(65:end);
RMSA = sqrt(mse(EA));

RMSE = sqrt(mse(EE));
% aa = 

fclose(file1);
%输出 到TXT文件 科学计算法
%save ResultRMS.txt RMS -ascii; 

% end


for i = 1 :64
    if data2(i,1) < 0 
        data2(i,1) = data2(i,1) + 360;
    end
end


subplot(2,1,1,'Position',[0.1 0.55 0.8 0.4]);
title('指向校正前高度角误差分布图',FontWeight='bold')
grid on;
box on; 
% axis equal;
xlim([-80 80]);
ylim([-100 200]);
xticks([-80 -60 -40 -20  0  20  40 60 80])
% yticks([0 30 60 90])
xlabel('校正前方位角偏差（角秒）');
ylabel('校正前高度角偏差（角秒）');


ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;
r = pdist([[0,0];[50,0]],'euclidean');
p1 = rectangle('position',[0-r,0-r,2*r,2*r],'curvature',[1,1],'LineStyle','--');

hold on;
r = pdist([[0,0];[150,0]],'euclidean');
p2 = rectangle('position',[0-r,0-r,2*r,2*r],'curvature',[1,1],'LineStyle','-.');


hold on;

scatter(data2(:,3), data2(:,5), 28)

% scatter(points(best_way,1), points(best_way,2), 20,'filled','MarkerFaceColor','#FF8000')
% xtickformat('usd')
% legend('All stars', 'Selected stars')
% xtickformat('degrees')
% ytickformat('degrees')
a = plot([[1000,2000],[1000,2000]],'--k');
hold on;
b = plot([[1000,2000],[1000,2000]],'-.k');
hold on
legend([a, b],{'50 arcsec','150 arcsec'})

subplot(2,1,2,'Position',[0.1 0.1 0.8 0.35]);
title('指向校正后高度角残差分布图',FontWeight='bold')
box on; 
% axis equal;
xlim([-60 60]);
ylim([-20 20]);
xticks([-80 -60 -40 -20  0  20  40 60 80])
% yticks([0 30 60 90])
xlabel('校正后方位角残差（角秒）');
ylabel('校正后高度角残差（角秒）');



ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;

r = pdist([[0,0];[10,0]],'euclidean');
p1 = rectangle('position',[0-r,0-r,2*r,2*r],'curvature',[1,1],'LineStyle','--');

hold on;
r = pdist([[0,0];[15,0]],'euclidean');
p2 = rectangle('position',[0-r,0-r,2*r,2*r],'curvature',[1,1],'LineStyle','-.');

hold on;



scatter(EA, EE, 28)


% scatter(points(best_way,1), points(best_way,2), 20,'filled','MarkerFaceColor','#FF8000')
% xtickformat('usd')
% legend('All stars', 'Selected stars')
% xtickformat('degrees')
% ytickformat('degrees')
hold on; 
a = plot([[100,200],[100,200]],'--k');
hold on;
b = plot([[100,200],[100,200]],'-.k');
hold on
legend([a, b],{'10 arcsec','15 arcsec'})

---

clear; 
% M = csvread('J2000_5hor.csv');
MM = csvread('J2000_5horzhengfangxing.csv');

% MM = sortrows(M,4);
% MMM = MM(22:85,:);

clf;
% box on;
axis equal;
xticks([0 60 120 180 240 300 360])
yticks([0 30 60 90])
xlim([-10 370]);
ylim([-10 100]);
% grid on;
ax = gca;
ax.XAxis.MinorTick = 'on';
ax.XAxis.MinorTickValues = ax.XAxis.Limits(1):10:ax.XAxis.Limits(2);
ax.YAxis.MinorTick = 'on';
ax.YAxis.MinorTickValues = ax.YAxis.Limits(1):10:ax.YAxis.Limits(2);
grid minor;
hold on;

xt=get(gca,'xtick');
for k=1:numel(xt);
	xt1{k}=sprintf('%d°',xt(k));
end
set(gca,'xticklabel',xt1);
yt=get(gca,'ytick');
for k=1:numel(yt);
	yt1{k}=sprintf('%d°',yt(k));
end
set(gca,'yticklabel',yt1);


xlabel('方位轴','FontWeight','bold');
ylabel('高度轴','FontWeight','bold');
line([0,360],[90,90],'color','k','linewidth',1.5);
line([0,0],[0,90],'color','k','linewidth',1.5);
line([0,360],[0,0],'color','k','linewidth',1.5);
line([360,360],[0,90],'color','k','linewidth',1.5);

scatter(MM(:,2), MM(:,3))
Copy_3_of_totalH(MM(:,2:3))