Refraction calculation excessive at negative horizons ?

[glenenglishgithub]

I found predicted moonrise times to be very early (2-3 minutes) compared to observations for negative horizons. In this case, altitude was 1270 m which has an ocean horizon ~ -1.14 deg. This was for 1013 Pa and 14 deg C.

I picked out the refract.c and ran refract.c alone and found that at -1.14deg Ta, refraction of 0.9 deg was being applied, resulting in *aa being set to ~ -2.04 deg ! Is this refraction calc invalid for negative horizons ? I see there are a few methods around.

I modified refract.c so that the refraction correction would never exceed 34/60 and now the moonrise times are very close (seconds!) with exception of a moon elevation question I posted on stackoverflow.

Great software Brandon, nice work. I am using it for very low angle ham radio moonbounce calculations for world records where the elevation of the moon is about 0.5 deg at each end , IE working at the antipodes, there is only a few minutes to make the communication.

regards glen english Canberra, Australia.

[brandon-rhodes]

Thanks very much for the kind comments, glen! Though most of my attention these days is focused on my Skyfield astronomy library, I'm happy to see PyEphem continue to get used.

So that I can be sure I can reproduce exactly your problem, could you provide a few-line Python file that, if I run it, will show me the output you're getting so I can compare it to what you'd like to receive instead?

[glenenglishgithub]

Hi Brandon, thanks very much for the reply. I guess the first Q is should I be using Skyfield ? OK so in a nutshell : the program , for a specified date and location pair looks for a moonrise at either location A (called HS- home station) and location B DX - the other station. From Moonrise derived using the rise_time class property. In runs through n days with m timesteps for our moon tracking and guesswork. Now.. For the example files provided, this was for the purposes of confirming the calcs. The HS site is Mt Wellington in Tasmania at 1270m , which you will see commented IN on line #52. My program also calculates the horizon which was this location is -1.14 degrees,, I calculate using simple geometric calc.

The other end, location B 'DX" which doesn't really matter for the moonrise check is out in the pacific ocean. At Mt Wellington, the ocean is visible on the rising AZ for that moonrise day.

For the dates used on line #20 PyEphem unmodified says HS moonrise is 2018/9/27 10:12:29 UTC that is produced in the STDOUT printout at line 20, and also the first line of the generate CSV file.

However! Actual moon rise (upper limb def) is at 10:14:24 UTC observed.. some 2 minutes later. comparison with some other calc options suggested pyEphem is over compensating for refraction.

I wrapped GCC around refract.c and found that pyEPhem was using 0.9 deg refraction correction. I modified refract.c to limit the refraction to 34/60 degrees and rebuild pyEphem and it came out within 10 seconds of observed !. Line #87 is modified, there is a #define for it. Interesting eh ? The reason why we are attempting this level of accuracy is that our window for stations at the antipodes where both have a viable moon is minutes ! And it takes minutes to complete a conversation (called a QSO) . We understand that even very minor variations of refractive index, possible over the 200km of atmosphere we look at the moon at the horizon of -1.14deg over ocean may affect moonrise times by minutes ....

The other question for me is that at all locations, sea level, whatever, pyEphem at the moonrise time reports the moon is as elevation (altitude in PyEphem) of the location horizon. IE when at sealevel with horizon 0, PyEphem reports 0.0 moon atltitude. . But since the upper limb visibility is the moonrise definition, I would have thought pyephem would have reported the moon at -0.25deg (IE the centre is a moon radius below the horizon) . Maybe I misunderstand defs.

cheers glen english

files attached. ephemtest.zip

[brandon-rhodes]

Wow — it's amazing to think of my astronomy library powering a real-life conversation between radios on opposite sides of the Earth! Thanks for sharing the details of your application.

You'll indeed find Skyfield a bit more accurate than PyEphem, but probably PyEphem is still accurate to within a second or so if properly managed. I'll take a look at your code — it's been a long time since I've fiddled with the old libastro diffraction routine to remember how much it fluctuates with the temperature and pressure settings it can be given.

I remember learning that the USNO doesn't even try to predict sunrise and moonrise with real simulations of diffraction; instead they just choose the moment when the center is 90.8333 degrees exactly from the zenith!

http://aa.usno.navy.mil/faq/docs/RST_defs.php#riseset

[glenenglishgithub]

Hi Brandon, cheers.

Is it worth a look at me looking at the refraction calc source used for Skyfield ? Happy to modify my program to use Skyfield also. I am a competent programmer /RF engineer.

