之前提到过FITS文件的投影问题（http://blog.sciencenet.cn/blog-117333-360909.html），FITS头文件里会有关键词描述投影的方式。需要投影的原因在于天文观测面对的是分布在天球上的天体，而实际使用数据的时候，将数据展在平面上比较方便，这和地图是一个原理。投影的方式很多，关于地图投影，专门有书。具体到FITS文件的投影问题，也有专门的文章描述（Calabretta & Greisen, 2002, A&A, 395, 1077）。
      原本投影的问题不需要专门写下来，看书和看文章都可以，但是其实我经常遇到的只有几种，每次看文章都得重新回忆一遍，费时费力。所以在这里记一下，主要为我自己。
      FITS头文件中CTYPE1和CTYPE2的值就表示了投影的方式。最常遇到的有GLS、SIN、CAR、TAN和NCP。现在分别来说一下。以下直接给出（赤经RA，赤纬DEC）到投影坐标（x, y）转换公式。
      按照前面提到文章的说法，GLS已经被重命名为SFL了，也就是Sanson-Flamsteed投影。按照文章中的公式直接计算会有问题，实际的转换应该为
begin{equation}
y=DEC
end{equation}
begin{equation}
x=(RA-RA0)cos(DEC*pi/180)+RA0
end{equation}
其中因为参考点的赤纬为0，所以上面没有写出来。
       SIN表示的是Slant orthographic投影。上述文章我看不明白，借鉴Wolfram Mathworld中正交投影的表达式（http://mathworld.wolfram.com/OrthographicProjection.html）可以得到
begin{equation}
x=cos(DEC*pi/180)sin[(RA-RA0)*pi/180]/pi*180+RA0
end{equation}
begin{equation}
y=cos(DEC0*pi/180)sin(DEC*pi/180)/pi*180-sin(DEC0*pi/180)cos(DEC*pi/180)cos[(RA-RA0)*pi/180]/pi*180+DEC0
end{equation}
       CAR表示的是Plate carrée投影，其实就是根本不作处理，$x=RA$，$y=DEC$。
       TAN表示的是Gnomonic投影。上述文章我也看不明白，在Wolfram Mathworld(http://140.177.205.23/GnomonicProjection.html )可以找到相应的表达式。
begin{equation}
x=frac{cos(DEC*pi/180)sin[(RA-RA0)*pi/180]}{cos c}/pi*180+RA0
end{equation}
begin{equation}
y=frac{cos(DEC0*pi/180)sin(DEC*pi/180)-sin(DEC0*pi/180)cos(DEC*pi/180)cos[(RA-RA0)*pi/180]}{cos c}/pi*180+DEC0
end{equation}
其中
begin{equation}
cos c=sin(DEC0*pi/180)sin(DEC*pi/180)+cos(DEC0*pi/180)cos(DEC*pi/180)cos[(RA-RA0)*pi/180]
end{equation}
        至于NCP，其实就是SIN。
   

2013年8月5日补记：修改了TAN投影的公式，感谢许铎指出错误。

用wcsxy2ad.pro相对人性化一点。

PRO projection2

fitsname='t13_new.fits';default file name

head=headfits(fitsname)

fxaddpar,head,'CTYPE1','RA---SFL'

fxaddpar,head,'CTYPE2','DEC--SFL'

;fxaddpar,head,'CTYPE1','RA---TAN'

;fxaddpar,head,'CTYPE2','DEC--TAN'

ctype1  =   fxpar(head,'CTYPE1');

ctype2  =   fxpar(head,'CTYPE2');

extast,head,astr,noparams

print, astr,noparams

convert=!pi/180.0

convert=1.0

ra=65.0*convert

dec=23.0*convert

xr=65.0*convert

yd=23.0*convert

extast,head,astr,noparams

print, astr,noparams

;adxy,head,ra,dec,xr,yd

;print,xr,yd

;xyad,head,xr,yd,ra,dec

;print,ra,dec

;ad2xy,ra,dec,astr,xr,yd

;print,xr,yd

;xy2ad,xr,yd,astr,ra,dec

;print,ra,dec

;ad2xy,ra,dec,astr,xr,yd

;print,xr,yd

;xy2ad,xr,yd,astr,ra,dec

;print,ra,dec

wcsxy2sph,xr,yd,ra,dec,CTYPE=[ctype1,ctype2]

print,ra,dec

END