CSPICE星历程序简介和其与python的交互使用

星表与星历
CSPICE星历程序
- CSPICE+Python
应用

星表与星历

作为天文爱好者，大家都会经常玩一些星空模拟软件，对于这些软件来说，可视化效果是其最重要的指标之一。而所有这些程序的基础是一些星表文件和星历文件，程序所要做的就是把这些文件中的数据漂亮的展示出来。这些程序为我们提供了非常丰富的功能，却无法帮助我们实现我们所要的一切。比如stellarium是一款非常优秀的星空模拟软件，可以给出（合理范围内星表中提供的…）任意天体在天球上的位置，但是比如我想要得到一年以内的“太阳8字曲线”，那么….我可以固定视角调整日期截好多屏，然后用PS叠加起来（好吧我赢了），也可以直接拿到星历数据自己在二维平面上做球极投影。说以说要想真正实现自己想要的功能，还是得会处理星历问题，自己做可视化。

通常要回答：“大哥哥，这颗星星在哪里呀？”这种问题，我们要借助于星表或者星历。

通常的星表提供的是（距离我们较远的）恒星相对地球的角位置、距离、自行等信息，而星历则更注重“历元”，提供的是太阳系内天体（行星，小行星，卫星，人造卫星等）随时间变化的位置与速度信息。本文向大家介绍CSPICE星历程序与python的交互使用，最后再做一些作为天爱，掌握了使用星历数据的能力后想做的一些有趣的事情。

CSPICE星历程序

SPICE是NASA旗下的NAIF(Navigationand Ancillary Information Facility) 开发的一套星历程序,可以直接读取NASA提供的星历文件来为科学工作者们提供太阳系内各种星际探测任务中卫星和科学仪器的位置、速度、姿态等信息。如下图所示

naif

SPICE是Spacecraft, Planet, Instrument, Camera-matrix和Events几个单词的首字母缩写，表明了SPICE要为我们大家提供的内容包括：

SPK文件：轨道飞行器星历数据行星、卫星、彗星、小行星星历数据自定义的其他天体的与以上天体的相对位置

PCK文件：行星、卫星、彗小行星的空间指向、大小和形状信息以上天体的一些其他的物理模型(如引力势模型)

IK文件：仪器的视场、指向、形状大小等信息仪器的一些其他附加信息

CK文件：仪器所在平台的飞行姿态信息一些其他指定目标的空间指向信息

EK文件： ESP:科学仪器所在平台事件 ESQ:卫星、科学仪器指令事件 ENB:实验笔记”与仪器资料日志

这些文件也被称作Kernels，需要从NAIF的FTP（需翻墙…）下载

比如我想得到各个行星相对地球的位置信息，需要下载de405.bsp星历文件，而如果想得到HST的轨道、姿态等数据，需要到http://naif.jpl.nasa.gov/pub/naif/HST/kernels下载对应的文件。

有了对应的星历文件后，SPICE可以输出星历文件中所包含的那些天体，在星历文件包含的时间范围内的，在指定坐标系下的相对位置和其他一些参数。

这里面的“指定坐标系”都有那些呢？我们看图

哈哈，只要有相应的星历文件，几乎每一个探测器，连火星车的萌大眼睛都有一个坐标系，你可以让星历程序直接返回火星车所看到的世界。星历程序中提供的常用坐标框架包括

惯性坐标系：J2000等

非惯性坐标系：地球固定坐标系（z轴北极，x轴赤道与子午线的交点，右手系），火星固定坐标系等

但是没有我们所期望的某些地面观测站的坐标系（北京等…）

所以说我们只能得到相对于地固坐标系的坐标，然后自己完成坐标转换。

=======非程序员天爱请跳过这一节，直接到最后的实战========

CSPICE+Python

CSPICE星历程序和其与python的交互使用：

SPICE程序包括4个语言版本（C，FORTRAN，IDL，MATLAB）和两个正在开发的语言版本（Java Native and Python）

由于本人使用matplotlib进行作图，而python版的又正在开发之中，所以使用python调用C语言版的CSPICE程序。

本工作的依赖环境是：

python
numpy
matplotlib>=1.2.0
CSPICE
pyastroCoords
GNU make
gcc