Yeah the antipode communication is really only possible with negative horizons at each end, and for our world record attempt , one station will be up in the mountains of the South Island of New Zealand looking at moonrise over the South Pacific, while the other station will be on a big hill in Portugal looking west out over the Atlantic ....

when I get this sorted, I will publish my program and it will probably become the standard for moon bounce schedule planning . Most of which usually work with maybe 1 to 2 hour common windows- that is stations that are far from the antipodes, with good common moon windows. cheers

[brandon-rhodes]

Come to think of it, the Skyfield moonrise code uses the USNO definition exactly; so you'd have to make your own little routine (see "almanac.py" for the pattern) to have a "realistic" mooonrise that uses temperature and pressure of your own devising. So maybe sticking with PyEphem makes more sense for now?

Here's the USNO form at http://aa.usno.navy.mil/data/docs/RS_OneYear.php filled out, I think correctly?, for Mount Wellington:

It gives:

http://aa.usno.navy.mil/cgi-bin/aa_rstablew.pl?ID=AA&year=2018&task=1&place=&lon_sign=1&lon_deg=147&lon_min=14&lat_sign=-1&lat_deg=42&lat_min=54&tz=0&tz_sign=-1

Which, if I'm reading correctly, gives a moonrise of 10:20, which is even later than your real observation.

The first thing that comes to mind is: when you set the "elevation" in PyEphem (and in Skyfield), all it does is assumes a flat horizon but up on higher ground — it does not, unless I'm forgetting, interpret "elevation" (height of the ground) as "altitude" (the observer's height above the ground) and so doesn't account for the fact that someone on a mountain — or even someone on a building several stories tall — will see sunrise and moonrise earlier because the light hits them earlier than someone down at the ground.

Your observation of moonrise 6 minutes earlier than the USNO's prediction might be because the mountain is above the surrounding terrain — which might have a much larger effect than temperature or pressure?

[glenenglishgithub]

Hi Brandon It is possible that the mountain , and the 200km of atmosphere are having an effect not fully appreciated understood. I think we will need more observations to get good data. The mountain is about 20km from the ocean. More datapoints/measurements required certainly to base some conclusion.

Wow Skyfield looks impressive ! I ran the refraction calc of SKYFIELD in a python scratchpad and I note that it returns 0 for alt_degrees less than -1.0 .

I calc the geometric horizon angle for the elevation and assume in this case, ocean, and I input the horizon angle into pyephem observer.horizon.

Given then PyEphem will provide the location centre of an object , if I am at sea level, why does pyephem return moon altitude = 0.0 for moonrise and not -0.25 given that the upper limb is just visible? maybe I need to reread the def. Anyway, I will take a good look at skyfield. Pure python makes it a bit easy to modify.

cheers and thanks for your assistance.

[drbitboy]

[Brandon wrote:]

Your observation of moonrise 6 minutes earlier than the USNO's prediction might be because the mountain is above the surrounding terrain — which might have a much larger effect than temperature or pressure?

