3
$\begingroup$

How do I get the days in this year when the sun will set exactly on top of Fuji as viewed from Enoshima.

I shamefully admit that I just want to take a nice picture of Fuji.

$\endgroup$
  • 1
    $\begingroup$ I just ran across this article and thought you might find its XKCD wisdom fun/interesting gizmodo.com/… $\endgroup$ – uhoh Sep 5 '19 at 23:26
2
$\begingroup$

You would like to do some Viewing “Diamond Fuji”!

Using Python and Skyfield and the GPS coordinates in your urls:

Mt_Fuji  = earth + Topos(latitude_degrees    =  +35.36304,
                         longitude_degrees   = +138.73040,
                         elevation_m         =  3776.0)

Enoshima = earth + Topos(latitude_degrees    =  +35.29875,
                        longitude_degrees    = +139.47457,
                        elevation_m          =  10.0)

I get the following. I make no guarantees!

Mt. Fuji peak from Enoshima
azimuth:     276.23 
altitude       2.86

Sunset over Mt. Fuji from Enoshima

Sunset over Mt. Fuji from Enoshima

import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Loader, Topos

load  = Loader('~/Documents/fishing/SkyData')  # single instance for big files
ts    = load.timescale()
de421 = load('de421.bsp')
sun   = de421['sun']
earth = de421['earth']
moon  = de421['moon']

Mt_Fuji  = earth + Topos(latitude_degrees    =  +35.36304,
                         longitude_degrees   = +138.73040,
                         elevation_m         =  3776.0)

Enoshima = earth + Topos(latitude_degrees    =  +35.29875,
                        longitude_degrees    = +139.47457,
                        elevation_m          =  10.0)

hours  = 17 - 9 + np.arange(91)/60.

days   = np.arange(1, 367)

alt_Fuji, az_Fuji, d_Fuji = Enoshima.at(ts.now()).observe(Mt_Fuji).apparent().altaz()
alt_Fuji, az_Fuji         = [thing.degrees for thing in (alt_Fuji, az_Fuji)]
d_Fuji_km                  = d_Fuji.km

Mt_Fuji_obs = Enoshima.at(ts.now()).observe(Mt_Fuji).apparent()

if True:
    seps = []
    for day in days:
        times = ts.utc(2019, 1, day, hours)
        sunpos  = Enoshima.at(times).observe(sun).apparent()
        Fujipos = Enoshima.at(times).observe(Mt_Fuji)
        sep     = Fujipos.separation_from(sunpos)
        seps.append(sep)

    sepz = [x.degrees for x in seps]
    SEP  = np.array(sepz)

    if True:
        plt.figure()
        plt.imshow(SEP, vmin=0, vmax=5)
        plt.colorbar()
        plt.xlabel('minutes after 17:00 JST', fontsize=14)
        plt.ylabel('day number in 2019 JST', fontsize=14)
        plt.title('Sun sep (deg) from Mt. Fuji from Enoshima', fontsize=14)
        plt.show()

# make a detailed plot
if True:
    days_1   = np.arange( 95, 100)   # april 5 thru 9
    days_2   = np.arange(246, 251)   # sept  3 thru 7
    both = []
    for days in (days_1, days_2):
        altazs = []
        for day in days:
            times = ts.utc(2019, 1, day, hours)
            alt, az, d = Enoshima.at(times).observe(sun).apparent().altaz()
            alt, az    = [thing.degrees for thing in (alt, az)]
            altazs.append((alt, az))
        both.append(altazs)

    if True:
        hw_deg         = 5.0
        altmin, altmax = alt_Fuji - hw_deg, alt_Fuji + hw_deg
        azmin,  azmax  = az_Fuji  - hw_deg, az_Fuji  + hw_deg
        xFuji = [az_Fuji - 2*alt_Fuji, az_Fuji,  az_Fuji + 2*alt_Fuji, az_Fuji - 2*alt_Fuji]
        yFuji = [0,                    alt_Fuji, 0,                    0                   ]
        plt.figure()
        for i, altazs in enumerate(both):
            plt.subplot(2, 1, i+1)
            for (alt, az) in altazs:
                plt.plot(az, alt)
            plt.plot(xFuji, yFuji, '-k', linewidth=2)
            plt.plot([azmin, azmax], [0, 0], '-k')
            plt.xlim(azmin, azmax)
            plt.ylim(altmin, altmax)
            plt.ylabel('altitude(deg)', fontsize=14)

        plt.xlabel('azimuth (deg)', fontsize=14)
        plt.suptitle('Sunset vs Mt. Fuji from Enoshima', fontsize=14)
        plt.show()
