地球上的椭圆与直线的交点问题
以地球上的某一点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); =ellipse1(y1,x1,,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);
=min(abs(slope2-slope1)); y3=elat1(I2(1)); x3=elon1(I2(1));
plot(x3,y3,'+m'); line(,)
[ 本帖最后由 ChaChing 于 2010-2-3 10:46 编辑 ] 从未使用过axes2ecc/ellipse1这些函数, 刚刚google下才学习到是Mapping Toolbox的东西!
office的版本无此工具箱, 无法试! 仅说说想法
楼主有先以x的关系筛选过, 为何不再加上y的关系筛选?
还有可能专业不了, 可否直接算出直线方程及椭圆方程, 再求解?
因都未试过, 或许不可行 1. 我感觉x和y筛选只能用一个,不然去掉的点太多了。
2.能够算出椭圆方程没有问题,但是就是不知道怎么考虑纬度效应。比如说:一个以50km为半径的圆,位于赤道和位于北纬60度,所跨越的经度是不一样的(不知道说清楚了没有)
谢谢 1.不解, 去掉的点愈多, 不是愈快速?
2.抱歉, 专业问题个人未研究过! 嗯,只用x有剩下500个点,用两个的话有200点剩余。
多谢,多谢。
页:
[1]