声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1996|回复: 4

[编程技巧] 地球上的椭圆与直线的交点问题

[复制链接]
发表于 2010-2-3 05:30 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
以地球上的某一点A(x1,y1)为椭圆的中心,半长轴是100km,半短轴是50km,倾斜角是315度。另外一个点是B(x2,y2),求过A(x1,y1)和B(x2,y2)的直线和该椭圆的交点(该直线的方向是从A到B,所以就一个交点)。注意:x1,y1,x2,y2是经度,纬度。

我的做法是:用ellipse1函数产生1000个点,用这1000个点代表这个椭圆。然后循环这1000个点与A的斜率 slope1,A与B的斜率slope2,这两个斜率最接近的那个点,就是交点。但是这样耗费很多时间,如果想精确的话,估计要产生更多的点(比如100000)。这样的好处就是ellipse1已经考虑了纬度效应,就是低纬度的100km和高纬度的100km代表不同度数。

怎么更快速和准确的找出这个交点?谢谢

我的代码:
clc; clear
%the ellipse
const=2*pi*6371/360; x1=35; y1=35; smaj=100/const; smin=50/const;
ecc=axes2ecc(smaj,smin); [elat,elon]=ellipse1(y1,x1,[smaj,ecc],360-45,[],[],[],1000);
figure; plot(elon,elat); hold on; plot(x1,y1,'+'); text(x1+0.03,y1,'A')
x2=35.5;y2=35.5; plot(x2,y2,'+'); text(x2+0.03,y2,'B')

index=sign(x2-x1)==sign(elon-x1); elon1=elon(index); elat1=elat(index);

%calculate the slope
slope1=(elat1-y1)./(elon1-x1); slope2=(y2-elat1)./(x2-elon1);
[I1,I2]=min(abs(slope2-slope1)); y3=elat1(I2(1)); x3=elon1(I2(1));
plot(x3,y3,'+m'); line([x1,x2],[y1,y2])

[ 本帖最后由 ChaChing 于 2010-2-3 10:46 编辑 ]
1.jpg
回复
分享到:

使用道具 举报

发表于 2010-2-3 10:43 | 显示全部楼层
从未使用过axes2ecc/ellipse1这些函数, 刚刚google下才学习到是Mapping Toolbox的东西!
office的版本无此工具箱, 无法试! 仅说说想法
楼主有先以x的关系筛选过, 为何不再加上y的关系筛选?
还有可能专业不了, 可否直接算出直线方程及椭圆方程, 再求解?
因都未试过, 或许不可行
 楼主| 发表于 2010-2-3 22:30 | 显示全部楼层
1. 我感觉x和y筛选只能用一个,不然去掉的点太多了。
2.能够算出椭圆方程没有问题,但是就是不知道怎么考虑纬度效应。比如说:一个以50km为半径的圆,位于赤道和位于北纬60度,所跨越的经度是不一样的(不知道说清楚了没有)

谢谢
发表于 2010-2-3 23:13 | 显示全部楼层
1.不解, 去掉的点愈多, 不是愈快速?
2.抱歉, 专业问题个人未研究过!
 楼主| 发表于 2010-2-5 02:55 | 显示全部楼层
嗯,只用x有剩下500个点,用两个的话有200点剩余。
多谢,多谢。
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-29 19:32 , Processed in 0.059818 second(s), 21 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表