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.
Astronomy Stack Exchange is a question and answer site for astronomers and astrophysicists. It only takes a minute to sign up.
Sign up to join this communityHow 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.
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()
.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$
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…
2019-04-07 17:49:00
. But now I'm worried that that site's scripts don't consider earths varying radius.
$\endgroup$