File:Twilight illumination by sun altitude 1 r 1.svg

-20,  0.000004,    0.0014,    0.0008, 564.3, 601.7
-19,  0.000005,    0.0017,    0.0009, 563.8, 605.7
-18,  0.000005,    0.0016,    0.0009, 563.5, 604.0
-17,  0.000005,    0.0016,    0.0008, 563.3, 603.0
-16,  0.000003,    0.0010,    0.0005, 563.0, 593.0
-15,  0.000004,    0.0017,    0.0008, 559.0, 570.7
-14,  0.000006,    0.0031,    0.0011, 552.8, 553.9
-13,  0.000011,    0.0058,    0.0016, 546.8, 541.0
-12,  0.000020,    0.0106,    0.0026, 542.2, 535.9
-11,  0.000047,    0.0247,    0.0060, 540.7, 540.4
-10,  0.000117,    0.0650,    0.0147, 537.6, 531.6
 -9,  0.000273,    0.1564,    0.0347, 536.7, 526.0
 -8,  0.000761,    0.4549,    0.1016, 536.6, 520.4
 -7,  0.002138,    1.3136,    0.2970, 536.8, 517.2
 -6,  0.006553,    4.0423,    0.9328, 537.5, 518.1
 -5,  0.018006,   11.1915,    2.6195, 537.9, 518.0
 -4,  0.052251,   32.5612,    7.7281, 538.4, 519.1
 -3,  0.146334,   90.3824,   22.2753, 539.8, 525.1
 -2,  0.352731,  218.6244,   55.4249, 540.7, 526.1
 -1,  0.763037,  474.1386,  125.0604, 542.1, 526.6
  0,  1.580555,  962.4795,  272.0179, 544.5, 532.2
  1,  2.772255, 1666.0518,  500.7170, 546.5, 537.4
  2,  4.411440, 2523.4746,  815.0906, 549.0, 548.1
  3,  6.702043, 3671.5632, 1254.1349, 551.0, 555.6
  4,  9.416510, 5084.9744, 1773.6052, 551.7, 557.2
  5, 12.865984, 7436.9430, 2507.0312, 550.2, 545.2
  6, 17.437111, 9649.4890, 3420.7512, 551.9, 553.3
  7, 20.121992,11071.6979, 3965.8271, 552.3, 553.7
  8, 26.330863,13829.1131, 5235.3776, 554.1, 562.6
  9, 45.112879,21760.3774, 9179.4651, 557.6, 576.3
 10, 36.358187,18189.3212, 7454.5144, 556.3, 572.2
 11, 39.137271,20176.6702, 8084.3033, 555.5, 567.0
 12, 61.641958,30708.4968,12903.1196, 557.0, 573.3
 13, 77.212419,37803.7065,16239.1084, 557.7, 576.0
 14, 90.655367,43652.7913,18982.5532, 558.3, 577.2
 15, 93.396421,45058.8386,19410.2871, 558.3, 575.2
 16,112.860479,53569.4518,23508.1913, 558.8, 578.0
 17,121.886868,60050.7316,26070.7039, 557.8, 575.6
 18,132.862405,65464.2749,28434.6294, 557.9, 575.4
 19,120.301796,59054.1295,25659.3754, 557.8, 576.5

1. print(df.to_string())

sunangle=np.array(df['sunangle']) photlx=np.array(df['photlx'])

logphotlx=np.log10(photlx)

1. print(sunangle)

x=sunangle y=logphotlx

x2 = np.linspace(-18, 0, 360) fit1=np.polyfit(x, y, 12) print(fit1) ## paste output of this to plot code to poly fit!

1. 1. illumination on twilight, sun 0-18 degrees below horizon
   2. python3 script
   3. print also moon illumination lux lx brightness
   4. uses ephem, pylunar,
3. fitted photlx from
   1. https://physics.stackexchange.com/questions/410633/what-are-the-average-wavelengths-and-brightnesses-of-sunlight-across-the-stages
5. 15.8.2023 0000.0001

import math import numpy as np import matplotlib.pyplot as plt import ephem import pylunar import re

def lx_from_frac_single(frac1):

   frac2=frac1*100
   frac_fit1=[-2.02393074e-07 , 5.52148721e-05, -5.23702354e-03,  2.19919440e-01, -5.00000220e+00]
   yn = np.poly1d(frac_fit1)
   y2=yn(frac2) 
   lux=np.power(10,y2)
   return(lux)
   

def lx_from_phase_deg_single(phasedeg1):

   phase_fit1=[-9.98825150e-09,  2.17143911e-06 ,-1.73248246e-04, -6.45885492e-03,-4.02728954e-01]
   yn = np.poly1d(phase_fit1)
   y2=yn(phasedeg1) 
   lux=np.power(10,y2)
   return(lux)

