地球上的椭圆与直线的交点问题

re-us

以地球上的某一点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 编辑 ]

ChaChing

从未使用过axes2ecc/ellipse1这些函数, 刚刚google下才学习到是Mapping Toolbox的东西!
office的版本无此工具箱, 无法试! 仅说说想法
楼主有先以x的关系筛选过, 为何不再加上y的关系筛选?
还有可能专业不了, 可否直接算出直线方程及椭圆方程, 再求解?
因都未试过, 或许不可行

re-us

1. 我感觉x和y筛选只能用一个，不然去掉的点太多了。
2.能够算出椭圆方程没有问题，但是就是不知道怎么考虑纬度效应。比如说：一个以50km为半径的圆，位于赤道和位于北纬60度，所跨越的经度是不一样的（不知道说清楚了没有）

谢谢

ChaChing

1.不解, 去掉的点愈多, 不是愈快速?
2.抱歉, 专业问题个人未研究过!

re-us

嗯，只用x有剩下500个点，用两个的话有200点剩余。
多谢，多谢。