其中pyAstroCoords是本人为实测天文学数据处理编写的一套Python天文坐标转换程序（还是自己写的用着舒服），它的DOC比程序长…请大家放心使用。

第一步：下载安装CSPICE

注意其安装脚本是用csh写的，请装csh…

第二步：编译C的动态链接库供Python调用

源码不太长，直接贴过来：

spiceSelflib.c

Make the height to be less than 400pixel



#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "SpiceUsr.h"
#define wordSize 40
#define fileSize 100
void loadFiles(char *setupFileName);
double utc2et(char* utcStr);
void utc2etArray(char* utcStart, char* utcEnd, int number, double* etArray);
void getPos(char* obv, char* target, char* frame, double et, double* posMat);
void getPosArray(char* obv, char* target, char* frame, double etArray[], int number, double posMat[][3]);
void et2utcArray(double* etArray, int num, char* format, int prec, char strOut[][40]);

void loadFiles(char *setupFileName) {
	furnsh_c(setupFileName);
	printf("\nLoad Kernel:%s OK!\n",setupFileName);
}

double utc2et(char* utcStr) {
	SpiceDouble et;
	utc2et_c(utcStr,&et);
	return et;
}

void et2utcArray(double* etArray, int num, char* format, int prec, char strOut[][40]) {
	int i;
	for (i = 0; i < num; i++) {
		et2utc_c(etArray[i], format, prec, 40, strOut[i]);
	}

}

void utc2etArray(char* utcStart, char* utcEnd, int number, double* etArray) {
	SpiceDouble etStart, etEnd, step;
	int i;

	utc2et_c(utcStart,&etStart);
	utc2et_c(utcEnd,&etEnd);
	step = (etEnd - etStart)/((double)(number-1));
	//etArray = (double*)malloc(number*sizeof(SpiceDouble));
	for (i = 0; i < number; i++) {
		etArray[i] = etStart + i*step;
	}
}

void getPos(char* obv, char* target, char* frame, double et, double* posMat) {
	SpiceDouble lt;
	SpiceChar abCorr[wordSize];

	sprintf(abCorr,"LT+S");
	spkpos_c(target, et, frame, abCorr, obv, posMat, &lt);
}

void getPosArray(char* obv, char* target, char* frame, double etArray[], int number, double posMat[][3]) {
	SpiceDouble lt;
	SpiceChar abCorr[wordSize];
	int i;

	sprintf(abCorr,"LT+S");
	for (i = 0; i < number; i++) {
		spkpos_c(target, etArray[i], frame, abCorr, obv, posMat[i], &lt);
	}
}

相对应的makefile如下，注意其中CSPICE头文件和静态链接库的位置

makefile

shell


#! /bin/bash
# makefile for the readKernels program
CSPICEPATH:=/home/wujinnnnn/softwares/cspice/# the root path you install your cspice softwares
INCLUDE:=$(CSPICEPATH)include/
LIBS:=$(CSPICEPATH)lib/
spiceSelfLib:  $(LIBS)cspice.a $(LIBS)csupport.a
	gcc -fPIC spiceSelfLib.c $(LIBS)cspice.a $(LIBS)csupport.a -o spiceSelfLib.so -shared -I/usr/include/python2.7 -I$(INCLUDE)

.PHONY: clean
clean:
	-rm *.so

第三步：编写Python接口程序，使用ctypes调用

spiceLib.py

python


import numpy as np
import pdb
import ctypes
import pyAstroCoords as ac
import matplotlib.pyplot as plt

# load kernels
spiceSelfLib = ctypes.cdll.LoadLibrary("./spiceSelfLib.so")
loadKernels = spiceSelfLib.loadFiles
loadKernels.argtypes = [ctypes.c_char_p]
loadKernels("setupFile.txt")

utc2et = spiceSelfLib.utc2et
utc2et.restype = ctypes.c_double
#test = utc2et('2010-07-27T12:23:43.23')

def utc2etArray(utcStart, utcEnd, number, rType='numpy'):
	utc2etArray_c = spiceSelfLib.utc2etArray
	etArray_c = (ctypes.c_double*number)()
	utc2etArray_c(utcStart, utcEnd, number, ctypes.byref(etArray_c))
	if rType=='numpy':
		etArray = np.array(etArray_c[:])
		return etArray
	elif rType=='c' or rType=='C':
		return etArray_c
	else:
		raise(Exception("Error return type!"))
