BTA_utils/Stellarium_control/angle_functions.c
2015-11-03 17:42:44 +03:00

144 lines
4.5 KiB
C

/*
* angle_functions.c - different functions for angles/times processing in different format
*
* Copyright 2015 Edward V. Emelianov <eddy@sao.ru, edward.emelianoff@gmail.com>
*
* This program is free software; you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation; either version 2 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
* MA 02110-1301, USA.
*/
#include <assert.h>
#include <math.h>
#include <stdlib.h>
#include <limits.h>
#include <time.h>
#include <string.h>
#include <slamac.h> // SLA macros
#include "bta_shdata.h"
#include "angle_functions.h"
#include "usefull_macros.h"
extern void sla_caldj(int*, int*, int*, double*, int*);
extern void sla_amp(double*, double*, double*, double*, double*, double*);
extern void sla_map(double*, double*, double*, double*, double*,double*, double*, double*, double*, double*);
void slacaldj(int y, int m, int d, double *djm, int *j){
int iy = y, im = m, id = d;
sla_caldj(&iy, &im, &id, djm, j);
}
void slaamp(double ra, double da, double date, double eq, double *rm, double *dm ){
double r = ra, d = da, mjd = date, equi = eq;
sla_amp(&r, &d, &mjd, &equi, rm, dm);
}
void slamap(double rm, double dm, double pr, double pd,
double px, double rv, double eq, double date,
double *ra, double *da){
double r = rm, d = dm, p1 = pr, p2 = pd, ppx = px, prv = rv, equi = eq, dd = date;
sla_map(&r, &d, &p1, &p2, &ppx, &prv, &equi, &dd, ra, da);
}
/**
* convert angle in seconds into degrees
* @return angle in range [0..360]
*/
double sec_to_degr(double sec){
double sign = 1.;
if(sec < 0.){
sign = -1.;
sec = -sec;
}
int d = ((int)sec / 3600.);
sec -= ((double)d) * 3600.;
d %= 360;
double angle = sign * (((double)d) + sec / 3600.);
if(angle < 0.) angle += 360.;
return (angle);
}
const double jd0 = 2400000.5; // JD for MJD==0
/**
* Calculate apparent place for given coordinates
* @param r,d (i) - RA/Decl for JD2000.0 (RA in hours, Decl in degrees)
* @param appRA, appDecl (o) - calculated apparent place (in seconds)
*/
void calc_AP(double r, double d, double *appRA, double *appDecl){
double mjd = 51544.5, pmra = 0., pmdecl = 0.; // mjd2000
// convert to radians
r *= DH2R;
d *= DD2R;
/* double ra2000, decl2000; // coordinates for 2000.0
DBG("slaamp(%g, %g, %g, 2000.0, ra, dec)", r,d,mjd);
slaamp(r, d, mjd, 2000.0, &ra2000, &decl2000);
DBG("2000: %g, %g", ra2000*DR2H, decl2000*DR2D);
// proper motion on R.A./Decl (mas/year)
double pmra = GP->pmra/1000.*DAS2R, pmdecl = GP->pmdecl/1000.*DAS2R;*/
mjd = JDate - jd0;
slamap(r, d, pmra, pmdecl, 0., 0., 2000.0, mjd, &r, &d);
DBG("APP: %g, %g", r*DR2H, d*DR2D);
r *= DR2S;
d *= DR2AS;
if(appRA) *appRA = r;
if(appDecl) *appDecl = d;
}
/**
* convert apparent coordinates (nowadays) to mean (JD2000)
* appRA, appDecl in seconds
* r, d in hours & degrees
*/
void calc_mean(double appRA, double appDecl, double *r, double *d){
double ra, dec;
appRA *= DS2R;
appDecl *= DAS2R;
DBG("appRa: %g, appDecl: %g", appRA, appDecl);
double mjd = JDate - jd0;
slaamp(appRA, appDecl, mjd, 2000.0, &ra, &dec);
ra *= DR2H;
dec *= DR2D;
if(r) *r = ra;
if(d) *d = dec;
}
/**
* calculate polar coordinates alpha, delta
* from horizontal az, zd for given siderial time stime
* ALL IN SECONDS!
*/
void calc_AD(double az, double zd, double stime, double *alpha, double *delta){
double sin_d, sin_a, cos_a, sin_z, cos_z;
double t, d, x, y;
const double cos_fi = 0.7235272793; /* Cos of SAO latitude */
const double sin_fi = 0.6902957888; /* Sin --- "" ----- */
DBG("AZ: %g, ZD: %g", az, zd);
az *= DAS2R;
zd *= DAS2R;
sincos(az, &sin_a, &cos_a);
sincos(zd, &sin_z, &cos_z);
y = sin_z * sin_a;
x = cos_a * sin_fi * sin_z + cos_fi * cos_z;
t = atan2(y, x);
if (t < 0.0)
t += 2.0*M_PI;
sin_d = sin_fi * cos_z - cos_fi * cos_a * sin_z;
d = asin(sin_d);
*delta = d * DR2AS;
*alpha = (stime - t * DR2S);
if (*alpha < 0.0)
*alpha += 86400.; // +24h
DBG("A: %g, Z: %g, alp: %g, del: %g", az, zd, *alpha, *delta);
}