Celestial Navigation Plugin shows incorrect altitude (Hc) on FIND celestial body form #24
Comments
|
If we find someone willing and able to fix this, are you available for testing? This is a great example. thanks. |
|
Absolutely!! Bob Bossert
…On Fri, Jan 17, 2025 at 8:34 AM Rick Gleason ***@***.***> wrote:
If we find someone willing and able to fix this, are you available for
testing?
—
Reply to this email directly, view it on GitHub
<#24 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/BOSWTTQKQBZMHPZMWLPQHAD2LEBFFAVCNFSM6AAAAABVL4P7PKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDKOJYGM3TSOBUGE>
.
You are receiving this because you authored the thread.Message ID:
***@***.***>
|
|
Yes, indeed, please submit them. We need experts such as you! I hope that some old programming sailors who are passionate about sextants and accuracy will focus on improving this tool, as it could become a great resource. Also if you have any suggestions or improvements to the manual, please advise. How to get to the manuals from OpenCPN Wiki Then https://opencpn-manuals.github.io/main/opencpn-plugins/index.html |
|
Sean Patrick has offered to take a look at the code. Bob If you can create a zip file and upload to here it would be very useful! |
|
I updated Github with the Code. I also updated the problem description
with some additional information. The Celestial Tools developer only send
high level code. The missing code we needed for the SD and HP. I asked
for specific missing functions that would be required to completly define
the SD/HP calculation. But the developer was reluctant to share. He
thinks by sharing the code, more questions will be generated and that will
require him to spend more time with us.
He said in his email, "don't pull the sting, you don't want to get the
ball".
Frankly, Celestial Tools clearly uses Jan Meeus book, Celestial
Algorithms. And, it made some changes....Probably simplifications.....And
the simplifications were done before the developer was involved and won't
know the answers. Jan Meeus did a good job of documenting his algorithm.
Testament is the spreadsheet I was able to create, test, and build unit
test cases.
I reviewed Celestial Tools code and read Meeus's book in detail. You'll
notice the spreadsheet is in the same order that Meeus described his
algorithm. The code also follows the algorithm's order. The spreadsheet
refers to the same variable names as in the algorithm and the code. If you
want the spreadsheet in Excel, let me know. Meeus's algorithm relies on 2
tables for 3 variables. They are completly defined in the spreadsheet.
This spreadsheet does more than calculate Horizontal Parallax, and
Semi-Diameter. That's because both are derived using the Distance of
Earth to Moon, I also included I also includled calculations for Moon
Longitude and latitude/declination in the Ecliptic Coordinate system and
Lat/Lon coordinate system. For implementation of the coding change in Open
CPN, only those variables needed for HP and SD. It starts with time, and
the algorithm starts with UT1 (open cpn plugin input), Julian Date (JD)
which I'm sure Open CPN already calculates, and Delta_T is added to JD to
create JDE. and then T (time) which is number of centures after the Jan
1, 2000 epoch (according to Meeus, carry T to at least 9 significant
digits....if not more. I'm sure Open CPN does all of these time
definitions already and changing this isn't necessary. Just use it. Then
the blue bold highlited text are the parts needed to calculate HP and SD.
Like I said, I think GHA and Dec for Moon are good.
So, please review the spreadsheet step by step. Inside the cells is the
formulas and equations used. I can send you Meeus's chapter 25 book which
the code and the spreadsheet closely follows. You can google a PDF version
of it. Let me know if you need me to send the Chapters (I primarily used
Chapter 25) I used.
The spreadsheet has test cases. you can easily build additional test
cases. just duplicate the last worksheet and input the UT1 date and time.
Pull off the corresponding Julian Date from Navy website I reference. then
compare results to a Nautical Almanac. There are free Nautical Almanacs on
the web..
I got close matches for SD and HP. I show output of Lat/Lon (RA and Dec
format) the algorithm calculates, but this part of the algorithm you don't
need because Open CPN does a good job with Lat/Lon.
Bob Bossert
PS - Numbers don't upload into gethub. If you can't read numbers, I can
send exel copy.
…On Sat, Jan 25, 2025 at 10:31 AM Rick Gleason ***@***.***> wrote:
@rbossert149419 <https://github.com/rbossert149419>
Sean Patrick has offered to take a look at the code.
https://www.cruisersforum.com/forums/f134/celestial-navigation-plugin-redux-98748.html#post3967665
Sean Patrick:
If you have any questions about celestial navigation, ask me!
Celestial Navigation Spreadsheet <https://www.cooknavigation.com/>
NavList Celestial Navigation Forum <https://navlist.net/>
Bob If you can create a zip file and upload to here it would be very
useful!
—
Reply to this email directly, view it on GitHub
<#24 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/BOSWTTW7XVAC3HSMZZW4TCL2MOU3ZAVCNFSM6AAAAABVL4P7PKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMMJUGAYDGMJTG4>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
|
Bob, I sent you a PM on CruiserForum with my email. artisthos on CF has agreed to do testing. |
|
Bob, would you mind uploading to here the converted "numbers" and the excel spreadsheets that you refer too? Also I am not so sure OpenCPN handles time in the way you think it might. It appears that we need to:
I was confused by this first statement "I updated Github with the Code. I also updated the problem description |
|
Julian Date is the standard method used in astronomy for tracking
celestial objects over a very long period of time. My C++ code takes the
Julian Date as input and converts it to the time variable TT. I just
assumed that open cpn used equations for GHA and Dec of the bodies based on
system of equations that are depend on the time variable TT or Julian Date.
But, if the code doesn't have TT / Julian Date, I can add to the C++ code
to convert UTC1 as input and then derive TT. Bob
…On Mon, Feb 10, 2025 at 8:16 PM Rick Gleason ***@***.***> wrote:
Bob, would you mind uploading the converted "numbers" and the excel
spreadsheet?
Also I am not so sure OpenCPN handles time in the way you think it might.
I may have to ask sebastien-rossert about that.
—
Reply to this email directly, view it on GitHub
<#24 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/BOSWTTTGTMZ4HPXI3M7BZRD2PFFPVAVCNFSM6AAAAABVL4P7PKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNBZGYYDGMRYHE>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
|
Bob, I think we better do that. I've written sebastien-rosset as he is in the process of implementing global time, and he'll probably have some thoughts about this. I will send him your note above. He may chime in here. |
|
I'm not successful uploading the .xls to github. Not sure why. So, I'll
email the spreadsheets and the C++ code here. Hope you can open these
files,
On Mon, Feb 10, 2025 at 9:30 PM Rick Gleason ***@***.***> wrote:
Bob, I think we better do that. I've written sebastien-rosset as he is in
the process of implementing global time, and he'll probably have some
thoughts about this. I will send him your note above. He may chime in here.
@sebastien-rosset <https://github.com/sebastien-rosset>
—
Reply to this email directly, view it on GitHub
<#24 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/BOSWTTXMFECNPQQYLYU6BMD2PFOFNAVCNFSM6AAAAABVL4P7PKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNBZGY3TKNBYGI>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
//
// main.cpp
// Moon SemiDiameter (arc min) and Moon Parallax (min)
// for OpenCPN Celestial Plug-in
// Created by Robert Bossert on 2/9/25.
#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;
double modulo_360(double big_degrees) {
//
// modulo_360
// Convert large Degree numbers (+ or - ) to 0 - 360 degrees.
// Calculate number of complete 360 degree rotations.
// Subtract or Add the rotations to BIG. Result is the remaining degrees
// If remaining Degrees is between 0 and -360 rotations were = 0, add 360
// BIG shouldn't be bigger than 2.13 billion degrees.
//
int z = 0;
double adj_degrees = big_degrees ;
z = abs(trunc(big_degrees/360));
if (big_degrees >= 360) {
adj_degrees = big_degrees - z * 360;
}
if (big_degrees < 0) {
adj_degrees = big_degrees + z * 360;
}
if (big_degrees < 0) {
adj_degrees = adj_degrees + 360;
}
cout << "input, rounds, adjdeg: "<<big_degrees<<", " <<z <<", "<< adj_degrees<< "\n";
return adj_degrees ;
}
double moon_distance(long double TT) {
//
// Calculates Geocentric Moon Distance from Earth in km at Julian Time TT.
// TT should be carried to at least 9 significant digits.
//
double BASE_MOON_DISTANCE = 385000.56; // Base Earth-Moon distance 385000.56 km
double m_distance = BASE_MOON_DISTANCE ; // Avg Earth-Moon distance 384400 km
double ml_prime, md, mm, mm_prime, mf, a1, a2, a3, e ;
//variables and equations Defined in J. Meeus Celestial Equations chapter 45....
// ml_prime Moon's Mean Longitude (degrees)
// md Mean elongation of the moon (degrees)
// mm Sun's mean anomaly (degrees)
// mm_prime Moon's mean anomaly (degrees)
// mf Moon's argument of latitude or
// Mean distance moon from it's ascending node (degrees)
// a1 due to action of Venus (degrees)
// a2 due to action of Jupiter (degrees)
// e Eccentricity of earth's orbit around sun. (ratio)
//
ml_prime = 218.3164591 + (481267.8813436*TT) - (.0013268*TT*TT)
+ (TT*TT*TT/538841) - (TT*TT*TT*TT/65194000) ;
md = 297.8502042 + (445267.1115168*TT) - (.00163*TT*TT)
+ (TT*TT*TT/545868) - (TT*TT*TT*TT/113065000) ;
mm = 357.5291092 + (35999.0502909*TT) - (.0001536*TT*TT)
+ (TT*TT*TT/24490000) ;
mm_prime = 134.9634114 + (477198.8676313*TT) + (.008997*TT*TT)
+ (TT*TT*TT/69699) - (TT*TT*TT*TT/14712000) ;
mf = 93.2720993 + (483202.0175273*TT) - (.0034029*TT*TT)
+ (TT*TT*TT/3526000) + (TT*TT*TT*TT/863310000) ;
a1 = 119.75 + (131.849*TT) ;
a2 = 53.09+(479264.29*TT) ;
a3 = 313.45 + (481266.484*TT) ;
e = (1 - (0.002516*TT) - (0.0000074*TT*TT) ) ; // orbit ecentricity.
ml_prime = modulo_360(ml_prime);
md = modulo_360(md);
mm = modulo_360(mm);
mm_prime = modulo_360(mm_prime);
mf = modulo_360(mf);
a1 = modulo_360(a1);
a2 = modulo_360(a2);
a3 = modulo_360(a3);
cout << "TT: " << std::setprecision(15) << std::fixed << TT ; cout << "\n";
cout << "ml_prime: " << ml_prime ; cout << "\n";
cout << "md: " << md; cout << "\n";
cout << "mm: " << mm; cout << "\n";
cout << "mm_prime: " << mm_prime; cout << "\n";
cout << "mf: " << mf; cout << "\n";
cout << "a1: " << a1; cout << "\n";
cout << "a2: " << a2; cout << "\n";
cout << "a3: " << a3; cout << "\n";
cout << "e: " << e; cout << "\n";
//
// Source J. Meeus, Celestial Equations. Table 45a page 309 to 310
// constant arrays. The numbers do not change.
//
int D[] = {0,2,2,0,0,0,2,2,2,2,0,1,0,2,0,0,4,0,4,2,2,1,1,2,2,4,2,0,2,2,1,2,0, 0,2,2,2,4,0,3,2,4,0,2,2,2,4,0,4,1,2,0,1,3,4,2,0,1,2,2};
int M[] = {0,0,0,0,1,0,0,-1,0,-1,1,0,1,0,0,0,0,0,0,1,1,0,1,-1,0,0,0,1,0,-1,0,-2,1, 2,-2,0,0,-1,0,0,1,-1,2,2,1,-1,0,0,-1,0,1,0,1,0,0,-1,2,1,0,0};
int M_PRIME[] = {1,-1,0,2,0,0,-2,-1,1,0,-1,0,1,0,1,1,-1,3,-2,-1,0,-1,0,1,2,0,-3,-2, -1,-2,1,0,2,0,-1,1,0,-1,2,-1,1,-2,-1,-1,-2,0,1,4,0,-2,0,2,1,-2,-3,2,1,-1,3,-1 };
int F[]= {0,0,0,0,0,2,0,0,0,0,0,0,0,-2,2,-2,0,0,0,0,0,0,0,0,0,0,0,0,2,0,0,0,0,0,0, -2,2,0,2,0,0,0,0,0,0,-2,0,0,0,0,-2,-2,0,0,0,0,0,0,0,-2};
int R_AMP[]={-20905355, -3699111, -2955968, -569925, 48888, -3149, 246158, -152138, -170733, -204586, -129620, 108743, 104755, 10321, 0, 79661, -34782, -23210, -21636, 24208, 30824, -8379, -16675, -12831, -10445, -11650, 14403, -7003, 0, 10056, 6322, -9884, 5751, 0, -4950, 4130, 0, -3958, 0, 3258, 2616, -1897, -2117, 2354, 0, 0, -1423, -1117, -1571, -1739, 0, -4421,0,0,0,0,1165,0,0,8752};
double sigma_r, e_mult, r;
//
// sigma_r - Periodic terms of moon motion (longitude, latitude, distance)
// for a given TT. Sigma_R is .001 Km to avoid calculation rounding.
// Table 45.A of Meeus Astronomical Algorithms (61 rows) contains the most
// important operiodic terms.
// then more adjustments to SigmaR for Flattening of Earth, motion of Venus.
// e_mult - Eccentricity multiplier. Depends on M. If M +/- 1, then multiply M
// by E. If M is +/- 2, then Multiply M by E*E. If M is 0, then Efactor is 1.
// r - Calculated for each period term's (calculation of rows in array).
// Temporary for testing and validation.
sigma_r = 0;
for (int i = 0; i<61; i++) {
e_mult = 1;
if (M[i] == 1 || M[i] == -1) {
e_mult=e;
}
if (M[i] == 2 || M[i] == -2) {
e_mult = e*e;
}
r = R_AMP[i] * cos( ( (D[i] * md) + (e_mult * M[i] * mm) +
(M_PRIME[i] * mm_prime) + (F[i] * mf)) * M_PI/180 ) ;
sigma_r = sigma_r + r;
// cout << "Table 45.A " <<i <<", " <<D[i] <<", " <<M[i] <<", " <<M_PRIME[i] <<", " <<F[i] <<", " <<e_mult <<", " <<R_AMP[i] <<", r:" <<r <<", Accum r: " <<sigma_r; cout << "\n";
}
//
// sigma_r - more periodic motion adjustements to sigma_r.
// terms are in .001 km for calculation precision.
// A1-mf are Actions of Venus
// ml_pime and mm_prime due to flattening of Earth
//
sigma_r = sigma_r + -2235 * sin( ml_prime * M_PI/180)
+ 382 * sin( a3 * M_PI/180)
+ 175 * sin( (a1 - mf) * M_PI/180)
+ 175 * sin( (a2 + mf) * M_PI/180)
+ 127 * sin( (ml_prime - mm_prime) * M_PI/180)
- 115 * sin( (ml_prime + mm_prime) * M_PI/180);
//
// m_distsmvr - After sigma_r and sigma_r adjustments,
// convert from .001 km units to whole km units.
// Then add a constant (representing the base moon distance to earth).
// Note that in moon measured to be 384,400km (2/2025).
// I used constant supplied by Meeus since sigma_r assumed the base number.
//
// m_distance is the Geocentric Distance of Moon from the Earth at time tt.
// Geocentric means, earth's center to moon's center.
m_distance = BASE_MOON_DISTANCE + sigma_r / 1000 ;
return m_distance;
}
int main() {
double JD, JDE;
double TT; // Terrestial Time in centuries.
//
// JD Julian Date
// JDE Julian Ephemerides Date. jde = jd found to be good enough in testing.
// TT Terrestial Julian Time in centuries starting from Jan 1 2000 Epoch
//
// For unit testing, used https://aa.usno.navy.mil/data/JulianDate.
// input UTC1 and the site returns Julian Date.
//
cout << "Enter Julian Date: ";
cin >> JD ;
JDE = JD ;
const double J2000 = 2451545.0 ; // Julian date 2451545.0 TT=2000 Jan 1, 12h TT
const double DAYS_IN_CENTURY = 36525 ; // numnber of days in century.
TT = (JDE - J2000)/DAYS_IN_CENTURY ; // Terestrial Time
//
// Function moon-distance determines Geocentric Moon Distance to Earth.
// based on Time.
// Moon semi-diameter and parallax is based on Moon Distance in km.
// moon_semidiameter is a simplified equation.
// Alt_Moon_Semidiameter is a more precise equation.
// Meeus used K=.27248. K = Moon Mean radius/Earth Radius at Equator
//
double moon_dist, moon_semidiameter, moon_parallax, alt_moon_semidiameter;
double EARTH_RADIUS = 6378.14 ; //Earth radius at equator in km.
double MOON_MEAN_RADIUS = 1737.5; // Moon mean radius in km
double K = MOON_MEAN_RADIUS / EARTH_RADIUS ;
//
// Function "moon_distance". Geocentric Moon to Earth distance in km at time TT.
//
moon_dist = moon_distance(TT) ; // geocentric distance in km
//
// SemiDiameter and Parallax is based on Distance, Earth Radius, Moon Radius
//
moon_semidiameter = 5974556.667/moon_dist;
// in arc-minutes (same units as Nautical Almanac.
moon_parallax = 60 * asin(EARTH_RADIUS/moon_dist) * 180/M_PI;
// in arc-minutes (same units as Nautical Almanac.
alt_moon_semidiameter=60 * asin(K*sin((moon_parallax/60)*(M_PI/180))) * 180/M_PI; // in arc minutes
cout << "Moon Distance (km): " << moon_dist; cout << "\n";
cout << "Moon Diameter (km): " << EARTH_RADIUS; cout << "\n";
cout << "MoonParallax (min) " << moon_parallax; cout << "\n";
cout << "MoonSemidiameter " << moon_semidiameter; cout << "\n";
cout << "Alt MoonSemidiameter " << alt_moon_semidiameter; cout << "\n";
return 0;
}
//
// main.cpp
// Altitude Computed
//
// Created by Robert & Terry Bossert on 2/9/25.
//
#include <iostream>
#include <cmath>
#include <iomanip>
using namespace std;
double modulo_360(double big_degrees) {
//
// modulo_360
// Convert large Degree numbers (+ or - ) to 0 - 360 degrees.
// Calculate number of complete 360 degree rotations.
// Subtract or Add the rotations to BIG. Result is the remaining degrees
// If remaining Degrees is between 0 and -360 rotations were = 0, add 360
// BIG shouldn't be bigger than 2.13 billion degrees.
//
int z = 0;
double adj_degrees = big_degrees ;
z = abs(trunc(big_degrees/360));
if (big_degrees >= 360) {
adj_degrees = big_degrees - z * 360;
}
if (big_degrees < 0) {
adj_degrees = big_degrees + z * 360;
}
if (big_degrees < 0) {
adj_degrees = adj_degrees + 360;
}
// cout << "input, rounds, adjdeg: "<<big_degrees<<", " <<z <<", "<< adj_degrees<< "\n";
return adj_degrees ;
}
long double Altitude_computed(double latitude, double longitude, double gha, double declination) {
//
// Calculated Altitude using Law of Cosines
//
// latitude - DR location. decimal degrees. Angle 0 to -90 deg. is South of Equator,
// 0 to +90 deg. if North of Equator
// longitude - DR location. decimal degrees. Angle is 0 to -180 deg. if West of
// Greenwich. 0 to +180 deg. if East of Greenwich.
// gha Greenwich Hour Angle of the Body. 0 to 360 decimal degrees.
// declination of the Body. decimal degrees. + is North of equator, - is South
// lha = gha + longitude. Local Hour Angle. decimal degrees.
//
double lha, adj_declination, adj_latitude;
//
// sin_hc, hc, z, zn need to be long double to get 4+ significant digits precision
// and to calculate intercept (ho minus hc)
//
long double sin_hc, hc;
lha = gha + longitude;
//
// lha Local Hour Angle must be between 0 and 360 degrees.
// Modulus operation if larger than 360 or negative
//
lha = modulo_360(lha);
cout << "lha (degrees between 0 to 360): "; cout<<lha<<"\n";
//
// adj_declination is negative when latitude and declination have contrary names.
//
if ( (latitude > 0 and declination < 0) or (latitude < 0 and declination > 0 ) ) {
adj_declination = -(abs(declination));
}
else {
adj_declination = (abs(declination));
}
//
// adj_latitude is always positive
//
adj_latitude = abs(latitude);
//
// Sin of Computed Altitude using Law of cosines
//
// double s, t, u, v, w, x, y, z;
// s =sin(adj_declination * M_PI/180);
// t =sin(adj_latitude * M_PI/180);
// u =cos(adj_latitude * M_PI/180);
// v =cos(adj_declination * M_PI/180);
// w =cos(lha * M_PI/180);
// x = s * t;
// y = u * v * w ;
// z = x + y;
// cout << "s "<<s << "\n";
// cout << "t "<<t << "\n";
// cout << "u "<<u << "\n";
// cout << "v "<<v << "\n";
// cout << "w "<<w << "\n";
// cout << "x "<<x << "\n";
// cout << "y "<<y << "\n";
// cout << "z "<<z << "\n";
sin_hc = ( ( sin(adj_declination * M_PI/180) * sin(adj_latitude * M_PI/180) )
+ ( cos(adj_latitude * M_PI/180)
* cos(adj_declination * M_PI/180)
* cos(lha * M_PI/180) ) );
//
// Computed Altitude (Hc)
//
hc = asin(sin_hc) * 180/M_PI ;
return hc;
}
long double Azimuth(double latitude, double longitude, double gha, double declination, long double hc) {
//
// Calculated Azimuth using Law of Cosines
//
// latitude - DR location. decimal degrees. Angle 0 to -90 deg. is South of Equator,
// 0 to +90 deg. if North of Equator
// longitude - DR location. decimal degrees. Angle is 0 to -180 deg. if West of
// Greenwich. 0 to +180 deg. if East of Greenwich.
// gha Greenwich Hour Angle of the Body. 0 to 360 decimal degrees.
// declination of the Body. decimal degrees.
// lha = gha + longitude. Local Hour Angle. decimal degrees.
//
double lha ;
double adj_declination, adj_latitude;
//
// since sin_hc, hc are long double (because asin(sin_hc)), make z, zn long double.
// and to calculate intercept (ho minus hc)
//
long double z, zn, cos_z;
//
// lha transformed to between 0 to 360 degrees by modulo_360.
// if not, then Zn will be wrong.
//
lha = gha + longitude;
lha = modulo_360(lha);
cout << "lha (degrees between 0 to 360): "; cout<<lha<<"\n";
//
// adj_declination is negative when latitude and declination have contrary names.
//
if ( (latitude > 0 and declination < 0) or (latitude < 0 and declination > 0 ) ) {
adj_declination = -(abs(declination));
}
else {
adj_declination = (abs(declination));
}
//
// adj_latitude is always positive
//
adj_latitude = abs(latitude);
//
// Law of cosines for deriving the internal angle of the Nautical Triangle.
//
// double s, t, u, v, w, x, y;
// s =sin(adj_declination * M_PI/180);
// t =sin(adj_latitude * M_PI/180);
// u =sin(hc*M_PI/180);
// v =cos(adj_latitude * M_PI/180);
// w =cos(hc * M_PI/180);
// x = s-t*u;
// y = v * w ;
// cout << "s "<<s << "\n";
// cout << "t "<<t << "\n";
// cout << "u "<<u << "\n";
// cout << "v "<<v << "\n";
// cout << "w "<<w << "\n";
// cout << "x "<<x << "\n";
// cout << "y "<<y << "\n";
// cout << "lha "<<lha << "\n";
cos_z = ( sin(adj_declination * M_PI/180) - sin(adj_latitude * M_PI/180)
* sin(hc*M_PI/180) )
/ ( cos(adj_latitude * M_PI/180) * cos(hc * M_PI/180) ) ;
//
// Internal Angle of the Nautical Triangle
//
z = acos(cos_z) * 180/M_PI ;
//
// lha must be between 0 to 360 for correct calculation of zn
//
if (lha <= 180 and latitude>0) {
zn=360-z; }
else if (lha>180 and latitude>0) {
zn = z; }
else if (lha<=180 and latitude < 0) {
zn = 180 + z; }
else { zn = 180-z; }
return zn;
}
double Find_Celestial_Bodies_form() {
double latitude, longitude, gha, declination, ho, intercept;
long double hc, zn;
char intercept_direction;
cout << "Enter latitude (degrees) (+ N, - S):";
cin>>latitude;
cout << "Enter longitude (degrees):";
cin>>longitude;
cout << "Enter declination (degrees) (+ N, - S):";
cin>>declination;
cout << "Enter gha (degrees):";
cin>>gha;
cout << "Enter ho (degrees)";
cin>>ho;
//
// Computed Altitude (Hc) computed using law of cosines
//
hc = Altitude_computed(latitude,longitude,gha,declination);
//
// Azimuth (Bearing) computed using law of cosines
//
zn = Azimuth(latitude, longitude, gha, declination, hc);
//
// Intercept in arc-minutes or nautical miles
//
intercept = abs((hc - ho)*60) ;
cout << "hc (degrees): "; cout<<hc<<"\n";
cout << "zn (degrees): "; cout<<zn<<"\n";
//
// Estimated Position is Towards the body, or Away from the Body
//
if (ho > hc) {
intercept_direction = 'T' ;
} else {
intercept_direction = 'A' ;
}
cout <<"intercept (nm): "; cout<<intercept<<", ";cout << "Toward or Away: "; cout<<intercept_direction<<"\n";
return 0;
}
int main() {
Find_Celestial_Bodies_form(); {
return 0;
}
}
|
This comment has been minimized.
This comment has been minimized.
This comment has been minimized.
This comment has been minimized.
This comment has been minimized.
This comment has been minimized.
|
I'm catching up on this thread. Wrt OpenCPN already calculates Julian Date, I'm assuming you meant the OpenCPN celestial navigation plugin. Because AFAIK, OpenCPN core does not deal with JD. The celestial navigation plugin does work with JD, and it leverages the celestial_navigation_pi/src/Sight.cpp Line 155 in 75a120a You also mention whether OpenCPN manages time calculations to 9+ significant digits, I would have to look into it, but again, it leverages the wxDateTime struct. |
|
Thank you Sebastien for the clarification. Below are the Hc files Bob has emailed: Open.CPN.Altitude.calculation.error.Hc.pdf |
Closed
|
I was able to reproduce (mostly) the first example in this test: https://github.com/rgleason/celestial_navigation_pi/pull/33/files#diff-538ede6a0b5e101ab74abef50e6f6364e7015154d64c21821f8084e56fc490b5R49 In particular for the given input parameters, I can confirm I see:
For the GHA, I see the same results with a different sign. ctest -V
UpdateCTestConfiguration from :/Users/serosset/git/celestial_navigation_pi/test/DartConfiguration.tcl
UpdateCTestConfiguration from :/Users/serosset/git/celestial_navigation_pi/test/DartConfiguration.tcl
Test project /Users/serosset/git/celestial_navigation_pi/test
Constructing a list of tests
Updating test list for fixtures
Added 0 tests to meet fixture requirements
Checking test dependency graph...
Checking test dependency graph end
No tests were found!!!
SEROSSET-M-R8PT:test serosset$ cd bui
SEROSSET-M-R8PT:test serosset$ pwd
/Users/serosset/git/celestial_navigation_pi/test
SEROSSET-M-R8PT:test serosset$ cd ..
SEROSSET-M-R8PT:celestial_navigation_pi serosset$ cd build
SEROSSET-M-R8PT:build serosset$ ctest -V
UpdateCTestConfiguration from :/Users/serosset/git/celestial_navigation_pi/build/DartConfiguration.tcl
UpdateCTestConfiguration from :/Users/serosset/git/celestial_navigation_pi/build/DartConfiguration.tcl
Test project /Users/serosset/git/celestial_navigation_pi/build
Constructing a list of tests
Done constructing a list of tests
Updating test list for fixtures
Added 0 tests to meet fixture requirements
Checking test dependency graph...
Checking test dependency graph end
test 1
Start 1: celestial_tests
1: Test command: /Users/serosset/git/celestial_navigation_pi/build/test/celestial_tests
1: Working Directory: /Users/serosset/git/celestial_navigation_pi/build/test
1: Test timeout computed to be: 10000000
1: Running main() from /tmp/googletest-20240731-4829-u2xmsk/googletest-1.15.2/googletest/src/gtest_main.cc
1: [==========] Running 1 test from 1 test suite.
1: [----------] Global test environment set-up.
1: [----------] 1 test from AltitudeTest
1: [ RUN ] AltitudeTest.SunLowerLimbExample
1: Using plugin data directory: /Users/serosset/git/celestial_navigation_pi/test/testdata
1: Plugin data directory requested for celestial_navigation_pi, returning: /Users/serosset/git/celestial_navigation_pi/test/testdata
1: Plugin data directory requested for celestial_navigation_pi, returning: /Users/serosset/git/celestial_navigation_pi/test/testdata
1: Plugin data directory requested for celestial_navigation_pi, returning: /Users/serosset/git/celestial_navigation_pi/test/testdata
1:
1: === Sun Lower Limb Example (Nov 13, 2024) ===
1: Main Correction Analysis:
1: Almanac: 11.6'
1: Actual : 5.78594'
1: /Users/serosset/git/celestial_navigation_pi/test/altitude_tests.cpp:95: Failure
1: The difference between actualCorrection and ALMANAC_MAIN_CORRECTION is 0.096901006898011816, which exceeds EPSILON_MIN, where
1: actualCorrection evaluates to 0.096432326435321514,
1: ALMANAC_MAIN_CORRECTION evaluates to 0.19333333333333333, and
1: EPSILON_MIN evaluates to 0.0016666666666666668.
1: Main correction differs from NA by -5.8140604138807088 minutes
1:
1:
1: Ho Analysis:
1: Almanac: 11° 30.6'
1: Actual : 11° 30.8'
1: /Users/serosset/git/celestial_navigation_pi/test/altitude_tests.cpp:105: Failure
1: The difference between sight.m_ObservedAltitude and ALMANAC_HO is 0.0030989931019878014, which exceeds EPSILON_MIN, where
1: sight.m_ObservedAltitude evaluates to 11.513098993101988,
1: ALMANAC_HO evaluates to 11.51, and
1: EPSILON_MIN evaluates to 0.0016666666666666668.
1: Ho differs from NA by 0.18593958611926809 minutes
1:
1: Plugin data directory requested for celestial_navigation_pi, returning: /Users/serosset/git/celestial_navigation_pi/test/testdata
1:
1: Body Position Analysis:
1: Almanac GHA: 128° 20.5', Dec: 18° 15.2' S
1: Actual GHA : -128° -20.5', Dec: 18° 15.2' S
1: /Users/serosset/git/celestial_navigation_pi/test/altitude_tests.cpp:121: Failure
1: The difference between gha and ALMANAC_GHA is 256.68261321350474, which exceeds EPSILON_MIN, where
1: gha evaluates to -128.3409465468381,
1: ALMANAC_GHA evaluates to 128.34166666666667, and
1: EPSILON_MIN evaluates to 0.0016666666666666668.
1: GHA differs from Almanac by -15400.956792810284 minutes
1:
1:
1: Detailed Calculation String:
1: Almanac Data For Sun
1: Geographical Position (lat, lon) = -18.2537 -128.3409
1: GHAAST = 357 52.4'
1: SHA = 130 28.1'
1: GHA = 128 20.5'
1: Dec = S 18 15.2'
1: SD = 16.2'
1: HP = 0.1'
1:
1: Formulas used to calculate sight
1:
1: Index Error is 0.0533 degrees
1:
1: Eye Height is 2.4000 meters
1: Height Correction Degrees = 1.758*sqrt(2.4000) / 60.0
1: Height Correction Degrees = 0.0454
1:
1: Apparent Altitude (Ha)
1: ApparentAltitude = Measurement - IndexCorrection - EyeHeightCorrection
1: ApparentAltitude = 11.4167 - 0.0533 - 0.0454
1: ApparentAltitude = 11.3179
1:
1: Refraction Correction
1: x = tan(Pi/180*ApparentAltitude + 4.848e-2*(Pi/180) / (tan(Pi/180*ApparentAltitude) + .028))
1: x = tan(Pi/180*11.3179 + 4.848e-2*(Pi/180) / (tan(Pi/180*11.3179) + .028))
1: x = 0.2040
1: RefractionCorrection = .267 * Pressure / (x*(Temperature + 273.15)) / 60.0
1: RefractionCorrection = .267 * 1013.0000 / (x*(15.0000 + 273.15)) / 60.0
1: RefractionCorrection = 0.0767
1:
1: Sun selected, Limb Correction
1: ra = 0.9894, lc = 0.266564/ra = 0.2694
1: Lower Limb
1: LimbCorrection = -0.2694
1:
1: Corrected Altitude
1: CorrectedAltitude = ApparentAltitude - RefractionCorrection - LimbCorrection
1: CorrectedAltitude = 11.3179 - 0.0767 - -0.2694
1: CorrectedAltitude = 11.5107
1:
1: Sun selected, parallax correction
1: rad = 0.9894, HP = 0.002442/rad = 0.0025
1: ParallaxCorrection = -180/Pi * asin( sin(Pi/180 * HP ) * cos(Pi/180 * CorrectedAltitude))
1: ParallaxCorrection = -180/Pi * asin( sin(Pi/180 * 0.0025 ) * cos(Pi/180 * 11.5107))
1: ParallaxCorrection = -0.0024
1:
1: Observed Altitude (Ho)
1: ObservedAltitude = CorrectedAltitude - ParallaxCorrection
1: ObservedAltitude = 11.5107 - -0.0024
1: ObservedAltitude = 11.5131
1:
1: [ FAILED ] AltitudeTest.SunLowerLimbExample (145 ms)
1: [----------] 1 test from AltitudeTest (145 ms total)
1:
1: [----------] Global test environment tear-down
1: [==========] 1 test from 1 test suite ran. (145 ms total)
1: [ PASSED ] 0 tests.
1: [ FAILED ] 1 test, listed below:
1: [ FAILED ] AltitudeTest.SunLowerLimbExample
1:
1: 1 FAILED TEST
1/1 Test #1: celestial_tests ..................***Failed 0.20 sec |
|
This is great. We will be able to quickly test the plugin as we make improvements. The plugin appears to be failing the standard set right now. Incidentally, at the bottom of the Manual, I did a series of tests and left this note after hitting "Fix"
|
|
Hello Sebastian. My name is Bob Bossert. I’m a big fan of open cpn and I
teach celestial navigation and taught my students how to use it. I’m the
one who spotted the problem with moon position to be because open cpn plug
in uses avg distance of moon distance to earth. And avg distance is not
good enough to calculate a moon position. And moon with Sun makes a good
2 body fix
I wrote up the issue and included examples in a PDF. The issue is only
with moon sights.
Semi diameter and parallax is mostly a function of distance and earth/moon
diameter. I wrote a C++ program to calculate moon distance based on time.
I uploaded the C++ so you can take a look at it. I tried to write the
program in such a way that when integrated into the plug in, nothing else
in the plug in would be disturbed.
The moon-distance function uses TT as input. TT is based on Julien date
but converted to centuries and y2000 epoch. Although y2000 epoch is part
of the definition, the algorithm works with negative and positive TT. I
believe the plug in already calculates this time nomenclature. The
function’s output is distance in km.
For integration into the existing plugin, You can call the function just
before you calculate moon semi diameter and moon parallax which is then
used to calculate Ho. In my code, SD and PARALLAX are in arc-minutes which
can be converted to degrees since I think the plug in carries all angles as
decimal degrees.
I also sent an xls spreadsheet for test cases from 1992, 200x, and 2024,
2025. It was uploaded as a zip file
Let me know if you haven’t seen them. Please call me if you have any
questions or doubts about what I did in the code. Rick Gleason has my
phone number
I’m the author of the other issues that Rick Gleason mentioned.
Bob
…On Tue, Feb 11, 2025 at 6:34 PM Sebastien Rosset ***@***.***> wrote:
I was able to reproduce (mostly) the first example in this test:
https://github.com/rgleason/celestial_navigation_pi/pull/33/files#diff-538ede6a0b5e101ab74abef50e6f6364e7015154d64c21821f8084e56fc490b5R49
In particular for the given input parameters, I can confirm I see:
Almanac OpenCPN
Ho 11° 30.6’ 11° 30.8’ (11.5131°)
GHA 128° 20.5', Dec: 18° 15.2' S -128° -20.5', Dec: 18° 15.2' S
For the GHA, I see the same results with a different sign.
—
Reply to this email directly, view it on GitHub
<#24 (comment)>,
or unsubscribe
<https://github.com/notifications/unsubscribe-auth/BOSWTTRF5CJMHC666VEVKCD2PKCIHAVCNFSM6AAAAABVL4P7PKVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDMNJSGI4TOMBTHA>
.
You are receiving this because you were mentioned.Message ID:
***@***.***>
|
|
@rbossert149419 I would like to help you, though realistically I don't think I'll be able to do that in the next several weeks. I'm currently making improvements to the VDR plugin, GRIB plugin and have plans to improve the weather routing plugin. |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Celestial Navigation Plugin shows Altitude of the Body on the Find celestial body form (as well as Azimuth). Altitude is known as Hc in the Nautical Almanac. The Altitude is wrong on the form. Attached 3 examples show a 3 to 6 minute error which resulted in 3-6 nm error in EP.
Steps to reproduce:
Drop a Mark to your DR location. "Move your Boat" to the Mark (as close as possible). Then enter a sight, select Sun, Planet, or a Star. Enter your sight data. Select the FIND (Find celestial body) button on the site tab. Altitude is incorrect and inconsistent to the Line of Position the plug in Plotted on the chart (Azimuth is correct). On the Celestial Navigation Sights form (shows all the sights), click on the View icon for that sight. The Line of Position plot which relies on Hc, Ho, Azimuth is shown. This plot is correct. The problem appears to be limited to the calculation of Hc using Law of Cosines formula on the FIND form. Attached are 3 examples that illustrate the Find Form Altitude value and the plugin's correct answer. In the examples, the error resulted in a 3-6 NM error in the EP or Fix. You can also confirm the error by Deducing Hc using the Open CPN measurement tool for the distance between the plotted LOP and DR mark (known as intercept). Intercept = Hc - Ho, thus Hc = Intercept + Ho.
Expected behavior
The Altitude should be within .1 minute of a degree compared to Commercial edition of Nautical Almanac using Sight Reduction formula pages 277-283 (attachment includes the Law of Cosines formula). Because of all of the combinations of sin, sin-1, cos, cos-1 going on in Law of Cosines equation, it has been recommended to use 5 decimal accuracy for LHA, DEC, and LAT to get .1 minute precision. Since inputs from Sextant (Hs), and Nautical Almanac daily page (GHA, Dec, SHA) are at 0.1 minute accuracy, Navigators don't expect Hc displayed beyond 0.1 minutes accuracy. Nautical Almanac examples shows Hc decimal degrees at 4 significant digits. I'd suggest 4 significant digits
Screen shots....
Open CPN Altitude calculation error PDF.pdf
OS: problem duplicated for macOS and Windows 10 (I've delayed my upgrade to 11). I believe all platforms have the problem..
plug in Version 2.4.41
Open CPN version 5.10.2
The text was updated successfully, but these errors were encountered: