笔记（二）几个matlab脚本

经常需要做一些不同坐标系中坐标的变换，已经折腾了不止一两次了。把脚本放这儿，以备不时之需。

1. 赤道坐标变换为银道坐标

%north pole R.A. 12h 51m 26.282s , Dec. +27° 07′ 42.01″
%zero point R.A. 17h 45m 37.224s , Dec. -28° 56′ 10.23″
%Galactic Coordinate <-> Equitorial Coordinate
%Galactic pole (b=90 deg) RA=12.8573 h=192.859508 deg Dec=+27.128336 deg
%zero point (b=0 deg, l=0 deg) RA=17.7603 h=266.405100 deg Dec=-28.936175 deg
function V=ec2gc(Ra,Dec)
% Theta=(90.0-27.128336)./180.0.*pi;
% Phi=(192.859508-270.0)./180.0.*pi;
% Theta1=(90.0-(-28.936175))./180.0.*pi;
% Phi1=266.405100./180.0.*pi;
% xs=[sin(Theta1).*cos(Phi1), sin(Theta1).*sin(Phi1), cos(Theta1)];
% M1=[1 0 0; 0 cos(Theta) -sin(Theta); 0 sin(Theta) cos(Theta)];
% M2=[cos(Phi) -sin(Phi) 0; sin(Phi) cos(Phi) 0; 0 0 1];
% RotMz=M2*M1;
% xp=[RotMz*[1 0 0]']';
% zt=[RotMz*[0 0 1]']';
% Phip=asin(cross(xp,xs)*zt');
% M3=[cos(Phip) -sin(Phip) 0; sin(Phip) cos(Phip) 0; 0 0 1];
% RotM=RotMz*M3;
% Rotation Matrix for coordinate system:
RotM=[ -0.054875545667945   0.494110771211000  -0.867665384961804;
       -0.873437545087187  -0.444828615859559  -0.198076645126421;
       -0.483834196104114   0.746981959812779   0.455985112031682];
ra=Ra./24.0.*2.0.*pi;
dec=Dec./180.0.*pi;
H = [cos(ra).*cos(dec),sin(ra).*cos(dec),sin(dec)];
Pos=[RotM'*H']';
xp=Pos(1);
yp=Pos(2);
zp=Pos(3);
if abs(xp)<1.e-8
    if yp>0
        b=pi./2.0;
    end
    if yp<0 || yp==0
        b=3.0.*pi./2.0;
    end
else
    if yp>0
        if xp>0
           b=atan(yp./xp);
        else
           b=atan(yp./xp)+pi;
        end
    else
        if xp>0
           b=atan(yp./xp)+2.0.*pi;
        else
           b=atan(yp./xp)+pi;
        end
    end
end
l=asin(zp);
b=b./2.0./pi.*24.0;
l=l./pi.*180.0;
V=[b l];

2. 银道坐标变换为赤道坐标

%north pole R.A. 12h 51m 26.282s , Dec. +27° 07′ 42.01″
%zero point R.A. 17h 45m 37.224s , Dec. -28° 56′ 10.23″
%Galactic Coordinate <-> Equitorial Coordinate
%Galactic pole (b=90 deg) RA=12.8573 h=192.859508 deg Dec=+27.128336 deg
%zero point (b=0 deg, l=0 deg) RA=17.7603 h=266.405100 deg Dec=-28.936175 deg
function V=gc2ec(b,l)
% Theta=(90.0-27.128336)./180.0.*pi;
% Phi=(192.859508-270.0)./180.0.*pi;
% Theta1=(90.0-(-28.936175))./180.0.*pi;
% Phi1=266.405100./180.0.*pi;
% xs=[sin(Theta1).*cos(Phi1), sin(Theta1).*sin(Phi1), cos(Theta1)];
% M1=[1 0 0; 0 cos(Theta) -sin(Theta); 0 sin(Theta) cos(Theta)];
% M2=[cos(Phi) -sin(Phi) 0; sin(Phi) cos(Phi) 0; 0 0 1];
% RotMz=M2*M1;
% xp=[RotMz*[1 0 0]']';
% zt=[RotMz*[0 0 1]']';
% Phip=asin(cross(xp,xs)*zt');
% M3=[cos(Phip) -sin(Phip) 0; sin(Phip) cos(Phip) 0; 0 0 1];
% RotM=RotMz*M3;
% Rotation Matrix for coordinate vector:
RotM=[ -0.054875545667945   0.494110771211000  -0.867665384961804;
       -0.873437545087187  -0.444828615859559  -0.198076645126421;
       -0.483834196104114   0.746981959812779   0.455985112031682];
b=b./24.0.*2.0.*pi;
l=l./180.0.*pi;
H = [cos(b).*cos(l),sin(b).*cos(l),sin(l)];
Pos=[RotM*H']';
xp=Pos(1);
yp=Pos(2);
zp=Pos(3);
if abs(xp)<1.e-8
    if yp>0
        ra=pi./2.0;
    end
    if yp<0 || yp==0
        ra=3.0.*pi./2.0;
    end
else
    if yp>0
        if xp>0
           ra=atan(yp./xp);
        else
           ra=atan(yp./xp)+pi;
        end
    else
        if xp>0
           ra=atan(yp./xp)+2.0.*pi;
        else
           ra=atan(yp./xp)+pi;
        end
    end
end
dec=asin(zp);
Ra=ra./2.0./pi.*24.0;
Dec=dec./pi.*180.0;
V=[Ra Dec];

3. 赤道坐标变换为超星系坐标

%north pole R.A. 18h 55m 01s, Dec. +15deg 42' 32"
%zero point R.A. 2h 49m 14s, Dec. +59deg 31' 42"
%Supergalactic Coordinate <-> Equitorial Coordinate
%Supergalactic pole (SGB=90 deg) RA=18.9169 h=283.7535 deg Dec=+15.7089 deg
%zero point (SGB=0 deg, SGL=0 deg) RA=2.8206 h=42.3090 deg Dec=+59.5283 deg
function V=ec2sgc(Ra,Dec)
% Theta=(90.0-15.7089)./180.0.*pi;
% Phi=(283.7535-270.0)./180.0.*pi;
% Theta1=(90.0-59.5283)./180.0.*pi;
% Phi1=42.3090./180.0.*pi;
% xs=[sin(Theta1).*cos(Phi1), sin(Theta1).*sin(Phi1), cos(Theta1)];
% M1=[1 0 0; 0 cos(Theta) -sin(Theta); 0 sin(Theta) cos(Theta)];
% M2=[cos(Phi) -sin(Phi) 0; sin(Phi) cos(Phi) 0; 0 0 1];
% RotMz=M2*M1;
% xp=[RotMz*[1 0 0]']';
% zt=[RotMz*[0 0 1]']';
% Phip=asin(cross(xp,xs)*zt');
% M3=[cos(Phip) -sin(Phip) 0; sin(Phip) cos(Phip) 0; 0 0 1];
% RotM=RotMz*M3;
% Rotation Matrix for coordinate system:
RotM=[0.375022041659904  -0.898320159816291   0.228865372515963;
      0.341354889828154  -0.095717033759227  -0.935048174501584;
      0.861878940141622   0.428787989472616   0.270749981762484];
ra=Ra./24.0.*2.0.*pi;
dec=Dec./180.0.*pi;
Pos = [cos(ra).*cos(dec),sin(ra).*cos(dec),sin(dec)];
V=[RotM'*Pos']';

4. 超星系坐标变换为赤道坐标

%north pole R.A. 18h 55m 01s, Dec. +15deg 42' 32"
%zero point R.A. 2h 49m 14s, Dec. +59deg 31' 42"
%Supergalactic Coordinate <-> Equitorial Coordinate
%Supergalactic pole (SGB=90 deg) RA=18.9169 h=283.7535 deg Dec=+15.7089 deg
%zero point (SGB=0 deg, SGL=0 deg) RA=2.8206 h=42.3090 deg Dec=+59.5283 deg
function V=sgc2ec(x,y,z)
% Theta=(90.0-15.7089)./180.0.*pi;
% Phi=(283.7535-270.0)./180.0.*pi;
% Theta1=(90.0-59.5283)./180.0.*pi;
% Phi1=42.3090./180.0.*pi;
% xs=[sin(Theta1).*cos(Phi1), sin(Theta1).*sin(Phi1), cos(Theta1)];
% M1=[1 0 0; 0 cos(Theta) -sin(Theta); 0 sin(Theta) cos(Theta)];
% M2=[cos(Phi) -sin(Phi) 0; sin(Phi) cos(Phi) 0; 0 0 1];
% RotMz=M2*M1;
% xp=[RotMz*[1 0 0]']';
% zt=[RotMz*[0 0 1]']';
% Phip=asin(cross(xp,xs)*zt');
% M3=[cos(Phip) -sin(Phip) 0; sin(Phip) cos(Phip) 0; 0 0 1];
% RotM=RotMz*M3;
% Rotation Matrix for coordinate vector:
RotM=[0.375022041659904  -0.898320159816291   0.228865372515963;
      0.341354889828154  -0.095717033759227  -0.935048174501584;
      0.861878940141622   0.428787989472616   0.270749981762484];
H=[x y z]./sqrt(x.^2+y.^2+z.^2);
Pos=[RotM*H']';
xp=Pos(1);
yp=Pos(2);
zp=Pos(3);
if abs(xp)<1.e-8
    if yp>0
        ra=pi./2.0;
    end
    if yp<0
        ra=3.0.*pi./2.0;
    end
else
    if yp>0
       ra=atan(yp./xp)+pi./2.0;
    else
       ra=atan(yp./xp)+3.0.*pi./2.0;
    end
end
dec=asin(zp);
Ra=ra./2.0./pi.*24.0;
Dec=dec./pi.*180.0;
V=[Ra Dec];

转载本文请联系原作者获取授权，同时请注明本文来自钱磊科学网博客。
链接地址：https://blog.sciencenet.cn/blog-117333-301628.html

上一篇：『读文献』（一）在Ia型超新星遗迹中寻找前身星的伴星
下一篇：笔记（四）