? ? ? ?大三的時(shí)候做的一個(gè)東西，當(dāng)時(shí)在上海天文臺(tái)實(shí)習(xí)了15天，在老師要求下寫了下面一大段代碼，現(xiàn)在想起來(lái)，都不記得這是干的什么了。

? ? ?? 人最可悲的，就是很多東西本來(lái)學(xué)會(huì)了，結(jié)果沒(méi)過(guò)多久就忘了。

?????? 我只記得，當(dāng)時(shí)是美國(guó)的一顆衛(wèi)星，我們登了它的網(wǎng)站，下載了它的觀測(cè)數(shù)據(jù)，經(jīng)過(guò)篩選，大概有4萬(wàn)多個(gè)點(diǎn)。然后用這些數(shù)據(jù)導(dǎo)入進(jìn)來(lái)，根據(jù)一些計(jì)算，求出絕對(duì)星等，以及距離，然后繪制南極星空分布圖，用的是極坐標(biāo)。最后證實(shí)了，恒星確實(shí)是呈網(wǎng)狀分布。

?????? 據(jù)說(shuō)，恒星的分布就像絲襪一樣：相近的恒星喜歡聚成團(tuán)，然后各個(gè)星系團(tuán)中間是巨大的空洞。

clc;
clear;
format long
Um=0.3;
Uk=0;
Uv=0.7;???????? %基本參數(shù)賦值
z=csvread('shixi_gliming.csv',0,6,'G1..G44103');%導(dǎo)入紅移z
for i=1:44103
DC(i)=quad('integrand',0,z(i),0.00001);?? %計(jì)算積分
??? DM(i)=DC(i);????? %因?yàn)閗=0
?DL(i)=(1+z(i))*DM(i);??????????? %計(jì)算光度距離
end
dered_r=csvread('shixi_gliming.csv',0,4,'E1..E44103');?? %導(dǎo)入視星等r
for i=1:44103
??? if z(i)<0.1
??????? p0=-0.194858;
??????? p1=0.912570;
??? else
??????? p0=-0.242659;
??????? p1=1.38731;
??? end
???? M(i)=dered_r(i)-5*log10(DL(i))-25-(p0+p1*z(i))-1.62*(z(i)-0.1);%計(jì)算絕對(duì)星等M
end
?csvwrite('M',M);
dered_g=csvread('shixi_gliming.csv',0,3,'D1..D44103');?? %導(dǎo)入視星等g
col=dered_g-dered_r;
figure(1);
plot(z,M,'r.','MarkerSize',5);
legend('絕對(duì)星等對(duì)紅移');xlabel('紅移z');ylabel('絕對(duì)星等M');?? %絕對(duì)星等對(duì)紅移
figure(2);
plot(M,col,'r.','MarkerSize',5);axis([-25 -10 -2 4]);
legend('顏色對(duì)絕對(duì)星等');xlabel('絕對(duì)星等M');ylabel('顏色c');?? %顏色對(duì)星等
pause;
i=1;n=1;p=1;q=1;s=1;u=1;v=1;l=1;
ra=csvread('shixi_gliming.csv',0,1,'B1..B44103');?? %導(dǎo)入經(jīng)度
for k=1:44103
??????? if M(k)>-24&M(k)<=-23
??????????? ra1(i)=ra(k);z1(i)=z(k);col1(i)=col(k);M1(i)=M(k);
??????????? i=i+1;
??????? end
end
??????? figure(3);
?? theta1=ra1*3.14159/180;
??? rho1=z1;
??? polar(theta1,rho1,'r.');title('-24<M<=-23');
pause;???????????? %畫第三張圖
for k=1:44103????????
??????? if M(k)>-23&M(k)<=-22
??????????? ra2(l)=ra(k);z2(l)=z(k);col2(l)=col(k);M2(l)=M(k);
??????????? l=l+1;
??????? end
end
figure(4);
?? theta2=ra2*3.14159/180;
??? rho2=z2;
??? polar(theta2,rho2,'m.');title('-23<M<=-22');
pause;??????? %畫第四張圖
for k=1:44103??
??????? if M(k)>-22&M(k)<=-21
??????????? ra3(n)=ra(k);z3(n)=z(k);col3(n)=col(k);M3(n)=M(k);
??????????? n=n+1;
??????? end
end
figure(5);
?? theta3=ra3*3.14159/180;
??? rho3=z3;
??? polar(theta3,rho3,'y.');title('-22<M<=-21');
pause;??????????? %畫第五張圖
for k=1:44103??
??????? if M(k)>-21&M(k)<=-20
??????????? ra4(p)=ra(k);z4(p)=z(k);col4(p)=col(k);M4(p)=M(k);
??????????? p=p+1;
??????? end
end
figure(6);
?? theta4=ra4*3.14159/180;
??? rho4=z4;
??? polar(theta4,rho4,'g.');title('-21<M<=-20');
pause;?????????? %畫第六張圖
for k=1:44103?
???????? if M(k)>-20&M(k)<=-19
??????????? ra5(q)=ra(k);z5(q)=z(k);col5(q)=col(k);M5(q)=M(k);
??????????? q=q+1;
???????? end
end
figure(7);
?? theta5=ra5*3.14159/180;
??? rho5=z5;
??? polar(theta5,rho5,'c.');title('-20<M<=-19');
pause;???????? %畫第七張圖
for k=1:44103?
??????? if M(k)>-19&M(k)<=-18
??????????? ra6(s)=ra(k);z6(s)=z(k);col6(s)=col(k);M6(s)=M(k);
??????????? s=s+1;
??????? end
end
figure(8);
?? theta6=ra6*3.14159/180;
??? rho6=z6;
??? polar(theta6,rho6,'b.');title('-19<M<=-18');
pause;?????????? %畫第八張圖
for k=1:44103?
??????? if M(k)>-18&M(k)<=-17
??????????? ra7(u)=ra(k);z7(u)=z(k);col7(u)=col(k);M7(u)=M(k);
??????????? u=u+1;
??????? end
end
figure(9);
?? theta7=ra7*3.14159/180;
??? rho7=z7;
??? polar(theta7,rho7,'k.');title('-18<M<=-17');
pause;?????????? %畫第九張圖
for k=1:44103?
??????? if M(k)>-17&M(k)<=-16
??????????? ra8(v)=ra(k);z8(v)=z(k);col8(v)=col(k);M8(v)=M(k);
??????????? v=v+1;???
??? end
end
figure(10);
?? theta8=ra8*3.14159/180;
??? rho8=z8;
??? polar(theta8,rho8,'r.');title('-17<M<=-16');
? pause;??????????? %畫第十張圖
figure(11)
?polar(theta1,rho1,'r.');
?hold on
?polar(theta2,rho2,'m.');
?hold on
?polar(theta3,rho3,'y.');
?hold on
?polar(theta4,rho4,'g.');
?hold on
?polar(theta5,rho5,'c.');
?hold on
?polar(theta6,rho6,'b.');
?hold on?
?polar(theta7,rho7,'k.');
?hold on
?polar(theta8,rho8,'w.');title('-24<M<-16');%八張圖繪在一起
?pause;
figure(12);
i=size(z1);w=1;t=1;
for j=1:i(1,2)
???? if col1(j)>-0.1328*M1(j)-1.991
??????? theta11(w)=theta1(j);rho11(w)=rho1(j);w=w+1;
???? else?
???????? theta12(t)=theta1(j);rho12(t)=rho1(j);t=t+1;
???? end
end
? polar(theta11,rho11,'r.');
?hold on
?polar(theta12,rho12,'b.');
?pause;??? %畫第一幅圖
?figure(13);
?l=size(z2);w=1;t=1;
?for j=1:l(1,2)
???? if col2(j)>-0.1328*M2(j)-1.991
???????? theta21(w)=theta2(j);rho21(w)=rho2(j);w=w+1;
???? else?
???????? theta22(t)=theta2(j);rho22(t)=rho2(j);t=t+1;
???? end
?end
? polar(theta21,rho21,'r.');
?hold on
?polar(theta22,rho22,'b.');
?pause;? %畫第二幅圖
?figure(14);
? n=size(z3);w=1;t=1;
?for j=1:n(1,2)??
???? if col3(j)>-0.1328*M3(j)-1.991
???????? theta31(w)=theta3(j);rho31(w)=rho3(j);w=w+1;
???? else?
??????? theta32(t)=theta3(j);rho32(t)=rho3(j);t=t+1;
???? end?????
?end
?polar(theta31,rho31,'r.');
?hold on
?polar(theta32,rho32,'b.');
?pause;? %畫第三幅圖
?figure(15);
p=size(z4);w=1;t=1;
?for j=1:p(1,2)
???? if col4(j)>-0.1328*M4(j)-1.991
???????? theta41(w)=theta4(j);rho41(w)=rho4(j);w=w+1;
???? else?
???????? theta42(t)=theta4(j);rho42(t)=rho4(j);t=t+1;
???? end
?end
?polar(theta41,rho41,'r.');
?hold on
?polar(theta42,rho42,'b.');
?pause;? %畫第四幅圖
?figure(16);
?q=size(z5);w=1;t=1;
?for j=1:q(1,2)
???? if col5(j)>-0.1328*M5(j)-1.991
??????? theta51(w)=theta5(j);rho51(w)=rho5(j);w=w+1;?
???? else?
????? theta52(t)=theta5(j);rho52(t)=rho5(j);t=t+1;???
???? end
?end
?polar(theta51,rho51,'r.');
?hold on
?polar(theta52,rho52,'b.');
?pause;? %畫第五幅圖
?figure(17);
?s=size(z6);w=1;t=1;
?for j=1:s(1,2)
???? if col6(j)>-0.1328*M6(j)-1.991
??????? theta61(w)=theta6(j);rho61(w)=rho6(j);w=w+1;??
???? else?
?????? theta62(t)=theta6(j);rho62(t)=rho6(j);t=t+1;?????
???? end
?end
? polar(theta61,rho61,'r.');
?hold on
?polar(theta62,rho62,'b.');
?pause;? %畫第六幅圖
?figure(18);
? u=size(z7);w=1;t=1;
?for j=1:u(1,2)
???? if col7(j)>-0.1328*M7(j)-1.991
???????? theta71(w)=theta7(j);rho71(w)=rho7(j);w=w+1;??
???? else?
?????? theta72(t)=theta7(j);rho72(t)=rho7(j);t=t+1;?????
???? end
?end
? polar(theta71,rho71,'r.');
?hold on
?polar(theta72,rho72,'b.');
?pause;? %畫第七幅圖
?figure(19);
?v=size(z8); w=1;t=1;
?for j=1:v(1,2)
???? if col8(j)>-0.1328*M8(j)-1.991
??????? theta81(w)=theta8(j);rho81(w)=rho8(j);w=w+1;??
???? else?
???????? theta82(t)=theta8(j);rho82(t)=rho8(j);t=t+1;?
???? end
?end
? polar(theta81,rho81,'r.');
?hold on
?polar(theta82,rho82,'b.');
?pause;? %畫第八幅圖
?figure(20);
?w=1;t=1;
?for j=1:44103
???? if col(j)>-0.1328*M(j)-1.991
???????? theta_x(w)=ra(j)*3.14159/180;
???????? rho_x(w)=z(j);w=w+1;
???? else?
???????? theta_y(t)=ra(j)*3.14159/180;
???????? rho_y(t)=z(j);t=t+1;
???? end
?end
?polar(theta_x,rho_x,'r.');
?hold on
?polar(theta_y,rho_y,'b.');

另附函數(shù)文件：integrand.M

function y=integrand(x)
DH=3000/0.71;
Um=0.3;
Uk=0;
Uv=0.7;
y=DH*1./sqrt(Um*(1+x).^3+Uk*(1+x).^2+Uv);

總結(jié)

以上是生活随笔為你收集整理的上海天文台实习的一个项目-根据卫星数据绘制南极星空分布图的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。