#test = utc2etArray('2010-07-27T12:23:43.23', '2010-07-27T12:33:43.23', 1000)

def getPos(obv, target, frame, et, rType='numpy'):
	getPos_c = spiceSelfLib.getPos
	pos_c = (ctypes.c_double*3)()
	#getPos_c.argtypes=[ctypes.c_char_p, ctypes.c_char_p, ctypes.c_char_p, ctypes.c_double, pos_c]
	getPos_c(obv, target, frame, ctypes.c_double(et), ctypes.byref(pos_c))
	if rType=='numpy':
		pos = np.array(pos_c[:])
		return pos
	elif rType=='c' or rType=='C':
		return pos_c
	else:
		raise(Exception("Error return type!"))
#test = getPos('SSB', 'EARTH', 'J2000', utc2et('2013-07-27T12:30:00.00') )

def getPosArray(obv, target, frame, etArray_c, rType='numpy'):
	getPosArray_c = spiceSelfLib.getPosArray
	number = len(etArray_c[:])
	posMat_c = (ctypes.c_double*3*number)()
	#pdb.set_trace()
	getPosArray_c(obv, target, frame, ctypes.pointer(etArray_c), number, ctypes.byref(posMat_c))
	if rType=='numpy':
		posMat = np.mat(np.ctypeslib.as_array(posMat_c).T)
		return posMat
	elif rType=='c' or rType=='C':
		return posMat_c
	else:
		raise(Exception("Error return type!"))
#test = getPosArray('SSB', 'EARTH', 'J2000', utc2etArray('2010-07-27T12:23:43.23', '2010-07-27T12:33:43.23', 10, 'C'))

def StarFixed2LocalHorizontalFrame(long_dec, lat_dec, oldCoors, r=6372.797):
	m, n = oldCoors.shape
	if m!=3:
		raise(Exception('oldCoors Error!'))

	deltaVec = -np.mat(ac.sphere2Rect(np.array([r, long_dec, lat_dec]))).T
	M = ac.starFixed2LocalHorizontalMat(long_dec, lat_dec)
	delta = np.kron(np.ones((1,n)), deltaVec)
	newCoors = M * (oldCoors + delta)
	return newCoors

def stationShift(long_dec, lat_dec, oldCoors, r=6372.797):
	m, n = oldCoors.shape
	if m!=3:
		raise(Exception('oldCoors Error!'))

	deltaVec = -np.mat(ac.sphere2Rect(np.array([r, long_dec, lat_dec]))).T
	delta = np.kron(np.ones((1,n)), deltaVec)
	newCoors = (oldCoors + delta)
	return newCoors

def incAngle(mat1, mat2):
	m, n = mat1.shape
	if m!=3:
		raise(Exception('oldCoors Error!'))

	array1 = mat1.A
	array2 = mat2.A
	sqrted = np.sqrt(np.sum(array1*array1, axis=0) * np.sum(array2*array2,axis=0))
	multed = np.sum(array1 * array2, axis=0)
	angles = np.arccos(multed/sqrted)/np.pi*180.
	return angles

def et2utcArray(etArray, format, prec):
	''' format
			"C"      "1986 APR 12 16:31:09.814"
			"D"      "1986-102 // 16:31:12.814"
			"J"      "JD 2446533.18834276"
			"ISOC"   "1987-04-12T16:31:12.814"
			"ISOD"   "1987-102T16:31:12.814"
	'''
	n = etArray.size
	et2utcArray_c = spiceSelfLib.et2utcArray
	strOut = (ctypes.c_char*40*n)()
	et2utcArray_c(etArray.ctypes.data, n, format, prec, ctypes.byref(strOut))
	utcList = [strOut[:][i][:].split()[0] for i in range(n)]
	return np.array(utcList)

现在一切OK，可以一战了~

===============恭喜空降成功==================

应用

下面就看看我们能用星历程序做些什么事情吧～

1、太阳8字

太阳8字就是在每天的同一时刻拍太阳，一年以来其轨迹长得像个8字，我们就看看北京一年以来在中午12点左右的太阳8字图吧

python


