木星大红斑位置

2009年9月12日晚天气晴好。使用FS78+30毫米目镜+ToUCam拍摄木星。恰逢大红斑位于视面之上。

利用2009年9月12日22时26分40秒的图片进行测量。

首先确定木星轮廓：

可以计算出扁率为： 5.54%

测量大红斑位置，转换成木星表面坐标（中央经线设为0度）：经度：-8.9448度， 纬度: -22.1141度

大红斑位于南纬20度，应用木星经度II系统。此系统的自转周期为9小时55分40秒。照片拍摄时刻的中央经度为147.9度，对应大红斑的经度为139度。大红斑的位置并不是固定的，以后再次测量时可以和这个数值作比较。

Matlab程序，放在这，以后用得着：

%木星坐标测量
clear
close all

%边界亮度阈值
cri=8;

%读入图像
[im,fn]=cloadimage;
img=sign(double(rgb2gray(im))-cri);
img=(img+abs(img))/2*255;
[h,w]=size(img);

%寻找边界
idx=1;
for i=(50:w-50)
lookforone=1;
for j=(50:h-50)
if(lookforone==1)
if(img(j,i)>128)
lookforone=0;
x(idx)=i;
y(idx)=j;
idx=idx+1;
continue;
end
else
if(img(j,i)<128)
x(idx)=i;
y(idx)=j;
idx=idx+1;
break;
end
end
end
end

%椭圆拟合: ax^2+bxy+cy^2+dx+ey+f=0
x=x';
y=y';
x2=x.*x;
y2=y.*y;
xy=x.*y;

A=[x2,xy,y2,x,y];
B=x*0-1;

v=A\B;

[hx,wx]=size(x);

%计算椭圆中心
A=[2*v(1),v(2);v(2),2*v(3)];
B=-[v(4);v(5)];

center=A\B;

%显示

figure(1)

image(im);
hold on
axis equal
plot(center(1),center(2),'*r');

%计算椭圆上的点，用于验证。其中一点坐标用于下步
for i=(1:hx)
x0=x(i);
a=v(3);
b=v(2)*x0+v(5);
c=v(1)*x0*x0+v(4)*x0+1;

if(b*b<4*a*c)
continue;
end
y0=(-b+sqrt(b*b-4*a*c))/2/a;
%plot(x0,y0,'b*');
y0=(-b-sqrt(b*b-4*a*c))/2/a;
%plot(x0,y0,'b*');
end

ptx=x0;
pty=y0;

if(1==0)
%横向和纵向极值点连线
y0=(1:h);
x0=(v(2)*y0+v(4))/(-2*v(1));
%plot(x0,y0);

x0=(1:w);
y0=(v(2)*x0+v(5))/(-2*v(3));
%plot(x0,y0);
end

%倾角，如果图像未反转，不加pi
ang=atan(v(2)/(v(1)-v(3)))/2+pi;

disp(['倾角：',num2str(ang*180/pi),' 度']);

%木星坐标到图片坐标的转换矩阵，转置为逆矩阵
T=[cos(ang),-sin(ang);sin(ang),cos(ang)];

%椭圆扁率
a=v(1)*(cos(ang))^2+v(3)*(sin(ang))^2+v(2)*sin(ang)*cos(ang);
b=v(1)*(sin(ang))^2+v(3)*(cos(ang))^2-v(2)*sin(ang)*cos(ang);
b=sqrt(a/b);
disp(['扁率: ',num2str(100-b*100),'%']);

x0=ptx-center(1);
y0=pty-center(2);

v0=T'*([x0;y0]);
%椭圆长短轴
a=sqrt(v0(1)^2+(v0(2)/b)^2);
b=a*b;

clear x y

%画椭圆
for i=(1:360)
th=i/180*pi;
x0=a*cos(th);
y0=b*sin(th);

v1=T*([x0;y0]);

x(i)=v1(1)+center(1);
y(i)=v1(2)+center(2);
end
plot(x,y,'y')

%测量位置坐标

[xr,yr] = GINPUT(1);

v0=T'*([xr-center(1);yr-center(2)]);

rx=sqrt(1-(v0(2)/b)^2)*a;

lat=-atan(v0(2)/rx)*180/pi;
lon=-asin(v0(1)/rx)*180/pi;

disp(['经度：',num2str(lon),'度 纬度: ',num2str(lat),'度']);

%展开图

for i=(1:181)
an=(i-91)*pi/180;
rx=sqrt((cos(an))^2/((cos(an)/a)^2+(sin(an)/b)^2));
y0=sqrt(1-(rx/a)^2)*b;

if(an<0)
y0=-y0;
end


for j=(1:180)

x0=rx*sin((j-90)*pi/180);

v0=T*([x0;y0]);

imprj(i,j,:)=im(round(v0(2)+center(2)),round(v0(1)+center(1)),:);

end
end

figure(2)
image(imprj)

imwrite(imprj,[fn,'_prj.bmp'],'bmp');