def dms_dd(dd):

   dd = f"""{dd}"""
   dd = re.sub('[^a-zA-Z0-9. ]', ' ', dd)
   dd = dd.split(" ")
   return round(float(dd[0])+(float(dd[1])/60)+(float(dd[2])/3600), 8)

def dms_list(dd):

   dd = f"""{dd}"""
   dd = re.sub('[^a-zA-Z0-9. ]', ' ', dd)
   dd2 = np.array(dd.split(" "))
   dd3=[]
   for n  in range(0,3):
       a=dd2[n]
       dd3.append(float(a))
       #print(a)
   
   if (dd2[3]=='S'): dd3[0]=dd3[0]*-1    
   if (dd2[3]=='W'): dd3[0]=dd3[0]*-1  
   return dd3

1. 1. tulilahti
2. 62°26'06.3"N 28°34'41.7"E
3. 62.435079, 28.578249
   1. tulilahti

latt1="62°26'06.3\"N" lonn1="28°34'41.7\"E" sdatetime1= (1959, 7, 27, 24-2, 0, 0) ##MUST be UTC! EET is UTC+2h

latt2=dms_dd(latt1) ## dec latt3=dms_list(latt1) ## dms list lonn2=dms_dd(lonn1) ## dec lonn3=dms_list(lonn1) ## dms list

print(latt2) print(latt3)

1. latt3=re.split('[°\'"]+', latt1)

1. print(latt3)

1. quit(-1)

1. Define the observer's location

observer = ephem.Observer()

obs_name="Tulilahti, Heinävesi"

1. observer.lat = '62.434967'
2. observer.long = '28.578223'

observer.lat = str(latt2) observer.long = str(lonn2)

1. mi = pylunar.MoonInfo((62, 25, 35), (28, 37, 50))

mi = pylunar.MoonInfo(latt3, lonn3)

1. mi.update((1959, 7, 27, 22, 0, 0)) ## needs utc

mi.update(sdatetime1) ## needs utc

## !!! date, time  must be in UTC, that is eg. helsinki time-3, in summertime helsinki time-3

1. observer.date=ephem.Date((1959, 7, 27, 22, 0, 0))

observer.date=ephem.Date(sdatetime1)

moon = ephem.Moon() moon.compute(observer) sun = ephem.Sun() sun.compute(observer)

print(obs_name) print("Paikka pituusaste leveysaste ", observer.long ,observer.lat ) print("Aika UTC ", observer.date )

print(" Auringon korkeuskulma ", sun.alt) print(" Auringon suuntakulma ", sun.az)

print(" Kuun korkeuskulma ", moon.alt) print(" Kuun suuntakulma " , moon.az)

1. Print the phase of the moon

print(" Kuun vaihe ", moon.moon_phase)

1. 1. altitude of sun, recommend 0 ... -18 degrees

alt1=sun.alt alt3=float(repr(sun.alt))*180/math.pi

1. print(alt2)

1. quit(-1)

1. 1. params for poly1d, fitted from data

1. 1. new angles of sun for plot

x2 = np.linspace(-20, 5, 400)

fit2=[ 1.14909168e-15, 5.05406644e-14, -2.31959774e-12 ,-5.29563006e-11,

 1.57147669e-09 , 2.25671280e-08, -5.52103489e-07, -5.24010954e-06,
 1.35247027e-04 , 3.17102665e-04, -2.32199880e-02 , 2.91689849e-01,
 2.42028605e+00]
 

yn = np.poly1d(fit2) y2=yn(x2)

photlx2=np.power(10,y2)

logphotlx3=yn(alt3) lux3=math.pow(10,logphotlx3 ) print(" Sun altitude deg ",alt3) print(" Estimated illumination lx ",lux3)

plt.xlabel("Sun altitude degrees", size=18) plt.ylabel("Log10 of illumination lx", size=18)

plt.yticks(fontsize=18) plt.xticks(fontsize=18)

plt.title(" Twilight illuminaltion ", fontsize=22)

1. 1. lat, lon tulilahti

mi.age()

mi.fractional_phase()

1. print(mi.is_visible())

1. print(mi.phase_name())

moonalt1=mi.altitude() moonaz1=mi.azimuth()

1. print(mi.magnitude())

mmoon1=mi.magnitude() frac1= mi.fractional_phase() alt1=mi.altitude() angsiz1=mi.angular_size() az1=mi.azimuth() distance1=mi.earth_distance() rise_set_ts1=mi.rise_set_times("EET")

rel_distance1= distance1/384399 magnitude1=mi.magnitude()

print(alt1) print(frac1) print(rel_distance1) print(magnitude1) print(rise_set_ts1)

mlux1=lx_from_frac_single(frac1)

1. 1. from hat