1.14deg * ( (1440min/d) / (360deg/d) ) => 4.56min  (not 5.56 d'Oh!)

Sounds like Brandon's conjecture is correct (gotta love good intuition; it makes debugging things so much easier).

I assume 1.14deg comes from

tan ⁻¹(√((2×1.27÷6378.1)+(1.27÷6378.1)^2))

with perhaps some other value between 6378.1 and 6371,0 (https://en.wikipedia.org/wiki/Earth)?

The obvious hack would be to find the azimuth to the moon and ask PyEphem for the moonrise at a point T, at sea level some

√((2×1.27*6378.1)+1.27^2)

= 127.29km along that azimuth from Mt. Wellington. Assuming that you make a good estimate of the pressure and temperature at T, and that PyEphem's refraction calculations moonward from that point are correct, unmodified PyEphem should give you a much better estimate of the moonrise; the remaining second order error should be the (slight?) refraction between T at sea level and Mt. Wellington at 1270m.

(This is very cool, I am surprised every other PyEphem/Github nut on the planet is not offering their two cents, as well;-)

P.S. Sorry, but cannot resist this plug: funniest moment in all of YouTube (IMNSHO): https://www.youtube.com/watch?v=Er5K_nR5lDQ&t=1840s

[glenenglishgithub]

Great suggestions thanks , on asking PyEphem to provide time at point T . I will write those in. Good point about the 2nd order effects from 1270 down to sea level (and ocean..). The reality is refraction can only be estimated short of using refractive index measuring radio or acoustic sounders...-glen

[brandon-rhodes]

I remember reading in an astronomy book about the fun trick of standing at the top of a 30-foot high stairway down to a beach, waiting til the very moment of sunrise, then running down the stairs and getting to watch the sun rise again a few moments later. Even a few dozen feet changes the moment of rising; the height of a mountain above the surrounding terrain, even more so.

[glenenglishgithub]

The refraction calc is probably the weak link. Maybe an integrated method is going to be best, given the varied atmospheric densities along a cold high altitude to warm humid low altitude path grazing over the ocean...

[drbitboy]

this quotes a nominal max of 35.4' of retraction i.e. 0.59deg, or roughly half of your 1.14deg, equivalent to 2.36min of earth rotation (or does Mt. Wellington not being at the equator affect that degrees-to-minutes calc?)

[drbitboy]

d'Oh! updated the conversion, several entries above, from 1.14deg to minutes of rotation.

[drbitboy]

You want to see some cool radio-atmospheric interactions, check out Figure 3, page 29 of this.

My attempt to summarize REX is on p. 110 here.

Yes, I know I am off topic.

[glenenglishgithub]

Original PyEphem computes refraction at 0.676deg for neg 1.14deg. However it computes 0.472 at zero deg for 1013 and 10 deg C.

The 64$ question is , what to use for temp and pressure for computing refractive index over a 127 km path to the horizon from a high location... and also the moon is of course going through some 127km of atmosphere on the other side of the horizon (thanks the Rex Moncur for pointing this out) , IE the moon is not a ship on the horizon).

[drbitboy]

The refraction calc is probably the weak link. Maybe an integrated method is going to be best, given the varied atmospheric densities along a cold high altitude to warm humid low altitude path grazing over the ocean...

You know, I'd guess the second-order effect of temperature and pressure are small; once you have several months or a year of data, I wouldn't be surprised if single number, somewhere not far from that 35.4' from that Wiki page, is close enough. Like you say, short of some kind of direct measurement, from a point in the ocean and 127km from where you want to be at the time no less, refraction can only be estimated anyway. 15 angular minutes of earth rotation is one time minute; so if you have a fixed number for refraction within 20% of the actual values' range, which we can assume is in that 35-40' range, then PyEphem should get you within 30s of the actual rise time.

[glenenglishgithub]

Agreed. When I limited PyEphem to 34/60 (0.5666 is was within seconds of the actual). but that was one datapoint... Any reason not to move the code into using Skyfield ? Or no reason ? I have no looked at skyfield yet except to check out what refraction calc it uses (and is invalid below -1.0 deg) -glen

[drbitboy]

Any reason not to move the code into using Skyfield ? Or no reason ? I have no looked at skyfield yet except to check out what refraction calc it uses (and is invalid below -1.0 deg) -glen

I don't see a reason to go either way, other than that SkyField is being actively developed. I think your best bet is to see how that 34/60 value holds up, or better yet make the observations and see how it varies, and whether you can correlate it with time of year, temperature, pressure, etc. 127km is a long way, and there are another wehavenoideahowmanymore km on the other side about which you know bupkis from your perch.

My guess is that you find the "typical" value and you will get down to 10-20s accuracy; is that good enough?

[glenenglishgithub]

yeah , 30 seconds is fine . A minute is desirable. For radio work, the majority of reflection if from the inner 1/3 of the radius of the disk, approx . So it is the centre of the moon we are interested in. The antipode attempt is a hard one- the moon's lower limb will be on the horizon moon-setting end simultaneously when the lower limb is on the horizon at the moon rising end ! It takes several minutes to exchange information because the information rate is very very low. I might point out, moonbounce is nothing new, but working with a common moon window of minutes is unusual and difficult.

[glenenglishgithub]

Hi Brian Remaining question I still don't understand- PyEphem reports moon rise as 0.0 elevation (for sea level). but the moon rise definition is for upper limb at horizon, implying centre of disk is at -0.25 Any (moon) light on that subject ? Is PyEphem reporting moon rise as centre of disk on horizon and not per navy upper limb ?

PS, I found this little gem "The horizon attribute defines your horizon, the altitude of the upper limb of a body at the moment you consider it to be rising and setting."

OK, so that will be different for different (angular) diameter bodies.

so that explains that ! PyEpehem reports it for MY horizon for THAT Body OK I understand now. regards, glen

[glenenglishgithub]

Readers of this forum that are interested in improving their refractive index calcs might find this reference useful :

http://www.dtic.mil/dtic/tr/fulltext/u2/a137095.pdf

Section 4, conclusions.and then the appendix are also very useful.

and this more recent one :

https://prod.sandia.gov/techlib-noauth/access-control.cgi/2012/1210690.pdf

[drbitboy]

Woo hoo, what a find. That's a keeper, thank you! Personal favorite extracts:

However, it must be emphasized that very accurate
propagation prediction still requires detailed ray trace calculations in the context of
accurate atmosnheric modI1.

Ray trace computations based on the concentric shell atmospheric
model have been performed using an HP-41CV system.

[SYSTEM, not calculator!]

Oh man, check out the HP41C Wiki page: there are emulators and simulators! Maybe we could get the original code.

Related topic: could you calibrate the surface refractivity Ns by measuring the depression angle to a known point e.g. the coast, or a distant buoy? Measuring accuracy would be the limitation, I suspect.

[glenenglishgithub]

Hi Brian, yeah there are quite a few ray tracing options around. I have been playing with Skyfield this morning, getting my head around it and managed to get it all agreeing with PyEphem. Skyfield is certainly a big hammer. As for the buoys- good idea. otherwise acoustic or baloon radio sondes are the gold standard I'll look at moving to Skyfield, though I might just use PyEphem and apply refraction elsewhere (setting pressure to 0 to force the internal compiled in refraction calculation to zero) For Skyfield , I'll need to generate my own moorise, moonset calc but that is an easy iterative job, and I am running my own refractive calcs ... Next step will be calculating topocentric librations... Cannot see any libration keywords in skyfield so guess can roll my own from all the other info. Librations cause a random doppler shift to be applied, spreading out the signals, proportional to frequency, in this case, 10368 MHz.. Best performance is when librations between two stations are a minima. regards,

[drbitboy]

Librations? As in some wobble in the moon's attitude, or earth's?

For high accuracy stuff you might want to go with NAIF/SPICE, an even bigger hammer than SkyField. It is an object-oriented-ish approach to dealing with 3D geometry.

It has a High accuracy EOP-based Earth rotation model, updated at least twice per week. You would need to add in refraction, but it sounds like you are doing that anyway.

Andrew Annex create a Python interface to SPICE here that is quite good.

The learning curve is steep, and kernel management (the data part of the object-iented approach) but I have been using it for nigh on to three decades, so I could help you with any problems.

For example, say you said that after your refraction calc, you know that moonrise will be when the center of the moon so (90+delta)degrees away from your local zenith; after the one-time setup of the data kernels, the time for that event would be a single call.

Right now I am using it to calculate the twist for the New Horizons flyby of KBO MU69 that results from an out-of-plane delivery with respect to the nominal trajectory, and how that affects the suite of detectors that will be trying to get a glimpse of that target.

[glenenglishgithub]

New Horizons eh? Inspirational mission for all engineers. Amazing work. Currently on my desk I have a pile of 2nd hand books I bought in the past couple of months, all on ephemeres, astrodynamics, and some astro physics, all with the aid of better understanding the background of the basic stuff. I am a RF and signal processing engineer in my professional and hobby work.

thanks for the tips on NAIF/ SPICE Charlie Suckling, a colleague of mine wrote up the modern moonbounce moon libration paper a while paper atteched below in next message

[drbitboy]

Link did not work.

Coolest site on the internet: https://eyes.nasa.gov/dsn/dsn.html

Brilliant display of info (IMNSHO;-). Don't forget to click on the information button to get the legend.

[glenenglishgithub]

here is the paper (attached) G3WDG_libration paper revised.pdf

great site !

[brandon-rhodes]

“Librations cause a random doppler shift to be applied” — wow, I had not even considered that effect! Am I imagining correctly that it's the derivative of the libration — not the instantaneous angle of libration — that spreads the signal?

[brandon-rhodes]

(And — very inspiring links from both of you, thanks!)

[glenenglishgithub]

Hi Brandon Yeah. While the mani application for my program was to produce very time dense high resolution moon calcs, I did a quick and dirty 1st derivative in my EMEcalc (EARTH MOON EARTH ) program using PyEphem to produce libration spreading values , using the geocentric librations from pyephem and doing a dirty low precision conversion to provide topocentric, but I need to improve on that. There is plenty of good solid literature on computing Topocentric librations. its on my list to do.

The bandwidth of the signals used is very narrow, IE information bandwidth , as there is not much signal, and up at 10368 MHz, the typical receiver bandwidth is ~ 1 Hz, and spreading varies from 2 to 5 Hz up to 100Hz ish. So knowing the spreading on a particular date is useful for contact planning. There are a couple of software packages out there that compute all this for EME but not to any great precision (which I have changed with PyEphem ) , nor are they maintained anymore and indeed one of the packages used the writer is deceased, so, it's time to provide the community something better,

Skyfield is certainly a full bottle ! I think PyEphem provides all the facilities required, however Skyfield is pure python so there is some customization advantage. I'll need to run some more computational requirement numbers, as we like to produce quite time dense almanacs, and PyEphem with the compiled C engine is very fast indeed, roughly 100x faster than the old programs used for EME calc.

Currently we are right at the limit of what can be done in narrow bandwidths. (other authors). I have a new method which uses a type of spread spectrum, where the message is spread over a wide bandwidth when it is transmitted, and then "de spread" at the receiver-decoded.....and so when I 'de-spread' the message, the libration spreading gets reduced also, reducing the degradation at high spreading values. This is necessary to enable EME communication at much higher frequencies with smaller stations- this means smaller dish sizes, lower power etc.

-glen

[brandon-rhodes]

Yes, Skyfield doesn't always have the zing of PyEphem's C code — alas! But I never get complaints about it not compiling, nor does it have odd bugs under 64-bit Windows like PyEphem does, so I'm for the moment content.

It's about time that Skyfield learn about lunar libration, so I just spent an hour on the problem. I now need to start the rest of my day, so here's where I'm at so far:

1. I did a Google search for SPICE and libration and got https://naif.jpl.nasa.gov/naif/lunar_kernels.txt which recommended the file moon_pa_de421_1900-2050.bpc which I found here: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/pck/
2. So I looked up "bpc" files and got the complete guide to how that file is laid out: ftp://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/pck.html
3. It looks very similar to SPK files so I cut and pasted and added a new experimental module to jplephem which you can get by doing pip install https://github.com/brandon-rhodes/python-jplephem/archive/master.zip
4. For real data to compare to (HORIZONS doesn't do libration, that I can see?), I'm going to use: http://aa.usno.navy.mil/imagery/disk?body=moon&year=2018&month=10&day=1&hour=12&minute=34
5. And just to make sure I get that Julian date correct: http://aa.usno.navy.mil/jdconverter?ID=AA&year=2018&month=10&day=1&era=1&hr=12&min=34&sec=0.0

Here's a program which asks the libration kernel for the Moon's attitude, and also asks Skyfield for the current position of the Moon relative to Earth:

from skyfield.api import load, tau

ts = load.timescale()
t = ts.utc(2018, 10, 1, 12, 34)

eph = load('de421.bsp')
e = eph['earth']
m = eph['moon']

p = e.at(t).observe(m).apparent()

from jplephem.binary_pck import BinaryPCK

m = BinaryPCK.open('moon_pa_de421_1900-2050.bpc')
print(m)
jd = 2458393.023611
for s in m.segments:
    for piece in s.compute(jd, 0):
        print([str(n) for n in (piece)])

ra, dec, distance = p.radec()
print(ra)
print(dec)
print(ra.radians)
print(dec.radians)

So there's a big mystery and a little one.

The little one is why the body # and frame # from the kernel don't make sense. The Required Reading doc that I linked above clearly says in "Binary PCK Kernel Format" that the body number should come first followed by the frame but the two numbers appear to be 31006 and 1 and that first number looks a lot like a SPICE Frame ID, not like a body number.

The big mystery is whether the Chebyshev that I cut and pasted from the SPK logic is entirely correct here; it ran without error, but the numbers coming out don't look like three angles:

File type DAF/PCK and format LTL-IEEE with 1 segments:
2415020.50..2470172.50 frame=1  Unknown Body (31006)
['-0.054470029238916196', '0.424337745960259', '4139.103878849478']
['0.0002157918637080179', '-0.0001307530018170278', '0.22976293055855648']
05h 45m 30.77s
+20deg 21' 11.1"
1.5075840638817268
0.3552282142495826

In particular 4139.103878849478 does not look with either a sensible number of radians nor of degrees (unless it's an angle that's "unwrapped" so it increases continuously as the body rotates? I guess we should graph it to see how it behaves through time).

In any case, if either of you gets the chance, I think the last step here is to take the 3 "angles" (if that's what I've really extracted from the ephemeris here), use the rotation matrix routines in Skyfield to turn them into polar angles (see functions.py in Skyfield), and then see if the RA and Dec direction of the Moon's central point when compared to the RA and Dec of us looking back produces the libration offsets given by the USNO.

And if they don't, we need to figure out what step I'm missing :)

[drbitboy]

1) the first number will almost certainly be a frame ID (possible typo in pck doc, I suspect), and the second number probably is the frame ID to that first frame ID is relative: frame ID 1 is the EME J2000 frame (here, EME means Earth Mean Equator-of-date; Glen, you can imagine how confused I was reading that paper - it's great, but it needs an acronym table for the non-cognoscenti, or anyone without a 6-character moniker ;-).

(I may be going on- and offline for a bit; I did something like [pip install -t /usr/ ...] as root, and after that there were only three files in /usr/bin/ and they all started with "pip" - hoo boy, did my butt cheeks tighten up.

[brandon-rhodes]

Wow. I always use conda environments or virtual environments these days, and try to avoid doing anything system-wide.

Your interpretation of the fields makes sense; based on that guess, I'll update the code tonight or tomorrow, and then email the SPICE folks to see if they'll update the docs.

[drbitboy]

Binary PCKs came into SPICE during NEAR (ca. 2000), but have been largely replaced by the FRAME infrastructure and CKs. That said, I think NAIF will have to continue supporting the high-accuracy Earth and Moon PCKs.

1.1. I forgot to mention in the previous post: you should be able to run the SPICE utility BRIEF on PCKs (or you could two decades ago). Use [brief -t -c] which should tell you the base-relative target frame and the base frame; the -n option will print them both as frame IDs i.e. numbers, but without an FK, SPICE may not know the frame name for frame ID 31006.

2. those three angles that you got from the Chebyshev evaluation are almost certainly some form of Euler angles wrt the base frame (J2000; frame ID 1) - the first two express small offsets wrt the J2000 frame, and the fourth is the around-Z rotation from that near-J2000 attiude. If you plot the value over 27days, I suspect it will change by about two PI. If they chose to limit that third value on output from the Cheby eval to a range of [0:2PI), then they would have to have one segment per synodic moon period. If you look at the SPICE code, it may reduce that number modulo 2PI, or not at all since it will be passed to trig sine and cosine functions anyway.

[drbitboy]

Wow. I always use conda environments or virtual environments these days, and try to avoid doing anything system-wide.

Yeah, I was in too much of a rush. I was trying to do that (update the installation under $HOME/.local/), but [pip] somehow got bitrot, and I was just trying to fix that so I could go back to saftey. Definitely just carelessness and stupidity. And arrogance: "yeah, I'm root, but this can't make it any worse"

Anyway, I was able to cut a LiveCD (16.04.5 LTS, same as my laptop - it's good to be smart, but way, way better to be lucky) on another system, boot up to [Live], then [mount /dev/sdb1 /mnt] [rsync -av /usr/bin/ /mnt/usr/bin/]. That put the wheels back on, but the carburetor still needs tuning, the timing is off, and there is a battery is in front of the passenger seat with jumper cables going out the window and under the hood. Went to bed at 5am. Phew.

[drbitboy]

yeah, 31006 and 1 are definitely the target and base frame ID; cf. from this FK:

   Frames Specified by this Kernel
   =====================================================================

   Frame Name       Relative to        Type   Frame ID
   --------------   -----------------  -----  --------
   MOON_PA          MOON_PA_DE421      FIXED  31000
   MOON_ME          MOON_ME_DE421      FIXED  31001
   MOON_PA_DE421    ICRF/J2000         PCK    31006
   MOON_ME_DE421    MOON_PA_DE421      FIXED  31007

[drbitboy]

And here is the brief SPICE utility output, as expected:

$ brief -t -c *.bpc

BRIEF -- Version 4.0.0, September 8, 2010 -- Toolkit Version N0066

Summary for: moon_pa_de421_1900-2050.bpc

Frames                Start of Interval (ET)          End of Interval (ET)
-------               -----------------------------   -----------------------------
31006 w.r.t. 1 J2000  1900 JAN 01 00:00:00.000        2051 JAN 01 00:00:00.000

[brandon-rhodes]

So it is an unwrapped angle — because otherwise they couldn't polynomialize it! That makes perfect sense but I hadn't even thought of it.

So both of my bare questions are answered, and I'll just need to jump in to using the resulting angles and see if I can get a relative libration angle. (If you look through the Required Reading, you'll find descriptions of the three angles — two polar adjustments, then the rotation term.)

[brandon-rhodes]

And I suppose at some point I'll have to learn about ck files. I'll add this to my to-read list:

https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/ck.html

[brandon-rhodes]

Or — or is that specific to spacecraft? And what would involve Moon orientation would be this one instead?

https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/frames.html

[drbitboy]

And, just to complete the picture, here is an abbreviate summary of that BPCK by the SPICE utility SPACIT:

$ spacit[...]
   ( S ) Summarize binary file[...]
   Option: s
   Binary file     : moon_pa_de421_1900-2050.bpc
   Leapseconds file: tmp/naif0012.tls
         PCK Summary Options
   ( F ) Summarize entire file.
   Option: f
********************************************************************************
Summary for PCK file: moon_pa_de421_1900-2050.bpc
Leapseconds File    : tmp/naif0012.tls
Summary Type        : Entire File
--------------------------------------------------------------------------------
   Segment ID     : de421.nio
   Body           : Body 31006
   Reference frame: Frame 1, J2000
   PCK Data Type  : 2
      Description : Fixed Width, Fixed Order Chebyshev Polynomials: Angles
   UTC Start time : 1899 DEC 31 23:59:18.816
   UTC Stop time  : 2050 DEC 31 23:58:50.816
   ET Start time  : 1900 JAN 01 00:00:00.000
   ET Stop time   : 2051 JAN 01 00:00:00.000
--------------------------------------------------------------------------------[...]

[drbitboy]

No, I think they dropped BPCKs because, once the FRAMES infrastructure was in place, BPCKs and CKs are basically doing the same thing, and so few people used BPCKs that it was better to only advertise one. In the meantime, BPCKs hang around, like stale gas, for backward compatibility. (this is all just conjecture, but they definitely push new missions toward CKs and away from BPCKs).

Funny story about BPCKs: at some point for NEAR, we had a rough shape model, and there was this navigator: brilliant; curmudgeon; thought of himself as an 800-pound gorilla and thought the mission revolved around him (hmm, a lot like me;-). Anyway, he sees the BPCKs as a way to model wobble based on a homogeneous-density rotation model using the shape we had (the text PCKs can do a little, but not much). And me, I'm working with the scientists, and I want to be able to debug and do simple things like read the text PCK in VI and say "the pole of Eros points to this RA and DEC," so I do not want a BPCK with all the hidden internals and extra (what, three?) step to get info (yes, I admit I am lazy). The curmudgeon sends me a binary PCK, and the first thing I want to do is validate it. Hmm, I know, plot the NEAR trajectory in the rotating, body-fixed NEAR frame: something like two lines of code and a print statement in a loop to get the data. Great. So I plot it up (gnuplot, probably - this was all FORTRAN and C in the day), and the trajectory looks like a conical spring, but the axis of the spring is the Eros X axis, not the Z axis. So I send it back, and he makes another, and it checks out. So, I check to see how much wobble there is in the Well-Made Super-Duper BPCK. So I create a simple POLE_RA, POLE_DEC, W, Wdot linear rotation model in a text PCK based on the mean info in the BPCK, and calculate the difference between that and the BPCK: less than a pixel in the camera during the closest orbital operations. Definitely the 800-pound tail wagging the dog.

[brandon-rhodes]

@drbitboy — Looking at:

ftp://naif.jpl.nasa.gov/pub/naif/misc/toolkit_docs_N0063/FORTRAN/req/spk.html

— it doesn't exactly describe PCK frames and CK frames as interchangeable. They sound like they have two different purposes: "PCK frames are bodyfixed frames" vs "CK frames are frames that are defined relative to a spacecraft structure".

So unless @glenenglishgithub gets to it before I do, I think I'm going to assume that the binary PCK module I created this morning is indeed how we'll access the data for Moon libration. But if you do run across a more recent Moon libration dataset for SPICE with a more recent extension, please let me know and I'll re-orient!

[glenenglishgithub]

Hi Brandon. I am going to leave that one for you , if you would oblige , many thanks. I have much reading to do before I would be close to being ready to manipulate such sources- glen

[drbitboy]

I know what it says, but a frame is a frame, and PCKs and CKs both express the transformation between two frames - an attitude ephemeris, so in that sense they are interchangeable. The FRAMES subsystem is completely flexibile, and there are five classes of frames: 1 inertial; 2 PCK; 3 CK; 4 Fixed; 5 Dynamic. Although that is important in defining the frames, it is irrelevant from a user perspective when accessing the transformation from one frame to another - SPICE finds a path the same way it finds a path for SPKs e.g. from the Moon to Earth to Earth Barycenter to Sun to Solar System Barycenter to Mars Barycenter to Mars to Phobos, if you ask for the state of Phobos wrt the Moon.

Actually, there is one difference, and that is that the internal time lookup in PCKs is ET (TDB/TDT/etc., linear s past the J2000 epoch); the internal time lookup in CKs is a "Spacecraft CLocK" or SCLK, which needs to be tied to the frame logicially, and tied to ET numerically.

For OSIRIS-REx we are starting with text PCKs, but we have been told that if Bennu is tumbling or otherwise has an attitude ephemeris that cannot be described with text PCKs, there will be a CK created for its frame.

[brandon-rhodes]

It looked, at first glance, like CK can only do simple time series (so good for a tumble?) but PCKs can do extended series of periods each defined by their own Chebyshev polynomial, good for long-term predictions of bodies that over those long periods are perturbed by other masses. So my guess is that a lunar ephemeris needs PCK for high precision?

But this is my first exposure to them so I'm still reading!

[drbitboy]

yes I agree PCKs definitely are more suited for bodies.

in the end both do the same thing: provide attitude ephemeris via the FRAMES subsystem.

So, once the kernels are loaded, the application is agnostic as to the type of kernel.

[glenenglishgithub]

EMECalc now is using Skyfield. That was easy. Will post in the next couple of days when I run some more numbers.
And like all Numpy etc based things, speed is no problem for Skyfield for the calcs I am doing ... Can't tell any speed difference...

Python is optimization the right way around in my opinion. IE "get it working first, optimize if required, later" glen.

[brandon-rhodes]

I'm glad that the conversion wasn't a trauma! Let me know about any sharp edges that I might have documented better.

I have a script which takes the Earth→Moon vector, rotates it by the three angles provided in moon_pa_de421_1900-2050.bpc, and very nearly gets the angles given by the USNO for the Moon's orientation (which I pasted in as a comment at the bottom of the script):

from numpy import einsum
from skyfield.api import load, tau
from skyfield.functions import rot_x, rot_y, rot_z, length_of, to_polar

ts = load.timescale()
t = ts.utc(2018, 10, 1, 12, 34)

eph = load('de421.bsp')
earth = eph['earth']
moon = eph['moon']

p = earth.at(t).observe(moon) #.apparent()?
ra, dec, distance = p.radec()
v = p.position.au / length_of(p.position.au)  # unit vector

from jplephem.binary_pck import BinaryPCK

m = BinaryPCK.open('moon_pa_de421_1900-2050.bpc')
jd = t.tdb
jd -= p.light_time
morient, morient_v = m.segments[0].compute(jd, 0)
a, b, c = morient

def rotate(rotation, vector):
    return einsum('ij...,j...->i...', rotation, vector)

print('Pole RA', a)
print('Pole dec', b)
print('Rotation angle', c)

print('Direction from Earth to Moon', v)

v = rotate(ts.J2000.M, v)  # adjust from ICRS to J2000?

v = rotate(rot_z((a)), v)
v = rotate(rot_x((-b)), v)
v = rotate(rot_z((c)), v)

print(v)

distance, lat, lon = to_polar(v)
print(lat * 360 / tau, 'or', lat * 360 / tau % 360)
print(lon * 360 / tau, 'or', lon * 360 / tau % 360)

# http://aa.usno.navy.mil/imagery/disk?body=moon&year=2018&month=10&day=1&hour=12&minute=34
# gives the numbers:
# Sub-Earth: Lat 3.9°, Lon 356.4° Sub-solar: Lat 1.4°, Lon 277.8° Phase: 0.600 Diameter: 31' 56.7"

The angles the script prints are around 3.8° and 356.9°. Is that because I'm doing the calculation nearly correctly but am neglecting some effect? Or am I doing the calculation all wrong and just happened to have fiddled until I by coincidence got the numbers close?

I guess I should pull down the Explanatory Supplement this evening and see if they outline how the USNO computes the value.

[glenenglishgithub]

Hi Brandon The skyfield doco is pretty good. The section "Positions and Coordinates¶" is probably the most important read. from python console : help(jplephem) was very useful !

with regard to :

from jplephem.binary_pck import BinaryPCK

I have package jplephem [2.8] but it does not contain the module binaryPCK. Is there a different flavor I need ? Using Python 2.7.x. My initial jplephem install was using pip. I download the tar files and installed that, same issue. any ideas ? with thanks.

[brandon-rhodes]

You should get the new experimental feature if you reinstall from GitHub:

pip install https://github.com/brandon-rhodes/python-jplephem/archive/master.zip