Merge lp:~georg-zotti/stellarium/gz_fix-ecliptic-obliquity into lp:stellarium

Proposed by gzotti
Status: Merged
Approved by: Alexander Wolf
Approved revision: 7625
Merged at revision: 7748
Proposed branch: lp:~georg-zotti/stellarium/gz_fix-ecliptic-obliquity
Merge into: lp:stellarium
Diff against target: 2869 lines (+1419/-206)
28 files modified
README (+5/-1)
data/default_config.ini.cmake (+22/-6)
src/CMakeLists.txt (+19/-0)
src/core/StelCore.cpp (+82/-15)
src/core/StelCore.hpp (+61/-36)
src/core/StelObject.cpp (+39/-19)
src/core/StelObject.hpp (+12/-2)
src/core/StelObserver.cpp (+28/-2)
src/core/VecMath.hpp (+42/-2)
src/core/modules/GridLinesMgr.cpp (+273/-24)
src/core/modules/GridLinesMgr.hpp (+140/-17)
src/core/modules/Landscape.cpp (+6/-0)
src/core/modules/Planet.cpp (+62/-18)
src/core/modules/Planet.hpp (+9/-3)
src/core/modules/SolarSystem.cpp (+4/-0)
src/core/planetsephems/elp82b.h (+1/-1)
src/core/planetsephems/gust86.c (+1/-1)
src/core/planetsephems/gust86.h (+3/-3)
src/core/planetsephems/precession.c (+260/-0)
src/core/planetsephems/precession.h (+59/-0)
src/core/planetsephems/sidereal_time.c (+35/-20)
src/core/planetsephems/vsop87.h (+1/-1)
src/gui/ConfigurationDialog.cpp (+4/-1)
src/gui/LocationDialog.cpp (+1/-1)
src/gui/ViewDialog.cpp (+5/-1)
src/gui/viewDialog.ui (+89/-32)
src/tests/testPrecession.cpp (+121/-0)
src/tests/testPrecession.hpp (+35/-0)
To merge this branch: bzr merge lp:~georg-zotti/stellarium/gz_fix-ecliptic-obliquity
Reviewer Review Type Date Requested Status
Alexander Wolf Approve
gzotti Needs Resubmitting
Review via email: mp+265992@code.launchpad.net

Description of the change

This brings an important change in accuracy, fixing several questions about accuracy of simulation in antiquity. Please see branch description for details.

There are a few issues to be fixed as discussed in the last few weeks. Precession (coming here), deltaT once again (I suggest doing this next) and Nutation (IAU-2000A should replace the currently used old model from the 1980s). Accuracy gains/testing against solutions known to be accurate) is best done in this sequence. I left a few remarks and TODO notes in the code in places that should soon be fixed (IMHO these are things that should be done before 0.14 is prepared).

To post a comment you must log in.
7615. By gzotti

Added reference to README. A few docfixes and removed some excess comments.

Revision history for this message
Alexander Wolf (alexwolf) wrote :

Georg, are you can add unit tests for this feature? It should be helpful for checking accuracy and checking for correct for the code.

review: Needs Information
Revision history for this message
Alexander Wolf (alexwolf) wrote :

Oh! One important note - pole precision line is not a circle, it's spirale. Current implementation can confuse many peoples. I think this line may be implemented as line with some range of the dates, like orbits in Satellites plugin.

review: Needs Information
Revision history for this message
Alexander Wolf (alexwolf) wrote :

Revision 7610 contains the code:

oss << q_("Ecliptic obliquity (of date): %1").arg(StelUtils::radToDmsStr(ecl, true)) << "<br>";

which will be displayed when any object is selected.

But (!):
1) this line give info about current planet, where observer is located, not about the selected object!
2) Ecliptic obliquity for Sun is not have physical sense.

Revision history for this message
Alexander Wolf (alexwolf) wrote :

Revision 7612: I think short keys for add/subtract Julian century should be added. :)

Revision history for this message
Alexander Wolf (alexwolf) wrote :

Revision 7613: same problem as in revision 7610 - the planetocentric distance shows info not for the selected planet, what will be confused users.

Revision history for this message
Alexander Wolf (alexwolf) wrote :

I'm sorry, but notes from my last 3 comments should be fixed IMHO.

review: Needs Fixing
Revision history for this message
gzotti (georg-zotti) wrote :

> Revision 7610 contains the code:
>
> oss << q_("Ecliptic obliquity (of date): %1").arg(StelUtils::radToDmsStr(ecl,
> true)) << "<br>";
>
> which will be displayed when any object is selected.
>
> But (!):
> 1) this line give info about current planet, where observer is located, not
> about the selected object!
> 2) Ecliptic obliquity for Sun is not have physical sense.

"Ecliptic" is only Earth's orbital plane, there is no source or reason for confusion. But given the fact that we also don't show ecliptical grid of date if located on other planets, I will change to not show it if location is not on Earth.

Revision history for this message
gzotti (georg-zotti) wrote :

> Revision 7612: I think short keys for add/subtract Julian century should be
> added. :)

Yes, only in comparing the merge submission I have seen that Sidereal Year has vanished (?)
I would even like to have Sidereal, Tropical, Julian, and Anomalistic year! But I am easily running out of hotkeys! Whatever time jumps we decide upon for 0.14 I think we can add them just before release.

Revision history for this message
Alexander Wolf (alexwolf) wrote :

>"Ecliptic" is only Earth's orbital plane, there is no source or reason for confusion. But given the fact that we also don't show ecliptical grid of date if located on other planets, I will change to not show it if location is not on Earth.

Why? When we say "ecliptic" then we saying about Earth, it's typical, but other planets has "path of the Sun" too and other planet has ecliptic too ;)

Revision history for this message
Alexander Wolf (alexwolf) wrote :

> Yes, only in comparing the merge submission I have seen that Sidereal Year has vanished (?)

No, it exists. (Hint: http://bazaar.launchpad.net/~stellarium/stellarium/trunk/revision/7734)

> I would even like to have Sidereal, Tropical, Julian, and Anomalistic year! But I am easily running out of hotkeys! Whatever time jumps we decide upon for 0.14 I think we can add them just before release.

What about filling feature request, because we can forget it?

Revision history for this message
Alexander Wolf (alexwolf) wrote :

Oops... Why "Ecliptic obliquity (of date)" string are ignores the option of settings for decimal format displaying? :)

Revision history for this message
Alexander Wolf (alexwolf) wrote :

Georg, please ignore my comment about precession :)

Revision history for this message
gzotti (georg-zotti) wrote :

> Revision 7613: same problem as in revision 7610 - the planetocentric distance
> shows info not for the selected planet, what will be confused users.

Oh yes, it does. It should display how many km observer is located away from his planet's center.

But ... argh - something is wrong for off-earth sites. Must fix that, sure. (Last-minute thing I wanted to put in here, only tested for Earth. *-(

Revision history for this message
gzotti (georg-zotti) wrote :

> >"Ecliptic" is only Earth's orbital plane, there is no source or reason for
> confusion. But given the fact that we also don't show ecliptical grid of date
> if located on other planets, I will change to not show it if location is not
> on Earth.
>
> Why? When we say "ecliptic" then we saying about Earth, it's typical, but
> other planets has "path of the Sun" too and other planet has ecliptic too ;)

No, "The Ecliptic" is specifically the projection of Earth's instantaneous orbital plane onto the celestial sphere, appearing as great circle.

Other planets' orbital planes are just called "orbital planes", and Stellarium so far never made an attempt to show other planets' orbital planes (although it would certainly be nice for astronaut training...). "Ecliptic" is particularly Earth-related, like "The Moon" is only our Earth's moon.

There is apparently some widespread idea that Earth's orbital plane is stable in VSOP coordinates. It is not the case, as we can see here and now! VSOP has the J2000.0 ecliptical plane as their XY-plane, which Stellarium can show as "Ecliptic of J2000.0". But Earth's orbit is wandering around a bit, which is called the Precession of the Ecliptic (Earlier name: "planetary precession", but in the post-2000 literature authors preferred "Precession of the Ecliptic").

Revision history for this message
gzotti (georg-zotti) wrote :

> Oops... Why "Ecliptic obliquity (of date)" string are ignores the option of
> settings for decimal format displaying? :)

fixed in r7616, thanks.

Revision history for this message
gzotti (georg-zotti) wrote :

> Oh! One important note - pole precision line is not a circle, it's spirale.
> Current implementation can confuse many peoples. I think this line may be
> implemented as line with some range of the dates, like orbits in Satellites
> plugin.

Yes it's a spiral. But what I show here is the "instantaneous circle" which at least visualizes the concept. Avid observers will note that stars orbiting the instantaneous pole of ecliptic of date next to this circle will intersect the circle, and that Polaris will have a different polar distance next time it is pole star. I think this is perfectly enough as visualisation for now.

When the feature request "orbit with labeled date marks" has been implemented, maybe someone can also implement the spiral with date labels then?

In these time ranges (tens of thousands of years) it's sad to see that VSOP87 only covers less than 1/2 cycle... But I have never seen the answer to "what happens outside -4000...+8000" as clearly.

DE431 should show a bit more than 1 cycle, at least.

Revision history for this message
gzotti (georg-zotti) wrote :

> Georg, please ignore my comment about precession :)

I hope you mean ignoring the spiral request? I have no plan to implement any spiral path before December (if at all), there are really more urgent things (e.g. moon rotation!) and I do also not have infinite time, sorry.

Revision history for this message
gzotti (georg-zotti) wrote :

> > Revision 7613: same problem as in revision 7610 - the planetocentric
> distance
> > shows info not for the selected planet, what will be confused users.
>
> Oh yes, it does. It should display how many km observer is located away from
> his planet's center.
>
> But ... argh - something is wrong for off-earth sites. Must fix that, sure.
> (Last-minute thing I wanted to put in here, only tested for Earth. *-(

OK, fixed that in r7617.

This value may actually only be useful while testing. Again, it does not give distance to observed object, but to the center of the planet you are standing on. (Needed for topocentric correction, which used purely-spherical planets until now.)

7616. By gzotti

show ecliptic obliquity also with non-decimal string

7617. By gzotti

ensure getDistanceFromCenter() does not cause problems while observing from SpaceShipObserver

Revision history for this message
Alexander Wolf (alexwolf) wrote :

OK, I think ecliptic obliquity and XYZ should display for Earth only and she should not mixed with other data for the chosen planets. What think about adding "For Earth:" line after info (after the one empty line) and put ecliptic obliquity and XYZ data after it.

The question for info with the planetocentric distance is more difficult here, because this infobox should not mixed with data about chosen planet too. But this data available for all planets in general. I think I can propose compromise for visualization of this data - add the planetocentric distance info as tooltip for height data on bottom toolbar.

Revision history for this message
Alexander Wolf (alexwolf) wrote :

I guess we can open question about implementation rendering line of precession after introducing DE4x parts.

7618. By gzotti

added unit tests (themselves untested, TBD on other platform)

7619. By gzotti

added test case files

7620. By gzotti

commented away XYZ and planetocentric observer coordinates

7621. By gzotti

harmonized visible (grid) and displayed ecliptical coordinates for non-earth locations. We cannot compute orbit-related coordinates for other planets yet :-(

Revision history for this message
gzotti (georg-zotti) wrote :

> OK, I think ecliptic obliquity and XYZ should display for Earth only and she
> should not mixed with other data for the chosen planets. What think about
> adding "For Earth:" line after info (after the one empty line) and put
> ecliptic obliquity and XYZ data after it.

Given that nobody outside our development team thinks that "Ecliptic" does not refer to the Earth only, I see no need to annotate "Ecliptic" with "Earth".

By the way, the "ecliptical coordinates" in the infostring are totally nonsense for observer locations outside earth. You would have to first find orbital plane for the planet, and then the "spring point" (ascending node between orbital plane and planet equator). The situation found here were in principle terrestrial ecliptical coordinates with different obliquity, that's nonsense. I now have changed this, for non-terrestrial observers, to

*) show "ecliptical coordinates J2000" referring to earth's orbital plane of J2000.
   The displayed coordinates fit to the displayed ecliptical grid of J2000.
*) Disable display of ecliptical coordinates of date because they were wrong.
*) Disable display of "Ecliptic obliquity" for non-earth observers to avoid confusion.

If we want to display "axis inclination to planet's orbit" or (90-Planet.getRotObliquity(JDE)), this may be a nice extra information.

Also XYZ is mostly here for testing, but it should show selected object's position in "parentocentric" XYZ coordinates (i.e., for moons we have the respective planets as zero). I have commented it away for your comfort. Please do not delete the lines, though, I think it can be helpful when integrating DE431 which have to be rotated into VSOP87 reference frame.

> The question for info with the planetocentric distance is more difficult here,
> because this infobox should not mixed with data about chosen planet too. But
> this data available for all planets in general. I think I can propose
> compromise for visualization of this data - add the planetocentric distance
> info as tooltip for height data on bottom toolbar.

Also here I just wanted a quick view while I move latitude with the location panel open, and not just qDebug() lines in the logfile. The current location was admittedly not good, but this is an intermediate state, not polished. If you can make a tooltip from that, it's fine.

r7621 has the latest changes. I now move to Linux for testing the unit test...

7622. By gzotti

avoid debug print also in "regular" debug builds.

7623. By gzotti

added toString() functions for matrices. Clarified (unused) upper3x3().

7624. By gzotti

First unit test for precession complete.

7625. By gzotti

added another test

Revision history for this message
gzotti (georg-zotti) wrote :

> Georg, are you can add unit tests for this feature? It should be helpful for
> checking accuracy and checking for correct for the code.

Alright, a few unit tests against the numerical example given in the paper are completed (r7625). There was even a minus missing in the (otherwise unused) reference solution :-P
Not sure about how the unit tests are used in an automatic workflow, are debug messages allowed? I leave them in but commented away, for further development.

Alex, I hope this completes the package of your requests. This branch does not care for other planets' orbital planes which were never properly handled by now, only applications in terrestrial history will see a benefit.

Again, this is Step 1 of 3 to be done in the next few weeks to solve several critical questions by historically inclined audience. #2 shall be correct application of DeltaT, and #3 shall bring in Nutation with the IAU-2000A model fitting to the 2006 precession of which this Vondrak model is a compatible, extended version.

Anything missing?

review: Needs Resubmitting
Revision history for this message
Alexander Wolf (alexwolf) wrote :

>Not sure about how the unit tests are used in an automatic workflow, are debug messages allowed?

No.

review: Approve

Preview Diff