import pyAstroCoords as ac
import matplotlib.pyplot as plt
from spiceLib import *

sunOrb = getPosArray('EARTH', 'SUN', 'IAU_EARTH', utc2etArray('2013-01-01T04:00:00', '2013-12-31T04:00:00', 365, 'C')) 
sunOrb_station = StarFixed2LocalHorizontalFrame(116.4, 39.9, sunOrb)
xx,yy = ac.StereographicProjection(180,90,1,"rect",sunOrb_station[0].A1, sunOrb_station[1].A1, sunOrb_station[2].A1)

ac.plotHorizon()
plt.scatter(xx,yy,s=2)
plt.show(0)

输出如下：

上北下南左东右西的天球投影图哦～

2、“预测”日食

我们既然已经有了太阳和月亮在我们眼中每个时刻的角位置，月亮视的大小大概0.5度，所以大致认为太阳和月亮的视位置夹角<0.25°且太阳在地平以上的时候发生了可见的日食。我们只需在时间线上扫描一下就行了…代码如下

python


y_c = utc2etArray('2000-01-01T04:00:00', '2013-12-31T04:00:00', 999999, 'C')
etArray = np.array(etArray_c[:])
sunOrb = getPosArray('EARTH', 'SUN', 'IAU_EARTH', etArray_c) 
sunOrb_station = stationShift(116.4, 39.9, sunOrb)
moonOrb = getPosArray('EARTH', 'MOON', 'IAU_EARTH', etArray_c)
moonOrb_station = stationShift(116.4, 39.9, moonOrb)
angles = incAngle(sunOrb_station, moonOrb_station)
inds = (angles<0.25)*(sunOrb_station[2,:].A1>0)
littleAngles = angles[inds]
times = et2utcArray(etArray[inds], 'ISOC', 3)
n = times.size
for i in range(n):
    print times[i], littleAngles[i]

输出结果如下：

python


2006-03-29T11:31:25.839 0.246793230066
2006-03-29T11:38:47.603 0.191481601497
2006-03-29T11:46:09.368 0.156433598549
2006-03-29T11:53:31.132 0.156905172331
2006-03-29T12:00:52.896 0.193391116388
2008-08-01T10:46:51.567 0.241370705366
2008-08-01T10:54:13.332 0.169313722133
2008-08-01T11:01:35.096 0.100649082644
2008-08-01T11:08:56.860 0.0541265150609
2008-08-01T11:16:18.624 0.0857005623547
2008-08-01T11:23:40.388 0.153792578973
2008-08-01T11:31:02.152 0.228313757749
2009-07-22T01:14:35.260 0.211695472114
2009-07-22T01:21:57.024 0.18195162594
2009-07-22T01:29:18.788 0.166216761301
2009-07-22T01:36:40.552 0.168020805425
2009-07-22T01:44:02.317 0.186484014454
2009-07-22T01:51:24.081 0.21710535425
2012-05-20T22:09:00.611 0.244895292214
2012-05-20T22:16:22.375 0.204921087569
2012-05-20T22:23:44.139 0.174888796591
2012-05-20T22:31:05.903 0.15979984084
2012-05-20T22:38:27.667 0.163175647615
2012-05-20T22:45:49.431 0.183421013565
2012-05-20T22:53:11.195 0.215350513449

3、可视化旅行者1、2号的航行轨迹

旅行者1、2号是美国在1977年发射的两颗“远征”的太空探测器，在一个宝贵的多行星引力加速的发射窗口时相继发射，目前两个探测器已行至太阳系边界，但他们究竟飞到哪了？就让我们来看一下

python


import numpy as np
import ctypes
import matplotlib as mpl
from mpl_toolkits.mplot3d import Axes3D
import mpl_toolkits as mplt
import  matplotlib.pyplot as plt
from spiceLib import *

loadKernels('voyagerSetup.dat')

