基于SOFA的空间坐标系变换

最新推荐文章于 2024-04-19 00:30:00 发布

eyepeak

最新推荐文章于 2024-04-19 00:30:00 发布

阅读量3k

点赞数 2

文章标签： cuda 算法矩阵

本文链接：https://blog.csdn.net/wangrongfeng/article/details/122236903

版权

基于SOFA的空间坐标系变换

2.1 基于SOFA的空间坐标系变换

空间坐标系变换过程中岁差和章动的计算比较繁琐，还涉及恒星时计算等问题，易于出错。值得庆幸的是，IAU以源代码形式提供了计算软件包SOFA，可从网上下载。SOFA提供了天文计算函数和辅助的矢量／矩阵计算函数。以此为基础可完成航天器轨道计算所需的时间计算和坐标系变换。基于SOFA库实现的IAU1976变换矩阵计算如下所示。

void		propagateTransMat(const timeOfSpace& t0, mat3x3& ECI2ECFmat, mat3x3& TEME2ECImat){
	double xp, yp, dut1, ddp80, dde80;
	find(t0._year, t0._month, t0._day, xp, yp, dut1, ddp80, dde80);//获取极移数据等
	iauCal2jd(t0._year, t0._month, t0._day, &djmjd0, &date);	//获得儒略日
	time = (60.0*(double)(60 * t0._hour + t0._minute) + t0._second + (double)t0._millisecond / 1000) / DAYSEC;
	utc = date + time;
	iauDat(t0._year, t0._month, t0._day, time, &dat);
	tai = utc + dat / DAYSEC;
	tt = tai + 32.184 / DAYSEC;
	tut = time + dut1 / DAYSEC;
	ut1 = date + tut;
	iauPmat76(djmjd0, tt, rp);	//计算岁差矩阵
	iauNut80(djmjd0, tt, &dp80, &de80);//计算章动矩阵
	dpsi = dp80 + ddp80;
	deps = de80 + dde80;
	epsa = iauObl80(djmjd0, tt); /* Mean obliquity. */
	iauNumat(epsa, dpsi, deps, rn);
	iauRxr(rn, rp, rnpb);
	ee = iauEqeq94(djmjd0, tt) + ddp80 * cos(epsa);
	double tobl = epsa + deps + dde80;
	double angle = dpsi * cos(tobl);
	iauIr(tmat);
	iauRz(angle, tmat);
	iauRxr(rnpb, tmat, rp);
	iauTr(rp, tmat);
	TEME2ECImat = tmat;
	gst = iauAnp(iauGmst82(djmjd0 + date, tut) + ee);
	iauCr(rnpb, rc2ti);
	iauRz(gst, rc2ti);
	iauIr(rpom);
	iauRx(-yp, rpom);
	iauRy(-xp, rpom);
	iauRxr(rpom, rc2ti, rc2it);
	ECI2ECFmat = rc2it ;}

上述代码是从完整代码中截取了部分，是在CPU上执行的代码，同时计算ECI到ECF的变换矩阵和TEME到ECI的变换矩阵，后者主要用于SGP4模型的计算。

在上述过程中，需要调用的SOFA函数（函数内部还会调用其他SOFA函数），均是采用标准C语言实现的，为使其可供CUDA使用，仅需要在每个函数定义加入__host__ __device__关键字即可。

其中的极移数据，来源于IERS所发布的Bulletin数据，其格式可在网上获得。如不考虑极移数据，会产生较大的位置误差。采用Bulletin B的数据，其数据包括：PM-x [arcsec]， PM-y [arcsec]，UT1-UTC [seconds]，Psi [milliarcsec]，dEps [milliarcsec]。在前述的代码中，开始就通过查找获得给定时间的上述数据，供计算使用。极移数据量比较大，按时间排序、二分查找。

完整的源代码（CUDA版本，结合到二体模型实现）已上传，链接为源代码