[H/L] Next/Prev Comment, [J/K] Next/Prev File, [N/P] Next/Prev Hunk
1=== modified file 'README'
2--- README 2015-03-23 00:02:31 +0000
3+++ README 2015-07-28 19:18:20 +0000
4@@ -207,6 +207,10 @@
5 [19] "Astronomy on the Personal Computer" by O. Montenbruck & T. Pfleger (4nd ed., 2000)
6 [20] "Calendrical Calculations" by E. M. Reingold & N. Dershowitz (2nd ed., 2001)
7 [21] DeltaT webpage by V. Reijs: http://www.iol.ie/~geniet/eng/DeltaTeval.htm
8+ 1.5 An accurate long-time precession model compatible with P03:
9+ J. Vondrak, N. Capitaine, P. Wallace: New precession expressions, valid for long time intervals.
10+ Astronomy&Astrophysics 534, A22 (2011); DOI: 10.1051/0004-6361/201117274
11+
12
13 2. Included source code
14 2.1 Some computation of the sideral time (sidereal_time.h/c) and pluto
15@@ -252,7 +256,7 @@
16 -) compute nr_of_measurements = the number of valid b,v,r values
17 -) throw away or keep stars (depending on magnitude,
18 nr_of_measurements, combination of flags, tycho2 number)
19- -) add all stars from hiparcos (incl. component solutions), and
20+ -) add all stars from Hipparcos (incl. component solutions), and
21 tycho2+1st supplement
22 -) reorganize the stars in several brigthness levels and
23 triangular zones according to position and magnitude
24
25=== modified file 'data/default_config.ini.cmake'
26--- data/default_config.ini.cmake 2015-06-19 18:02:55 +0000
27+++ data/default_config.ini.cmake 2015-07-28 19:18:20 +0000
28@@ -113,7 +113,6 @@
29 #base_font_name = DejaVu Sans
30 mouse_cursor_timeout = 10
31 flag_mouse_cursor_timeout = false
32-day_key_mode = calendar
33 selected_object_info = all
34 auto_hide_horizontal_toolbar = true
35 auto_hide_vertical_toolbar = true
36@@ -124,13 +123,27 @@
37 gui_base_color = 0.5,0.5,0.7
38 gui_text_color = 0.8,0.9,0.9
39 azimuthal_color = 0.3,0.2,0.1
40-equatorial_color = 0.1,0.2,0.3
41-equatorial_J2000_color = 0.1,0.3,0.4
42+daylight_text_color = 0.0,0.0,0.0
43+#
44+#equatorial_color = 0.1,0.2,0.3
45+#equatorial_J2000_color = 0.1,0.3,0.4
46+#equator_color = 0.2,0.2,0.6
47+#ecliptic_color = 0.6,0.2,0.2
48+#ecliptic_J2000_color = 0.3,0.3,0.1
49+# GZ instead: Note -al colors go to the grids.
50+ecliptic_J2000_color = 0.7,0.2,0.2
51+ecliptic_color = 0.9,0.6,0.2
52+ecliptical_J2000_color = 0.4,0.1,0.1
53+ecliptical_color = 0.6,0.3,0.1
54+# default: same as ecliptic_color
55+precession_circle_color = 0.9,0.6,0.2
56+equator_J2000_color = 0.2,0.2,0.6
57+equator_color = 0.3,0.5,1.0
58+equatorial_J2000_color = 0.1,0.1,0.5
59+equatorial_color = 0.2,0.3,0.8
60+#
61 galactic_color = 0.0,0.3,0.2
62 galactic_equator_color = 0.5,0.3,0.1
63-equator_color = 0.2,0.2,0.6
64-ecliptic_color = 0.6,0.2,0.2
65-ecliptic_J2000_color = 0.3,0.3,0.1
66 meridian_color = 0.2,0.6,0.2
67 longitude_color = 0.2,0.4,0.4
68 horizon_color = 0.2,0.6,0.2
69@@ -207,11 +220,14 @@
70 flag_azimuthal_grid = false
71 flag_equatorial_grid = false
72 flag_equatorial_J2000_grid = false
73+flag_ecliptic_grid = false
74 flag_ecliptic_J2000_grid = false
75 flag_galactic_grid = false
76 flag_galactic_equator_line = false
77 flag_equator_line = false
78+flag_equator_J2000_line = false
79 flag_ecliptic_line = false
80+flag_ecliptic__J2000_line = false
81 flag_meridian_line = false
82 flag_longitude_line = false
83 flag_horizon_line = false
84
85=== modified file 'src/CMakeLists.txt'
86--- src/CMakeLists.txt 2015-06-13 20:19:16 +0000
87+++ src/CMakeLists.txt 2015-07-28 19:18:20 +0000
88@@ -174,6 +174,8 @@
89 core/planetsephems/marssat.h
90 core/planetsephems/pluto.c
91 core/planetsephems/pluto.h
92+ core/planetsephems/precession.c
93+ core/planetsephems/precession.h
94 core/planetsephems/sidereal_time.c
95 core/planetsephems/sidereal_time.h
96 core/planetsephems/stellplanet.c
97@@ -620,6 +622,23 @@
98 ADD_DEPENDENCIES(buildTests testRefraction)
99 ADD_TEST(testRefraction)
100
101+SET(tests_testPrecession_SRCS
102+ tests/testPrecession.hpp
103+ tests/testPrecession.cpp
104+ core/planetsephems/precession.h
105+ core/planetsephems/precession.c
106+ core/StelUtils.hpp
107+ core/StelUtils.cpp)
108+IF(WIN32)
109+ # StelUtils required zlib sources
110+ SET(tests_testPrecession_SRCS ${tests_testPrecession_SRCS} ${zlib_SRCS})
111+ENDIF()
112+ADD_EXECUTABLE(testPrecession EXCLUDE_FROM_ALL ${tests_testPrecession_SRCS})
113+QT5_USE_MODULES(testPrecession Core Gui Widgets Script Declarative Test)
114+TARGET_LINK_LIBRARIES(testPrecession ${extLinkerOptionTest})
115+ADD_DEPENDENCIES(buildTests testPrecession)
116+ADD_TEST(testPrecession)
117+
118 ADD_CUSTOM_TARGET(tests COMMENT "Run the Stellarium unit tests")
119 FOREACH(NAME ${STELLARIUM_TESTS})
120 IF(MSVC)
121
122=== modified file 'src/core/StelCore.cpp'
123--- src/core/StelCore.cpp 2015-07-25 17:05:44 +0000
124+++ src/core/StelCore.cpp 2015-07-28 19:18:20 +0000
125@@ -37,6 +37,7 @@
126 #include "LandscapeMgr.hpp"
127 #include "StelTranslator.hpp"
128 #include "StelActionMgr.hpp"
129+#include "precession.h"
130
131 #include <QSettings>
132 #include <QDebug>
133@@ -311,10 +312,12 @@
134 {
135 case FrameAltAz:
136 return getProjection(getAltAzModelViewTransform(refractionMode));
137- case FrameHeliocentricEcliptic:
138+ case FrameHeliocentricEclipticJ2000:
139 return getProjection(getHeliocentricEclipticModelViewTransform(refractionMode));
140- case FrameObservercentricEcliptic:
141- return getProjection(getObservercentricEclipticModelViewTransform(refractionMode));
142+ case FrameObservercentricEclipticJ2000:
143+ return getProjection(getObservercentricEclipticJ2000ModelViewTransform(refractionMode));
144+ case FrameObservercentricEclipticOfDate:
145+ return getProjection(getObservercentricEclipticOfDateModelViewTransform(refractionMode));
146 case FrameEquinoxEqu:
147 return getProjection(getEquinoxEquModelViewTransform(refractionMode));
148 case FrameJ2000:
149@@ -635,9 +638,9 @@
150 Vec3d StelCore::heliocentricEclipticToAltAz(const Vec3d& v, RefractionMode refMode) const
151 {
152 if (refMode==RefractionOff || skyDrawer==NULL || (refMode==RefractionAuto && skyDrawer->getFlagHasAtmosphere()==false))
153- return matHeliocentricEclipticToAltAz*v;
154+ return matHeliocentricEclipticJ2000ToAltAz*v;
155 Vec3d r(v);
156- r.transfo4d(matHeliocentricEclipticToAltAz);
157+ r.transfo4d(matHeliocentricEclipticJ2000ToAltAz);
158 skyDrawer->getRefraction().forward(r);
159 return r;
160 }
161@@ -648,26 +651,29 @@
162 return matHeliocentricEclipticToEquinoxEqu*v;
163 }
164
165+/*
166 //! Transform vector from heliocentric coordinate to false equatorial : equatorial
167-//! coordinate but centered on the observer position (usefull for objects close to earth)
168+//! coordinate but centered on the observer position (useful for objects close to earth)
169+//! Unused as of V0.13
170 Vec3d StelCore::heliocentricEclipticToEarthPosEquinoxEqu(const Vec3d& v) const
171 {
172 return matAltAzToEquinoxEqu*matHeliocentricEclipticToAltAz*v;
173 }
174+*/
175
176 StelProjector::ModelViewTranformP StelCore::getHeliocentricEclipticModelViewTransform(RefractionMode refMode) const
177 {
178 if (refMode==RefractionOff || skyDrawer==NULL || (refMode==RefractionAuto && skyDrawer->getFlagHasAtmosphere()==false))
179- return StelProjector::ModelViewTranformP(new StelProjector::Mat4dTransform(matAltAzModelView*matHeliocentricEclipticToAltAz));
180+ return StelProjector::ModelViewTranformP(new StelProjector::Mat4dTransform(matAltAzModelView*matHeliocentricEclipticJ2000ToAltAz));
181 Refraction* refr = new Refraction(skyDrawer->getRefraction());
182 // The pretransform matrix will convert from input coordinates to AltAz needed by the refraction function.
183- refr->setPreTransfoMat(matHeliocentricEclipticToAltAz);
184+ refr->setPreTransfoMat(matHeliocentricEclipticJ2000ToAltAz);
185 refr->setPostTransfoMat(matAltAzModelView);
186 return StelProjector::ModelViewTranformP(refr);
187 }
188
189-//! Get the modelview matrix for observer-centric ecliptic (Vsop87) drawing
190-StelProjector::ModelViewTranformP StelCore::getObservercentricEclipticModelViewTransform(RefractionMode refMode) const
191+//! Get the modelview matrix for observer-centric ecliptic J2000 (Vsop87A) drawing
192+StelProjector::ModelViewTranformP StelCore::getObservercentricEclipticJ2000ModelViewTransform(RefractionMode refMode) const
193 {
194 if (refMode==RefractionOff || skyDrawer==NULL || (refMode==RefractionAuto && skyDrawer->getFlagHasAtmosphere()==false))
195 return StelProjector::ModelViewTranformP(new StelProjector::Mat4dTransform(matAltAzModelView*matJ2000ToAltAz*matVsop87ToJ2000));
196@@ -678,6 +684,19 @@
197 return StelProjector::ModelViewTranformP(refr);
198 }
199
200+//! Get the modelview matrix for observer-centric ecliptic-of-date drawing
201+StelProjector::ModelViewTranformP StelCore::getObservercentricEclipticOfDateModelViewTransform(RefractionMode refMode) const
202+{
203+ double eps_A=getPrecessionAngleVondrakCurrentEpsilonA();
204+ if (refMode==RefractionOff || skyDrawer==NULL || (refMode==RefractionAuto && skyDrawer->getFlagHasAtmosphere()==false))
205+ return StelProjector::ModelViewTranformP(new StelProjector::Mat4dTransform(matAltAzModelView*matEquinoxEquToAltAz* Mat4d::xrotation(eps_A)));
206+ Refraction* refr = new Refraction(skyDrawer->getRefraction());
207+ // The pretransform matrix will convert from input coordinates to AltAz needed by the refraction function.
208+ refr->setPreTransfoMat(matEquinoxEquToAltAz* Mat4d::xrotation(eps_A));
209+ refr->setPostTransfoMat(matAltAzModelView);
210+ return StelProjector::ModelViewTranformP(refr);
211+}
212+
213 //! Get the modelview matrix for observer-centric equatorial at equinox drawing
214 StelProjector::ModelViewTranformP StelCore::getEquinoxEquModelViewTransform(RefractionMode refMode) const
215 {
216@@ -729,31 +748,52 @@
217 return StelProjector::ModelViewTranformP(refr);
218 }
219
220+// GZ: One of the most important functions, totally void of doc. :-(
221+// called in update() (for every frame)
222 void StelCore::updateTransformMatrices()
223 {
224 matAltAzToEquinoxEqu = position->getRotAltAzToEquatorial(JDay);
225 matEquinoxEquToAltAz = matAltAzToEquinoxEqu.transpose();
226
227+ // multiply static J2000 earth axis tilt (eclipticalJ2000<->equatorialJ2000)
228+ // in effect, this matrix transforms from VSOP87 ecliptical J2000 to planet-based equatorial coordinates.
229+ // For earth, matJ2000ToEquinoxEqu is the precession matrix.
230+ // TODO: rename matEquinoxEquToJ2000 to matEquinoxEquDateToJ2000
231 matEquinoxEquToJ2000 = matVsop87ToJ2000 * position->getRotEquatorialToVsop87();
232 matJ2000ToEquinoxEqu = matEquinoxEquToJ2000.transpose();
233 matJ2000ToAltAz = matEquinoxEquToAltAz*matJ2000ToEquinoxEqu;
234
235 matHeliocentricEclipticToEquinoxEqu = matJ2000ToEquinoxEqu * matVsop87ToJ2000 * Mat4d::translation(-position->getCenterVsop87Pos());
236
237- // These two next have to take into account the position of the observer on the earth
238+ // These two next have to take into account the position of the observer on the earth/planet of observation.
239+ // GZ tmp could be called matAltAzToVsop87
240 Mat4d tmp = matJ2000ToVsop87 * matEquinoxEquToJ2000 * matAltAzToEquinoxEqu;
241
242- matAltAzToHeliocentricEcliptic = Mat4d::translation(position->getCenterVsop87Pos()) * tmp *
243+ // TODO: getDistanceFromCenter assumes spherical planets. Use rectangular coordinates for observer!
244+ matAltAzToHeliocentricEclipticJ2000 = Mat4d::translation(position->getCenterVsop87Pos()) * tmp *
245 Mat4d::translation(Vec3d(0.,0., position->getDistanceFromCenter()));
246
247- matHeliocentricEclipticToAltAz = Mat4d::translation(Vec3d(0.,0.,-position->getDistanceFromCenter())) * tmp.transpose() *
248+ matHeliocentricEclipticJ2000ToAltAz = Mat4d::translation(Vec3d(0.,0.,-position->getDistanceFromCenter())) * tmp.transpose() *
249 Mat4d::translation(-position->getCenterVsop87Pos());
250+
251+
252+// Given that ecliptic of date is apparently not used any longer as reference plane for astronomical computations, it is enough to display it.
253+// This half-finished block of comments can the likely go away for V0.14. Else, do we have use for carrying the extra matrices with us? Then complete this:
254+// // GZ Finally, new: rotation of the ecliptic of date. Not sure which direction, though...
255+// double epsilon_A, chi_A, omega_A, psi_A;
256+// // MAKE SURE JDay IS TT!!!
257+// getPrecessionAnglesVondrak(JDay, &epsilon_A, &chi_A, &omega_A, &psi_A);
258+// //matVsop87ToEclOfDate= Mat4d::xrotation(epsilon_A)*Mat4d::zrotation(chi_A)*Mat4d::xrotation(-omega_A)*Mat4d::zrotation(-psi_A);
259+// matVsop87ToEclOfDate= Mat4d::xrotation(epsilon_A)*Mat4d::zrotation(chi_A)*Mat4d::xrotation(-23.4392803055555555556*(M_PI/180));
260+// matEclOfDateToVsop87=matVsop87ToEclOfDate.transpose();
261+// //matEclOfDateToVsop87= Mat4d::xrotation(epsilon_A)*Mat4d::zrotation(chi_A)*Mat4d::xrotation(-omega_A)*Mat4d::zrotation(-psi_A);
262+// //matVsop87ToEclOfDate=matVsop87ToEclOfDate.transpose();
263 }
264
265 // Return the observer heliocentric position
266 Vec3d StelCore::getObserverHeliocentricEclipticPos() const
267 {
268- return Vec3d(matAltAzToHeliocentricEcliptic[12], matAltAzToHeliocentricEcliptic[13], matAltAzToHeliocentricEcliptic[14]);
269+ return Vec3d(matAltAzToHeliocentricEclipticJ2000[12], matAltAzToHeliocentricEclipticJ2000[13], matAltAzToHeliocentricEclipticJ2000[14]);
270 }
271
272 // Set the location to use by default at startup
273@@ -905,6 +945,11 @@
274 return position->getHomePlanet();
275 }
276
277+const StelObserver *StelCore::getCurrentObserver() const
278+{
279+ return position;
280+}
281+
282 // Smoothly move the observer to the given location
283 void StelCore::moveObserverTo(const StelLocation& target, double duration, double durationIfPlanetChange)
284 {
285@@ -1045,6 +1090,16 @@
286 addSolarDays(365.242190419);
287 }
288
289+void StelCore::addTropicalCentury()
290+{
291+ addSolarDays(36524.21897);
292+}
293+
294+void StelCore::addJulianCentury()
295+{
296+ addSolarDays(36525.0);
297+}
298+
299 void StelCore::subtractHour()
300 {
301 addSolarDays(-JD_HOUR);
302@@ -1105,6 +1160,16 @@
303 addSolarDays(-365.242190419);
304 }
305
306+void StelCore::subtractTropicalCentury()
307+{
308+ addSolarDays(-36524.21897);
309+}
310+
311+void StelCore::subtractJulianCentury()
312+{
313+ addSolarDays(-36525.0);
314+}
315+
316 void StelCore::addSolarDays(double d)
317 {
318 const PlanetP& home = position->getHomePlanet();
319@@ -1717,10 +1782,12 @@
320 return QString(" %1").arg(validRange);
321 }
322
323+// return if sky plus atmosphere is bright enough from sunlight so that e.g. screen labels should be rendered dark.
324+// Could be renamed to isBrightDaylight()
325 bool StelCore::isDay() const
326 {
327 const Vec3d& sunPos = GETSTELMODULE(SolarSystem)->getSun()->getAltAzPosGeometric(this);
328- return sunPos[2] > -0.12; // Nautical twilight
329+ return sunPos[2] > -0.12; // Nautical twilight (0.12 > sin (6 deg),
330 }
331
332 double StelCore::getCurrentEpoch() const
333
334=== modified file 'src/core/StelCore.hpp'
335--- src/core/StelCore.hpp 2015-07-08 10:37:11 +0000
336+++ src/core/StelCore.hpp 2015-07-28 19:18:20 +0000
337@@ -1,6 +1,7 @@
338 /*
339 * Copyright (C) 2003 Fabien Chereau
340 * Copyright (C) 2012 Matthew Gates
341+ * Copyright (C) 2015 Georg Zotti (Precession fixes)
342 *
343 * This program is free software; you can redistribute it and/or
344 * modify it under the terms of the GNU General Public License
345@@ -38,11 +39,11 @@
346 //! Main class for Stellarium core processing.
347 //! This class provides services like management of sky projections,
348 //! tone conversion, or reference frame conversion. It is used by the
349-//! various StelModules to update and display themself.
350+//! various StelModules to update and display themselves.
351 //! There is currently only one StelCore instance in Stellarium, but
352 //! in the future they may be more, allowing for example to display
353 //! several independent views of the sky at the same time.
354-//! @author Fabien Chereau
355+//! @author Fabien Chereau, Matthew Gates, Georg Zotti
356 class StelCore : public QObject
357 {
358 Q_OBJECT
359@@ -56,15 +57,17 @@
360 //! Supported reference frame types
361 enum FrameType
362 {
363- FrameUninitialized, //!< Reference frame is not set (FMajerech: Added to avoid condition on uninitialized value in StelSkyLayerMgr::draw())
364- FrameAltAz, //!< Altazimuthal reference frame centered on observer.
365- FrameHeliocentricEcliptic, //!< Ecliptic reference frame centered on the Sun
366- FrameObservercentricEcliptic, //!< Ecliptic reference frame centered on the Observer
367- FrameEquinoxEqu, //!< Equatorial reference frame at the current equinox centered on the observer.
368- //! The north pole follows the precession of the planet on which the observer is located.
369- FrameJ2000, //!< Equatorial reference frame at the J2000 equinox centered on the observer.
370- //! This is also the ICRS reference frame.
371- FrameGalactic //! Galactic reference frame centered on observer.
372+ FrameUninitialized, //!< Reference frame is not set (FMajerech: Added to avoid condition on uninitialized value in StelSkyLayerMgr::draw())
373+ FrameAltAz, //!< Altazimuthal reference frame centered on observer.
374+ FrameHeliocentricEclipticJ2000, //!< Fixed-ecliptic reference frame centered on the Sun. GZ: This is J2000 ecliptical / almost VSOP87.
375+ FrameObservercentricEclipticJ2000, //!< Fixed-ecliptic reference frame centered on the Observer. GZ: was ObservercentricEcliptic, but renamed because it is Ecliptic of J2000!
376+ FrameObservercentricEclipticOfDate, //!< Moving ecliptic reference frame centered on the Observer. GZ new for V0.14: Ecliptic of date, i.e. includes the precession of the ecliptic.
377+ FrameEquinoxEqu, //!< Equatorial reference frame at the current equinox centered on the observer.
378+ //!< The north pole follows the precession of the planet on which the observer is located.
379+ //!< TBD: To be Corrected for V0.14 to really properly reflect ecliptical motion and precession (Vondrak 2011 model)
380+ FrameJ2000, //!< Equatorial reference frame at the J2000 equinox centered on the observer.
381+ //!< This is also the ICRS reference frame.
382+ FrameGalactic //!< Galactic reference frame centered on observer.
383 };
384
385 //! @enum ProjectionType
386@@ -150,11 +153,11 @@
387 //! only for 2d painting
388 StelProjectorP getProjection2d() const;
389
390- //! Get a new instance of projector using a modelview transformation corresponding to the the given frame.
391+ //! Get a new instance of projector using a modelview transformation corresponding to the given frame.
392 //! If not specified the refraction effect is included if atmosphere is on.
393 StelProjectorP getProjection(FrameType frameType, RefractionMode refractionMode=RefractionAuto) const;
394
395- //! Get a new instance of projector using the given modelview transformatione.
396+ //! Get a new instance of projector using the given modelview transformation.
397 //! If not specified the projection used is the one currently used as default.
398 StelProjectorP getProjection(StelProjector::ModelViewTranformP modelViewTransform, ProjectionType projType=(ProjectionType)1000) const;
399
400@@ -210,31 +213,35 @@
401
402 //! Transform from heliocentric coordinate to equatorial at current equinox (for the planet where the observer stands)
403 Vec3d heliocentricEclipticToEquinoxEqu(const Vec3d& v) const;
404- //! Transform vector from heliocentric coordinate to false equatorial : equatorial
405- //! coordinate but centered on the observer position (usefull for objects close to earth)
406- Vec3d heliocentricEclipticToEarthPosEquinoxEqu(const Vec3d& v) const;
407+// //! Transform vector from heliocentric coordinate to false equatorial : equatorial
408+// //! coordinate but centered on the observer position (useful for objects close to earth)
409+// //! Unused as of V0.13
410+// Vec3d heliocentricEclipticToEarthPosEquinoxEqu(const Vec3d& v) const;
411
412- //! Get the modelview matrix for heliocentric ecliptic (Vsop87) drawing
413+ //! Get the modelview matrix for heliocentric ecliptic (Vsop87) drawing.
414 StelProjector::ModelViewTranformP getHeliocentricEclipticModelViewTransform(RefractionMode refMode=RefractionAuto) const;
415
416- //! Get the modelview matrix for observer-centric ecliptic (Vsop87) drawing
417- StelProjector::ModelViewTranformP getObservercentricEclipticModelViewTransform(RefractionMode refMode=RefractionAuto) const;
418-
419- //! Get the modelview matrix for observer-centric equatorial at equinox drawing
420+ //! Get the modelview matrix for observer-centric ecliptic (Vsop87) drawing.
421+ StelProjector::ModelViewTranformP getObservercentricEclipticJ2000ModelViewTransform(RefractionMode refMode=RefractionAuto) const;
422+
423+ //! Get the modelview matrix for observer-centric ecliptic of Date drawing.
424+ StelProjector::ModelViewTranformP getObservercentricEclipticOfDateModelViewTransform(RefractionMode refMode=RefractionAuto) const;
425+
426+ //! Get the modelview matrix for observer-centric equatorial at equinox drawing.
427 StelProjector::ModelViewTranformP getEquinoxEquModelViewTransform(RefractionMode refMode=RefractionAuto) const;
428
429- //! Get the modelview matrix for observer-centric altazimuthal drawing
430+ //! Get the modelview matrix for observer-centric altazimuthal drawing.
431 StelProjector::ModelViewTranformP getAltAzModelViewTransform(RefractionMode refMode=RefractionAuto) const;
432
433- //! Get the modelview matrix for observer-centric J2000 equatorial drawing
434+ //! Get the modelview matrix for observer-centric J2000 equatorial drawing.
435 StelProjector::ModelViewTranformP getJ2000ModelViewTransform(RefractionMode refMode=RefractionAuto) const;
436
437- //! Get the modelview matrix for observer-centric Galactic equatorial drawing
438+ //! Get the modelview matrix for observer-centric Galactic equatorial drawing.
439 StelProjector::ModelViewTranformP getGalacticModelViewTransform(RefractionMode refMode=RefractionAuto) const;
440
441- //! Rotation matrix from equatorial J2000 to ecliptic (Vsop87)
442+ //! Rotation matrix from equatorial J2000 to ecliptic (VSOP87A).
443 static const Mat4d matJ2000ToVsop87;
444- //! Rotation matrix from ecliptic (Vsop87) to equatorial J2000
445+ //! Rotation matrix from ecliptic (VSOP87A) to equatorial J2000.
446 static const Mat4d matVsop87ToJ2000;
447 //! Rotation matrix from J2000 to Galactic reference frame, using FITS convention.
448 static const Mat4d matJ2000ToGalactic;
449@@ -249,6 +256,10 @@
450
451 const QSharedPointer<class Planet> getCurrentPlanet() const;
452
453+ //! Unfortunately we also need this.
454+ const StelObserver* getCurrentObserver() const;
455+
456+
457 SphericalCap getVisibleSkyArea() const;
458
459 //! Smoothly move the observer to the given location
460@@ -291,10 +302,12 @@
461 //! @return valid range
462 QString getCurrentDeltaTAlgorithmValidRange(double jDay, QString* marker) const;
463
464- //! Checks for current time of day - it's night or day?
465+ //! Checks for altitude of sun - is it night or day?
466+ //! @return true if sun higher than about -6 degrees, i.e. "day" includes civil twilight.
467+ //! @note Useful mostly for brightness-controlled GUI decisions like font colors.
468 bool isDay() const;
469
470- //! Get value of the current Julian epoch
471+ //! Get value of the current Julian epoch (i.e. current year with decimal fraction, e.g. 2012.34567)
472 double getCurrentEpoch() const;
473
474 //! Get the default Mapping used by the Projection
475@@ -451,6 +464,10 @@
476 void addTropicalMonth();
477 //! Add one mean tropical year to the simulation time.
478 void addTropicalYear();
479+ //! Add one mean tropical century to the simulation time.
480+ void addTropicalCentury();
481+ //! Add one Julian century to the simulation time.
482+ void addJulianCentury();
483
484 //! Subtract one synodic month to the simulation time.
485 void subtractSynodicMonth();
486@@ -466,7 +483,11 @@
487 //! Subtract one mean tropical month to the simulation time.
488 void subtractTropicalMonth();
489 //! Subtract one mean tropical year to the simulation time.
490- void subtractTropicalYear();
491+ void subtractTropicalYear();
492+ //! Subtract one mean tropical century to the simulation time.
493+ void subtractTropicalCentury();
494+ //! Subtract one Julian century to the simulation time.
495+ void subtractJulianCentury();
496
497 //! Add a number of Earth Solar days to the current simulation time
498 //! @param d the decimal number of days to add (use negative values to subtract)
499@@ -525,13 +546,17 @@
500 void resetSync();
501
502 // Matrices used for every coordinate transfo
503- Mat4d matHeliocentricEclipticToAltAz; // Transform from heliocentric ecliptic (Vsop87) to observer-centric altazimuthal coordinate
504- Mat4d matAltAzToHeliocentricEcliptic; // Transform from observer-centric altazimuthal coordinate to heliocentric ecliptic (Vsop87)
505- Mat4d matAltAzToEquinoxEqu; // Transform from observer-centric altazimuthal coordinate to Earth Equatorial
506- Mat4d matEquinoxEquToAltAz; // Transform from observer-centric altazimuthal coordinate to Earth Equatorial
507- Mat4d matHeliocentricEclipticToEquinoxEqu; // Transform from heliocentric ecliptic (Vsop87) to earth equatorial coordinate
508- Mat4d matEquinoxEquToJ2000;
509- Mat4d matJ2000ToEquinoxEqu;
510+ Mat4d matHeliocentricEclipticJ2000ToAltAz; // Transform from heliocentric ecliptic Cartesian (VSOP87A) to topocentric (StelObserver) altazimuthal coordinate
511+ Mat4d matAltAzToHeliocentricEclipticJ2000; // Transform from topocentric (StelObserver) altazimuthal coordinate to heliocentric ecliptic Cartesian (VSOP87A)
512+ Mat4d matAltAzToEquinoxEqu; // Transform from topocentric altazimuthal coordinate to Earth Equatorial
513+ Mat4d matEquinoxEquToAltAz; // Transform from Earth Equatorial to topocentric (StelObserver) altazimuthal coordinate
514+ Mat4d matHeliocentricEclipticToEquinoxEqu; // Transform from heliocentric ecliptic Cartesian (VSOP87A) to earth equatorial coordinate
515+ Mat4d matEquinoxEquToJ2000; // GZ For Earth, this should be inverse precession matrix. Yes, =Rz(VSOPoffset)Rx(eps0)Rz(-psiA)Rx(-omA)Rz(chiA)
516+ Mat4d matJ2000ToEquinoxEqu; // GZ: Should be precession matrix?
517+// Mat4d matEclOfDateToVsop87; // GZ NEW: precession of the ecliptic
518+// Mat4d matVsop87ToEclOfDate; // GZ NEW: precession of the ecliptic
519+
520+
521 Mat4d matJ2000ToAltAz;
522
523 Mat4d matAltAzModelView; // Modelview matrix for observer-centric altazimuthal drawing
524
525=== modified file 'src/core/StelObject.cpp'
526--- src/core/StelObject.cpp 2015-06-12 18:45:39 +0000
527+++ src/core/StelObject.cpp 2015-07-28 19:18:20 +0000
528@@ -28,6 +28,7 @@
529 #include "RefractionExtinction.hpp"
530 #include "StelLocation.hpp"
531 #include "SolarSystem.hpp"
532+#include "StelModuleMgr.hpp"
533
534 #include <QRegExp>
535 #include <QDebug>
536@@ -42,6 +43,7 @@
537 Vec3d StelObject::getSiderealPosGeometric(const StelCore* core) const
538 {
539 // Hour Angle corrected to Delta-T value
540+ // TODO: make code readable by calling siderealTime(JD_UT), this should not contain a deltaT in its algorithm.
541 double dt = (core->getDeltaT(core->getJDay())/240.)*M_PI/180.;
542 return Mat4d::zrotation(-core->getLocalSiderealTime()+dt)* getEquinoxEquatorialPos(core);
543 }
544@@ -52,6 +54,7 @@
545 Vec3d v=getAltAzPosApparent(core);
546 v = core->altAzToEquinoxEqu(v, StelCore::RefractionOff);
547 // Hour Angle corrected to Delta-T value
548+ // TODO: make code readable by calling siderealTime(JD_UT), this should not contain a deltaT in its algorithm.
549 double dt = (core->getDeltaT(core->getJDay())/240.)*M_PI/180.;
550 return Mat4d::zrotation(-core->getLocalSiderealTime()+dt)*v;
551 }
552@@ -109,6 +112,9 @@
553 {
554 bool withAtmosphere = core->getSkyDrawer()->getFlagHasAtmosphere();
555 bool withDecimalDegree = StelApp::getInstance().getFlagShowDecimalDegrees();;
556+ double az_app, alt_app;
557+ StelUtils::rectToSphe(&az_app,&alt_app,getAltAzPosApparent(core));
558+ Q_UNUSED(az_app);
559 QString cepoch = qc_("on date", "coordinates for current epoch");
560 QString res;
561 if (flags&RaDecJ2000)
562@@ -137,7 +143,7 @@
563 QString hadec;
564 StelUtils::rectToSphe(&ra_sidereal,&dec_sidereal,getSiderealPosGeometric(core));
565 ra_sidereal = 2.*M_PI-ra_sidereal;
566- if (withAtmosphere)
567+ if (withAtmosphere && (alt_app>-3.0*M_PI/180.0)) // Don't show refracted values much below horizon where model is meaningless.
568 {
569 StelUtils::rectToSphe(&ra_sidereal,&dec_sidereal,getSiderealPosApparent(core));
570 ra_sidereal = 2.*M_PI-ra_sidereal;
571@@ -175,16 +181,12 @@
572 az = 3.*M_PI - az; // N is zero, E is 90 degrees
573 if (az > M_PI*2)
574 az -= M_PI*2;
575- if (withAtmosphere)
576+ if (withAtmosphere && (alt_app>-3.0*M_PI/180.0)) // Don't show refracted altitude much below horizon where model is meaningless.
577 {
578- StelUtils::rectToSphe(&az,&alt,getAltAzPosApparent(core));
579- az = 3.*M_PI - az; // N is zero, E is 90 degrees
580- if (az > M_PI*2)
581- az -= M_PI*2;
582 if (withDecimalDegree)
583- res += q_("Az/Alt: %1/%2").arg(StelUtils::radToDecDegStr(az), StelUtils::radToDecDegStr(alt)) + " " + q_("(apparent)") + "<br>";
584+ res += q_("Az/Alt: %1/%2").arg(StelUtils::radToDecDegStr(az), StelUtils::radToDecDegStr(alt_app)) + " " + q_("(apparent)") + "<br>";
585 else
586- res += q_("Az/Alt: %1/%2").arg(StelUtils::radToDmsStr(az,true), StelUtils::radToDmsStr(alt,true)) + " " + q_("(apparent)") + "<br>";
587+ res += q_("Az/Alt: %1/%2").arg(StelUtils::radToDmsStr(az,true), StelUtils::radToDmsStr(alt_app,true)) + " " + q_("(apparent)") + "<br>";
588 }
589 else
590 {
591@@ -197,7 +199,15 @@
592
593 if (flags&EclipticCoord)
594 {
595- double ecl = core->getCurrentPlanet()->getRotObliquity(2451545.0);
596+ // N.B. Ecliptical coordinates are particularly earth-bound.
597+ // It may be OK to have terrestrial ecliptical coordinates of J2000.0 (standard epoch) because those are in practice linked with VSOP XY plane,
598+ // and because the ecliptical grid of J2000 is also shown for observers on other planets.
599+ // The formulation here has never computed the true position of any observer planet's orbital plane except for Earth,
600+ // or ever displayed the coordinates in the observer planet's equivalent to Earth's ecliptical coordinates.
601+ // As quick test you can observe if in any "Ecliptic coordinate" as seen from e.g. Mars or Jupiter the Sun was ever close to beta=0.
602+
603+ //double ecl = core->getCurrentPlanet()->getRotObliquity(2451545.0);
604+ double ecl=GETSTELMODULE(SolarSystem)->getEarth()->getRotObliquity(2451545.0);
605 double ra_equ, dec_equ, lambda, beta;
606 StelUtils::rectToSphe(&ra_equ,&dec_equ,getJ2000EquatorialPos(core));
607 StelUtils::equToEcl(ra_equ, dec_equ, ecl, &lambda, &beta);
608@@ -206,14 +216,24 @@
609 res += q_("Ecliptic longitude/latitude") + QString(" (J%1): %2/%3").arg(QString::number(2000.f, 'f', 1), StelUtils::radToDecDegStr(lambda), StelUtils::radToDecDegStr(beta)) + "<br>";
610 else
611 res += q_("Ecliptic longitude/latitude") + QString(" (J%1): %2/%3").arg(QString::number(2000.f, 'f', 1), StelUtils::radToDmsStr(lambda, true), StelUtils::radToDmsStr(beta, true)) + "<br>";
612- ecl = core->getCurrentPlanet()->getRotObliquity(core->getJDay());
613- StelUtils::rectToSphe(&ra_equ,&dec_equ,getEquinoxEquatorialPos(core));
614- StelUtils::equToEcl(ra_equ, dec_equ, ecl, &lambda, &beta);
615- if (lambda<0) lambda+=2.0*M_PI;
616- if (withDecimalDegree)
617- res += q_("Ecliptic longitude/latitude") + QString(" (%1): %2/%3").arg(cepoch, StelUtils::radToDecDegStr(lambda), StelUtils::radToDecDegStr(beta)) + "<br>";
618- else
619- res += q_("Ecliptic longitude/latitude") + QString(" (%1): %2/%3").arg(cepoch, StelUtils::radToDmsStr(lambda, true), StelUtils::radToDmsStr(beta, true)) + "<br>";
620+
621+ if (core->getCurrentPlanet()->getEnglishName()=="Earth")
622+ {
623+ ecl = GETSTELMODULE(SolarSystem)->getEarth()->getRotObliquity(2451545.0);
624+
625+ StelUtils::rectToSphe(&ra_equ,&dec_equ,getEquinoxEquatorialPos(core));
626+ StelUtils::equToEcl(ra_equ, dec_equ, ecl, &lambda, &beta);
627+ if (lambda<0) lambda+=2.0*M_PI;
628+ if (withDecimalDegree)
629+ res += q_("Ecliptic longitude/latitude") + QString(" (%1): %2/%3").arg(cepoch, StelUtils::radToDecDegStr(lambda), StelUtils::radToDecDegStr(beta)) + "<br>";
630+ else
631+ res += q_("Ecliptic longitude/latitude") + QString(" (%1): %2/%3").arg(cepoch, StelUtils::radToDmsStr(lambda, true), StelUtils::radToDmsStr(beta, true)) + "<br>";
632+ // GZ Only for now: display epsilon_A, angle between Earth's Axis and ecl. of date.
633+ if (withDecimalDegree)
634+ res += q_("Ecliptic obliquity") + QString(" (%1): %2").arg(cepoch, StelUtils::radToDecDegStr(ecl)) + "<br>";
635+ else
636+ res += q_("Ecliptic obliquity") + QString(" (%1): %2").arg(cepoch, StelUtils::radToDmsStr(ecl)) + "<br>";
637+ }
638 }
639
640 if (flags&GalacticCoord)
641@@ -249,8 +269,8 @@
642 StelCore* core = StelApp::getInstance().getCore();
643 if (core->isDay() && core->getSkyDrawer()->getFlagHasAtmosphere()==true)
644 {
645- // Let's make info text is more readable when atmosphere enabled at daylight.
646- color = StelUtils::strToVec3f(StelApp::getInstance().getSettings()->value("color/daylight_color", "0.0,0.0,0.0").toString());
647+ // make info text more readable when atmosphere enabled at daylight.
648+ color = StelUtils::strToVec3f(StelApp::getInstance().getSettings()->value("color/daylight_text_color", "0.0,0.0,0.0").toString());
649 }
650 str.prepend(QString("<font color=%1>").arg(StelUtils::vec3fToHtmlColor(color)));
651 str.append(QString("</font>"));
652
653=== modified file 'src/core/StelObject.hpp'
654--- src/core/StelObject.hpp 2014-10-17 22:54:15 +0000
655+++ src/core/StelObject.hpp 2015-07-28 19:18:20 +0000
656@@ -57,13 +57,23 @@
657 GalacticCoord = 0x00000800, //!< The galactic position
658 ObjectType = 0x00001000, //!< The type of the object (star, planet, etc.)
659 EclipticCoord = 0x00002000, //!< The ecliptic position
660- PlainText = 0x00004000 //!< Strip HTML tags from output
661+ EclipticCoordXYZ = 0x00004000, //!< The ecliptic position, XYZ of VSOP87A (used mainly for debugging, not public)
662+ PlainText = 0x00010000, //!< Strip HTML tags from output
663+// TODO GZ
664+// RaDecJ2000Planetocentric = 0x00020000, //!< The planetocentric equatorial position (J2000 ref) [Mostly to compare with almanacs]
665+// RaDecOfDatePlanetocentric = 0x00040000 //!< The planetocentric equatorial position (of date)
666+// // and split Ecliptical into
667+// EclipticCoordJ2000 = 0x00002000, //!< The ecliptic position w.r.t. ecliptic of eq.J2000.0
668+// EclipticCoordOfDate = 0x00002000, //!< The ecliptic position w.r.t. ecliptic of eq. of date
669+
670+
671 };
672 typedef QFlags<InfoStringGroupFlags> InfoStringGroup;
673 Q_FLAGS(InfoStringGroup)
674
675 //! A pre-defined set of specifiers for the getInfoString flags argument to getInfoString
676- static const InfoStringGroupFlags AllInfo = (InfoStringGroupFlags)(Name|CatalogNumber|Magnitude|RaDecJ2000|RaDecOfDate|AltAzi|Distance|Size|Extra|HourAngle|AbsoluteMagnitude|GalacticCoord|ObjectType|EclipticCoord);
677+ static const InfoStringGroupFlags AllInfo = (InfoStringGroupFlags)(Name|CatalogNumber|Magnitude|RaDecJ2000|RaDecOfDate|AltAzi|Distance|Size|Extra|HourAngle|
678+ AbsoluteMagnitude|GalacticCoord|ObjectType|EclipticCoord|EclipticCoordXYZ);
679 //! A pre-defined set of specifiers for the getInfoString flags argument to getInfoString
680 static const InfoStringGroupFlags ShortInfo = (InfoStringGroupFlags)(Name|CatalogNumber|Magnitude|RaDecJ2000);
681
682
683=== modified file 'src/core/StelObserver.cpp'
684--- src/core/StelObserver.cpp 2015-02-01 20:44:54 +0000
685+++ src/core/StelObserver.cpp 2015-07-28 19:18:20 +0000
686@@ -183,6 +183,7 @@
687
688 const QSharedPointer<Planet> StelObserver::getHomePlanet(void) const
689 {
690+ Q_ASSERT(planet);
691 return planet;
692 }
693
694@@ -191,22 +192,47 @@
695 return getHomePlanet()->getHeliocentricEclipticPos();
696 }
697
698+// Used to approximate solution with assuming a spherical planet.
699+// Since V0.14, following Meeus, Astr. Alg. 2nd ed, Ch.11.
700 double StelObserver::getDistanceFromCenter(void) const
701 {
702- return getHomePlanet()->getRadius() + (currentLocation.altitude/(1000*AU));
703+ if (getHomePlanet()->getRadius()==0.0) // the transitional ArtificialPlanet od SpaceShipObserver has this
704+ return currentLocation.altitude/(1000*AU);
705+
706+ double a=getHomePlanet()->getRadius();
707+ double bByA = getHomePlanet()->getOneMinusOblateness(); // b/a;
708+
709+ if (fabs(currentLocation.latitude)>=89.9) // avoid tan(90) issues.
710+ return a * bByA;
711+
712+ double latRad=currentLocation.latitude*(M_PI/180.0);
713+ double u = atan( bByA * tan(latRad));
714+ // qDebug() << "getDistanceFromCenter: a=" << a*AU << "b/a=" << bByA << "b=" << bByA*a *AU << "latRad=" << latRad << "u=" << u;
715+ Q_ASSERT(fabs(u)<= fabs(latRad));
716+ double altFix=(currentLocation.altitude/(1000.0*AU)) / a;
717+
718+ double rhoSinPhiPrime= bByA * sin(u) + altFix*sin(latRad);
719+ double rhoCosPhiPrime= cos(u) + altFix*cos(latRad);
720+
721+ double rho = sqrt(rhoSinPhiPrime*rhoSinPhiPrime+rhoCosPhiPrime*rhoCosPhiPrime);
722+ return rho*a;
723 }
724
725 Mat4d StelObserver::getRotAltAzToEquatorial(double jd) const
726 {
727 double lat = currentLocation.latitude;
728- // TODO: Figure out how to keep continuity in sky as reach poles
729+ // TODO: Figure out how to keep continuity in sky as we reach poles
730 // otherwise sky jumps in rotation when reach poles in equatorial mode
731 // This is a kludge
732+ // GZ: Actually, why would that be? Lat should be clamped elsewhere. Added tests.
733+ Q_ASSERT(lat <= 90.0);
734+ Q_ASSERT(lat >= -90.0);
735 if( lat > 90.0 ) lat = 90.0;
736 if( lat < -90.0 ) lat = -90.0;
737 // Include a DeltaT correction. Sidereal time and longitude here are both in degrees, but DeltaT in seconds of time.
738 // 360 degrees = 24hrs; 15 degrees = 1hr = 3600s; 1 degree = 240s
739 // Apply DeltaT correction only for Earth
740+ // TODO: make code readable by calling siderealTime(JD_UT), this should not contain a deltaT in its algorithm.
741 double deltaT = 0.;
742 if (getHomePlanet()->getEnglishName()=="Earth")
743 deltaT = StelApp::getInstance().getCore()->getDeltaT(jd)/240.;
744
745=== modified file 'src/core/VecMath.hpp'
746--- src/core/VecMath.hpp 2015-04-01 14:36:44 +0000
747+++ src/core/VecMath.hpp 2015-07-28 19:18:20 +0000
748@@ -241,6 +241,7 @@
749 inline void normalize();
750
751 inline void transfo4d(const Mat4d&);
752+ QString toString() const {return QString("[%1, %2, %3, %4]").arg(v[0]).arg(v[1]).arg(v[2]).arg(v[3]);}
753
754 T v[4]; // The 4 values
755 };
756@@ -273,8 +274,22 @@
757
758 Matrix3<T> transpose() const;
759 Matrix3<T> inverse() const;
760+ //! return trace (sum of diagonal elements).
761+ inline T trace() const {return r[0]+r[4]+r[8];}
762+ //! return rotational angle
763+ inline T angle() const {return acos(0.5*(this->trace()-1.0));}
764
765 inline void print(void) const;
766+ QString toString(int fieldWidth=0, char format='g', int precision=-1) const {return QString("[[%1, %2, %3], [%4, %5, %6], [%7, %8, %9]]")
767+ .arg(r[0], fieldWidth, format, precision)
768+ .arg(r[1], fieldWidth, format, precision)
769+ .arg(r[2], fieldWidth, format, precision)
770+ .arg(r[3], fieldWidth, format, precision)
771+ .arg(r[4], fieldWidth, format, precision)
772+ .arg(r[5], fieldWidth, format, precision)
773+ .arg(r[6], fieldWidth, format, precision)
774+ .arg(r[7], fieldWidth, format, precision)
775+ .arg(r[8], fieldWidth, format, precision);}
776
777 T r[9];
778 };
779@@ -322,8 +337,9 @@
780 Matrix4<T> transpose() const;
781 Matrix4<T> inverse() const;
782
783- //Returns the upper 3x3 Matrix of the current 4x4 Matrix
784- Matrix3<T> upper3x3() const;
785+ //Returns the upper 3x3 Matrix of the current 4x4 Matrix
786+ Matrix3<T> upper3x3() const;
787+ Matrix3<T> upper3x3Transposed() const;
788
789 inline Vector4<T> getRow(const int row) const;
790 inline Vector4<T> getColumn(const int column) const;
791@@ -332,6 +348,23 @@
792 inline QMatrix4x4 convertToQMatrix() const;
793
794 inline void print(void) const;
795+ QString toString(int fieldWidth=0, char format='g', int precision=-1) const {return QString("[[%1, %2, %3, %4], [%5, %6, %7, %8], [%9, %10, %11, %12], [%13, %14, %15, %16]]")
796+ .arg(r[0], fieldWidth, format, precision)
797+ .arg(r[1], fieldWidth, format, precision)
798+ .arg(r[2], fieldWidth, format, precision)
799+ .arg(r[3], fieldWidth, format, precision)
800+ .arg(r[4], fieldWidth, format, precision)
801+ .arg(r[5], fieldWidth, format, precision)
802+ .arg(r[6], fieldWidth, format, precision)
803+ .arg(r[7], fieldWidth, format, precision)
804+ .arg(r[8], fieldWidth, format, precision)
805+ .arg(r[9], fieldWidth, format, precision)
806+ .arg(r[10], fieldWidth, format, precision)
807+ .arg(r[11], fieldWidth, format, precision)
808+ .arg(r[12], fieldWidth, format, precision)
809+ .arg(r[13], fieldWidth, format, precision)
810+ .arg(r[14], fieldWidth, format, precision)
811+ .arg(r[15], fieldWidth, format, precision);}
812
813 T r[16];
814 };
815@@ -1467,8 +1500,15 @@
816 #undef SWAP_ROWS
817 }
818
819+
820 template<class T> Matrix3<T> Matrix4<T>::upper3x3() const
821 {
822+ return Matrix3<T>(r[0], r[1], r[2],
823+ r[4], r[5], r[6],
824+ r[8], r[9], r[10]);
825+}
826+template<class T> Matrix3<T> Matrix4<T>::upper3x3Transposed() const
827+{
828 return Matrix3<T>(r[0], r[4], r[8],
829 r[1], r[5], r[9],
830 r[2], r[6], r[10]);
831
832=== modified file 'src/core/modules/GridLinesMgr.cpp'
833--- src/core/modules/GridLinesMgr.cpp 2015-03-26 15:44:34 +0000
834+++ src/core/modules/GridLinesMgr.cpp 2015-07-28 19:18:20 +0000
835@@ -29,6 +29,7 @@
836 #include "StelCore.hpp"
837 #include "StelPainter.hpp"
838 #include "StelSkyDrawer.hpp"
839+#include "precession.h"
840
841 #include <set>
842 #include <QSettings>
843@@ -66,15 +67,19 @@
844 public:
845 enum SKY_LINE_TYPE
846 {
847- EQUATOR,
848- ECLIPTIC,
849+ EQUATOR_J2000,
850+ EQUATOR_OF_DATE,
851+ ECLIPTIC_J2000,
852+ ECLIPTIC_OF_DATE,
853+ PRECESSIONCIRCLE_N,
854+ PRECESSIONCIRCLE_S,
855 MERIDIAN,
856 HORIZON,
857 GALACTICEQUATOR,
858 LONGITUDE
859 };
860 // Create and precompute positions of a SkyGrid
861- SkyLine(SKY_LINE_TYPE _line_type = EQUATOR);
862+ SkyLine(SKY_LINE_TYPE _line_type = EQUATOR_J2000);
863 virtual ~SkyLine();
864 void draw(StelCore* core) const;
865 void setColor(const Vec3f& c) {color = c;}
866@@ -202,7 +207,8 @@
867
868 break;
869 }
870- case StelCore::FrameObservercentricEcliptic:
871+ case StelCore::FrameObservercentricEclipticJ2000:
872+ case StelCore::FrameObservercentricEclipticOfDate:
873 {
874 double raAngle = d->raAngle;
875 if (raAngle<0.)
876@@ -355,12 +361,11 @@
877 if (QOpenGLContext::currentContext()->format().renderableType()==QSurfaceFormat::OpenGL)
878 glEnable(GL_LINE_SMOOTH);
879 #endif
880- Vec4f textColor(color[0], color[1], color[2], 0);
881+
882+ // make text colors just a bit brighter. (But if >1, QColor::setRgb fails and makes text invisible.)
883+ Vec4f textColor(qMin(1.0f, 1.25f*color[0]), qMin(1.0f, 1.25f*color[1]), qMin(1.0f, 1.25f*color[2]), fader.getInterstate());
884 sPainter.setColor(color[0],color[1],color[2], fader.getInterstate());
885
886- textColor*=2;
887- textColor[3]=fader.getInterstate();
888-
889 sPainter.setFont(font);
890 ViewportEdgeIntersectCallbackData userData(&sPainter);
891 userData.textColor = textColor;
892@@ -586,14 +591,27 @@
893 frameType = StelCore::FrameAltAz;
894 label = q_("Meridian");
895 break;
896- case ECLIPTIC:
897- frameType = StelCore::FrameObservercentricEcliptic;
898- label = q_("Ecliptic");
899- break;
900- case EQUATOR:
901+ case ECLIPTIC_J2000:
902+ frameType = StelCore::FrameObservercentricEclipticJ2000;
903+ label = q_("Ecliptic of J2000.0");
904+ break;
905+ case ECLIPTIC_OF_DATE:
906+ frameType = StelCore::FrameObservercentricEclipticOfDate;
907+ label = q_("Ecliptic of Date");
908+ break;
909+ case EQUATOR_J2000:
910+ frameType = StelCore::FrameJ2000;
911+ label = q_("Equator of J2000.0");
912+ break;
913+ case EQUATOR_OF_DATE:
914 frameType = StelCore::FrameEquinoxEqu;
915 label = q_("Equator");
916 break;
917+ case PRECESSIONCIRCLE_N:
918+ case PRECESSIONCIRCLE_S:
919+ frameType = StelCore::FrameObservercentricEclipticOfDate;
920+ label = q_("Precession Circle");
921+ break;
922 case HORIZON:
923 frameType = StelCore::FrameAltAz;
924 label = q_("Horizon");
925@@ -603,10 +621,12 @@
926 label = q_("Galactic Equator");
927 break;
928 case LONGITUDE:
929- frameType = StelCore::FrameObservercentricEcliptic;
930+ frameType = StelCore::FrameObservercentricEclipticJ2000; // TODO: Switch to FrameObservercentriEclipticOfDate
931 // TRANSLATORS: Full term is "opposition/conjunction longitude"
932 label = q_("O./C. longitude");
933 break;
934+ default:
935+ Q_ASSERT(0);
936 }
937 }
938
939@@ -625,7 +645,10 @@
940 sPainter.setColor(color[0], color[1], color[2], fader.getInterstate());
941 glEnable(GL_BLEND);
942 glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA); // Normal transparency mode
943-
944+ #ifdef GL_LINE_SMOOTH
945+ if (QOpenGLContext::currentContext()->format().renderableType()==QSurfaceFormat::OpenGL)
946+ glEnable(GL_LINE_SMOOTH);
947+ #endif
948 Vec4f textColor(color[0], color[1], color[2], 0);
949 textColor[3]=fader.getInterstate();
950
951@@ -635,6 +658,79 @@
952 userData.text = label;
953 /////////////////////////////////////////////////
954 // Draw the line
955+
956+ // Precession circles are Small Circles, all others are Great Circles.
957+ if (line_type==PRECESSIONCIRCLE_N || line_type==PRECESSIONCIRCLE_S)
958+ {
959+ const double lat=(line_type==PRECESSIONCIRCLE_S ? -1.0 : 1.0) * (M_PI/2.0-getPrecessionAngleVondrakCurrentEpsilonA());
960+ SphericalCap declinationCap(Vec3d(0,0,1), std::sin(lat));
961+ const Vec3d rotCenter(0,0,declinationCap.d);
962+
963+ Vec3d p1, p2;
964+ if (!SphericalCap::intersectionPoints(viewPortSphericalCap, declinationCap, p1, p2))
965+ {
966+ if ((viewPortSphericalCap.d<declinationCap.d && viewPortSphericalCap.contains(declinationCap.n))
967+ || (viewPortSphericalCap.d<-declinationCap.d && viewPortSphericalCap.contains(-declinationCap.n)))
968+ {
969+ // The line is fully included in the viewport, draw it in 3 sub-arcs to avoid length > 180.
970+ Vec3d pt1;
971+ Vec3d pt2;
972+ Vec3d pt3;
973+ const double lon1=0.0;
974+ const double lon2=120.0*M_PI/180.0;
975+ const double lon3=240.0*M_PI/180.0;
976+ StelUtils::spheToRect(lon1, lat, pt1); pt1.normalize();
977+ StelUtils::spheToRect(lon2, lat, pt2); pt2.normalize();
978+ StelUtils::spheToRect(lon3, lat, pt3); pt3.normalize();
979+
980+ sPainter.drawSmallCircleArc(pt1, pt2, rotCenter, viewportEdgeIntersectCallback, &userData);
981+ sPainter.drawSmallCircleArc(pt2, pt3, rotCenter, viewportEdgeIntersectCallback, &userData);
982+ sPainter.drawSmallCircleArc(pt3, pt1, rotCenter, viewportEdgeIntersectCallback, &userData);
983+ #ifdef GL_LINE_SMOOTH
984+ if (QOpenGLContext::currentContext()->format().renderableType()==QSurfaceFormat::OpenGL)
985+ glDisable(GL_LINE_SMOOTH);
986+ #endif
987+ glDisable(GL_BLEND);
988+ return;
989+ }
990+ else
991+ {
992+ #ifdef GL_LINE_SMOOTH
993+ if (QOpenGLContext::currentContext()->format().renderableType()==QSurfaceFormat::OpenGL)
994+ glDisable(GL_LINE_SMOOTH);
995+ #endif
996+ glDisable(GL_BLEND);
997+ return;
998+ }
999+ }
1000+ // Draw the arc in 2 sub-arcs to avoid lengths > 180 deg
1001+ Vec3d middlePoint = p1-rotCenter+p2-rotCenter;
1002+ middlePoint.normalize();
1003+ middlePoint*=(p1-rotCenter).length();
1004+ middlePoint+=rotCenter;
1005+ if (!viewPortSphericalCap.contains(middlePoint))
1006+ {
1007+ middlePoint-=rotCenter;
1008+ middlePoint*=-1.;
1009+ middlePoint+=rotCenter;
1010+ }
1011+
1012+ sPainter.drawSmallCircleArc(p1, middlePoint, rotCenter,viewportEdgeIntersectCallback, &userData);
1013+ sPainter.drawSmallCircleArc(p2, middlePoint, rotCenter, viewportEdgeIntersectCallback, &userData);
1014+
1015+ // OpenGL ES 2.0 doesn't have GL_LINE_SMOOTH
1016+ #ifdef GL_LINE_SMOOTH
1017+ if (QOpenGLContext::currentContext()->format().renderableType()==QSurfaceFormat::OpenGL)
1018+ glDisable(GL_LINE_SMOOTH);
1019+ #endif
1020+
1021+ glDisable(GL_BLEND);
1022+
1023+
1024+ return;
1025+ }
1026+
1027+ // All the other "lines" are Great Circles
1028 SphericalCap meridianSphericalCap(Vec3d(0,0,1), 0);
1029 Vec3d fpt(1,0,0);
1030 if (line_type==MERIDIAN)
1031@@ -683,6 +779,14 @@
1032 sPainter.drawGreatCircleArc(p1, middlePoint, NULL, viewportEdgeIntersectCallback, &userData);
1033 sPainter.drawGreatCircleArc(p2, middlePoint, NULL, viewportEdgeIntersectCallback, &userData);
1034
1035+ // OpenGL ES 2.0 doesn't have GL_LINE_SMOOTH
1036+ #ifdef GL_LINE_SMOOTH
1037+ if (QOpenGLContext::currentContext()->format().renderableType()==QSurfaceFormat::OpenGL)
1038+ glDisable(GL_LINE_SMOOTH);
1039+ #endif
1040+
1041+ glDisable(GL_BLEND);
1042+
1043 // // Johannes: use a big radius as a dirty workaround for the bug that the
1044 // // ecliptic line is not drawn around the observer, but around the sun:
1045 // const Vec3d vv(1000000,0,0);
1046@@ -694,11 +798,16 @@
1047 setObjectName("GridLinesMgr");
1048 equGrid = new SkyGrid(StelCore::FrameEquinoxEqu);
1049 equJ2000Grid = new SkyGrid(StelCore::FrameJ2000);
1050- eclJ2000Grid = new SkyGrid(StelCore::FrameObservercentricEcliptic);
1051+ eclJ2000Grid = new SkyGrid(StelCore::FrameObservercentricEclipticJ2000);
1052+ eclGrid = new SkyGrid(StelCore::FrameObservercentricEclipticOfDate);
1053 galacticGrid = new SkyGrid(StelCore::FrameGalactic);
1054 aziGrid = new SkyGrid(StelCore::FrameAltAz);
1055- equatorLine = new SkyLine(SkyLine::EQUATOR);
1056- eclipticLine = new SkyLine(SkyLine::ECLIPTIC);
1057+ equatorLine = new SkyLine(SkyLine::EQUATOR_OF_DATE);
1058+ equatorJ2000Line = new SkyLine(SkyLine::EQUATOR_J2000);
1059+ eclipticJ2000Line = new SkyLine(SkyLine::ECLIPTIC_J2000); // previous eclipticLine
1060+ eclipticLine = new SkyLine(SkyLine::ECLIPTIC_OF_DATE); // GZ NEW
1061+ precessionCircleN = new SkyLine(SkyLine::PRECESSIONCIRCLE_N);
1062+ precessionCircleS = new SkyLine(SkyLine::PRECESSIONCIRCLE_S);
1063 meridianLine = new SkyLine(SkyLine::MERIDIAN);
1064 horizonLine = new SkyLine(SkyLine::HORIZON);
1065 galacticEquatorLine = new SkyLine(SkyLine::GALACTICEQUATOR);
1066@@ -710,10 +819,15 @@
1067 delete equGrid;
1068 delete equJ2000Grid;
1069 delete eclJ2000Grid;
1070+ delete eclGrid;
1071 delete galacticGrid;
1072 delete aziGrid;
1073 delete equatorLine;
1074+ delete equatorJ2000Line;
1075 delete eclipticLine;
1076+ delete eclipticJ2000Line;
1077+ delete precessionCircleN;
1078+ delete precessionCircleS;
1079 delete meridianLine;
1080 delete horizonLine;
1081 delete galacticEquatorLine;
1082@@ -739,9 +853,13 @@
1083 setFlagEquatorGrid(conf->value("viewing/flag_equatorial_grid").toBool());
1084 setFlagEquatorJ2000Grid(conf->value("viewing/flag_equatorial_J2000_grid").toBool());
1085 setFlagEclipticJ2000Grid(conf->value("viewing/flag_ecliptic_J2000_grid").toBool());
1086+ setFlagEclipticGrid(conf->value("viewing/flag_ecliptic_grid").toBool());
1087 setFlagGalacticGrid(conf->value("viewing/flag_galactic_grid").toBool());
1088 setFlagEquatorLine(conf->value("viewing/flag_equator_line").toBool());
1089+ setFlagEquatorJ2000Line(conf->value("viewing/flag_equator_J2000_line").toBool());
1090 setFlagEclipticLine(conf->value("viewing/flag_ecliptic_line").toBool());
1091+ setFlagEclipticJ2000Line(conf->value("viewing/flag_ecliptic_J2000_line").toBool());
1092+ setFlagPrecessionCircles(conf->value("viewing/flag_precession_circles").toBool());
1093 setFlagMeridianLine(conf->value("viewing/flag_meridian_line").toBool());
1094 setFlagHorizonLine(conf->value("viewing/flag_horizon_line").toBool());
1095 setFlagGalacticEquatorLine(conf->value("viewing/flag_galactic_equator_line").toBool());
1096@@ -755,14 +873,18 @@
1097 addAction("actionShow_Equatorial_Grid", displayGroup, N_("Equatorial grid"), "equatorGridDisplayed", "E");
1098 addAction("actionShow_Azimuthal_Grid", displayGroup, N_("Azimuthal grid"), "azimuthalGridDisplayed", "Z");
1099 addAction("actionShow_Ecliptic_Line", displayGroup, N_("Ecliptic line"), "eclipticLineDisplayed", ",");
1100+ addAction("actionShow_Ecliptic_J2000_Line", displayGroup, N_("Ecliptic J2000 line"), "eclipticJ2000LineDisplayed");
1101 addAction("actionShow_Equator_Line", displayGroup, N_("Equator line"), "equatorLineDisplayed", ".");
1102+ addAction("actionShow_Equator_J2000_Line", displayGroup, N_("Equator J2000 line"), "equatorJ2000LineDisplayed"); // or with Hotkey??
1103 addAction("actionShow_Meridian_Line", displayGroup, N_("Meridian line"), "meridianLineDisplayed", ";");
1104 addAction("actionShow_Horizon_Line", displayGroup, N_("Horizon line"), "horizonLineDisplayed", "H");
1105 addAction("actionShow_Equatorial_J2000_Grid", displayGroup, N_("Equatorial J2000 grid"), "equatorJ2000GridDisplayed");
1106 addAction("actionShow_Ecliptic_J2000_Grid", displayGroup, N_("Ecliptic J2000 grid"), "eclipticJ2000GridDisplayed");
1107+ addAction("actionShow_Ecliptic_Grid", displayGroup, N_("Ecliptic grid"), "eclipticGridDisplayed");
1108 addAction("actionShow_Galactic_Grid", displayGroup, N_("Galactic grid"), "galacticGridDisplayed");
1109 addAction("actionShow_Galactic_Equator_Line", displayGroup, N_("Galactic equator"), "galacticEquatorLineDisplayed");
1110 addAction("actionShow_Longitude_Line", displayGroup, N_("Opposition/conjunction longitude line"), "longitudeLineDisplayed");
1111+ addAction("actionShow_Precession_Circles", displayGroup, N_("Precession Circles"), "precessionCirclesDisplayed");
1112 }
1113
1114 void GridLinesMgr::update(double deltaTime)
1115@@ -771,10 +893,15 @@
1116 equGrid->update(deltaTime);
1117 equJ2000Grid->update(deltaTime);
1118 eclJ2000Grid->update(deltaTime);
1119+ eclGrid->update(deltaTime);
1120 galacticGrid->update(deltaTime);
1121 aziGrid->update(deltaTime);
1122 equatorLine->update(deltaTime);
1123+ equatorJ2000Line->update(deltaTime);
1124 eclipticLine->update(deltaTime);
1125+ eclipticJ2000Line->update(deltaTime);
1126+ precessionCircleN->update(deltaTime);
1127+ precessionCircleS->update(deltaTime);
1128 meridianLine->update(deltaTime);
1129 horizonLine->update(deltaTime);
1130 galacticEquatorLine->update(deltaTime);
1131@@ -783,17 +910,30 @@
1132
1133 void GridLinesMgr::draw(StelCore* core)
1134 {
1135- equGrid->draw(core);
1136 galacticGrid->draw(core);
1137+ eclJ2000Grid->draw(core);
1138+ // While ecliptic of J2000 may be helpful to get a feeling of the Z=0 plane of VSOP87,
1139+ // ecliptic of date is related to Earth and does not make much sense for the other planets.
1140+ // Of course, orbital plane of respective planet would be better, but is not implemented.
1141+ if (core->getCurrentPlanet()->getEnglishName()=="Earth")
1142+ eclGrid->draw(core);
1143 equJ2000Grid->draw(core);
1144- eclJ2000Grid->draw(core);
1145+ equGrid->draw(core);
1146 aziGrid->draw(core);
1147+ // Lines after grids, to be able to e.g. draw equators in different color!
1148+ galacticEquatorLine->draw(core);
1149+ eclipticJ2000Line->draw(core);
1150+ if (core->getCurrentPlanet()->getEnglishName()=="Earth")
1151+ {
1152+ eclipticLine->draw(core);
1153+ precessionCircleN->draw(core);
1154+ precessionCircleS->draw(core);
1155+ }
1156+ longitudeLine->draw(core);
1157+ equatorJ2000Line->draw(core);
1158 equatorLine->draw(core);
1159- eclipticLine->draw(core);
1160 meridianLine->draw(core);
1161 horizonLine->draw(core);
1162- galacticEquatorLine->draw(core);
1163- longitudeLine->draw(core);
1164 }
1165
1166 void GridLinesMgr::setStelStyle(const QString& section)
1167@@ -804,11 +944,14 @@
1168 QString defaultColor = conf->value(section+"/default_color").toString();
1169 setColorEquatorGrid(StelUtils::strToVec3f(conf->value(section+"/equatorial_color", defaultColor).toString()));
1170 setColorEquatorJ2000Grid(StelUtils::strToVec3f(conf->value(section+"/equatorial_J2000_color", defaultColor).toString()));
1171- setColorEclipticJ2000Grid(StelUtils::strToVec3f(conf->value(section+"/ecliptic_J2000_color", defaultColor).toString()));
1172+ setColorEclipticJ2000Grid(StelUtils::strToVec3f(conf->value(section+"/ecliptical_J2000_color", defaultColor).toString()));
1173+ setColorEclipticGrid(StelUtils::strToVec3f(conf->value(section+"/ecliptical_color", defaultColor).toString()));
1174 setColorGalacticGrid(StelUtils::strToVec3f(conf->value(section+"/galactic_color", defaultColor).toString()));
1175 setColorAzimuthalGrid(StelUtils::strToVec3f(conf->value(section+"/azimuthal_color", defaultColor).toString()));
1176 setColorEquatorLine(StelUtils::strToVec3f(conf->value(section+"/equator_color", defaultColor).toString()));
1177 setColorEclipticLine(StelUtils::strToVec3f(conf->value(section+"/ecliptic_color", defaultColor).toString()));
1178+ setColorEclipticJ2000Line(StelUtils::strToVec3f(conf->value(section+"/ecliptic_J2000_color", defaultColor).toString()));
1179+ setColorPrecessionCircles(StelUtils::strToVec3f(conf->value(section+"/precession_circles_color", defaultColor).toString()));
1180 setColorMeridianLine(StelUtils::strToVec3f(conf->value(section+"/meridian_color", defaultColor).toString()));
1181 setColorHorizonLine(StelUtils::strToVec3f(conf->value(section+"/horizon_color", defaultColor).toString()));
1182 setColorGalacticEquatorLine(StelUtils::strToVec3f(conf->value(section+"/galactic_equator_color", defaultColor).toString()));
1183@@ -817,8 +960,12 @@
1184
1185 void GridLinesMgr::updateLineLabels()
1186 {
1187+ equatorJ2000Line->updateLabel();
1188 equatorLine->updateLabel();
1189 eclipticLine->updateLabel();
1190+ eclipticJ2000Line->updateLabel();
1191+ precessionCircleN->updateLabel();
1192+ precessionCircleS->updateLabel();
1193 meridianLine->updateLabel();
1194 horizonLine->updateLabel();
1195 galacticEquatorLine->updateLabel();
1196@@ -925,6 +1072,31 @@
1197 }
1198 }
1199
1200+//! Set flag for displaying Ecliptic of Date Grid
1201+void GridLinesMgr::setFlagEclipticGrid(const bool displayed)
1202+{
1203+ if(displayed != eclGrid->isDisplayed()) {
1204+ eclGrid->setDisplayed(displayed);
1205+ emit eclipticGridDisplayedChanged(displayed);
1206+ }
1207+}
1208+//! Get flag for displaying Ecliptic of Date Grid
1209+bool GridLinesMgr::getFlagEclipticGrid(void) const
1210+{
1211+ return eclGrid->isDisplayed();
1212+}
1213+Vec3f GridLinesMgr::getColorEclipticGrid(void) const
1214+{
1215+ return eclGrid->getColor();
1216+}
1217+void GridLinesMgr::setColorEclipticGrid(const Vec3f& newColor)
1218+{
1219+ if(newColor != eclGrid->getColor()) {
1220+ eclGrid->setColor(newColor);
1221+ emit eclipticGridColorChanged(newColor);
1222+ }
1223+}
1224+
1225 //! Set flag for displaying Galactic Grid
1226 void GridLinesMgr::setFlagGalacticGrid(const bool displayed)
1227 {
1228@@ -975,6 +1147,31 @@
1229 }
1230 }
1231
1232+//! Set flag for displaying J2000 Equatorial Line
1233+void GridLinesMgr::setFlagEquatorJ2000Line(const bool displayed)
1234+{
1235+ if(displayed != equatorJ2000Line->isDisplayed()) {
1236+ equatorJ2000Line->setDisplayed(displayed);
1237+ emit equatorJ2000LineDisplayedChanged(displayed);
1238+ }
1239+}
1240+//! Get flag for displaying J2000 Equatorial Line
1241+bool GridLinesMgr::getFlagEquatorJ2000Line(void) const
1242+{
1243+ return equatorJ2000Line->isDisplayed();
1244+}
1245+Vec3f GridLinesMgr::getColorEquatorJ2000Line(void) const
1246+{
1247+ return equatorJ2000Line->getColor();
1248+}
1249+void GridLinesMgr::setColorEquatorJ2000Line(const Vec3f& newColor)
1250+{
1251+ if(newColor != equatorJ2000Line->getColor()) {
1252+ equatorJ2000Line->setColor(newColor);
1253+ emit equatorJ2000LineColorChanged(newColor);
1254+ }
1255+}
1256+
1257 //! Set flag for displaying Ecliptic Line
1258 void GridLinesMgr::setFlagEclipticLine(const bool displayed)
1259 {
1260@@ -1000,6 +1197,58 @@
1261 }
1262 }
1263
1264+//! Set flag for displaying Ecliptic J2000 Line
1265+void GridLinesMgr::setFlagEclipticJ2000Line(const bool displayed)
1266+{
1267+ if(displayed != eclipticJ2000Line->isDisplayed()) {
1268+ eclipticJ2000Line->setDisplayed(displayed);
1269+ emit eclipticJ2000LineDisplayedChanged(displayed);
1270+ }
1271+}
1272+//! Get flag for displaying Ecliptic J2000 Line
1273+bool GridLinesMgr::getFlagEclipticJ2000Line(void) const
1274+{
1275+ return eclipticJ2000Line->isDisplayed();
1276+}
1277+Vec3f GridLinesMgr::getColorEclipticJ2000Line(void) const
1278+{
1279+ return eclipticJ2000Line->getColor();
1280+}
1281+void GridLinesMgr::setColorEclipticJ2000Line(const Vec3f& newColor)
1282+{
1283+ if(newColor != eclipticJ2000Line->getColor()) {
1284+ eclipticJ2000Line->setColor(newColor);
1285+ emit eclipticJ2000LineColorChanged(newColor);
1286+ }
1287+}
1288+
1289+//! Set flag for displaying Precession Circles
1290+void GridLinesMgr::setFlagPrecessionCircles(const bool displayed)
1291+{
1292+ if(displayed != precessionCircleN->isDisplayed()) {
1293+ precessionCircleN->setDisplayed(displayed);
1294+ precessionCircleS->setDisplayed(displayed);
1295+ emit precessionCirclesDisplayedChanged(displayed);
1296+ }
1297+}
1298+//! Get flag for displaying Precession Circles
1299+bool GridLinesMgr::getFlagPrecessionCircles(void) const
1300+{
1301+ // precessionCircleS is always synchronous, no separate queries.
1302+ return precessionCircleN->isDisplayed();
1303+}
1304+Vec3f GridLinesMgr::getColorPrecessionCircles(void) const
1305+{
1306+ return precessionCircleN->getColor();
1307+}
1308+void GridLinesMgr::setColorPrecessionCircles(const Vec3f& newColor)
1309+{
1310+ if(newColor != precessionCircleN->getColor()) {
1311+ precessionCircleN->setColor(newColor);
1312+ precessionCircleS->setColor(newColor);
1313+ emit precessionCirclesColorChanged(newColor);
1314+ }
1315+}
1316
1317 //! Set flag for displaying Meridian Line
1318 void GridLinesMgr::setFlagMeridianLine(const bool displayed)
1319
1320=== modified file 'src/core/modules/GridLinesMgr.hpp'
1321--- src/core/modules/GridLinesMgr.hpp 2015-02-24 16:53:57 +0000
1322+++ src/core/modules/GridLinesMgr.hpp 2015-07-28 19:18:20 +0000
1323@@ -1,6 +1,7 @@
1324 /*
1325 * Stellarium
1326 * Copyright (C) 2007 Fabien Chereau
1327+ * Copyright (C) 2015 Georg Zotti (more grids to illustrate precession issues)
1328 *
1329 * This program is free software; you can redistribute it and/or
1330 * modify it under the terms of the GNU General Public License
1331@@ -27,8 +28,9 @@
1332 class SkyLine;
1333
1334 //! @class GridLinesMgr
1335-//! The GridLinesMgr controls the drawing of the Azimuthal and Equatorial Grids,
1336-//! as well as the great circles: Meridian Line, Ecliptic Line and Equator Line.
1337+//! The GridLinesMgr controls the drawing of the Azimuthal, Equatorial, Ecliptical and Galactic Grids,
1338+//! as well as the great circles: Meridian Line, Ecliptic Lines of J2000.0 and date, Equator Line (of J2000.0 and date),
1339+//! Precession Circles, and a special line marking conjunction or opposition in ecliptical longitude (of date).
1340 class GridLinesMgr : public StelModule
1341 {
1342 Q_OBJECT
1343@@ -40,6 +42,7 @@
1344 READ getColorAzimuthalGrid
1345 WRITE setColorAzimuthalGrid
1346 NOTIFY azimuthalGridColorChanged)
1347+
1348 Q_PROPERTY(bool equatorGridDisplayed
1349 READ getFlagEquatorGrid
1350 WRITE setFlagEquatorGrid
1351@@ -48,18 +51,34 @@
1352 READ getColorEquatorGrid
1353 WRITE setColorEquatorGrid
1354 NOTIFY equatorGridColorChanged)
1355+
1356 Q_PROPERTY(bool equatorJ2000GridDisplayed
1357 READ getFlagEquatorJ2000Grid
1358 WRITE setFlagEquatorJ2000Grid
1359 NOTIFY equatorJ2000GridDisplayedChanged)
1360+ Q_PROPERTY(Vec3f equatorJ2000GridColor
1361+ READ getColorEquatorJ2000Grid
1362+ WRITE setColorEquatorJ2000Grid
1363+ NOTIFY equatorJ2000GridColorChanged)
1364+
1365 Q_PROPERTY(bool eclipticJ2000GridDisplayed
1366 READ getFlagEclipticJ2000Grid
1367 WRITE setFlagEclipticJ2000Grid
1368 NOTIFY eclipticJ2000GridDisplayedChanged)
1369- Q_PROPERTY(Vec3f equatorJ2000GridColor
1370- READ getColorEquatorJ2000Grid
1371- WRITE setColorEquatorJ2000Grid
1372- NOTIFY equatorJ2000GridColorChanged)
1373+ Q_PROPERTY(Vec3f eclipticJ2000GridColor
1374+ READ getColorEclipticJ2000Grid
1375+ WRITE setColorEclipticJ2000Grid
1376+ NOTIFY eclipticJ2000GridColorChanged)
1377+ // NEW
1378+ Q_PROPERTY(bool eclipticGridDisplayed
1379+ READ getFlagEclipticGrid
1380+ WRITE setFlagEclipticGrid
1381+ NOTIFY eclipticGridDisplayedChanged)
1382+ Q_PROPERTY(Vec3f eclipticGridColor
1383+ READ getColorEclipticGrid
1384+ WRITE setColorEclipticGrid
1385+ NOTIFY eclipticGridColorChanged)
1386+
1387 Q_PROPERTY(bool galacticGridDisplayed
1388 READ getFlagGalacticGrid
1389 WRITE setFlagGalacticGrid
1390@@ -68,6 +87,7 @@
1391 READ getColorGalacticGrid
1392 WRITE setColorGalacticGrid
1393 NOTIFY galacticGridColorChanged)
1394+
1395 Q_PROPERTY(bool equatorLineDisplayed
1396 READ getFlagEquatorLine
1397 WRITE setFlagEquatorLine
1398@@ -76,6 +96,17 @@
1399 READ getColorEquatorLine
1400 WRITE setColorEquatorLine
1401 NOTIFY equatorLineColorChanged)
1402+
1403+ Q_PROPERTY(bool equatorJ2000LineDisplayed
1404+ READ getFlagEquatorJ2000Line
1405+ WRITE setFlagEquatorJ2000Line
1406+ NOTIFY equatorJ2000LineDisplayedChanged)
1407+ Q_PROPERTY(Vec3f equatorJ2000LineColor
1408+ READ getColorEquatorJ2000Line
1409+ WRITE setColorEquatorJ2000Line
1410+ NOTIFY equatorJ2000LineColorChanged)
1411+
1412+ // This is now ecl. of date.
1413 Q_PROPERTY(bool eclipticLineDisplayed
1414 READ getFlagEclipticLine
1415 WRITE setFlagEclipticLine
1416@@ -84,6 +115,26 @@
1417 READ getColorEclipticLine
1418 WRITE setColorEclipticLine
1419 NOTIFY eclipticLineColorChanged)
1420+
1421+ // new name, but replaces old ecliptic.
1422+ Q_PROPERTY(bool eclipticJ2000LineDisplayed
1423+ READ getFlagEclipticJ2000Line
1424+ WRITE setFlagEclipticJ2000Line
1425+ NOTIFY eclipticJ2000LineDisplayedChanged)
1426+ Q_PROPERTY(Vec3f eclipticJ2000LineColor
1427+ READ getColorEclipticJ2000Line
1428+ WRITE setColorEclipticJ2000Line
1429+ NOTIFY eclipticJ2000LineColorChanged)
1430+
1431+ Q_PROPERTY(bool precessionCirclesDisplayed
1432+ READ getFlagPrecessionCircles
1433+ WRITE setFlagPrecessionCircles
1434+ NOTIFY precessionCirclesDisplayedChanged)
1435+ Q_PROPERTY(Vec3f precessionCirclesColor
1436+ READ getColorPrecessionCircles
1437+ WRITE setColorPrecessionCircles
1438+ NOTIFY precessionCirclesColorChanged)
1439+
1440 Q_PROPERTY(bool meridianLineDisplayed
1441 READ getFlagMeridianLine
1442 WRITE setFlagMeridianLine
1443@@ -92,6 +143,7 @@
1444 READ getColorMeridianLine
1445 WRITE setColorMeridianLine
1446 NOTIFY meridianLineColorChanged)
1447+
1448 Q_PROPERTY(bool longitudeLineDisplayed
1449 READ getFlagLongitudeLine
1450 WRITE setFlagLongitudeLine
1451@@ -100,6 +152,7 @@
1452 READ getColorLongitudeLine
1453 WRITE setColorLongitudeLine
1454 NOTIFY longitudeLineColorChanged)
1455+
1456 Q_PROPERTY(bool horizonLineDisplayed
1457 READ getFlagHorizonLine
1458 WRITE setFlagHorizonLine
1459@@ -108,6 +161,7 @@
1460 READ getColorHorizonLine
1461 WRITE setColorHorizonLine
1462 NOTIFY horizonLineColorChanged)
1463+
1464 Q_PROPERTY(bool galacticEquatorLineDisplayed
1465 READ getFlagGalacticEquatorLine
1466 WRITE setFlagGalacticEquatorLine
1467@@ -125,17 +179,17 @@
1468 // Methods defined in the StelModule class
1469 //! Initialize the GridLinesMgr. This process checks the values in the
1470 //! application settings, and according to the settings there
1471- //! sets the visibility of the Equatorial Grid, Azimuthal Grid, Meridian Line,
1472- //! Equator Line and Ecliptic Line.
1473+ //! sets the visibility of the Equatorial Grids, Ecliptical Grids, Azimuthal Grid, Meridian Line,
1474+ //! Equator Line and Ecliptic Lines.
1475 virtual void init();
1476
1477- //! Get the module ID, returns, "gridlines".
1478+ //! Get the module ID, returns "GridLinesMgr".
1479 virtual QString getModuleID() const {return "GridLinesMgr";}
1480
1481 //! Draw the grids and great circle lines.
1482- //! Draws the Equatorial Grid, Azimuthal Grid, Meridian Line, Equator Line
1483- //! and Ecliptic Line according to the various flags which control their
1484- //! visibility.
1485+ //! Draws the Equatorial Grids, Ecliptical Grids, Azimuthal Grid, Meridian Line, Equator Line,
1486+ //! Ecliptic Lines, Precession Circles, and Conjunction-Opposition Line according to the
1487+ //! various flags which control their visibility.
1488 virtual void draw(StelCore* core);
1489
1490 //! Update time-dependent features.
1491@@ -176,9 +230,9 @@
1492 //! @endcode
1493 void setColorEquatorGrid(const Vec3f& newColor);
1494
1495- //! Setter for displaying Equatorial Grid.
1496+ //! Setter for displaying Equatorial J2000 Grid.
1497 void setFlagEquatorJ2000Grid(const bool displayed);
1498- //! Accessor for displaying Equatorial Grid.
1499+ //! Accessor for displaying Equatorial J2000 Grid.
1500 bool getFlagEquatorJ2000Grid(void) const;
1501 //! Get the current color of the Equatorial J2000 Grid.
1502 Vec3f getColorEquatorJ2000Grid(void) const;
1503@@ -190,7 +244,7 @@
1504 //! @endcode
1505 void setColorEquatorJ2000Grid(const Vec3f& newColor);
1506
1507- //! Setter for displaying Ecliptic Grid.
1508+ //! Setter for displaying Ecliptic Grid of J2000.0.
1509 void setFlagEclipticJ2000Grid(const bool displayed);
1510 //! Accessor for displaying Ecliptic Grid.
1511 bool getFlagEclipticJ2000Grid(void) const;
1512@@ -204,6 +258,20 @@
1513 //! @endcode
1514 void setColorEclipticJ2000Grid(const Vec3f& newColor);
1515
1516+ //! Setter for displaying Ecliptic Grid of Date.
1517+ void setFlagEclipticGrid(const bool displayed);
1518+ //! Accessor for displaying Ecliptic Grid.
1519+ bool getFlagEclipticGrid(void) const;
1520+ //! Get the current color of the Ecliptic of Date Grid.
1521+ Vec3f getColorEclipticGrid(void) const;
1522+ //! Set the color of the Ecliptic Grid.
1523+ //! @param newColor The color of Ecliptic of Date grid
1524+ //! @code
1525+ //! // example of usage in scripts
1526+ //! GridLinesMgr.setColorEclipticGrid(Vec3f(1.0,0.0,0.0));
1527+ //! @endcode
1528+ void setColorEclipticGrid(const Vec3f& newColor);
1529+
1530 //! Setter for displaying Galactic Grid.
1531 void setFlagGalacticGrid(const bool displayed);
1532 //! Accessor for displaying Galactic Grid.
1533@@ -232,6 +300,34 @@
1534 //! @endcode
1535 void setColorEquatorLine(const Vec3f& newColor);
1536
1537+ //! Setter for displaying J2000 Equatorial Line.
1538+ void setFlagEquatorJ2000Line(const bool displayed);
1539+ //! Accessor for displaying J2000 Equatorial Line.
1540+ bool getFlagEquatorJ2000Line(void) const;
1541+ //! Get the current color of the J2000 Equatorial Line.
1542+ Vec3f getColorEquatorJ2000Line(void) const;
1543+ //! Set the color of the J2000 Equator Line.
1544+ //! @param newColor The color of J2000 equator line
1545+ //! @code
1546+ //! // example of usage in scripts
1547+ //! GridLinesMgr.setColorEquatorJ2000Line(Vec3f(1.0,0.0,0.0));
1548+ //! @endcode
1549+ void setColorEquatorJ2000Line(const Vec3f& newColor);
1550+
1551+ //! Setter for displaying Ecliptic of J2000 Line.
1552+ void setFlagEclipticJ2000Line(const bool displayed);
1553+ //! Accessor for displaying Ecliptic of J2000 Line.
1554+ bool getFlagEclipticJ2000Line(void) const;
1555+ //! Get the current color of the Ecliptic of J2000 Line.
1556+ Vec3f getColorEclipticJ2000Line(void) const;
1557+ //! Set the color of the Ecliptic of J2000 Line.
1558+ //! @param newColor The color of ecliptic 2000 line
1559+ //! @code
1560+ //! // example of usage in scripts
1561+ //! GridLinesMgr.setColorEcliptic2000Line(Vec3f(1.0,0.0,0.0));
1562+ //! @endcode
1563+ void setColorEclipticJ2000Line(const Vec3f& newColor);
1564+
1565 //! Setter for displaying Ecliptic Line.
1566 void setFlagEclipticLine(const bool displayed);
1567 //! Accessor for displaying Ecliptic Line.
1568@@ -246,6 +342,20 @@
1569 //! @endcode
1570 void setColorEclipticLine(const Vec3f& newColor);
1571
1572+ //! Setter for displaying precession circles.
1573+ void setFlagPrecessionCircles(const bool displayed);
1574+ //! Accessor for displaying precession circles.
1575+ bool getFlagPrecessionCircles(void) const;
1576+ //! Get the current color of the precession circles.
1577+ Vec3f getColorPrecessionCircles(void) const;
1578+ //! Set the color of the precession circles.
1579+ //! @param newColor The color of precession circles
1580+ //! @code
1581+ //! // example of usage in scripts
1582+ //! GridLinesMgr.setColorPrecessionCircles(Vec3f(1.0,0.0,0.0));
1583+ //! @endcode
1584+ void setColorPrecessionCircles(const Vec3f& newColor);
1585+
1586 //! Setter for displaying Meridian Line.
1587 void setFlagMeridianLine(const bool displayed);
1588 //! Accessor for displaying Meridian Line.
1589@@ -291,12 +401,12 @@
1590 //! Setter for displaying Galactic Equator Line.
1591 void setFlagGalacticEquatorLine(const bool displayed);
1592 //! @deprecated Setter for displaying Galactic "Plane" (i.e., Equator) Line. Left here for compatibility with older scripts.
1593- //! @note will be delete in version 0.14
1594+ //! @note will be deleted in version 0.14
1595 void setFlagGalacticPlaneLine(const bool displayed) { setFlagGalacticEquatorLine(displayed); }
1596 //! Accessor for displaying Galactic Equator Line.
1597 bool getFlagGalacticEquatorLine(void) const;
1598 //! @deprecated Accessor for displaying Galactic "Plane" (i.e., Equator) Line. Left here for compatibility with older scripts.
1599- //! @note will be delete in version 0.14
1600+ //! @note will be deleted in version 0.14
1601 bool getFlagGalacticPlaneLine(void) const { return getFlagGalacticEquatorLine(); }
1602 //! Get the current color of the Galactic Equator Line.
1603 Vec3f getColorGalacticEquatorLine(void) const;
1604@@ -314,14 +424,22 @@
1605 void equatorGridColorChanged(const Vec3f & newColor) const;
1606 void equatorJ2000GridDisplayedChanged(const bool displayed) const;
1607 void equatorJ2000GridColorChanged(const Vec3f & newColor) const;
1608+ void eclipticGridDisplayedChanged(const bool displayed) const;
1609+ void eclipticGridColorChanged(const Vec3f & newColor) const;
1610 void eclipticJ2000GridDisplayedChanged(const bool displayed) const;
1611 void eclipticJ2000GridColorChanged(const Vec3f & newColor) const;
1612 void galacticGridDisplayedChanged(const bool displayed) const;
1613 void galacticGridColorChanged(const Vec3f & newColor) const;
1614 void equatorLineDisplayedChanged(const bool displayed) const;
1615 void equatorLineColorChanged(const Vec3f & newColor) const;
1616+ void equatorJ2000LineDisplayedChanged(const bool displayed) const;
1617+ void equatorJ2000LineColorChanged(const Vec3f & newColor) const;
1618 void eclipticLineDisplayedChanged(const bool displayed) const;
1619 void eclipticLineColorChanged(const Vec3f & newColor) const;
1620+ void eclipticJ2000LineDisplayedChanged(const bool displayed) const;
1621+ void eclipticJ2000LineColorChanged(const Vec3f & newColor) const;
1622+ void precessionCirclesDisplayedChanged(const bool displayed) const;
1623+ void precessionCirclesColorChanged(const Vec3f & newColor) const;
1624 void meridianLineDisplayedChanged(const bool displayed) const;
1625 void meridianLineColorChanged(const Vec3f & newColor) const;
1626 void longitudeLineDisplayedChanged(const bool displayed) const;
1627@@ -344,10 +462,15 @@
1628 SkyGrid * equGrid; // Equatorial grid
1629 SkyGrid * equJ2000Grid; // Equatorial J2000 grid
1630 SkyGrid * galacticGrid; // Galactic grid
1631+ SkyGrid * eclGrid; // Ecliptic of Date grid
1632 SkyGrid * eclJ2000Grid; // Ecliptic J2000 grid
1633 SkyGrid * aziGrid; // Azimuthal grid
1634 SkyLine * equatorLine; // Celestial Equator line
1635+ SkyLine * equatorJ2000Line; // Celestial Equator of J2000 line
1636 SkyLine * eclipticLine; // Ecliptic line
1637+ SkyLine * eclipticJ2000Line; // Ecliptic line of J2000
1638+ SkyLine * precessionCircleN; // Northern precession circle
1639+ SkyLine * precessionCircleS; // Southern precession circle
1640 SkyLine * meridianLine; // Meridian line
1641 SkyLine * longitudeLine; // Opposition/conjunction longitude line
1642 SkyLine * horizonLine; // Horizon line
1643
1644=== modified file 'src/core/modules/Landscape.cpp'
1645--- src/core/modules/Landscape.cpp 2015-06-30 12:08:33 +0000
1646+++ src/core/modules/Landscape.cpp 2015-07-28 19:18:20 +0000
1647@@ -803,6 +803,7 @@
1648 float y_baseImg_1 = sides[currentSide].texCoords[1]+ y_img_1*(sides[currentSide].texCoords[3]-sides[currentSide].texCoords[1]);
1649 int y=(1.0-y_baseImg_1)*sidesImages[currentSide]->height(); // pixel Y from top.
1650 QRgb pixVal=sidesImages[currentSide]->pixel(x, y);
1651+/*
1652 #ifndef NDEBUG
1653 // GZ: please leave the comment available for further development!
1654 qDebug() << "Oldstyle Landscape sampling: az=" << az*180.0 << "° alt=" << alt_rad*180.0f/M_PI
1655@@ -811,6 +812,7 @@
1656 << ", w=" << sidesImages[currentSide]->width() << " h=" << sidesImages[currentSide]->height()
1657 << " --> x:" << x << " y:" << y << " alpha:" << qAlpha(pixVal)/255.0f;
1658 #endif
1659+*/
1660 return qAlpha(pixVal)/255.0f;
1661 }
1662
1663@@ -1022,12 +1024,14 @@
1664 int y= mapImage->height()/2*(1 + radius*std::cos(az));
1665
1666 QRgb pixVal=mapImage->pixel(x, y);
1667+/*
1668 #ifndef NDEBUG
1669 // GZ: please leave the comment available for further development!
1670 qDebug() << "Landscape sampling: az=" << (az+angleRotateZ)/M_PI*180.0f << "° alt=" << alt_rad/M_PI*180.f
1671 << "°, w=" << mapImage->width() << " h=" << mapImage->height()
1672 << " --> x:" << x << " y:" << y << " alpha:" << qAlpha(pixVal)/255.0f;
1673 #endif
1674+*/
1675 return qAlpha(pixVal)/255.0f;
1676
1677
1678@@ -1210,6 +1214,7 @@
1679 int x=(az_phot/2.0f) * mapImage->width(); // pixel X from left.
1680
1681 QRgb pixVal=mapImage->pixel(x, y);
1682+/*
1683 #ifndef NDEBUG
1684 // GZ: please leave the comment available for further development!
1685 qDebug() << "Landscape sampling: az=" << az*180.0 << "° alt=" << alt_pm1*90.0f
1686@@ -1217,6 +1222,7 @@
1687 << ", w=" << mapImage->width() << " h=" << mapImage->height()
1688 << " --> x:" << x << " y:" << y << " alpha:" << qAlpha(pixVal)/255.0f;
1689 #endif
1690+*/
1691 return qAlpha(pixVal)/255.0f;
1692
1693 }
1694
1695=== modified file 'src/core/modules/Planet.cpp'
1696--- src/core/modules/Planet.cpp 2015-06-21 11:42:04 +0000
1697+++ src/core/modules/Planet.cpp 2015-07-28 19:18:20 +0000
1698@@ -24,7 +24,8 @@
1699 #include "StelSkyDrawer.hpp"
1700 #include "SolarSystem.hpp"
1701 #include "Planet.hpp"
1702-
1703+#include "planetsephems/precession.h"
1704+#include "StelObserver.hpp"
1705 #include "StelProjector.hpp"
1706 #include "sidereal_time.h"
1707 #include "StelTextureMgr.hpp"
1708@@ -225,15 +226,21 @@
1709
1710 oss << getPositionInfoString(core, flags);
1711
1712- if (flags&Extra)
1713- {
1714- static SolarSystem *ssystem=GETSTELMODULE(SolarSystem);
1715- double ecl= ssystem->getEarth()->getRotObliquity(core->getJDay());
1716- if (core->getCurrentLocation().planetName=="Earth")
1717- oss << q_("Obliquity (of date, for Earth): %1").arg(StelUtils::radToDmsStr(ecl, true)) << "<br>";
1718- //if (englishName!="Sun")
1719- // oss << q_("Obliquity (of date): %1").arg(StelUtils::radToDmsStr(getRotObliquity(core->getJDay()), true)) << "<br>";
1720- }
1721+ // GZ This is mostly for debugging. Maybe also useful for letting people use our results to cross-check theirs, but we should not act as reference, currently...
1722+ // TODO: maybe separate this out into:
1723+ //if (flags&EclipticCoordXYZ)
1724+ // For now: add to EclipticCoord
1725+ //if (flags&EclipticCoord)
1726+ // oss << q_("Ecliptical XYZ (VSOP87A): %1/%2/%3").arg(QString::number(eclipticPos[0], 'f', 3), QString::number(eclipticPos[1], 'f', 3), QString::number(eclipticPos[2], 'f', 3)) << "<br>";
1727+// if (flags&Extra)
1728+// {
1729+// static SolarSystem *ssystem=GETSTELMODULE(SolarSystem);
1730+// double ecl= ssystem->getEarth()->getRotObliquity(core->getJDay());
1731+// if (core->getCurrentLocation().planetName=="Earth")
1732+// oss << q_("Ecliptic obliquity (of date): %1").arg(StelUtils::radToDmsStr(ecl, true)) << "<br>";
1733+// //if (englishName!="Sun")
1734+// // oss << q_("Obliquity (of date): %1").arg(StelUtils::radToDmsStr(getRotObliquity(core->getJDay()), true)) << "<br>";
1735+// }
1736
1737 if (flags&Distance)
1738 {
1739@@ -277,6 +284,10 @@
1740 double siderealDay = getSiderealDay();
1741 if (flags&Extra)
1742 {
1743+ // This is a string you can activate for debugging. It shows the distance between observer and center of the body you are standing on.
1744+ // May be helpful for debugging critical parallax corrections for eclipses.
1745+ // For general use, find a better location first.
1746+ // oss << q_("Planetocentric distance &rho;: %1 (km)").arg(core->getCurrentObserver()->getDistanceFromCenter() * AU) <<"<br>";
1747 if (siderealPeriod>0)
1748 {
1749 // TRANSLATORS: Sidereal (orbital) period for solar system bodies in days and in Julian years (symbol: a)
1750@@ -291,7 +302,7 @@
1751 oss << q_("The period of rotation is chaotic") << "<br>";
1752 }
1753 }
1754- if (englishName.compare("Sun")!=0)
1755+ if (englishName != "Sun")
1756 {
1757 const Vec3d& observerHelioPos = core->getObserverHeliocentricEclipticPos();
1758 oss << QString(q_("Phase Angle: %1")).arg(StelUtils::radToDmsStr(getPhaseAngle(observerHelioPos))) << "<br>";
1759@@ -391,6 +402,10 @@
1760 }
1761 }
1762
1763+// return value in radians!
1764+// TODO: For earth, decide whether this should be Capitaine's omega_A (angle eclipticPoleJ2000//axis) or epsilon_A (angle axis//current ecliptic pole)
1765+// For now, it is epsilon_A, the angle between earth's rotational axis and mean ecliptic of date.
1766+// Details: e.g. Hilton etal, Report on Precession and the Ecliptic, Cel.Mech.Dyn.Astr.94:351-67 (2006), Fig1.
1767 double Planet::getRotObliquity(double JDay) const
1768 {
1769 // JDay=2451545.0 for J2000.0
1770@@ -575,16 +590,36 @@
1771
1772 }
1773
1774-// Compute the transformation matrix from the local Planet coordinate to the parent Planet coordinate
1775+// Compute the transformation matrix from the local Planet coordinate system to the parent Planet coordinate system.
1776+// In case of the planets, this makes the axis point to the respective celestial pole.
1777+// TODO: Verify for the other planets if their axes are relative to J2000 ecliptic (VSOP87A XY plane) or relative to (precessed) ecliptic of date?
1778 void Planet::computeTransMatrix(double jd)
1779 {
1780+ // TODO: correct this for earth with the new model.
1781 axisRotation = getSiderealTime(jd);
1782
1783- // Special case - heliocentric coordinates are on ecliptic,
1784+ // Special case - heliocentric coordinates are relative to eclipticJ2000 (VSOP87A XY plane),
1785 // not solar equator...
1786+
1787 if (parent)
1788 {
1789- rotLocalToParent = Mat4d::zrotation(re.ascendingNode - re.precessionRate*(jd-re.epoch)) * Mat4d::xrotation(re.obliquity);
1790+ // We can inject a proper precession plus even nutation matrix in this stage, if available.
1791+ if (englishName=="Earth")
1792+ {
1793+ // rotLocalToParent = Mat4d::zrotation(re.ascendingNode - re.precessionRate*(jd-re.epoch)) * Mat4d::xrotation(-getRotObliquity(jd));
1794+ // GZ: We follow Capitaine's (2003) formulation P=Rz(Chi_A)*Rx(-omega_A)*Rz(-psi_A)*Rx(eps_o).
1795+ // ADS: 2011A&A...534A..22V = A&A 534, A22 (2011): Vondrak, Capitane, Wallace: New Precession Expressions, valid for long time intervals:
1796+ // See also Hilton et al, Report on Precession and the Ecliptic. Cel.Mech.Dyn.Astr. 94:351-367 (2006).
1797+ double eps_A, chi_A, omega_A, psi_A;
1798+ getPrecessionAnglesVondrak(jd, &eps_A, &chi_A, &omega_A, &psi_A);
1799+ // GZ This is the right combination for precession of the equator: Nodal rotation psi_A,
1800+ // then rotation by omega_A, the angle between EclPoleJ2000 and EarthPoleOfDate.
1801+ // The final rotation by chi_A rotates the equinox (zero degree).
1802+ // To achieve ecliptical coords of date, you just have now to add a rotX by epsilon_A (obliquity of date).
1803+ rotLocalToParent= Mat4d::zrotation(-psi_A) * Mat4d::xrotation(-omega_A) * Mat4d::zrotation(chi_A);
1804+ }
1805+ else
1806+ rotLocalToParent = Mat4d::zrotation(re.ascendingNode - re.precessionRate*(jd-re.epoch)) * Mat4d::xrotation(re.obliquity);
1807 }
1808 }
1809
1810@@ -615,8 +650,17 @@
1811 double Planet::getSiderealTime(double jd) const
1812 {
1813 if (englishName=="Earth")
1814- {
1815- return get_apparent_sidereal_time(jd);
1816+ { // GZ I wanted to be sure that nutation is just those ignorable few arcseconds.
1817+ //qDebug() << "Difference apparent-mean sidereal times (s): " << (get_apparent_sidereal_time(jd)- get_mean_sidereal_time(jd))* 240.0; // 1degree=4min=240s.
1818+ // TODO: Reactivate Nutation (but following IAU-2000A) after fixing JD_UT/JD_ET issues, then change this call back to apparent_sidereal_time.
1819+ //return get_apparent_sidereal_time(jd);
1820+ return get_mean_sidereal_time(jd);
1821+
1822+ // In the newer precession/nutation literature (starting around 2006) there are two sets of algorithms:
1823+ // "Classical" sidereal time (Greenwich hour angle GHA) and a solution based on Earth Rotational Angle ERA.
1824+ // We keep using the classical set, but make sure we are correct!
1825+
1826+
1827 }
1828
1829 double t = jd - re.epoch;
1830@@ -1801,7 +1845,7 @@
1831 if (!re.siderealPeriod)
1832 return;
1833
1834- const StelProjectorP prj = core->getProjection(StelCore::FrameHeliocentricEcliptic);
1835+ const StelProjectorP prj = core->getProjection(StelCore::FrameHeliocentricEclipticJ2000);
1836
1837 StelPainter sPainter(prj);
1838
1839@@ -1812,7 +1856,7 @@
1840 sPainter.setColor(orbitColor[0], orbitColor[1], orbitColor[2], orbitFader.getInterstate());
1841 Vec3d onscreen;
1842 // special case - use current Planet position as center vertex so that draws
1843- // on it's orbit all the time (since segmented rather than smooth curve)
1844+ // on its orbit all the time (since segmented rather than smooth curve)
1845 Vec3d savePos = orbit[ORBIT_SEGMENTS/2];
1846 orbit[ORBIT_SEGMENTS/2]=getHeliocentricEclipticPos();
1847 orbit[ORBIT_SEGMENTS]=orbit[0];
1848
1849=== modified file 'src/core/modules/Planet.hpp'
1850--- src/core/modules/Planet.hpp 2015-04-17 21:30:03 +0000
1851+++ src/core/modules/Planet.hpp 2015-07-28 19:18:20 +0000
1852@@ -172,6 +172,8 @@
1853 //! Get the radius of the planet in AU.
1854 //! @return the radius of the planet in astronomical units.
1855 double getRadius(void) const {return radius;}
1856+ //! Get the value (1-f) for oblateness f.
1857+ double getOneMinusOblateness(void) const {return oneMinusOblateness;}
1858 //! Get duration of sidereal day
1859 double getSiderealDay(void) const {return re.period;}
1860 //! Get duration of sidereal year
1861@@ -219,6 +221,8 @@
1862 float _obliquity, float _ascendingNode,
1863 float _precessionRate, double _siderealPeriod);
1864 double getRotAscendingnode(void) const {return re.ascendingNode;}
1865+ // return angle between axis and normal of ecliptic plane (or, for a moon, equatorial/reference plane defined by parent).
1866+ // TODO: decide if this is always angle between axis and J2000 ecliptic, or should be axis//current ecliptic!
1867 double getRotObliquity(double JDay) const;
1868
1869 //! Get the Planet position in the parent Planet ecliptic coordinate in AU
1870@@ -321,8 +325,10 @@
1871 Vec3f color; // exclusively used for drawing the planet halo
1872
1873 float albedo; // Planet albedo. Used for magnitude computation (but formula dubious!)
1874- Mat4d rotLocalToParent;
1875- float axisRotation; // Rotation angle of the Planet on it's axis
1876+ Mat4d rotLocalToParent; // GZ2015: was undocumented.
1877+ // Apparently this is the axis orientation with respect to the parent body. For planets, this is axis orientation w.r.t. VSOP87A/J2000 ecliptical system.
1878+ float axisRotation; // Rotation angle of the Planet on its axis.
1879+ // For Earth, this should be Greenwich Mean Sidereal Time GMST.
1880 StelTextureSP texMap; // Planet map texture
1881 StelTextureSP normalMap; // Planet normal map texture
1882
1883@@ -333,7 +339,7 @@
1884 double lastJD; // caches JD of last positional computation
1885 // The callback for the calculation of the equatorial rect heliocentric position at time JD.
1886 posFuncType coordFunc;
1887- void* userDataPtr;
1888+ void* userDataPtr; // this is always used with an Orbit object.
1889
1890 OsculatingFunctType *const osculatingFunc;
1891 QSharedPointer<Planet> parent; // Planet parent i.e. sun for earth
1892
1893=== modified file 'src/core/modules/SolarSystem.cpp'
1894--- src/core/modules/SolarSystem.cpp 2015-07-05 19:48:33 +0000
1895+++ src/core/modules/SolarSystem.cpp 2015-07-28 19:18:20 +0000
1896@@ -1691,6 +1691,10 @@
1897 systemPlanets.clear();
1898 // Memory leak? What's the proper way of cleaning shared pointers?
1899
1900+ // Also delete Comet textures (loaded in loadPlanets()
1901+ Comet::tailTexture.clear();
1902+ Comet::comaTexture.clear();
1903+
1904 // Re-load the ssystem.ini file
1905 loadPlanets();
1906 computePositions(StelUtils::getJDFromSystem());
1907
1908=== modified file 'src/core/planetsephems/elp82b.h'
1909--- src/core/planetsephems/elp82b.h 2010-04-15 12:28:34 +0000
1910+++ src/core/planetsephems/elp82b.h 2015-07-28 19:18:20 +0000
1911@@ -82,7 +82,7 @@
1912
1913
1914 #ifdef __cplusplus
1915-};
1916+}
1917 #endif
1918
1919 #endif
1920
1921=== modified file 'src/core/planetsephems/gust86.c'
1922--- src/core/planetsephems/gust86.c 2010-04-15 12:28:34 +0000
1923+++ src/core/planetsephems/gust86.c 2015-07-28 19:18:20 +0000
1924@@ -438,7 +438,7 @@
1925 static double gust86_jd0 = -1e100;
1926 static double gust86_elem[GUST86_DIM];
1927
1928-void GetGust86Coor(double jd,int body,double *xyz) {
1929+void GetGust86Coor(const double jd,const int body,double *xyz) {
1930 GetGust86OsculatingCoor(jd,jd,body,xyz);
1931 }
1932
1933
1934=== modified file 'src/core/planetsephems/gust86.h'
1935--- src/core/planetsephems/gust86.h 2013-06-02 03:39:47 +0000
1936+++ src/core/planetsephems/gust86.h 2015-07-28 19:18:20 +0000
1937@@ -59,7 +59,7 @@
1938 #define GUST86_TITANIA 3
1939 #define GUST86_OBERON 4
1940
1941-void GetGust86Coor(double jd,int body,double *xyz);
1942+void GetGust86Coor(const double jd, const int body, double *xyz);
1943 /* Return the rectangular coordinates of the given satellite
1944 and the given julian date jd expressed in dynamical time (TAI+32.184s).
1945 The origin of the xyz-coordinates is the center of the planet.
1946@@ -88,8 +88,8 @@
1947 ICRF <-> VSOP87 must be done with the matrix given above.
1948 */
1949
1950-void GetGust86OsculatingCoor(const double jd0, const double jd, const int body,double *xyz);
1951- /* The oculating orbit of epoch jd0, evatuated at jd, is returned.
1952+void GetGust86OsculatingCoor(const double jd0, const double jd, const int body, double *xyz);
1953+ /* The oculating orbit of epoch jd0, evaluated at jd, is returned.
1954 */
1955
1956 #ifdef __cplusplus
1957
1958=== added file 'src/core/planetsephems/precession.c'
1959--- src/core/planetsephems/precession.c 1970-01-01 00:00:00 +0000
1960+++ src/core/planetsephems/precession.c 2015-07-28 19:18:20 +0000
1961@@ -0,0 +1,260 @@
1962+/*
1963+Copyright (C) 2015 Georg Zotti
1964+
1965+This program is free software; you can redistribute it and/or modify
1966+it under the terms of the GNU Library General Public License as published by
1967+the Free Software Foundation; either version 2 of the License, or
1968+(at your option) any later version.
1969+
1970+This program is distributed in the hope that it will be useful,
1971+but WITHOUT ANY WARRANTY; without even the implied warranty of
1972+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
1973+GNU General Public License for more details.
1974+
1975+You should have received a copy of the GNU General Public License
1976+along with this program; if not, write to the Free Software
1977+Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
1978+*/
1979+
1980+/*
1981+ * Precession solution following methods from:
1982+ * J. Vondrak, N. Capitaine, P. Wallace
1983+ * New precession expressions, valid for long time intervals
1984+ * Astronomy&Astrophysics 534, A22 (2011)
1985+ * DOI: 10.1051/0004-6361/201117274
1986+ * The data are only applicable for a time range of 200.000 years around J2000.
1987+ * This is by far enough for Stellarium as of 2015, but just to make sure I added a few asserts.
1988+ */
1989+
1990+#include <math.h>
1991+#include <assert.h>
1992+
1993+/* Interval threshold (days) for re-computing these values. with 1, compute only 1/day: */
1994+#define PRECESSION_EPOCH_THRESHOLD 1.0
1995+
1996+/* cache results for retrieval if recomputation is not required */
1997+
1998+static double c_psi_A=0.0, c_omega_A=0.0, c_chi_A=0.0, /*c_p_A=0.0, */ c_epsilon_A=0.0,
1999+ c_Y_A=0.0, c_X_A=0.0, c_Q_A=0.0, c_P_A=0.0,
2000+ c_lastJDE=-1e100;
2001+
2002+static const double arcSec2Rad=M_PI*2.0/(360.0*3600.0);
2003+
2004+static const double PQvals[8][5]=
2005+{ // 1/Pn P_A:Cn Q_A:Cn P_A:Sn Q_A:Sn
2006+ { 1.0/ 708.15, -5486.751211, -684.661560, 667.666730, -5523.863691 },
2007+ { 1.0/2309.00, -17.127623, 2446.283880, -2354.886252, -549.747450 },
2008+ { 1.0/1620.00, -617.517403, 399.671049, -428.152441, -310.998056 },
2009+ { 1.0/ 492.20, 413.442940, -356.652376, 376.202861, 421.535876 },
2010+ { 1.0/1183.00, 78.614193, -186.387003, 184.778874, -36.776172 },
2011+ { 1.0/ 622.00, -180.732815, -316.800070, 335.321713, -145.278396 },
2012+ { 1.0/ 882.00, -87.676083, 198.296071, -185.138669, -34.744450 },
2013+ { 1.0/ 547.00, 46.140315, 101.135679, -120.972830, 22.885731 }};
2014+
2015+static const double XYvals[14][5]=
2016+{ // 1/Pn Xa:Cn Ya:Cn Xa:Sn Ya:Sn
2017+ { 1.0/ 256.75, -819.940624, 75004.344875, 81491.287984, 1558.515853 },
2018+ { 1.0/ 708.15, -8444.676815, 624.033993, 787.163481, 7774.939698 },
2019+ { 1.0/ 274.20, 2600.009459, 1251.136893, 1251.296102, -2219.534038 },
2020+ { 1.0/ 241.45, 2755.175630, -1102.212834, -1257.950837, -2523.969396 },
2021+ { 1.0/2309.00, -167.659835, -2660.664980, -2966.799730, 247.850422 },
2022+ { 1.0/ 492.20, 871.855056, 699.291817, 639.744522, -846.485643 },
2023+ { 1.0/ 396.10, 44.769698, 153.167220, 131.600209, -1393.124055 },
2024+ { 1.0/ 288.90, -512.313065, -950.865637, -445.040117, 368.526116 },
2025+ { 1.0/ 231.10, -819.415595, 499.754645, 584.522874, 749.045012 },
2026+ { 1.0/1610.00, -538.071099, -145.188210, -89.756563, 444.704518 },
2027+ { 1.0/ 620.00, -189.793622, 558.116553, 524.429630, 235.934465 },
2028+ { 1.0/ 157.87, -402.922932, -23.923029, -13.549067, 374.049623 },
2029+ { 1.0/ 220.30, 179.516345, -165.405086, -210.157124, -171.330180 },
2030+ { 1.0/1200.00, -9.814756, 9.344131, -44.919798, -22.899655 }};
2031+
2032+
2033+static const double precVals[18][7]=
2034+{ // 1/Pn psi_A:Cn om_A:Cn chi_A:Cn psi_A:Sn om_A:Sn chi_A:Sn
2035+
2036+ { 1.0/402.90, -22206.325946, 1267.727824, -13765.924050, -3243.236469, -8571.476251, -2206.967126 },
2037+ { 1.0/256.75, 12236.649447, 1702.324248, 13511.858383, -3969.723769, 5309.796459, -4186.752711 },
2038+ { 1.0/292.00, -1589.008343, -2970.553839, -1455.229106, 7099.207893, -610.393953, 6737.949677 },
2039+ { 1.0/537.22, 2482.103195, 693.790312, 1054.394467, -1903.696711, 923.201931, -856.922846 },
2040+ { 1.0/241.45, 150.322920, -14.724451, 0.0 , 146.435014, 3.759055, 0.0 },
2041+ { 1.0/375.22, -13.632066, -516.649401, -112.300144, 1300.630106, -40.691114, 957.149088 },
2042+ { 1.0/157.87, 389.437420, -356.794454, 202.769908, 1727.498039, 80.437484, 1709.440735 },
2043+ { 1.0/274.20, 2031.433792, -129.552058, 1936.050095, 299.854055, 807.300668, 154.425505 },
2044+ { 1.0/203.00, 363.748303, 256.129314, 0.0 , -1217.125982, 83.712326, 0.0 },
2045+ { 1.0/440.00, -896.747562, 190.266114, -655.484214, -471.367487, -368.654854, -243.520976 },
2046+ { 1.0/170.72, -926.995700, 95.103991, -891.898637, -441.682145, -191.881064, -406.539008 },
2047+ { 1.0/713.37, 37.070667, -332.907067, 0.0 , -86.169171, -4.263770, 0.0 },
2048+ { 1.0/313.00, -597.682468, 131.337633, 0.0 , -308.320429, -270.353691, 0.0 },
2049+ { 1.0/128.38, 66.282812, 82.731919, -333.322021, -422.815629, 11.602861, -446.656435 },
2050+ { 1.0/202.00, 0.0 , 0.0 , 327.517465, 0.0 , 0.0 , -1049.071786 },
2051+ { 1.0/315.00, 0.0 , 0.0 , -494.780332, 0.0 , 0.0 , -301.504189 },
2052+ { 1.0/136.32, 0.0 , 0.0 , 585.492621, 0.0 , 0.0 , 41.348740 },
2053+ { 1.0/490.00, 0.0 , 0.0 , 110.512834, 0.0 , 0.0 , 142.525186 }};
2054+
2055+static const double p_epsVals[10][5]=
2056+{ // 1/Pn p_A:Cn eps_A:Cn p_A:Sn eps_A:Sn
2057+ { 1.0/ 409.90, -6908.287473, 753.872780, -2845.175469, -1704.720302},
2058+ { 1.0/ 396.15, -3198.706291, -247.805823, 449.844989, -862.308358},
2059+ { 1.0/ 537.22, 1453.674527, 379.471484, -1255.915323, 447.832178},
2060+ { 1.0/ 402.90, -857.748557, -53.880558, 886.736783, -889.571909},
2061+ { 1.0/ 417.15, 1173.231614, -90.109153, 418.887514, 190.402846},
2062+ { 1.0/ 288.92, -156.981465, -353.600190, 997.912441, -56.564991},
2063+ { 1.0/4043.00, 371.836550, -63.115353, -240.979710, -296.222622},
2064+ { 1.0/ 306.00, -216.619040, -28.248187, 76.541307, -75.859952},
2065+ { 1.0/ 277.00, 193.691479, 17.703387, -36.788069, 67.473503},
2066+ { 1.0/ 203.00, 11.891524, 38.911307, -170.964086, 3.014055}};
2067+
2068+// compute angles for the series we are in fact using.
2069+// jde: date JD_TT
2070+//
2071+void getPrecessionAnglesVondrak(const double jde, double *epsilon_A, double *chi_A, double *omega_A, double *psi_A)
2072+{
2073+ if (fabs(jde-c_lastJDE) > PRECESSION_EPOCH_THRESHOLD)
2074+ {
2075+ double T=(jde-2451545.0)* (1.0/36525.0); // Julian centuries from J2000.0
2076+ assert(fabs(T)<=2000); // MAKES SURE YOU NEVER OVERSTRETCH THIS!
2077+ double T2pi= T*(2.0*M_PI); // Julian centuries from J2000.0, premultiplied by 2Pi
2078+ // these are actually small greek letters in the papers.
2079+ double Psi_A=0.0;
2080+ double Omega_A=0.0;
2081+ double Chi_A=0.0;
2082+ double Epsilon_A=0.0;
2083+ //double p_A=0.0; // currently unused. The data don't disturb.
2084+ int i;
2085+ for (i=0; i<18; ++i)
2086+ {
2087+ double invP=precVals[i][0];
2088+ double sin2piT_P, cos2piT_P;
2089+#ifdef _GNU_SOURCE
2090+ sincos(T2pi*invP, &sin2piT_P, &cos2piT_P);
2091+#else
2092+ double phase=T2pi*invP;
2093+ sin2piT_P= sin(phase);
2094+ cos2piT_P= cos(phase);
2095+#endif
2096+ Psi_A += precVals[i][1]*cos2piT_P + precVals[i][4]*sin2piT_P;
2097+ Omega_A += precVals[i][2]*cos2piT_P + precVals[i][5]*sin2piT_P;
2098+ Chi_A += precVals[i][3]*cos2piT_P + precVals[i][6]*sin2piT_P;
2099+ }
2100+
2101+ for (i=0; i<10; ++i)
2102+ {
2103+ double invP=p_epsVals[i][0];
2104+ double sin2piT_P, cos2piT_P;
2105+#ifdef _GNU_SOURCE
2106+ sincos(T2pi*invP, &sin2piT_P, &cos2piT_P);
2107+#else
2108+ double phase=T2pi*invP;
2109+ sin2piT_P= sin(phase);
2110+ cos2piT_P= cos(phase);
2111+#endif
2112+ //p_A += p_epsVals[i][1]*cos2piT_P + p_epsVals[i][3]*sin2piT_P;
2113+ Epsilon_A += p_epsVals[i][2]*cos2piT_P + p_epsVals[i][4]*sin2piT_P;
2114+ }
2115+
2116+ Psi_A += (( 289.e-9*T - 0.00740913)*T + 5042.7980307)*T + 8473.343527;
2117+ Omega_A += (( 151.e-9*T + 0.00000146)*T - 0.4436568)*T + 84283.175915;
2118+ Chi_A += (( -61.e-9*T + 0.00001472)*T + 0.0790159)*T - 19.657270;
2119+ //p_A += ((271.e-9*T - 0.00710733)*T + 5043.0520035)*T + 8134.017132;
2120+ Epsilon_A += ((-110.e-9*T - 0.00004039)*T + 0.3624445)*T + 84028.206305;
2121+ c_psi_A = arcSec2Rad*Psi_A;
2122+ c_omega_A = arcSec2Rad*Omega_A;
2123+ c_chi_A = arcSec2Rad*Chi_A;
2124+ // c_p_A = arcSec2Rad*p_A;
2125+ c_epsilon_A = arcSec2Rad*Epsilon_A;
2126+ }
2127+ *psi_A = c_psi_A;
2128+ *omega_A = c_omega_A;
2129+ *chi_A = c_chi_A;
2130+ *epsilon_A = c_epsilon_A;
2131+}
2132+
2133+void getPrecessionAnglesVondrakPQXYe(const double jde, double *vP_A, double *vQ_A, double *vX_A, double *vY_A, double *vepsilon_A)
2134+{
2135+ if (fabs(jde-c_lastJDE) > PRECESSION_EPOCH_THRESHOLD)
2136+ {
2137+ double T=(jde-2451545.0)* (1.0/36525.0);
2138+ assert(fabs(T)<=2000); // MAKES SURE YOU NEVER OVERSTRETCH THIS!
2139+ double T2pi= T*(2.0*M_PI); // Julian centuries from J2000.0, premultiplied by 2Pi
2140+ // these are actually small greek letters in the papers.
2141+ double P_A=0.0;
2142+ double Q_A=0.0;
2143+ double X_A=0.0;
2144+ double Y_A=0.0;
2145+ double Epsilon_A=0.0;
2146+ int i;
2147+ for (i=0; i<8; ++i)
2148+ {
2149+ double invP=PQvals[i][0];
2150+ double sin2piT_P, cos2piT_P;
2151+#ifdef _GNU_SOURCE
2152+ sincos(T2pi*invP, &sin2piT_P, &cos2piT_P);
2153+#else
2154+ double phase=T2pi*invP;
2155+ sin2piT_P= sin(phase);
2156+ cos2piT_P= cos(phase);
2157+#endif
2158+ P_A += PQvals[i][1]*cos2piT_P + PQvals[i][3]*sin2piT_P;
2159+ Q_A += PQvals[i][2]*cos2piT_P + PQvals[i][4]*sin2piT_P;
2160+ }
2161+ for (i=0; i<14; ++i)
2162+ {
2163+ double invP=XYvals[i][0];
2164+ double sin2piT_P, cos2piT_P;
2165+#ifdef _GNU_SOURCE
2166+ sincos(T2pi*invP, &sin2piT_P, &cos2piT_P);
2167+#else
2168+ double phase=T2pi*invP;
2169+ sin2piT_P= sin(phase);
2170+ cos2piT_P= cos(phase);
2171+#endif
2172+ X_A += XYvals[i][1]*cos2piT_P + XYvals[i][3]*sin2piT_P;
2173+ Y_A += XYvals[i][2]*cos2piT_P + XYvals[i][4]*sin2piT_P;
2174+ }
2175+ for (i=0; i<10; ++i)
2176+ {
2177+ double invP=p_epsVals[i][0];
2178+ double sin2piT_P, cos2piT_P;
2179+#ifdef _GNU_SOURCE
2180+ sincos(T2pi*invP, &sin2piT_P, &cos2piT_P);
2181+#else
2182+ double phase=T2pi*invP;
2183+ sin2piT_P= sin(phase);
2184+ cos2piT_P= cos(phase);
2185+#endif
2186+ //p_A += p_epsVals[i][1]*cos2piT_P + p_epsVals[i][3]*sin2piT_P;
2187+ Epsilon_A += p_epsVals[i][2]*cos2piT_P + p_epsVals[i][4]*sin2piT_P;
2188+ }
2189+
2190+ // Now the polynomial terms in T. Horner's scheme is best again.
2191+ P_A += (( 110.e-9*T - 0.00028913)*T - 0.1189000)*T + 5851.607687;
2192+ Q_A += ((-437.e-9*T - 0.00000020)*T + 1.1689818)*T - 1600.886300;
2193+ X_A += ((-152.e-9*T - 0.00037173)*T + 0.4252841)*T + 5453.282155;
2194+ Y_A += ((+231.e-9*T - 0.00018725)*T - 0.7675452)*T - 73750.930350;
2195+ Epsilon_A += (( 110.e-9*T - 0.00004039)*T + 0.3624445)*T + 84028.206305;
2196+ c_P_A = arcSec2Rad*P_A;
2197+ c_Q_A = arcSec2Rad*Q_A;
2198+ c_X_A = arcSec2Rad*X_A;
2199+ c_Y_A = arcSec2Rad*Y_A;
2200+ c_epsilon_A = arcSec2Rad*Epsilon_A;
2201+ }
2202+ *vP_A = c_P_A;
2203+ *vQ_A = c_Q_A;
2204+ *vX_A = c_X_A;
2205+ *vY_A = c_Y_A;
2206+ *vepsilon_A = c_epsilon_A;
2207+
2208+}
2209+
2210+//! Just return (presumably precomputed) ecliptic obliquity.
2211+double getPrecessionAngleVondrakEpsilon(const double jde)
2212+{
2213+ double epsilon_A, dummy_chi_A, dummy_omega_A, dummy_psi_A;
2214+ getPrecessionAnglesVondrak(jde, &epsilon_A, &dummy_chi_A, &dummy_omega_A, &dummy_psi_A);
2215+ return epsilon_A;
2216+}
2217+//! Just return (presumably precomputed) ecliptic obliquity.
2218+double getPrecessionAngleVondrakCurrentEpsilonA(void)
2219+{
2220+ return c_epsilon_A;
2221+}
2222
2223=== added file 'src/core/planetsephems/precession.h'
2224--- src/core/planetsephems/precession.h 1970-01-01 00:00:00 +0000
2225+++ src/core/planetsephems/precession.h 2015-07-28 19:18:20 +0000
2226@@ -0,0 +1,59 @@
2227+/*
2228+Copyright (C) 2015 Georg Zotti
2229+
2230+This program is free software; you can redistribute it and/or modify
2231+it under the terms of the GNU Library General Public License as published by
2232+the Free Software Foundation; either version 2 of the License, or
2233+(at your option) any later version.
2234+
2235+This program is distributed in the hope that it will be useful,
2236+but WITHOUT ANY WARRANTY; without even the implied warranty of
2237+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2238+GNU General Public License for more details.
2239+
2240+You should have received a copy of the GNU General Public License
2241+along with this program; if not, write to the Free Software
2242+Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
2243+*/
2244+#ifndef _PRECESSION_H_
2245+#define _PRECESSION_H_
2246+
2247+#ifdef __cplusplus
2248+extern "C" {
2249+#endif
2250+
2251+//! Precession modelled from:
2252+//! J. Vondrák, N. Capitaine, and P. Wallace: New precession expressions, valid for long time intervals
2253+//! A&A (Astronomy&Astrophysics) 534, A22 (2011)
2254+//! DOI: http://dx.doi.org/10.1051/0004-6361/201117274
2255+//! This paper describes a precession model valid for +/-200.000 years from J2000.0 and consistent with P03 precession accepted as IAU2006 Precession.
2256+//! Some better understanding of the angles can be found in:
2257+//! + 1994AJ____108__711W J.G.Williams: Contributions to the Earth's Obliquity Rate, Precession and Nutation (Angles of eq. (35))
2258+//! + A&A 459, 981–985 (2006) DOI: 10.1051/0004-6361:20065897: Wallace&Capitaine: Precession-nutation procedures consistent with IAU 2006 resolutions
2259+//!
2260+//! The angles computed therein are used to rotate the planet Earth's axis, and also to rotate an "Ecliptic of Date", i.e. the current orbital plane of Earth.
2261+//! Currently this is without Nutation.
2262+//! Return values are in radians
2263+void getPrecessionAnglesVondrak(const double jde, double *epsilon_A, double *chi_A, double *omega_A, double *psi_A);
2264+
2265+//! Alternative solution, the one also implemented in the paper,
2266+//! combining matrix P from P_A, Q_A, X_A, Y_A and, for the ecliptic of date, rotate back by epsilon_A.
2267+//! Return values are in radians.
2268+//! This solution is currently unused, it seems easier to use the Capitaine sequence above.
2269+void getPrecessionAnglesVondrakPQXYe(const double jde, double *vP_A, double *vQ_A, double *vX_A, double *vY_A, double *vepsilon_A);
2270+
2271+//! Return ecliptic obliquity. [radians]
2272+double getPrecessionAngleVondrakEpsilon(const double jde);
2273+
2274+//! Just return (previously computed) ecliptic obliquity. [radians]
2275+double getPrecessionAngleVondrakCurrentEpsilonA(void);
2276+
2277+//! TODO: To complete the task of correct&accurate precession-nutation handling, find IAU-2000A Nutation and how this fits in here.
2278+//! E.g. A&A 459, 981–985 (2006) P. T. Wallace and N. Capitaine: Precession-nutation procedures consistent with IAU 2006 resolutions. DOI: 10.1051/0004-6361:20065897
2279+
2280+
2281+#ifdef __cplusplus
2282+}
2283+#endif
2284+
2285+#endif
2286
2287=== modified file 'src/core/planetsephems/sidereal_time.c'
2288--- src/core/planetsephems/sidereal_time.c 2014-02-05 16:51:48 +0000
2289+++ src/core/planetsephems/sidereal_time.c 2015-07-28 19:18:20 +0000
2290@@ -18,6 +18,8 @@
2291 */
2292
2293 #include <math.h>
2294+#include <assert.h>
2295+#include "precession.h"
2296
2297 #ifndef M_PI
2298 #define M_PI 3.14159265358979323846
2299@@ -44,7 +46,7 @@
2300 #define LN_NUTATION_EPOCH_THRESHOLD 0.1
2301
2302
2303-/* Nutation is a period oscillation of the Earths rotational axis around it's mean position.*/
2304+/* Nutation is a periodic oscillation of the Earth's rotational axis around its mean position.*/
2305
2306 /*
2307 Contains Nutation in longitude, obliquity and ecliptic obliquity.
2308@@ -52,9 +54,9 @@
2309 */
2310 struct ln_nutation
2311 {
2312- double longitude; /*!< Nutation in longitude */
2313- double obliquity; /*!< Nutation in obliquity */
2314- double ecliptic; /*!< Obliquity of the ecliptic */
2315+ double deltaPsi; /*!< DeltaPsi, Nutation in longitude */
2316+ double deltaEps; /*!< DeltaEps, Nutation in obliquity */
2317+ double ecliptic; /*!< epsilon, Mean Obliquity of the ecliptic */
2318 };
2319
2320 struct nutation_arguments
2321@@ -215,6 +217,8 @@
2322 * Meeus, Astr. Alg. (1st ed., 1994), Chapter 21 pg 131-134 Using Table 21A */
2323 /* GZ: Changed: ecliptic obliquity used to be constant J2000.0.
2324 * If you don't compute this, you may as well forget about nutation!
2325+ * 2015: Laskar's 1986 formula replaced by Vondrak 2011. c_ecliptic is now epsilon_A, Vondrak's obliquity of date.
2326+ * TODO: replace this whole function with Nutation IAU2000A or compatible version.
2327 */
2328 void get_nutation (double JD, struct ln_nutation * nutation)
2329 {
2330@@ -245,7 +249,7 @@
2331 F = 93.2719100 + 483202.017538 * T - 0.0036825 * T2 + T3 / 327270.0;
2332 O = 125.04452 - 1934.136261 * T + 0.0020708 * T2 + T3 / 450000.0;
2333
2334- / * GZotti: Laskar's formula, only valid for J2000+/-10000 years! */
2335+ / * GZotti: Laskar's formula (1986), only valid for J2000+/-10000 years! */
2336 if (fabs(T)<100)
2337 {
2338 double U=T/100.0;
2339@@ -322,25 +326,32 @@
2340 }
2341
2342 /* return results */
2343- nutation->longitude = c_longitude;
2344- nutation->obliquity = c_obliquity;
2345+ nutation->deltaPsi = c_longitude;
2346+ nutation->deltaEps = c_obliquity;
2347 nutation->ecliptic = c_ecliptic;
2348 }
2349
2350-/* Calculate the mean sidereal time at the meridian of Greenwich of a given date.
2351- * returns apparent sidereal time (degree).
2352- * Formula 11.1, 11.4 pg 83 */
2353+/* Calculate the mean sidereal time at the meridian of Greenwich (GMST) of a given date.
2354+ * returns mean sidereal time (degrees).
2355+ * Formula 11.1, 11.4 pg 83
2356+ * MAKE SURE argument JD is UT, not TT!
2357+ */
2358 double get_mean_sidereal_time (double JD)
2359 {
2360 double sidereal;
2361- double T;
2362+ double T, T1;
2363
2364- T = (JD - 2451545.0) / 36525.0;
2365+
2366+ T1 = (JD - 2451545.0);
2367+ T= T1 * (1.0/ 36525.0);
2368
2369 /* calc mean angle */
2370- sidereal = 280.46061837 + (360.98564736629 * (JD - 2451545.0)) + (0.000387933 * T * T) - (T * T * T / 38710000.0);
2371+ //sidereal = 280.46061837 + (360.98564736629 * (JD - 2451545.0)) + (0.000387933 * T * T) - (T * T * T / 38710000.0);
2372
2373- /* add a convenient multiple of 360 degrees */
2374+ sidereal = range_degrees(T1*360.98564736629); // Number gets large, better bring that to interval soon.
2375+ sidereal += (-1.0/38710000.0 * T + 0.000387933)*T*T+280.46061837;
2376+
2377+ /* add again a convenient multiple of 360 degrees */
2378 sidereal = range_degrees (sidereal);
2379
2380 return sidereal;
2381@@ -348,13 +359,15 @@
2382
2383
2384 /* Calculate the apparent sidereal time at the meridian of Greenwich of a given date.
2385- * returns apparent sidereal time (degree).
2386+ * returns apparent sidereal time (degrees).
2387 * Formula 11.1, 11.4 pg 83 */
2388 double get_apparent_sidereal_time (double JD)
2389 {
2390 double correction, sidereal;
2391 struct ln_nutation nutation;
2392
2393+ // GZ For now this is forbidden!
2394+ //assert(0);
2395 /* get the mean sidereal time */
2396 sidereal = get_mean_sidereal_time (JD);
2397
2398@@ -363,23 +376,25 @@
2399 get_nutation (JD, &nutation);
2400
2401 /* GZ: This was the only place where this was used. I added the summation here. */
2402- correction = (nutation.longitude * cos ((nutation.ecliptic+nutation.obliquity)*M_PI/180.));
2403+ correction = (nutation.deltaPsi * cos ((nutation.ecliptic+nutation.deltaEps)*M_PI/180.));
2404
2405 sidereal += correction;
2406
2407 return (sidereal);
2408 }
2409
2410+// return value in degrees
2411 double get_mean_ecliptical_obliquity(double JDE)
2412 {
2413- struct ln_nutation nutation;
2414- get_nutation(JDE, &nutation);
2415- return nutation.ecliptic;
2416+// struct ln_nutation nutation;
2417+// get_nutation(JDE, &nutation);
2418+// return nutation.ecliptic;
2419+ return getPrecessionAngleVondrakEpsilon(JDE)*180.0/M_PI;
2420 }
2421
2422 double get_nutation_longitude(double JDE)
2423 {
2424 struct ln_nutation nutation;
2425 get_nutation(JDE, &nutation);
2426- return nutation.longitude;
2427+ return nutation.deltaPsi;
2428 }
2429
2430=== modified file 'src/core/planetsephems/vsop87.h'
2431--- src/core/planetsephems/vsop87.h 2013-06-02 03:39:47 +0000
2432+++ src/core/planetsephems/vsop87.h 2015-07-28 19:18:20 +0000
2433@@ -69,7 +69,7 @@
2434 */
2435
2436 void GetVsop87OsculatingCoor(const double jd0,const double jd, const int body,double *xyz);
2437- /* The oculating orbit of epoch jd0, evatuated at jd, is returned.
2438+ /* The oculating orbit of epoch jd0, evaluated at jd, is returned.
2439 */
2440
2441 #ifdef __cplusplus
2442
2443=== modified file 'src/gui/ConfigurationDialog.cpp'
2444--- src/gui/ConfigurationDialog.cpp 2015-07-08 13:31:50 +0000
2445+++ src/gui/ConfigurationDialog.cpp 2015-07-28 19:18:20 +0000
2446@@ -563,13 +563,16 @@
2447 // view dialog / markings tab settings
2448 conf->setValue("viewing/flag_azimuthal_grid", glmgr->getFlagAzimuthalGrid());
2449 conf->setValue("viewing/flag_equatorial_grid", glmgr->getFlagEquatorGrid());
2450+ conf->setValue("viewing/flag_equatorial_J2000_grid", glmgr->getFlagEquatorJ2000Grid());
2451 conf->setValue("viewing/flag_equator_line", glmgr->getFlagEquatorLine());
2452+ conf->setValue("viewing/flag_equator_J2000_line", glmgr->getFlagEquatorJ2000Line());
2453 conf->setValue("viewing/flag_ecliptic_line", glmgr->getFlagEclipticLine());
2454+ conf->setValue("viewing/flag_ecliptic_J2000_line", glmgr->getFlagEclipticJ2000Line());
2455+ conf->setValue("viewing/flag_ecliptic_grid", glmgr->getFlagEclipticGrid());
2456 conf->setValue("viewing/flag_ecliptic_J2000_grid", glmgr->getFlagEclipticJ2000Grid());
2457 conf->setValue("viewing/flag_meridian_line", glmgr->getFlagMeridianLine());
2458 conf->setValue("viewing/flag_longitude_line", glmgr->getFlagLongitudeLine());
2459 conf->setValue("viewing/flag_horizon_line", glmgr->getFlagHorizonLine());
2460- conf->setValue("viewing/flag_equatorial_J2000_grid", glmgr->getFlagEquatorJ2000Grid());
2461 conf->setValue("viewing/flag_galactic_grid", glmgr->getFlagGalacticGrid());
2462 conf->setValue("viewing/flag_galactic_equator_line", glmgr->getFlagGalacticEquatorLine());
2463 conf->setValue("viewing/flag_cardinal_points", lmgr->getFlagCardinalsPoints());
2464
2465=== modified file 'src/gui/LocationDialog.cpp'
2466--- src/gui/LocationDialog.cpp 2014-11-25 13:02:40 +0000
2467+++ src/gui/LocationDialog.cpp 2015-07-28 19:18:20 +0000
2468@@ -340,7 +340,7 @@
2469 else
2470 loc.planetName = ui->planetNameComboBox->itemData(index).toString();
2471 loc.name = ui->cityNameLineEdit->text();
2472- loc.latitude = ui->latitudeSpinBox->valueDegrees();
2473+ loc.latitude = qMin(90.0, qMax(-90.0, ui->latitudeSpinBox->valueDegrees()));
2474 loc.longitude = ui->longitudeSpinBox->valueDegrees();
2475 loc.altitude = ui->altitudeSpinBox->value();
2476 index = ui->countryNameComboBox->currentIndex();
2477
2478=== modified file 'src/gui/ViewDialog.cpp'
2479--- src/gui/ViewDialog.cpp 2015-07-10 09:38:01 +0000
2480+++ src/gui/ViewDialog.cpp 2015-07-28 19:18:20 +0000
2481@@ -271,7 +271,9 @@
2482
2483 // Grid and lines
2484 connectCheckBox(ui->showEquatorLineCheckBox, "actionShow_Equator_Line");
2485- connectCheckBox(ui->showEclipticLineCheckBox, "actionShow_Ecliptic_Line");
2486+ connectCheckBox(ui->showEquatorJ2000LineCheckBox, "actionShow_Equator_J2000_Line");
2487+ connectCheckBox(ui->showEclipticLineJ2000CheckBox, "actionShow_Ecliptic_J2000_Line");
2488+ connectCheckBox(ui->showEclipticLineOfDateCheckBox, "actionShow_Ecliptic_Line");
2489 connectCheckBox(ui->showMeridianLineCheckBox, "actionShow_Meridian_Line");
2490 connectCheckBox(ui->showLongitudeLineCheckBox, "actionShow_Longitude_Line");
2491 connectCheckBox(ui->showHorizonLineCheckBox, "actionShow_Horizon_Line");
2492@@ -281,7 +283,9 @@
2493 connectCheckBox(ui->showAzimuthalGridCheckBox, "actionShow_Azimuthal_Grid");
2494 connectCheckBox(ui->showEquatorialJ2000GridCheckBox, "actionShow_Equatorial_J2000_Grid");
2495 connectCheckBox(ui->showEclipticGridJ2000CheckBox, "actionShow_Ecliptic_J2000_Grid");
2496+ connectCheckBox(ui->showEclipticGridOfDateCheckBox, "actionShow_Ecliptic_Grid");
2497 connectCheckBox(ui->showCardinalPointsCheckBox, "actionShow_Cardinal_Points");
2498+ connectCheckBox(ui->showPrecessionCirclesCheckBox, "actionShow_Precession_Circles");
2499
2500 // Constellations
2501 ConstellationMgr* cmgr = GETSTELMODULE(ConstellationMgr);
2502
2503=== modified file 'src/gui/viewDialog.ui'
2504--- src/gui/viewDialog.ui 2015-07-10 09:38:01 +0000
2505+++ src/gui/viewDialog.ui 2015-07-28 19:18:20 +0000
2506@@ -1453,27 +1453,49 @@
2507 <layout class="QVBoxLayout" name="verticalLayout_13">
2508 <item>
2509 <widget class="QCheckBox" name="showEquatorialJ2000GridCheckBox">
2510+ <property name="toolTip">
2511+ <string>Equatorial coordinates of J2000.0.</string>
2512+ </property>
2513 <property name="text">
2514- <string>Equatorial grid (on J2000)</string>
2515+ <string>Equatorial grid (J2000)</string>
2516 </property>
2517 </widget>
2518 </item>
2519 <item>
2520 <widget class="QCheckBox" name="showEquatorialGridCheckBox">
2521+ <property name="toolTip">
2522+ <string>Equatorial coordinates of current date and planet.</string>
2523+ </property>
2524 <property name="text">
2525- <string>Equatorial grid (on date)</string>
2526+ <string>Equatorial grid (of date)</string>
2527 </property>
2528 </widget>
2529 </item>
2530 <item>
2531 <widget class="QCheckBox" name="showEclipticGridJ2000CheckBox">
2532- <property name="text">
2533- <string>Ecliptic grid (on J2000)</string>
2534+ <property name="toolTip">
2535+ <string>Ecliptical coordinates for J2000.0.</string>
2536+ </property>
2537+ <property name="text">
2538+ <string>Ecliptic grid (J2000)</string>
2539+ </property>
2540+ </widget>
2541+ </item>
2542+ <item>
2543+ <widget class="QCheckBox" name="showEclipticGridOfDateCheckBox">
2544+ <property name="toolTip">
2545+ <string>Ecliptical coordinates for current date. Displayed on Earth only.</string>
2546+ </property>
2547+ <property name="text">
2548+ <string>Ecliptic grid (of date)</string>
2549 </property>
2550 </widget>
2551 </item>
2552 <item>
2553 <widget class="QCheckBox" name="showAzimuthalGridCheckBox">
2554+ <property name="toolTip">
2555+ <string>Altitudes and azimuth (counted from North towards East).</string>
2556+ </property>
2557 <property name="text">
2558 <string>Azimuthal grid</string>
2559 </property>
2560@@ -1481,6 +1503,9 @@
2561 </item>
2562 <item>
2563 <widget class="QCheckBox" name="showGalacticGridCheckBox">
2564+ <property name="toolTip">
2565+ <string>Galactic Coordinates, System II (IAU 1958).</string>
2566+ </property>
2567 <property name="text">
2568 <string>Galactic grid</string>
2569 </property>
2570@@ -1488,49 +1513,72 @@
2571 </item>
2572 <item>
2573 <widget class="QCheckBox" name="showCardinalPointsCheckBox">
2574+ <property name="toolTip">
2575+ <string>Labels for the Cardinal directions.</string>
2576+ </property>
2577 <property name="text">
2578 <string>Cardinal points</string>
2579 </property>
2580 </widget>
2581 </item>
2582+ <item>
2583+ <widget class="QCheckBox" name="showPrecessionCirclesCheckBox">
2584+ <property name="toolTip">
2585+ <string>Instantaneous circles of earth's axis on its motion around ecliptical poles. Displayed on Earth only.</string>
2586+ </property>
2587+ <property name="text">
2588+ <string>Precession circles</string>
2589+ </property>
2590+ </widget>
2591+ </item>
2592 </layout>
2593 </item>
2594 <item>
2595 <layout class="QVBoxLayout" name="verticalLayout_15">
2596 <item>
2597+ <widget class="QCheckBox" name="showEquatorJ2000LineCheckBox">
2598+ <property name="toolTip">
2599+ <string>Show celestial equator of J2000.0.</string>
2600+ </property>
2601+ <property name="text">
2602+ <string>Equator (J2000)</string>
2603+ </property>
2604+ </widget>
2605+ </item>
2606+ <item>
2607 <widget class="QCheckBox" name="showEquatorLineCheckBox">
2608 <property name="toolTip">
2609- <string>Show equator line</string>
2610- </property>
2611- <property name="text">
2612- <string>Equator</string>
2613- </property>
2614- </widget>
2615- </item>
2616- <item>
2617- <widget class="QCheckBox" name="showMeridianLineCheckBox">
2618- <property name="toolTip">
2619- <string>Show meridian line</string>
2620- </property>
2621- <property name="text">
2622- <string>Meridian</string>
2623- </property>
2624- </widget>
2625- </item>
2626- <item>
2627- <widget class="QCheckBox" name="showEclipticLineCheckBox">
2628- <property name="toolTip">
2629- <string>Show ecliptic line</string>
2630- </property>
2631- <property name="text">
2632- <string>Ecliptic</string>
2633+ <string>Show celestial equator of current planet and date.</string>
2634+ </property>
2635+ <property name="text">
2636+ <string>Equator (of date)</string>
2637+ </property>
2638+ </widget>
2639+ </item>
2640+ <item>
2641+ <widget class="QCheckBox" name="showEclipticLineJ2000CheckBox">
2642+ <property name="toolTip">
2643+ <string>Show ecliptic line of J2000.0 (VSOP87A fundamental plane).</string>
2644+ </property>
2645+ <property name="text">
2646+ <string>Ecliptic (J2000)</string>
2647+ </property>
2648+ </widget>
2649+ </item>
2650+ <item>
2651+ <widget class="QCheckBox" name="showEclipticLineOfDateCheckBox">
2652+ <property name="toolTip">
2653+ <string>Show ecliptic line of current date.</string>
2654+ </property>
2655+ <property name="text">
2656+ <string>Ecliptic (of date)</string>
2657 </property>
2658 </widget>
2659 </item>
2660 <item>
2661 <widget class="QCheckBox" name="showHorizonLineCheckBox">
2662 <property name="toolTip">
2663- <string>Show horizon line</string>
2664+ <string>Show horizon line.</string>
2665 </property>
2666 <property name="text">
2667 <string>Horizon</string>
2668@@ -1540,7 +1588,7 @@
2669 <item>
2670 <widget class="QCheckBox" name="showGalacticEquatorLineCheckBox">
2671 <property name="toolTip">
2672- <string>Show Galactic equator line</string>
2673+ <string>Show Galactic equator line.</string>
2674 </property>
2675 <property name="text">
2676 <string>Galactic equator</string>
2677@@ -1548,6 +1596,16 @@
2678 </widget>
2679 </item>
2680 <item>
2681+ <widget class="QCheckBox" name="showMeridianLineCheckBox">
2682+ <property name="toolTip">
2683+ <string>Show meridian line.</string>
2684+ </property>
2685+ <property name="text">
2686+ <string>Meridian</string>
2687+ </property>
2688+ </widget>
2689+ </item>
2690+ <item>
2691 <widget class="QCheckBox" name="showLongitudeLineCheckBox">
2692 <property name="toolTip">
2693 <string>Opposition/conjunction longitude line - the line of ecliptic longitude which passes through both ecliptic poles, the Sun and opposition point.</string>
2694@@ -2219,8 +2277,7 @@
2695 <tabstop>showAzimuthalGridCheckBox</tabstop>
2696 <tabstop>showGalacticGridCheckBox</tabstop>
2697 <tabstop>showEquatorLineCheckBox</tabstop>
2698- <tabstop>showMeridianLineCheckBox</tabstop>
2699- <tabstop>showEclipticLineCheckBox</tabstop>
2700+ <tabstop>showEclipticLineOfDateCheckBox</tabstop>
2701 <tabstop>showHorizonLineCheckBox</tabstop>
2702 <tabstop>showGalacticEquatorLineCheckBox</tabstop>
2703 <tabstop>showConstellationLinesCheckBox</tabstop>
2704
2705=== added file 'src/tests/testPrecession.cpp'
2706--- src/tests/testPrecession.cpp 1970-01-01 00:00:00 +0000
2707+++ src/tests/testPrecession.cpp 2015-07-28 19:18:20 +0000
2708@@ -0,0 +1,121 @@
2709+/*
2710+ * Stellarium
2711+ * Copyright (C) 2015 Georg Zotti
2712+ *
2713+ * This program is free software; you can redistribute it and/or
2714+ * modify it under the terms of the GNU General Public License
2715+ * as published by the Free Software Foundation; either version 2
2716+ * of the License, or (at your option) any later version.
2717+ *
2718+ * This program is distributed in the hope that it will be useful,
2719+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
2720+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2721+ * GNU General Public License for more details.
2722+ *
2723+ * You should have received a copy of the GNU General Public License
2724+ * along with this program; if not, write to the Free Software
2725+ * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
2726+ */
2727+
2728+#include <QObject>
2729+#include <QtDebug>
2730+#include <QVariantList>
2731+#include <QtTest>
2732+
2733+#include "tests/testPrecession.hpp"
2734+#include "StelUtils.hpp"
2735+#include "VecMath.hpp"
2736+
2737+QTEST_MAIN(TestPrecession)
2738+
2739+static const double arcSec2Rad=M_PI*2.0/(360.0*3600.0);
2740+static const double eps0=84381.406*arcSec2Rad;
2741+
2742+void TestPrecession::initTestCase()
2743+{
2744+}
2745+
2746+// example data found in Vondrak et al. 2011: New Precession Expressions, Valid for Long Time Intervals. A&A534, A22.
2747+// Section A.5
2748+// pecl = ( +0.00041724785764001342 -0.40495491104576162693 +0.91433656053126552350 )
2749+// = ( P_A -Q_A*C-Z*S -Q_A*S+Z*C )
2750+// pequ = ( -0.29437643797369031532 -0.11719098023370257855 +0.94847708824082091796 )
2751+// = ( X Y sqrt(1-x^2-y^2) || 0 )
2752+// This provides only the axes formulation. But the final Precession matrix should then be used as reference
2753+// comparable to the matrix that can be reached via the Capitaine parametrisation that we are using.
2754+
2755+void TestPrecession::testPrecessionAnglesVondrak()
2756+{
2757+ const double S=sin(eps0);
2758+ const double C=cos(eps0);
2759+ double Z, W;
2760+ double JulianDay=1219339.078000; // TT
2761+ double epsilon_A, chi_A, omega_A, psi_A;
2762+ double P_A, Q_A, X_A, Y_A;
2763+ double VEC2, VEC3, VEQ3;
2764+
2765+ // get angles for Capitaine parameterisation
2766+ getPrecessionAnglesVondrak(JulianDay, &epsilon_A, &chi_A, &omega_A, &psi_A);
2767+ // Get reference angles. We have to call this twice with different dates to avoid returning the cached (zero) angles.
2768+ getPrecessionAnglesVondrakPQXYe(-1234.567, &P_A, &Q_A, &X_A, &Y_A, &epsilon_A);
2769+ getPrecessionAnglesVondrakPQXYe(JulianDay, &P_A, &Q_A, &X_A, &Y_A, &epsilon_A);
2770+ Z=sqrt(qMax(1.0-P_A*P_A-Q_A*Q_A, 0.0));
2771+ W=X_A*X_A+Y_A*Y_A;
2772+ VEC2=-Q_A*C -Z*S;
2773+ VEC3=-Q_A*S +Z*C;
2774+ VEQ3=(W<1.0 ? sqrt(1.0-W) : 0.0);
2775+
2776+ QVERIFY2(fabs(P_A -0.00041724785764001342)<=0.000001, QString("JD %1: Pecl,x: %2 Difference: %3").arg(JulianDay).arg(P_A ).arg(P_A -0.00041724785764001342).toUtf8());
2777+ QVERIFY2(fabs(VEC2+0.40495491104576162693)<=0.000001, QString("JD %1: Pecl,y: %2 Difference: %3").arg(JulianDay).arg(VEC2).arg(VEC2+0.40495491104576162693).toUtf8());
2778+ QVERIFY2(fabs(VEC3-0.91433656053126552350)<=0.000001, QString("JD %1: Pecl,z: %2 Difference: %3").arg(JulianDay).arg(VEC3).arg(VEC3-0.91433656053126552350).toUtf8());
2779+ QVERIFY2(fabs(X_A +0.29437643797369031532)<=0.000001, QString("JD %1: Pequ,x: %2 Difference: %3").arg(JulianDay).arg(X_A ).arg(X_A +0.29437643797369031532).toUtf8());
2780+ QVERIFY2(fabs(Y_A +0.11719098023370257855)<=0.000001, QString("JD %1: Pequ,y: %2 Difference: %3").arg(JulianDay).arg(Y_A ).arg(Y_A +0.11719098023370257855).toUtf8());
2781+ QVERIFY2(fabs(VEQ3-0.94847708824082091796)<=0.000001, QString("JD %1: Pequ,z: %2 Difference: %3").arg(JulianDay).arg(VEQ3).arg(VEQ3-0.94847708824082091796).toUtf8());
2782+ // the same, to be seen ...
2783+// qDebug() << QString("JD %1: Pecl,x: %2 Difference: %3").arg(JulianDay, 8, 'f', 5).arg(P_A , 15, 'f', 12).arg(P_A -0.00041724785764001342, 15, 'f', 12);
2784+// qDebug() << QString("JD %1: Pecl,y: %2 Difference: %3").arg(JulianDay, 8, 'f', 5).arg(VEC2, 15, 'f', 12).arg(VEC2+0.40495491104576162693, 15, 'f', 12);
2785+// qDebug() << QString("JD %1: Pecl,z: %2 Difference: %3").arg(JulianDay, 8, 'f', 5).arg(VEC3, 15, 'f', 12).arg(VEC3-0.91433656053126552350, 15, 'f', 12);
2786+// qDebug() << QString("JD %1: Pequ,x: %2 Difference: %3").arg(JulianDay, 8, 'f', 5).arg(X_A , 15, 'f', 12).arg(X_A +0.29437643797369031532, 15, 'f', 12);
2787+// qDebug() << QString("JD %1: Pequ,y: %2 Difference: %3").arg(JulianDay, 8, 'f', 5).arg(Y_A , 15, 'f', 12).arg(Y_A +0.11719098023370257855, 15, 'f', 12);
2788+// qDebug() << QString("JD %1: Pequ,z: %2 Difference: %3").arg(JulianDay, 8, 'f', 5).arg(VEQ3, 15, 'f', 12).arg(VEQ3-0.94847708824082091796, 15, 'f', 12);
2789+
2790+
2791+ Vec3d PECL(P_A, VEC2, VEC3);
2792+ Vec3d PEQR(X_A, Y_A, VEQ3);
2793+ Vec3d EQX=PEQR^PECL;
2794+ EQX.normalize();
2795+ Vec3d V=PEQR^EQX;
2796+ Mat3d RP(EQX[0], EQX[1], EQX[2], V[0], V[1], V[2], PEQR[0], PEQR[1], PEQR[2]); // result from paper, section A.3
2797+ // Now we create the (hopefully) same rotation matrix from the Capitaine angles.
2798+ // Mat4d Rrot=Mat4d::zrotation(chi_A) * Mat4d::xrotation(-omega_A) * Mat4d::zrotation(-psi_A) * Mat4d::xrotation(eps0);
2799+
2800+ // Uh-oh - it seems the A&A paper has the matrix operations in the other direction.
2801+ // So the matrix to be compared must be composed in reverse.
2802+ Mat4d RRot=Mat4d::xrotation(eps0)*Mat4d::zrotation(-psi_A) * Mat4d::xrotation(-omega_A) * Mat4d::zrotation(chi_A);
2803+
2804+ //qDebug() << "RP : " << RP.toString(15, 'f', 12);
2805+ //qDebug() << "RcapTr : " << RRot.upper3x3().toString(15, 'f', 12);
2806+
2807+ Mat4d RP4(EQX[0], EQX[1], EQX[2], 0, V[0], V[1], V[2], 0, PEQR[0], PEQR[1], PEQR[2], 0, 0, 0, 0, 1);
2808+ //QMatrix4x4 matDiff=((RRot-RP4)).convertToQMatrix();
2809+ //qDebug() << matDiff;
2810+ Mat3d matDiff3x3=(RRot-RP4).upper3x3();
2811+ double max=0.0;
2812+ for (int i=0; i<9; ++i)
2813+ {
2814+ if (fabs(matDiff3x3[i]) > max)
2815+ max=fabs(matDiff3x3[i]);
2816+ }
2817+ //qDebug() << "largest value in difference of matrices:" << max;
2818+ QVERIFY2(max<2e-5, QString("Some values in the precession matrices differ by too much.").toUtf8());
2819+
2820+ double angleRef=RP.angle()*180.0/M_PI;
2821+ double angleCap=RRot.upper3x3().angle()*180.0/M_PI;
2822+ //qDebug() << "Rotation angle of reference matrix:" << angleRef;
2823+ //qDebug() << "Rotation angle of Capitaine matrix:" << angleCap;
2824+ //qDebug() << "Difference (arcseconds):" << (angleCap-angleRef)*3600.0;
2825+ // according to Fig12 in the paper, a few arcseconds of difference between PQXY and Capitaine parametrisation are allowed.
2826+ QVERIFY2((angleCap-angleRef)*3600.0<6.0, QString("Angle between rotation matrices too different!").toUtf8());
2827+
2828+ //TODO: Add more dates and verify this angle difference is limited to what we can see in Fig.12
2829+}
2830
2831=== added file 'src/tests/testPrecession.hpp'
2832--- src/tests/testPrecession.hpp 1970-01-01 00:00:00 +0000
2833+++ src/tests/testPrecession.hpp 2015-07-28 19:18:20 +0000
2834@@ -0,0 +1,35 @@
2835+/*
2836+ * Stellarium
2837+ * Copyright (C) 2015 Georg Zotti
2838+ *
2839+ * This program is free software; you can redistribute it and/or
2840+ * modify it under the terms of the GNU General Public License
2841+ * as published by the Free Software Foundation; either version 2
2842+ * of the License, or (at your option) any later version.
2843+ *
2844+ * This program is distributed in the hope that it will be useful,
2845+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
2846+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
2847+ * GNU General Public License for more details.
2848+ *
2849+ * You should have received a copy of the GNU General Public License
2850+ * along with this program; if not, write to the Free Software
2851+ * Foundation, Inc., 51 Franklin Street, Suite 500, Boston, MA 02110-1335, USA.
2852+ */
2853+
2854+#ifndef _TESTPRECESSION_HPP_
2855+#define _TESTPRECESSION_HPP_
2856+
2857+#include <QObject>
2858+#include <QtTest>
2859+#include "precession.h"
2860+
2861+class TestPrecession : public QObject
2862+{
2863+ Q_OBJECT
2864+private slots:
2865+ void initTestCase();
2866+ void testPrecessionAnglesVondrak();
2867+};
2868+
2869+#endif // _TESTPRECESSION_HPP_