Skip to content

Instantly share code, notes, and snippets.

@matwey

matwey/ADefine.h

Last active April 23, 2023 07:11
Show Gist options

  • Select an option

    Learn more about clone URLs
  • Save matwey/27a4aadbbab99908d78afb493bbda6fc to your computer and use it in GitHub Desktop.

Select an option

Learn more about clone URLs
Save matwey/27a4aadbbab99908d78afb493bbda6fc to your computer and use it in GitHub Desktop.
Download ZIP
LMST test
Raw
ADefine.h
/**
Name : ADefine.h 常量及宏定义声明文件
Author : Xiaomeng Lu
Version : 0.1
Date : Oct 13, 2012
Last Date : Oct 13, 2012
Description : 天文数字图像处理中常用的常量及宏定义
**/
#ifndef _ADEFINE_H_
#define _ADEFINE_H_
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <ctype.h>
#include <stdio.h>
#include <assert.h>
#include <time.h>
#include <vector>
#include <string>
#include <algorithm>
#if !defined(WIN32) && !defined(_WINDOWS_) // Linux平台下使用该类库
#include <unistd.h>
#else
#include <io.h> // Windows平台使用该类库
#endif
namespace AstroUtil
{
///////////////////////////////////////////////////////////////////////////////
/********************************** 宏定义 **********************************/
#if !defined(BZERO) // 内存清零
#define BZERO(ptr, szbyte) memset(ptr, 0, szbyte)
#endif
#if !defined(BYTE)
typedef unsigned char BYTE;
#endif
#if !defined(USHORT)
typedef unsigned short USHORT;
#endif
#if !defined(UINT)
typedef unsigned int UINT;
#endif
#ifndef WORD
typedef unsigned short WORD;
#endif
// 圆周率相关常量
#if !defined(API) // 圆周率
#define API 3.141592653589793
#endif
#define PI45 0.785398163397448
#define PI90 1.570796326794897
#define PI180 API
#define PI270 4.712388980384690
#define PI360 6.283185307179586
#define GtoR 0.017453292519943 // 使用乘法, 角度转换为弧度的系数
#define RtoG 57.295779513082323 // 使用乘法, 弧度转换为角度的系数
#define StoR 4.848136811095360e-06 // 使用乘法, 角秒转换为弧度的系数
#define RtoS 2.062648062470964e+05 // 使用乘法, 弧度转换为角秒的系数
#define ASEC StoR // 角秒对应的弧度
// 时间相关常量
#define MILLISEC 3600000 // 1度对应的毫角秒
#define MILLISEC5 18000000 // 5度对应的毫角秒
#define MILLISEC90 324000000 // 90度对应的毫角秒
#define MILLISEC180 648000000 // 180度对应的毫角秒
#define MILLISEC360 1296000000 // 360度对应的毫角秒
// 软件中限定最大/最小数
#define AMAX 1E30
#define AEPS 1E-4
#define cot(x) (tan(PI90 - (x))) // 计算弧度对应的余切
///////////////////////////////////////////////////////////////////////////////
// 定义与测站位置有关信息
typedef struct param_mount_site{
double lgt; // 地理经度, 量纲: 弧度
double lat; // 地理纬度, 量纲: 弧度
double alt; // 海拔高度, 量纲: 米
double airp; // 大气压, 量纲: 百帕
double temp; // 环境温度, 量纲: 摄氏度
int tmc; // 时区修正, 所使用本地时的时间与地理位置所在时区的差
}* ptr_param_mntst;
// 定义实数型'点'
typedef struct point_2d {// 两值实数点
double x;
double y;
public:
point_2d() {
x = y = 0;
}
point_2d(double x, double y) {
this->x = x; this->y = y;
}
point_2d& operator=(const point_2d &pt) {
if (this != &pt) memcpy(this, &pt, sizeof(point_2d));
return *this;
}
}* ptr_point_2d;
typedef struct point_3d {// 三值实数点
double x;
double y;
double z;
public:
point_3d() {
x = y = z = 0;
}
point_3d(point_3d &pt) {
this->x = pt.x;
this->y = pt.y;
this->z = pt.z;
}
point_3d(double x, double y, double z) {
this->x = x; this->y = y; this->z = z;
}
point_3d& operator=(const point_2d &pt) {
this->x = pt.x; this->y = pt.y; this->z = 0;
return *this;
}
point_3d& operator=(const point_3d &pt) {
if (this != &pt) memcpy(this, &pt, sizeof(point_3d));
return *this;
}
}* ptr_point_3d;
///////////////////////////////////////////////////////////////////////////////
};
#endif
Raw
AMath.cpp
/**
Name : AMath.cpp 常用数学函数定义文件
Author : Xiaomeng Lu
Version : 0.1
Date : Oct 13, 2012
Last Date : Oct 13, 2012
Description : 天文数字图像处理中常用的数学函数
**/
#include "ADefine.h"
#include "AMath.h"
namespace AstroUtil
{
/*---------------------------------------------------------------------------*/
///////////////////////////////////////////////////////////////////////////////
/*------------------------------- 部分球坐标转换 -------------------------------*/
void Sphere2Cart(double r, double alpha, double beta, double& x, double& y, double& z)
{
x = r * cos(beta) * cos(alpha);
y = r * cos(beta) * sin(alpha);
z = r * sin(beta);
}
void Cart2Sphere(double x, double y, double z, double& r, double& alpha, double& beta)
{
r = sqrt(x * x + y * y + z * z);
if (fabs(y) < AEPS && fabs(x) < AEPS) alpha = 0;
else if ((alpha = atan2(y, x)) < 0) alpha += PI360;
beta = atan2(z, sqrt(x * x + y * y));
}
/*---------------------------------------------------------------------------*/
}
Raw
AMath.h
/**
Name : AMath.h 常用数学函数声明文件
Author : Xiaomeng Lu
Version : 0.1
Date : Oct 13, 2012
Last Date : Oct 13, 2012
Description : 天文数字图像处理中常用的数学函数
**/
#ifndef _AMATH_H_
#define _AMATH_H_
#include <math.h>
namespace AstroUtil
{
// 计算实数的小数部分
#define frac(x) ((x) - floor(x))
// 将具有周期性的实数调整到一个周期内, 该周期的最小值为０
#define reduce(x, period) ((x) - floor((x) / (period)) * (period))
/*!
* \brief 将球坐标转换为笛卡尔坐标. 默认为右手坐标系, 且XY对应球坐标基平面, Z为极轴.
* alpha从X轴正向逆时针增加, X轴对应0点; beta为与XY平面的夹角, Z轴正向为正.
* 此转换中球坐标与赤道坐标系相同.
* 若球坐标为左手坐标系, 则alpha应取坐标系数值的负值
* \param[in] r 位置矢量的模. 对于天球坐标系取为1
* \param[in] alpha alpha角, 量纲: 弧度
* \param[in] beta beta角, 量纲: 弧度
* \param[out] x X轴坐标
* \param[out] y Y轴坐标
* \param[out] z Z轴坐标
**/
extern void Sphere2Cart(double r, double alpha, double beta, double& x, double& y, double& z);
/*!
* \fn void Cart2Sphere(double r, double alpha, double beta, double& x, double& y, double& z)
* \brief 将笛卡尔坐标转换为球坐标. 默认为右手坐标系, 且XY对应球坐标基平面, Z为极轴.
* alpha从X轴正向逆时针增加, X轴对应0点; beta为与XY平面的夹角, Z轴正向为正.
* 此转换中球坐标与赤道坐标系相同.
* 若球坐标为左手坐标系, 则alpha为坐标系数值的负值
* \param[in] x X轴坐标
* \param[in] y Y轴坐标
* \param[in] z Z轴坐标
* \param[out] r 位置矢量的模. 对于天球坐标系取为1
* \param[out] alpha alpha角, 量纲: 弧度
* \param[out] beta beta角, 量纲: 弧度
**/
extern void Cart2Sphere(double x, double y, double z, double& r, double& alpha, double& beta);
}
#endif
Raw
ATimeSpace.cpp
/*
* ATimeSpace.cpp
*
*/
#include "ADefine.h"
#include "ATimeSpace.h"
#include "AMath.h"
namespace AstroUtil
{
///////////////////////////////////////////////////////////////////////////////
ATimeSpace::ATimeSpace()
{
}
ATimeSpace::~ATimeSpace()
{
}
void ATimeSpace::Hour2HMS(const double hour, int &hh, int &mm, double &ss)
{
double t = hour;
hh = (int) t;
t = (t - hh) * 60;
mm = (int) t;
ss = (t - mm) * 60;
}
double ATimeSpace::HMS2Hour(const int hh, const int mm, const double ss)
{
return (hh + mm / 60.0 + ss/ 3600.0);
}
void ATimeSpace::Deg2DMS(const double deg, int &dd, int &mm, double &ss, int &sign)
{
double t;
sign = deg < 0 ? -1 : 1;
dd = (int) (t = deg * sign);
t = (t - dd) * 60;
mm = (int) t;
ss = (t - mm) * 60;
}
double ATimeSpace::DMS2Deg(const int sign, const int dd, const int mm, const double ss)
{
return (dd + mm / 60.0 + ss / 3600.0) * sign;
}
bool ATimeSpace::RAStr2Dbl(const char *pszVal, double &val)
{
if (pszVal == NULL || strlen(pszVal) == 0) return false;
int len = strlen(pszVal);
char *buff = new char[len + 1];
char ch;
int i, ii;
double temp;
bool bDot = false;
int part = 0;
for (i = 0, ii = 0, val = 0; i < len; ++i) {// scan input string
if ((ch = pszVal[i]) >= '0' && ch <= '9') {// is digit
if (ii == 2 && !bDot) {// scan two sequential digits
if (part < 2) {
buff[ii] = 0;
ii = 0;
temp = atof(buff);
if (part == 0 && temp >= 0 && temp < 24) // hour complete
val = temp;
else if (part == 1 && temp >= 0) // minute complete
val += temp / 60.0;
else return false;
part++;
}
else {
bDot = true;
buff[ii++] = '.';
}
}
buff[ii++] = ch;
}
else if (ch == ':' || ch ==' ') {// scan separator
if (bDot || part == 2) // more separator in second section
return false;
if (ii != 0) {
buff[ii] = 0;
ii = 0;
temp = atoi(buff);
if (part == 0 && temp >= 0 && temp < 24) val = temp;
else if (part == 1 && temp >= 0) val += temp / 60.0;
else return false;
}
part++;
}
else if (ch == '.') {
if(bDot) break;
bDot = true;
buff[ii++] = '.';
}
else break;
}
buff[ii] = 0;
if (ii > 1 || (ii == 1 && buff[0] != '.')) {
temp = atof(buff);
if (part == 0 && temp >= 0 && temp < 24) val = temp;
else if (part == 1 && temp >= 0) val += temp / 60;
else if (part == 2 && temp >= 0) val += temp / 3600;
}
delete []buff;
return (part <= 2);
}
bool ATimeSpace::DECStr2Dbl(const char *pszVal, double &val)
{
if(pszVal == NULL || strlen(pszVal) == 0) return false;
int len = strlen(pszVal);
char *buff;
int i, ii;
char ch;
double temp;
bool bDot = false;
bool bFlag = false;
int part = 0;
i = 0;
switch (ch = pszVal[0])
{
case '+':
i = 1;
break;
case '-':
bFlag = true;
i = 1;
break;
default:
if(!isdigit(ch) && ch != '.' && ch != ':' && ch != ' ') return false;
break;
}
buff = new char[len + 1];
for (ii = 0, temp = 0; i < len; ++i) {
if ((ch = pszVal[i]) >= '0' && ch <= '9') {
if (ii == 2 && !bDot) {
if (part < 2) {
buff[ii] = 0;
ii = 0;
temp = atof(buff);
if(part == 0) val = temp;
else if (part == 1 && val >= 0) val += temp / 60.0;
else return false;
part++;
}
else {
bDot = true;
buff[ii++] = '.';
}
}
buff[ii++] = ch;
}
else if (ch == ':' || ch ==' ') {
if (bDot || part == 2) // more separator in second section
return false;
if (ii != 0) {
buff[ii] = 0;
ii = 0;
temp = atoi(buff);
if (part == 0) val = temp;
else if (part == 1 && val >= 0) val += temp / 60.0;
else return false;
}
part++;
}
else if (ch == '.') {
if (bDot) return false;
bDot = true;
buff[ii++] = '.';
}
else return false;
}
buff[ii] = 0;
if (ii > 1 || (ii == 1 && buff[0] != '.')) {
temp = atof(buff);
if(part == 0) val = temp;
else if (part == 1 && val >= 0) val += temp / 60;
else if (part == 2 && val >= 0) val += temp / 3600;
}
if(bFlag) val = -val;
delete []buff;
return (part <= 2);
}
void ATimeSpace::RADbl2Str(const double val, char *pszVal)
{
int hh, mm;
double ss;
Hour2HMS(val, hh, mm, ss);
sprintf(pszVal, "%02d:%02d:%06.3f", hh, mm, ss);
}
void ATimeSpace::DECDbl2Str(const double val, char *pszVal)
{
int dd, mm, sign;
double ss;
Deg2DMS(val, dd, mm, ss, sign);
if (sign == 1) sprintf(pszVal, "+%02d:%02d:%05.2f", dd, mm, ss);
else sprintf(pszVal, "-%02d:%02d:%05.2f", dd, mm, ss);
}
double ATimeSpace::GetEpoch(const int year, const int month, const int day, const double hour)
{
return 2000 + (JulianDay(year, month, day, hour) - 2451545) / 365.25;
}
double ATimeSpace::GetEpoch(double mjd)
{
return 2000 + (mjd - 51544.5) / 365.25;
}
double ATimeSpace::Epoch2JC(double epoch)
{
return (epoch - 2000) * 0.01;
}
double ATimeSpace::JulianDay(int year, int month, int day, double hour)
{
int a = (14 - month) / 12;
int y = year + 4800 - a;
int m = month + 12 * a - 3;
int JDN = day + (153 * m + 2) / 5 + 365 * y + y / 4 - y / 100 + y / 400 - 32045;
return (JDN + hour / 24 - 0.5);
}
double ATimeSpace::ModifiedJulianDay(int nYear, int nMonth, int nDay, double vHour)
{
return (JulianDay(nYear, nMonth, nDay, vHour) - 2400000.5);
}
void ATimeSpace::JD2YMDH(double jd, int &year, int &month, int &day, double &hour)
{
int jdn = (int) (jd + 0.5);
int f = jdn + 1401 + (((4 * jdn + 274277) / 146097) * 3) / 4 - 38;
int e = 4 * f + 3;
int g = (e % 1461) / 4;
int h = 5 * g + 2;
day = (h % 153) / 5 + 1;
month = (h / 153 + 2) % 12 + 1;
year = e / 1461 - 4716 + (14 - month) / 12;
hour = (jd - jdn) * 24;
}
///////////////////////////////////////////////////////////////////////////////
void ATimeSpace::PMatEqu(double T1, double T2)
{
double dt = T2 - T1;
double Z, Zeta, Theta;
const double sec = 3600;
double c1, c2, c3; // Z, Zeta, Theta对应的余弦项
double s1, s2, s3; // Z, Zeta, Theta对应的正弦项
Zeta = ((2306.2181 + (1.39656 - 0.000139 * T1) * T1)
+ ((0.30188 - 0.000345 * T1) + 0.017998 * dt) * dt) * dt / sec;
Z = Zeta + ((0.7928 + 0.000411 * T1) + 0.000205 * dt) * dt * dt / sec;
Theta = ((2004.3109 - (0.8533 + 0.000217 * T1) * T1)
- ((0.42665 + 0.000217 * T1) + 0.041833 * dt) * dt) * dt / sec;
Z *= GtoR;
Zeta *= GtoR;
Theta*= GtoR;
c1 = cos(Z);
c2 = cos(Theta);
c3 = cos(Zeta);
s1 = sin(Z);
s2 = sin(Theta);
s3 = sin(Zeta);
m_matPN[0][0] = -s1 * s3 + c1 * c2 * c3;
m_matPN[0][1] = -s1 * c3 - c1 * c2 * s3;
m_matPN[0][2] = -c1 * s2;
m_matPN[1][0] = c1 * s3 + s1 * c2 * c3;
m_matPN[1][1] = c1 * c3 - s1 * c2 * s3;
m_matPN[1][2] = -s1 * s2;
m_matPN[2][0] = s2 * c3;
m_matPN[2][1] = -s2 * s3;
m_matPN[2][2] = c2;
}
void ATimeSpace::PN_Matrix(double T1, double T2)
{
PMatEqu(T1, T2);
EqMean2True(T2, m_matPN[0][0], m_matPN[1][0], m_matPN[2][0]);
EqMean2True(T2, m_matPN[0][1], m_matPN[1][1], m_matPN[2][1]);
EqMean2True(T2, m_matPN[0][2], m_matPN[1][2], m_matPN[2][2]);
}
void ATimeSpace::PError(double &X, double &Y, double &Z)
{
double U, V, W;
U = m_matPN[0][0] * X + m_matPN[0][1] * Y + m_matPN[0][2] * Z;
V = m_matPN[1][0] * X + m_matPN[1][1] * Y + m_matPN[1][2] * Z;
W = m_matPN[2][0] * X + m_matPN[2][1] * Y + m_matPN[2][2] * Z;
X = U;
Y = V;
Z = W;
}
void ATimeSpace::Aberrat(double JC, double &VX, double &VY, double &VZ)
{
double l;
double CL;
l = PI360 * frac(0.27908 + 100.00214 * JC);
CL = cos(l);
VX = -0.0000994 * sin(l);
VY = 0.0000912 * CL;
VZ = 0.0000395 * CL;
}
void ATimeSpace::EqTrue2Mean(double JC, double &X, double &Y, double &Z)
{
double dx, dy, dz;
DeltaMean2True(JC, X, Y, Z, dx, dy, dz);
X -= dx;
Y -= dy;
Z -= dz;
}
void ATimeSpace::EqMean2True(double JC, double &X, double &Y, double &Z)
{
double dx, dy, dz;
DeltaMean2True(JC, X, Y, Z, dx, dy, dz);
X += dx;
Y += dy;
Z += dz;
}
void ATimeSpace::DeltaMean2True(double JC, double X, double Y, double Z, double &dx, double &dy, double &dz)
{
double LS;
double D;
double F;
double N;
double EPS;
double DPSI;
double DEPS;
double c;
double s;
LS = PI360 * frac(0.993133 + 99.997306 * JC); // 平近点角
D = PI360 * frac(0.827362 + 1236.853087 * JC); //
F = PI360 * frac(0.259089 + 1342.227826 * JC);
N = PI360 * frac(0.347346 - 5.372447 * JC);
EPS = 0.4090928 - 0.00022696 * JC;
DPSI = (-17.2 * sin(N) - 1.319 * sin(2 * (F - D + N))
- 0.227 * sin(2 * (F + N)) + 0.206 * sin(2 * N) + 0.143 * sin(LS)) / 3600 * GtoR;
DEPS = (9.203 * cos(N) + 0.574 * cos(2 * (F - D + N))
+ 0.098 * cos(2 * (F + N)) - 0.09 * cos(2 * N)) / 3600 * GtoR;
c = DPSI * cos(EPS);
s = DPSI * sin(EPS);
dx = -(c * Y + s * Z);
dy = (c * X - DEPS * Z);
dz = (s * X + DEPS * Y);
}
void ATimeSpace::Actual2Epoch(double E0, double ra_a, double de_a, double E1, double &ra_e, double &de_e)
{
double r;
double X, Y, Z;
double VX, VY, VZ;
double T0 = Epoch2JC(E0);
double T1 = Epoch2JC(E1);
Aberrat(T0, VX, VY, VZ);
Sphere2Cart(1.0, ra_a, de_a, X, Y, Z);
X -= VX;
Y -= VY;
Z -= VZ;
EqTrue2Mean(T0, X, Y, Z);
PMatEqu(T0, T1);
PError(X, Y, Z);
Cart2Sphere(X, Y, Z, r, ra_e, de_e);
}
void ATimeSpace::Epoch2Actual(double E0, double ra_e, double de_e, double E1, double &ra_a, double &de_a)
{
double r;
double X, Y, Z;
double VX, VY, VZ;
double T0 = Epoch2JC(E0);
double T1 = Epoch2JC(E1);
PN_Matrix(T0, T1);
Aberrat(T1, VX, VY, VZ);
Sphere2Cart(1.0, ra_e, de_e, X, Y, Z);
PError(X, Y, Z);
X += VX;
Y += VY;
Z += VZ;
Cart2Sphere(X, Y, Z, r, ra_a, de_a);
}
double ATimeSpace::GetGMST0(const double mjd)
{
double t = (mjd - 51544.5) / 36525;
double gmst0;
gmst0 = reduce(100.46061837 + (36000.770053608 + (0.000387933 - t / 38710000) * t) * t, 360.0);
return (gmst0 * GtoR);
}
double ATimeSpace::LocalMeanSiderialTime(const double mjd, const double longitude)
{
double t = (mjd - 51544.5) / 36525;
double gmst; // 格林威治平恒星时
double lmst; // 本地平恒星时
gmst = 280.46061837 + (13185000.77005374225 + (0.000387933 - t / 38710000) * t) * t;
lmst = gmst * GtoR + longitude;
return reduce(lmst, PI360);
}
void ATimeSpace::Eq2AltAzi(double ha, double dec, double latitude, double &azi, double &alt)
{
double sinz, cosz;
double sina, cosa;
sinz = sin(latitude) * sin(dec) + cos(latitude) * cos(dec) * cos(ha);
cosz = sqrt(1 - sinz * sinz);
alt = atan2(sinz, cosz);
sina = cos(dec) * sin(ha);
cosa = sin(latitude) * cos(dec) * cos(ha) - cos(latitude) * sin(dec);
azi = reduce(atan2(sina, cosa) + API, PI360);
}
void ATimeSpace::AltAzi2Eq(double azi, double alt, double latitude, double &ha, double &dec)
{
double z, sinz, cosz;
double sint, cost;
z = API * 0.5 - alt;
sinz = sin(latitude) * cos(z) + cos(latitude) * sin(z) * cos(azi);
cosz = sqrt(1 - sinz * sinz);
dec = atan2(sinz, cosz);
sint = -1 * sin(z) * sin(azi);
cost = cos(latitude) * cos(z) - sin(latitude) * sin(z) * cos(azi);
ha = reduce(atan2(sint, cost), PI360);
}
double ATimeSpace::ParAngle(double ha, double dec, double latitude, char mode)
{
double y = sin(ha);
double x = tan(latitude) * cos(dec) - sin(dec) * cos(ha);
double q = atan2(y, x);
double azi(0), alt(0);
if(mode != 'N')
{
Eq2AltAzi(ha, dec, latitude, azi, alt);
if(mode == '+')
{
q += alt;
}
else if(mode == '-')
{
q -= alt;
//q=q/2.0; // optical derotator case
}
}
return reduce(q, PI360);
}
double ATimeSpace::RefractReal(const double vAlt, const double pressure, const double temperature)
{
double arcdeg = vAlt + 10.3 * GtoR / (vAlt * RtoG + 5.11);
double R = 1.02 * cot(arcdeg) * GtoR / 60;
double scale = pressure * 283 / 1010 / (273 + temperature);
return scale * R;
}
double ATimeSpace::RefractVisual(const double vAlt, const double pressure, const double temperature)
{
double arcdeg = vAlt + 7.31 * GtoR / (vAlt * RtoG + 4.4);
double R = cot(arcdeg) * GtoR / 60;
double scale = pressure * 283 / 1010 / (273 + temperature);
return scale * R;
}
double ATimeSpace::AirMass(const double vAlt)
{
double x = vAlt + 244.0 * GtoR / (165 + 47 * pow(vAlt * RtoG, 1.1));
return 1.0 / sin(x);
}
///////////////////////////////////////////////////////////////////////////////
}
Raw
ATimeSpace.h
/*!
* class ATimeSpace 常用天文时空相关转换接口
* Version：0.2
*
* 功能列表:
* 1. 赤经/赤纬格式转换: 字符串<->实数
* 2. 时间/角度格式转换: 实数<->时/度, 分, 秒
*
* 创建时间: 2012年3月1日
* 完成时间: 2012年3月10日
* 作者: 卢晓猛, lxm@nao.cas.cn
*/
#ifndef _ATIMESPACE_H_
#define _ATIMESPACE_H_
//#include "QtGui"
namespace AstroUtil {
///////////////////////////////////////////////////////////////////////////////
class ATimeSpace
{
public:
ATimeSpace();
virtual ~ATimeSpace();
protected:
// 赤道坐标系岁差与章动系数, 3x3矩阵
double m_matPN[3][3];
// Function
public:
/*!
* \fn bool RAStr2Dbl(const char *pszVal, double &val)
* \brief Transfer R.A. from string to double
* \param[in] pszVal R.A. in string style
* \param[out] val R.A. in double, in hour
* \return
* if pszVal format is right then transfer it to be double and return true, or else return false
* \note
* pszVal could be the following style:
* HH:MM:SS.SS, colon could be omitted or replaced by space, and both MM and SS could be omitted too.
*/
bool RAStr2Dbl(const char *pszVal, double &val);
/*!
* \fn bool DECStr2Dbl(const char *pszVal, double &val)
* \brief Transfer DEC. from string to double
* \param[in] pszVal DEC. in string style
* \param[out] val DEC. in double, in degree
* \return
* if pszVal format is right then transfer it to be double and return true, or else return false
* \note
* pszVal could be the following style:
* sDD:MM:SS.SS, colon could be omitted or replaced by space, and both MM and SS could be omitted too.
* 's' means plus or minus sign
*/
bool DECStr2Dbl(const char *pszVal, double &val);
/*!
* \fn bool RADbl2Str(const double val, char *pszVal)
* \brief Transfer R.A. from double to string
* \param[in] val R.A. in double, in hour
* \param[out] pszVal R.A. in string style
* \return
* regulate val to be between 0 and 24 and transfer it to be string style
* \note
* pszVal could be the following style:
* HH:MM:SS.SS, colon could be omitted or replaced by space, and both MM and SS could be omitted too.
*/
void RADbl2Str(const double val, char *pszVal);
/*!
* \fn bool DECDbl2Str(const double val, char *pszVal)
* \brief Transfer DEC. from double to string
* \param[in] val DEC. in double, in hour
* \param[out] pszVal DEC in string style
* \return
* regulate val to be between -90 and +90 and transfer it to be string style
* \note
* pszVal could be the following style:
* sDD:MM:SS.SS, colon could be omitted or replaced by space, and both MM and SS could be omitted too.
* 's' means plus or minus sign
*/
void DECDbl2Str(const double val, char *pszVal);
/*!
* \fn void Hour2HMS(const double hour, int &hh, int &mm, double &ss)
* \brief resolve double hour to hour, minute and second
* \param[in] hour double time, in hour
* \param[out] hh hour
* \param[out] mm minute
* \param[out] ss second
*/
static void Hour2HMS(const double hour, int &hh, int &mm, double &ss);
/*!
* \fn double HMS2Hour(const int hh, const int mm, const double ss)
* \brief construct hour, minute and second to double hour
* \param[in] hh hour
* \param[in] mm minute
* \param[in] ss second
* \return
* double hour
*/
double HMS2Hour(const int hh, const int mm, const double ss);
/*!
* \fn void Deg2DMS(const double deg, int &dd, int &mm, double &ss, int &sign)
* \brief resolve double degree to degree, minute and second
* \param[in] deg double degree, in degree
* \param[out] dd degree
* \param[out] mm minute
* \param[out] ss second
* \param[out] sign plus or minus sign. if deg is less than 0 then sign is -1, or else 1
*/
void Deg2DMS(const double deg, int &dd, int &mm, double &ss, int &sign);
/*!
* \fn double DMS2Deg(const int sign, const int dd, const int mm, const double ss)
* \brief construct degree, minute and second to double degree
* \param[in] sign plus or minus sign. if deg is less than 0 then sign is -1, or else 1
* \param[in] dd degree
* \param[in] mm minute
* \param[in] ss second
* \return
* double degree
*/
double DMS2Deg(const int sign, const int dd, const int mm, const double ss);
public:
/*!
* \brief calculate epoch referred by year\month\day\hour
*/
static double GetEpoch(const int year, const int month, const int day, const double hour = 0.0);
/*!
* @brief calculate epoch from modified julian day
*/
static double GetEpoch(double mjd);
/*! \fn double Epoch2JC(const double epoch)
* \brief 根据历元计算对应的儒略纪元.
* \param[in] epoch 由函数GetEpoch()计算或者用户输入的历元
* \return epoch对应的儒略纪元. 从J2000开始计算
*/
double Epoch2JC(double epoch);
/*!
* @brief 由格里高利历计算儒略日
* \param[in] year 年
* \param[in] month 月
* \param[in] day 日
* \param[in] hour 当日小时数
* @return
* 公元历对应的儒略日
*/
static double JulianDay(int year, int month, int day, double hour);
/*! \fn double ModifiedJulianDay(int nYear, int nMonth, int nDay, double vHour)
* \brief 计算修正儒略日
* \param[in] nYear 年
* \param[in] nMonth 月
* \param[in] nDay 日
* \param[in] vHour 时
* \return 年月日时对应的修正儒略日
**/
static double ModifiedJulianDay(int nYear, int nMonth, int nDay, double vHour);
/*!
* @brief 由儒略日计算对应的公元历
* @param jd 儒略日
* @param year 年
* @param month 月
* @param day 日
* @param hour 当日小时数
*/
void JD2YMDH(double jd, int &year, int &month, int &day, double &hour);
protected:
/*! \fn void PMatEqu(double T1, double T2)
* \brief 根据儒略世纪计算岁差系数.
* \param[in] T1 儒略世纪
* \param[in] T2 儒略世纪
*
* \note
* 如果已经知道T1历元对应的赤道坐标位置, 需要计算T2历元时的坐标位置. 使用本函数
* 计算历元变化对目标位置的影响.<br>
* 本函数计算的结果, 并不是已经分解到赤经/赤纬轴的坐标修正量, 而是岁差系数.<br>
* 本函数计算的结果, 存储在类的矩阵m_matPN中。
**/
void PMatEqu(double T1, double T2);
/*! \fn void PN_Matrix(double T1, double T2)
* \brief 计算从平历元T1到真历元T2的岁差系数.
* \param[in] T1 儒略世纪
* \param[in] T2 儒略世纪
**/
void PN_Matrix(double T1, double T2);
/*! \fn void PError(double& X, double& Y, double &Z)
* \brief 计算岁差带来的偏差
* \param X X轴坐标
* \param Y Y轴坐标
* \param Z Z轴坐标
**/
void PError(double &X, double &Y, double &Z);
/*! \fn void Aberrat(double JC, double& VX, double& VY, double& VZ)
* \brief 赤道坐标中, 地球运动速度矢量
* \param[in] JC 儒略世纪
* \param[out] VX X轴速度, 单位：光速
* \param[out] VY Y轴速度, 单位: 光速
* \param[out] VZ Z轴速度, 单位: 光速
**/
void Aberrat(double JC, double &VX, double &VY, double &VZ);
/*! \fn void EqTrue2Mean(double JC, double& X, double& Y, double& Z)
* \brief 从真坐标转换为平坐标
* \param[in] JC 儒略世纪
* \param X X轴坐标
* \param Y Y轴坐标
* \param Z Z轴坐标
**/
void EqTrue2Mean(double JC, double &X, double &Y, double &Z);
/*! \fn void EqMean2True(double JC, double& X, double& Y, double& Z)
* \brief 从平坐标转换为真坐标
* \param[in] JC 儒略世纪
* \param X X轴坐标
* \param Y Y轴坐标
* \param Z Z轴坐标
**/
void EqMean2True(double JC, double &X, double &Y, double &Z);
/*! \fn void DeltaMean2True(double JC, double X, double Y, double Z, double& dx, double& dy, double& dz)
* \brief 真坐标与平坐标的差值
* \param[in] JC 儒略世纪
* \param[in] X 笛卡尔(直角)坐标系X轴坐标
* \param[in] Y 笛卡尔(直角)坐标系Y轴坐标
* \param[in] Z 笛卡尔(直角)坐标系Z轴坐标
* \param[out] dx X轴差值
* \param[out] dy Y轴差值
* \param[out] dz Z轴差值
**/
void DeltaMean2True(double JC, double X, double Y, double Z, double &dx, double &dy, double &dz);
public:
/*! \fn void Actual2Epoch(double E0, double ra_a, double de_a, double E1, double &ra_e, double &de_e)
* \brief 从当前历元坐标转换为指定历元坐标
* \param[in] E0 当前西元历元
* \param[in] ra_a 当前历元下的赤经, 量纲: 弧度
* \param[in] de_a 当前历元下的赤纬, 量纲: 弧度
* \param[in] E1 指定西元历元
* \param[out] ra_e 指定历元下的赤经, 量纲: 弧度
* \param[out] de_e 指定历元下的赤纬, 量纲: 弧度
**/
void Actual2Epoch(double E0, double ra_a, double de_a, double E1, double &ra_e, double &de_e);
/*! \fn void Epoch2Actual(double T0, double ra_e, double de_e, double T1, double &ra_a, double &de_a)
* \brief 从指定历元坐标转换为当前历元坐标
* \param[in] E0 指定西元历元
* \param[in] ra_a 指定历元下的赤经, 量纲: 弧度
* \param[in] de_a 指定历元下的赤纬, 量纲: 弧度
* \param[in] E1 西元历元
* \param[out] ra_e 当前指定历元下的赤经, 量纲: 弧度
* \param[out] de_e 当前指定历元下的赤纬, 量纲: 弧度
**/
void Epoch2Actual(double E0, double ra_e, double de_e, double E1, double &ra_a, double &de_a);
public:
/*!
* \fn double GetGMST0(const double mjd)
* \brief 根据修正儒略日计算格林威治0时的平恒星时.
* \param[in] mjd 修正儒略日
* \return 以弧度为单位的平恒星时
*/
double GetGMST0(const double mjd);
/*!
* \fn double LocalMeanSiderialTime(const double mjd)
* \brief 根据修正儒略日计算本地平恒星时.
* \param[in] mjd 由ModifiedJulianDay计算得到的修正儒略日
* \param[in] longitude 地理经度, 量纲: 弧度
* \return 以弧度为单位的本地平恒星时
*/
static double LocalMeanSiderialTime(const double mjd, const double longitude);
public:
/*!
* \fn void Eq2AltAzi(double ha, double dec, double latitude, double &azi, double &alt)
* \brief 赤道坐标转换为地平坐标
* \param[in] ha 时角, 量纲: 弧度
* \param[in] dec 赤纬, 量纲: 弧度
* \param[in] latitude 地理纬度, 量纲: 弧度
* \param[out] azi 方位角, 量纲: 弧度, 北零点
* \param[out] alt 俯仰角, 量纲: 弧度
**/
static void Eq2AltAzi(double ha, double dec, double latitude, double &azi, double &alt);
/*!
* \fn void AltAzi2Eq(double azi, double alt, double latitude, double &ha, double &dec)
* \brief 地平坐标转换为赤道坐标
* \param[in] azi 方位角, 量纲: 弧度, 北零点
* \param[in] alt 俯仰角, 量纲: 弧度
* \param[in] latitude 地理纬度, 量纲: 弧度
* \param[out] ha 时角, 量纲: 弧度
* \param[out] dec 赤纬, 量纲: 弧度
**/
void AltAzi2Eq(double azi, double alt, double latitude, double &ha, double &dec);
/**
* \fn double ParAngle(double ha, double dec, double latitude, char mode)
* \brief 计算地平式望远镜抵消地平自转, 消转器的转动弧度
* \param[in] ha 时角, 量纲: 弧度
* \param[in] dec 赤纬, 量纲: 弧度
* \param[in] latitude 地理纬度, 量纲: 弧度
* \param[in] mode 卡塞格林焦点, mode='N'
* 耐氏焦点, 东边焦点mode='+'
* 西边焦点mode='-'
* \return 抵消地球自转需要的弧度
**/
double ParAngle(double ha, double dec, double latitude, char mode);
/*! \fn double RefractReal(const double vAlt, const double pressure, const double temperature)
* \brief 计算大气折射
* \param[in] vAlt 真高度角, 量纲: 弧度
* \param[in] presure 气压, 量纲: 毫巴
* \param[in] temperature 温度, 量纲: 摄氏度
* \return
* 大气折射, 量纲: 弧度
* \note
* 毫巴: 100帕
**/
static double RefractReal(const double vAlt, const double pressure, const double temperature);
/*! \fn double RefractVisual(const double vAlt, const double pressure, const double temperature)
* \brief 计算大气折射
* \param[in] vAlt 视高度角, 量纲: 弧度
* \param[in] presure 气压, 量纲: 毫巴
* \param[in] temperature 温度, 量纲: 摄氏度
* \return
* 大气折射, 量纲: 弧度
* \note
* 毫巴: 100帕
**/
double RefractVisual(const double vAlt, const double pressure, const double temperature);
/*! \fn double AirMass(const double vAlt)
* \brief 计算大气质量
* \param[in] vAlt 视高度角, 量纲: 弧度
* \return 大气质量
* \note
* 大气质量计算公式采用Pickering(2002)模型, 其公式为:<br>
* \f$X = \frac{1}{sin\left(h + \frac{244}{165 + 47 \times \mbox{h}^{1.1}}\right)} \f$<br>
* 上式中, h为视高度角, 量纲为角度
**/
double AirMass(const double vAlt);
};
///////////////////////////////////////////////////////////////////////////////
}
#endif /* _ATIMESPACE_H_ */
Raw
Makefile
all: test
test: test.cpp
$(CXX) -m32 -o test test.cpp ATimeSpace.cpp AMath.cpp libtimeandcoordinate.a
Raw
test.cpp
#include <iostream>
#include "ATimeSpace.h"
#include "TimeAndCoordinate.h"
int main(int argc, char** argv) {
int year,mon,mday,hour,min;
double sec;
double lgt = 42.5;
double lat = 42.5;
TimeAndCoordinate t;
AstroUtil::ATimeSpace s;
while (std::cin >> year >> mon >> mday >> hour >> min >> sec) {
double fhour = hour + static_cast<double>(min)/60 + static_cast<double>(sec)/3600;
int nsec = sec;
int msec = (sec - nsec) * 1000;
std::cout << t.GetLMST(lgt, lat, year, mon, mday, hour, min, sec, msec) << " " << s.LocalMeanSiderialTime(s.ModifiedJulianDay(year, mon, mday, fhour), lgt) << std::endl;
}
return 0;
}
Raw
TimeAndCoordinate.h
#pragma once
class TimeAndCoordinate
{
public:
/**
*程序功能：将J2000历元下的赤经赤纬转换至当前历元下赤经赤纬
*参数：
*inRA - 输入赤经，量纲：角度
*inDEC - 输入赤纬，量纲：角度
*outEpoch - 输入当前历元。量纲：简约儒略日（MJD）
*outRA - 输出赤经，量纲：角度
*outDEC - 输出赤纬，量纲：角度
*****************************************************************************************
* inRA - input Right ascension(decimal degree)
* inDEC - input Declination(decimal degree)
* outEpoch - output Epoch。(real epoch)
* outRA - output Right ascension(decimal degree)
* outDEC - output Declination(decimal degree)
*Return:void
*Notes：Convert the Right ascension and Declination from J2000 to real epoch
**/
void Epoch2Actual(double inRA,double inDEC,double outEpoch,double &outRA, double &outDEC);
/**
*程序功能：将当前历元下的赤经赤纬转换至J2000历元下赤经赤纬
*inEpoch - 输入当前历元。量纲：简约儒略日
*inRA - 输入赤经，量纲：角度
*inDEC - 输入赤纬，量纲：角度
*outRA - 输出赤经，量纲：角度
*outDEC - 输出赤纬，量纲：角度
*****************************************************************************************
* inEpoch - input Epoch。(real epoch)
* inRA - input Right ascension(decimal degree)
* inDEC - input Declination(decimal degree)
* outRA - output Right ascension(decimal degree)
* outDEC - output Declination(decimal degree)
*Return: void
*Notes: Convert the Right ascension and Declination from real epoch to J2000
**/
void Actual2Epoch(const double inEpoch, double inRA, double inDEC,double& outRA, double& outDEC) ;
/**
*程序功能：将年月日时分秒形式的日期转换为以简约儒略日表示的历元
*Year - 年份，UTC时间
*Month - 月份，UTC时间
*Day – 日数，UTC时间
*Hour - 小时，UTC时间
*Minute - 分钟，UTC时间
*Second - 秒，UTC时间
*MicroSec - 微秒，UTC时间
*返回值：Epoch - 简约儒略日，UTC时间
*****************************************************************************************
* Year - Year, UTC(coordinated universal time)
* Month - Month, UTC(coordinated universal time)
* Day - Day, UTC(coordinated universal time)
* Hour - Hour, UTC(coordinated universal time)
* Minute - Minute, UTC(coordinated universal time)
* Second - Second, UTC(coordinated universal time)
* MicroSec - Microsecond, UTC(coordinated universal time)
*Return：modified Julian date **/
double GetEpoch(int Year,int Month,int Day,int Hour,int Minute,int Second,int MicroSec);
/**
*程序功能：获取输入地理经纬度当地的地方恒星时计算
*Lgt - 地理经度，量纲：角度，东经为正
*Lat - 地理纬度，量纲：角度，北纬为正
*Year - 年份，UTC时间
*Month - 月份，UTC时间
*Day - 日数，UTC时间
*Hour - 小时，UTC时间
*Minute - 分钟，UTC时间
*Second - 秒，UTC时间
*MicroSec - 微秒，UTC时间
*返回值: LMST - 地方恒星时,(角度)
*****************************************************************************************
* Lgt - longitude (decimal degree, east longitude is positive)
* Lat - latitude (decimal degree, northern latitude is positive)
* Year - Year, UTC(coordinated universal time)
* Month - Month, UTC(coordinated universal time)
* Day - Day, UTC(coordinated universal time)
* Hour - Hour, UTC(coordinated universal time)
* Minute - Minute, UTC(coordinated universal time)
* Second - Second, UTC(coordinated universal time)
* MicroSec - Microsecond, UTC(coordinated universal time)
*Return：local sidereal time (decimal degree)
**/
double GetLMST(double Lgt,double Lat,int Year,int Month,int Day,int Hour,int Minute,int Second,int MicroSec);
/**
*程序功能：将测站地平坐标系下方位角、高度角转换至赤道坐标系下赤经赤纬
*Azi - 输入方位角，量纲：角度，北零点
*Alt - 输入俯仰角，量纲：角度，水平零点
*LMST - 输入本地恒星时，量纲：角度
*Lgt - 输入地理经度，量纲：角度，东经为正
*Lat - 输入地理纬度，量纲：角度，北纬为正
*RA - 输出赤经，量纲：角度，坐标系：当前历元
*DEC - 输出赤纬，量纲：角度，坐标系：当前历元
*****************************************************************************************
* Azi - input Azimuth (decimal degree, North zero)
* Alt - input Altitude (decimal degree, horizon zero)
* LMST - input Local Mean sidereal time (decimal degree)
* Lgt - input longitude(decimal degree, east longitude is positive)
* Lat - input latitude(decimal degree, northern latitude is positive)
* RA - output Right ascension(decimal degree,(Coordinate system: real epoch)
* DEC - output Declination(decimal degree Coordinate system: real epoch)
*Return：void
*notes：Transfer function from horizontal coordinate system to equatorial coordinate system.
**/
void AltAzi2Eq(double Azi, double Alt, double LMST, double Lgt, double Lat, double &RA, double &DEC);
/**
*程序功能：将测站赤道坐标系下赤经赤纬转换至地平坐标系下方位角、高度角,同时根据焦点类型的选择确定像场旋转角
*RA - 输入赤经，量纲：角度，坐标系：当前历元
*DEC - 输入赤纬，量纲：角度，坐标系：当前历元
*LMST - 输入本地恒星时，量纲：角度
*Lgt - 输入地理经度，量纲：角度，东经为正
*Lat - 输入地理纬度，量纲：角度，北纬为正
*Focus - 输入焦点类型。1：卡焦；2：东耐焦；3：西耐焦
*Azi - 输出方位角，量纲：角度，北零点
*Alt - 输出俯仰角，量纲：角度，水平零点
*Rot - 输出像场旋转角，量纲：角度
*****************************************************************************************
* RA - input Right ascension(decimal degree,(Coordinate system: real epoch)
* DEC - input Declination(decimal degree Coordinate system: real epoch)
* LMST - input Local Mean sidereal time (decimal degree)
* Lgt - input longitude(decimal degree, east longitude is positive)
* Lat - input latitude(decimal degree, northern latitude is positive)
* Focus - input Focus type. 1: Cassegrain focus; 2: East Nasmyth focus 3: West Nasmyth focus
* Azi - output Azimuth (decimal degree, North zero)
* Alt - output Altitude (decimal degree, horizon zero)
* Rot - output Rotation angle of image field(dimension: decimal degree)
*Return：void
*Notes：Transfer function from equatorial coordinate system to horizontal coordinate system.
**/
void Eq2AltAzi(double RA,double DEC,double LMST,double Lgt,double Lat,int Focus,double &Azi,double &Alt,double &Rot);
};
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment