ExternalVendorCode/k3ng_rotator_controller
1
0
Fork 0
mirror of https://github.com/k3ng/k3ng_rotator_controller.git synced 2024-12-22 06:27:50 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity
044b768c41
k3ng_rotator_controller/libraries/moon2.cpp
Anthony Good 090304d481 restructured with library file directory
2015-03-06 19:18:53 -05:00

331 lines
11 KiB
C++
Executable File
Raw Blame History

#include <Arduino.h>
#include <math.h>
#include <string.h>
#include <ctype.h>
// Translated from the WSJT Fortran code by Pete VE5VA
/////////////////////////////////////////////////////////
///////////// G R I D 2 D E G /////////////
/////////////////////////////////////////////////////////
// Choose which version of the grid conversion to use.
// Defining WSJT uses the one from WSJT
// Removing the define of WSJT uses the code based on
// the PERL script at wikipedia which seems to be
// slightly more accurate.
//#define WSJT
#ifdef WSJT
// grid = DO62QC <-> lat = 52.104167 lon = -106.658332
// grid = JO62MM <-> lat = 52.520832 lon = 13.008334
// The WSJT code returns West as positive but moon2 uses
// West is negative
// CHANGE this so that it returns West longitude as NEGATIVE
void grid2deg(char *grid0,double *dlong,double *dlat)
{
char grid[8];
int nlong,n20d;
int nlat;
double xminlong,xminlat;
// Initialize the grid with a default string
strncpy(grid,"AA00MM",6);
// copy in the grid
strncpy(grid,grid0,6);
// Convert the grid to upper case
for(int i = 0;i < 6;i++)grid[i] = toupper(grid[i]);
// Fix any errors
if(!isalpha(grid[0]))grid[0] = 'A';
if(!isalpha(grid[1]))grid[1] = 'A';
if(!isdigit(grid[2]))grid[2] = '0';
if(!isdigit(grid[3]))grid[3] = '0';
if(!isalpha(grid[4]))grid[4] = 'M';
if(!isalpha(grid[5]))grid[5] = 'M';
nlong = 180 - 20*(grid[0] - 'A');
n20d = 2*(grid[2] - '0');
xminlong = 5*(grid[4]-'A')+ 0.5;
// Make west longitude negative
*dlong = -(nlong - n20d -xminlong/60.);
nlat = -90 + 10*(grid[1] - 'A') + grid[3] - '0';
xminlat = 2.5*(grid[5] - 'A' + 0.5);
*dlat = nlat + xminlat/60.;
}
#else
// grid = DO62QC <-> lat = 52.104164 lon = -106.625000
// grid = JO62MM <-> lat = 52.520832 lon = 13.041667
// From: http://en.wikipedia.org/wiki/Maidenhead_Locator_System
void grid2deg(char *grid0,double *dlong,double *dlat)
{
char grid[8];
// Initialize the grid with a default string
strncpy(grid,"AA00MM",6);
// copy in the grid
strncpy(grid,grid0,6);
// Convert the grid to upper case
for(int i = 0;i < 6;i++)grid[i] = toupper(grid[i]);
// Fix any errors
if(!isalpha(grid[0]))grid[0] = 'A';
if(!isalpha(grid[1]))grid[1] = 'A';
if(!isdigit(grid[2]))grid[2] = '0';
if(!isdigit(grid[3]))grid[3] = '0';
if(!isalpha(grid[4]))grid[4] = 'M';
if(!isalpha(grid[5]))grid[5] = 'M';
*dlong = 20*(grid[0] - 'A') - 180;
*dlat = 10*(grid[1] - 'A') - 90;
*dlong += (grid[2] - '0') * 2;
*dlat += (grid[3] - '0');
// subsquares
*dlong += (grid[4] - 'A') * 5/60.;
*dlat += (grid[5] - 'A') * 2.5/60.;
*dlong += 2.5/60.;
*dlat += 1.25/60;
}
#endif
/////////////////////////////////////////////////////////
////////////// D C O O R D ///////////////
/////////////////////////////////////////////////////////
// In WSJT this is used in various places to do coordinate
// system conversions but moon2 only uses it once.
void DCOORD(double xA0,double xB0,double AP,double BP,
double xA1,double xB1,double *xA2,double *B2)
{
double TA2O2;
double SB0,CB0,SBP,CBP,SB1,CB1,SB2,CB2;
double SAA,CAA,SBB,CBB,CA2,SA2;
SB0 = sin(xB0);
CB0 = cos(xB0);
SBP = sin(BP);
CBP = cos(BP);
SB1 = sin(xB1);
CB1 = cos(xB1);
SB2 = SBP*SB1 + CBP*CB1*cos(AP-xA1);
CB2 = sqrt(1.-(SB2*SB2));
*B2 = atan(SB2/CB2);
SAA = sin(AP-xA1)*CB1/CB2;
CAA = (SB1-SB2*SBP)/(CB2*CBP);
CBB = SB0/CBP;
SBB = sin(AP-xA0)*CB0;
SA2 = SAA*CBB-CAA*SBB;
CA2 = CAA*CBB+SAA*SBB;
TA2O2 = 0.0;
if(CA2 <= 0.) TA2O2 = (1.-CA2)/SA2;
if(CA2 > 0.) TA2O2 = SA2/(1.+CA2);
*xA2=2.*atan(TA2O2);
if(*xA2 < 0.) *xA2 = *xA2+6.2831853071795864;
}
/////////////////////////////////////////////////////////
//////////////// M O O N 2 ////////////////
/////////////////////////////////////////////////////////
// You can derive the lat/long from the grid square using
// the grid2deg function to translate from grid to lat/long
// Example call to this function for 2014/01/04 1709Z:
// moon2(2014,1,4,17+9/60.,-106.625,52.104168,
// &RA, &Dec, &topRA, &topDec, &LST, &HA, &Az, &El, &dist);
// I have not used or checked any of the outputs other than Az and El
void moon2(int y,int m,int Day,
double UT,
double lon,double lat,
double *RA,double *Dec,
double *topRA,double *topDec,
double *LST,double *HA,
double *Az,double *El,double *dist)
{
// The strange position of some of the semicolons is because some
// of this was translated using a TCL script that I wrote - it isn't
// too smart but it saved some typing
double NN ;//Longitude of ascending node
double i ;//Inclination to the ecliptic
double w ;//Argument of perigee
double a ;//Semi-major axis
double e ;//Eccentricity
double MM ;//Mean anomaly
double v ;//True anomaly
double EE ;//Eccentric anomaly
double ecl ;//Obliquity of the ecliptic
double d ;//Ephemeris time argument in days
double r ;//Distance to sun, AU
double xv,yv ;//x and y coords in ecliptic
double lonecl,latecl ;//Ecliptic long and lat of moon
double xg,yg,zg ;//Ecliptic rectangular coords
double Ms ;//Mean anomaly of sun
double ws ;//Argument of perihelion of sun
double Ls ;//Mean longitude of sun (Ns=0)
double Lm ;//Mean longitude of moon
double DD ;//Mean elongation of moon
double FF ;//Argument of latitude for moon
double xe,ye,ze ;//Equatorial geocentric coords of moon
double mpar ;//Parallax of moon (r_E / d)
// double lat,lon ;//Station coordinates on earth
double gclat ;//Geocentric latitude
double rho ;//Earth radius factor
double GMST0; //,LST,HA;
double g;
// double topRA,topDec ;//Topocentric coordinates of Moon
// double Az,El;
// double dist;
double rad = 57.2957795131,twopi = 6.283185307,pi,pio2;
// data rad/57.2957795131d0/,twopi/6.283185307d0/
//Note the use of 'L' to force 32-bit integer arithmetic here.
d=367L*y - 7L*(y+(m+9L)/12L)/4L + 275L*m/9L + Day - 730530L + UT/24.;
ecl = 23.4393 - 3.563e-7 * d;
//Serial.print("d = ");
//Serial.println(d,3);
// Orbital elements for Moon:
NN = 125.1228 - 0.0529538083 * d;
i = 5.1454;
w = fmod(318.0634 + 0.1643573223 * d + 360000.,360.);
a = 60.2666;
e = 0.054900;
MM = fmod(115.3654 + 13.0649929509 * d + 360000.,360.);
// Orbital elements for Sun:
/*
NN = 0.0;
i = 0.0;
w = fmod(282.9404 + 4.70935e-5 * d + 360000.,360.);
a = 1.000000;
e = 0.016709 - 1.151e-9 * d;
MM = fmod(356.0470 + 0.99856002585 * d + 360000.,360.);
*/
/*
Serial.print("\nmoon2: w=");
Serial.print(w,3);
Serial.print(" e=");
Serial.print(e,3);
Serial.print(" MM=");
Serial.println(MM,3);
*/
EE = MM + e*rad*sin(MM/rad) * (1. + e*cos(MM/rad));
EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1. - e*cos(EE/rad));
EE = EE - (EE - e*rad*sin(EE/rad)-MM) / (1. - e*cos(EE/rad));
xv = a * (cos(EE/rad) - e);
yv = a * (sqrt(1.-e*e) * sin(EE/rad));
v = fmod(rad*atan2(yv,xv)+720.,360.);
r = sqrt(xv*xv + yv*yv);
// Get geocentric position in ecliptic rectangular coordinates:
xg = r * (cos(NN/rad)*cos((v+w)/rad) - sin(NN/rad)*sin((v+w)/rad)*cos(i/rad));
yg = r * (sin(NN/rad)*cos((v+w)/rad) + cos(NN/rad)*sin((v+w)/rad)*cos(i/rad));
zg = r * (sin((v+w)/rad)*sin(i/rad));
// Ecliptic longitude and latitude of moon:
lonecl = fmod(rad*atan2(yg/rad,xg/rad)+720.,360.);
latecl = rad * atan2(zg/rad,sqrt(xg*xg + yg*yg)/rad);
// Now include orbital perturbations:
Ms = fmod(356.0470 + 0.9856002585 * d + 3600000.,360.);
ws = 282.9404 + 4.70935e-5*d;
Ls = fmod(Ms + ws + 720.,360.);
Lm = fmod(MM + w + NN+720.,360.);
DD = fmod(Lm - Ls + 360.,360.);
FF = fmod(Lm - NN + 360.,360.);
lonecl = lonecl
-1.274 * sin((MM-2.*DD)/rad)
+0.658 * sin(2.*DD/rad)
-0.186 * sin(Ms/rad)
-0.059 * sin((2.*MM-2.*DD)/rad)
-0.057 * sin((MM-2.*DD+Ms)/rad)
+0.053 * sin((MM+2.*DD)/rad)
+0.046 * sin((2.*DD-Ms)/rad)
+0.041 * sin((MM-Ms)/rad)
-0.035 * sin(DD/rad)
-0.031 * sin((MM+Ms)/rad)
-0.015 * sin((2.*FF-2.*DD)/rad)
+0.011 * sin((MM-4.*DD)/rad);
latecl = latecl
-0.173 * sin((FF-2.*DD)/rad)
-0.055 * sin((MM-FF-2.*DD)/rad)
-0.046 * sin((MM+FF-2.*DD)/rad)
+0.033 * sin((FF+2.*DD)/rad)
+0.017 * sin((2.*MM+FF)/rad);
r = 60.36298
- 3.27746*cos(MM/rad)
- 0.57994*cos((MM-2.*DD)/rad)
- 0.46357*cos(2.*DD/rad)
- 0.08904*cos(2.*MM/rad)
+ 0.03865*cos((2.*MM-2.*DD)/rad)
- 0.03237*cos((2.*DD-Ms)/rad)
- 0.02688*cos((MM+2.*DD)/rad)
- 0.02358*cos((MM-2.*DD+Ms)/rad)
- 0.02030*cos((MM-Ms)/rad)
+ 0.01719*cos(DD/rad)
+ 0.01671*cos((MM+Ms)/rad);
*dist = r * 6378.140;
// Geocentric coordinates:
// Rectangular ecliptic coordinates of the moon:
xg = r * cos(lonecl/rad)*cos(latecl/rad);
yg = r * sin(lonecl/rad)*cos(latecl/rad);
zg = r * sin(latecl/rad);
// Rectangular equatorial coordinates of the moon:
xe = xg;
ye = yg * cos(ecl/rad) - zg*sin(ecl/rad);
ze = yg * sin(ecl/rad) + zg*cos(ecl/rad);
// Right Ascension, Declination:
*RA = fmod(rad*atan2(ye,xe)+360.,360.);
*Dec = rad * atan2(ze,sqrt(xe*xe + ye*ye));
// Now convert to topocentric system:
mpar=rad * asin(1./r);
// alt_topoc = alt_geoc - mpar*cos(alt_geoc);
gclat = lat - 0.1924 * sin(2.*lat/rad);
rho = 0.99883 + 0.00167 * cos(2.*lat/rad);
GMST0 = (Ls + 180.)/15.;
*LST = fmod(GMST0+UT+lon/15.+48.,24.) ;//LST in hours
*HA = 15. * *LST - *RA ;//HA in degrees
g = rad*atan(tan(gclat/rad)/cos(*HA/rad));
*topRA = *RA - mpar*rho*cos(gclat/rad)*sin(*HA/rad)/cos(*Dec/rad);
*topDec = *Dec - mpar*rho*sin(gclat/rad)*sin((g-*Dec)/rad)/sin(g/rad);
*HA = 15. * *LST - *topRA ;//HA in degrees
if(*HA > 180.) *HA=*HA-360.;
if(*HA < -180.) *HA=*HA+360.;
pi = 0.5 * twopi;
pio2 = 0.5 * pi;
DCOORD(pi,pio2-lat/rad,0.,lat/rad,*HA*twopi/360,*topDec/rad,Az,El);
*Az = *Az * rad;
*El = *El * rad;
return;
}
Reference in New Issue View Git Blame Copy Permalink