马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
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 编辑 ] |