$\endgroup$
  • 1
    $\begingroup$ Pointing me to skyfield would have been enough for an accepted answer, so… thanks for doing all the footwork. Sadly, it was pretty hazy, so I didn't even go to to Enoshima but took pictures from a bit further away/more north. i.stack.imgur.com/JIiLp.jpg One question though: wouldn't it be better to call .altaz('standard') to have at least some correction for atmospheric refraction? Not that it makes a big difference. $\endgroup$ – Caesar Apr 6 '19 at 13:36
  • 1
    $\begingroup$ @Caesar Oh that is beautiful, thank you so much for sharing the photo! :-) The idea behind Stack Exchange is to write answers that not only help the OP, but will possibly be helpful to future readers as well. Perhaps someone else decides to learn python, or adopt Skyfield by seeng what it can do. It would be great if you posted your photo in a new answer, calculated for the new lat-lon, and showed the comparison, of course in your "spare time" ;-) $\endgroup$ – uhoh Apr 6 '19 at 23:36
  • $\begingroup$ @Caesar About .atlaz(), I started using Skyfield from an early version, and have never read through the documentation. I'd thought that standard was the default, even going so far as to use .altaz(pressure_mbar=0) to turn off /refraction here. So posting the script not only keeps me in practice but brings in helpful comments like yours. Thanks! $\endgroup$ – uhoh Apr 6 '19 at 23:38
3
$\begingroup$

My attempt at an answer of my own question. I have no confidence that it is correct, and I'd like to know whether there's an easier way to do this…

https://www.movable-type.co.uk/scripts/latlong.html tells me that the bearing is 90.3°.

A quick application of middle school trigonometry tells me that the visible elevation of Fuji is about 3413 m or 2.90°.

Next, I went to Stellarium, fixed my position to Enoshima, and exported the Ephemeris (F10 menu) of the Sun in a 5 minute intervals with horizontal coordinates from February till May.

Finally, I used a quick python script to find the line in the Ephimeris export where the sun is closest to my desired coordinates (90.3°, 2.9°)

import re
from pprint import pprint

height = "Height"
azi = "Azimut"

# Quick and dirty parser
pre = re.compile("([+-])(\\d+)°(\\d+)'([0-9.]+)\"")
def ph(v):
    m = pre.match(v)
    if m:
        sign, *nums = m.groups()
        return sum([float(s) * 60 ** (-e) for e,s in enumerate(nums)]) * float(sign + "1")
    try:
        return float(v)
    except:
        pass
    return v

res = []
with open('ephemeris2.csv') as f:
    head = next(f).split(", ")
    for line in f:
        php = [ph(v) for v in line.split(", ")]
        res.append(dict(zip(head, php)))

# This is where the magic happens!
filt = [k for k in res if k is not None and height in k and azi in k]
filt.sort(key = lambda v: (v[height] - 2.9) ** 2 + (v[azi] - 270.3) ** 2)
pprint(filt[:5])

This yielded me

{ 'Azimut': 270.13194444444446,
  'Datum and Time': '2019-03-25 17:39:00',
  'Height': 3.0792222222222225 }

I'm unsure whether I can just convert the Bearing to Azimuth like that, but I suspect that a simple mathematical mistake is even more likely…

$\endgroup$
  • 1
    $\begingroup$ When I checked the azimuth of Mount Fuji as seen from Enoshima, I got a value of 276.2°, a bit different from yours. Which coordinates did you use for Enoshima and Mount Fuji? The elevation angle seems OK, I got 2.86° $\endgroup$ – FSimardGIS Mar 24 '19 at 17:53
  • $\begingroup$ Hm, you're right. I thought I checked the positions on the map… Anyway, the date changes to 2019-04-07 17:49:00. But now I'm worried that that site's scripts don't consider earths varying radius. $\endgroup$ – Caesar Mar 25 '19 at 14:23
  • 1
    $\begingroup$ I'm not sure I understand what you mean, do you refer to the fact that Earth is a spheroid (squashed sphere)? I have done the calculation using the WGS84 spheroid and 2.9° is the correct elevation angle for the line of sight. The bearing is good as well. I think April 7 would be a good bet, too. The biggest uncertainty will be refraction on that day, as well as the weather :-) $\endgroup$ – FSimardGIS Mar 25 '19 at 16:33
  • $\begingroup$ Hi @Caesar happy to see an other Python advocate! Try this out if you get a chance. $\endgroup$ – uhoh Mar 27 '19 at 12:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.