Sunset on [Mountain] from [Viewing Point]

Question

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.

Answer 1

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

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()

Answer 2

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…