1. tau=.5
   1. alt 10 0,4519
   2. alt 30 0.801

tau=.1828

z=(math.pi/180)*(90-alt1) zdeg=90-alt1

1. 1. air mass
2. X=1/math.cos(z)
3. https://en.wikipedia.org/wiki/Air_mass_(solar_energy)

X=1/(math.cos(z)+0.50572*math.pow((96.07995-zdeg), -1.6364))

1. print(X)

E=math.exp(-tau*X*math.cos(1/z))

print(alt1, E)

mlux2=mlux1*E mlux3=mlux2*1/math.pow(rel_distance1,2)

print(" Moon lux from phase ", round(mlux1,6)) print(" Moon lux from altitude ", round(mlux2,6)) print(" Moon lux from distance ", round(mlux3,6))

1. 1. sun -26.9

msun1=-26.832 insun1=-26.832*0.4

1. 1. mean full moon -12.74
   2. first quarter 1/11
   3. last quarter 1/12 #3 die to dark mare dirrerent position on moon surface

1. 1. first quarter -11.15

1. 1. full moon 0.25 lx

inmoon1=mmoon1*0.4

inmoon2=insun1-inmoon1

inmoon3=math.pow(10, inmoon2)*100000

1. moon params

print("Kuun korkeuskulma ", moonalt1) print("Kuun suuntakulma ", moonaz1) print("Kuun suuruusluokka ", mmoon1) print("Kuun valo ", round(mlux3,6)) ## estimeted moon illum lux

1. print("Kuunvalo lx ", inmoon3)

illum_total_est_1=inmoon3+lux3

1. print("Valaistus yhteensä ", illum_total_est_1)

plt.xlim([-20,5]) plt.xlim([-20,5])

1. plt.plot(x2, photlx2)

plt.plot(x2, np.log10(photlx2) , lw=2, color="#0000ff")

1. plt.scatter(alt3,logphotlx3 )

plt.axvline(x = 0, color = '#FFFF00', label = 'Sunset') plt.axvline(x = -6, color = '#7f7f00', label = 'Civil twilight ends') plt.axvline(x = -12, color = '#3f3f00', label = 'Nautical twilight ends') plt.axvline(x = -18, color = '#1f1f00', label = 'Astronomical twilight ends')

1. plt.plot(x2,y2)

1. 1. horizon illumination
   2. https://fi.wikipedia.org/wiki/H%C3%A4m%C3%A4r%C3%A4
2. plt.axhline(y = math.log10(442.5), color = '#FFFF00', label = 'Sunset') ## CT horiz illum lux
3. plt.axhline(y = math.log10(2.75), color = '#FFFF00', label = 'CT') ## CT horiz illum lux
4. plt.axhline(y = math.log10(0.008), color = '#FFFF00', label = 'NT') ## CT horiz illum lux
5. plt.axhline(y = math.log10(6e-4), color = '#FFFF00', label = 'AT') ## CT horiz illum lux
   1. full moon teor max 0.4lx, ofren 0.1 ... 0.25 lx
6. plt.axhline(y = math.log10(0.25), color = '#FFFF00', linestyle=":", label = 'FullMoon') ## Full Moon illum lux
7. plt.axhline(y = math.log10(0.1), color = '#FFFF00', linestyle=":", label = 'FullMoon') ## Full Moon illum lux
8. https://travislongcore.net/2017/08/06/how-bright-the-moon-correcting-a-propagated-figure-error-in-the-literature/

1. 1. half moon max 0.043, if alt 90 deg, 0.004 if alt 10 deg

1. 1. half moon
2. plt.axhline(y = math.log10(0.01), color = '#FFFF00', linestyle=":", label = 'FullMoon') ## Full Moon illum lux
   1. quarter crescent moon, if full moon 0.345, altitude 90 deg
3. https://travislongcore.net/2017/08/06/how-bright-the-moon-correcting-a-propagated-figure-error-in-the-literature/
4. brown 1952
5. plt.axhline(y = math.log10(0.0124), color = '#FFFF00', linestyle=":", label = 'FullMoon') ## Full Moon illum lux
   1. Krisciunas & Schaefer 1991
6. plt.axhline(y = math.log10(0.008), color = '#FFFF00', linestyle=":", label = 'FullMoon') ## Full Moon illum lux

1. 1. quarter crescent moon, if full moon 0.1
2. plt.axhline(y = math.log10(0.003), color = '#FFFF00', linestyle=":", label = 'FullMoon') ## Full Moon illum lux

1. 1. if moon under 10 above horiz, lx ca lx*0.1

1. 1. starlight
2. plt.axhline(y = math.log10(0.001), color = '#FFFF00', linestyle=":", label = 'FullMoon') ## Full Moon illum lux

plt.show()