mercuryOrb = getPosArray('SSB', 'MERCURY', 'ECLIPJ2000', utc2etArray('1977-10-01T00:00:00', '2012-12-31T23:59:59.99', 10000, 'C'))
venusOrb   = getPosArray('SSB', 'VENUS'  , 'ECLIPJ2000', utc2etArray('1977-10-01T00:00:00', '2012-12-31T23:59:59.99', 10000, 'C'))
earthOrb   = getPosArray('SSB', 'EARTH'  , 'ECLIPJ2000', utc2etArray('1977-10-01T00:00:00', '2012-12-31T23:59:59.99', 10000, 'C'))
marsOrb    = getPosArray('SSB', 'MARS'   , 'ECLIPJ2000', utc2etArray('1977-10-01T00:00:00', '2012-12-31T23:59:59.99', 10000, 'C'))
jupiterOrb = getPosArray('SSB', '5'      , 'ECLIPJ2000', utc2etArray('1977-10-01T00:00:00', '2012-12-31T23:59:59.99', 10000, 'C'))
saturnOrb  = getPosArray('SSB', '6'      , 'ECLIPJ2000', utc2etArray('1977-10-01T00:00:00', '2012-12-31T23:59:59.99', 10000, 'C'))
uranusOrb  = getPosArray('SSB', '7'      , 'ECLIPJ2000', utc2etArray('1977-10-01T00:00:00', '2012-12-31T23:59:59.99', 10000, 'C'))
neptuneOrb = getPosArray('SSB', '8'      , 'ECLIPJ2000', utc2etArray('1977-10-01T00:00:00', '2012-12-31T23:59:59.99', 10000, 'C'))
plotoOrb   = getPosArray('SSB', '9'      , 'ECLIPJ2000', utc2etArray('1977-10-01T00:00:00', '2012-12-31T23:59:59.99', 10000, 'C'))
voyager1Orb= getPosArray('SSB', '-31'    , 'ECLIPJ2000', utc2etArray('1977-10-01T00:00:00', '2012-12-31T23:59:59.99', 10000, 'C'))
voyager2Orb= getPosArray('SSB', '-32'    , 'ECLIPJ2000', utc2etArray('1977-10-01T00:00:00', '2012-12-31T23:59:59.99', 10000, 'C'))

mpl.rcParams['legend.fontsize'] = 10
fig = plt.figure()
ax = fig.gca(projection='3d')
plt.axis('equal')
ax.scatter(0,0,0,label='center')
ax.plot(mercuryOrb[0].A1, mercuryOrb[1].A1, mercuryOrb[2].A1, label='mercury orbit')
ax.plot(venusOrb[0].A1  , venusOrb[1].A1  , venusOrb[2].A1  , label='  venus orbit')
ax.plot(earthOrb[0].A1  , earthOrb[1].A1  , earthOrb[2].A1  , label='  earth orbit')
ax.plot(marsOrb[0].A1   , marsOrb[1].A1   , marsOrb[2].A1   , label='   mars orbit')
ax.plot(jupiterOrb[0].A1, jupiterOrb[1].A1, jupiterOrb[2].A1, label='jupiter orbit')
ax.plot(saturnOrb[0].A1 , saturnOrb[1].A1 , saturnOrb[2].A1 , label=' saturn orbit')
ax.plot(uranusOrb[0].A1 , uranusOrb[1].A1 , uranusOrb[2].A1 , label=' uranus orbit')
ax.plot(neptuneOrb[0].A1, neptuneOrb[1].A1, neptuneOrb[2].A1, label='neptune orbit')
ax.plot(plotoOrb[0].A1  , plotoOrb[1].A1  , plotoOrb[2].A1  , label='  ploto orbit')
ax.plot(voyager1Orb[0].A1, voyager1Orb[1].A1, voyager1Orb[2].A1, label='voyager1 orbit')
ax.plot(voyager2Orb[0].A1, voyager2Orb[1].A1, voyager2Orb[2].A1, label='voyager2 orbit')
ax.legend()
ax.set_xlim(-2e10, 2e10)
ax.set_ylim(-0.5e10,0e10)
ax.set_zlim(-2e10, 2e10)
ax.set_axis_off()
plt.show(0)

输出图像:

naif

有了数据，妈妈再也不用担心小朋友们问我问题了，大家have fun啦～

2013.07.29

最后给出题图的完整版～

naif

CSPICE星历程序简介和其与python的交互使用

星表与星历

CSPICE星历程序

CSPICE+Python

应用

1、太阳8字

2、“预测”日食

3、可视化旅行者1、2号的航行轨迹

标签

分类目录

存档

近期文章

友链

致谢

版权声明