initial commit

.gitignore

@@ -0,0 +1,8 @@
+.venv/
+build/
+*.png
+*.so
+*.c
+skydome.c
+*.html
+__pycache__/

app.py

@@ -0,0 +1,73 @@
+import tkinter
+import tkintermapview
+import requests
+import skydome
+from PIL import Image, ImageTk
+
+
+#api altitude endpoint
+URL="https://api.open-elevation.com/api/v1/lookup"
+
+def get_api_params(coords):
+    return {"locations":"{},{}".format(coords[0],coords[1])}
+
+# create tkinter window
+root_tk = tkinter.Tk()
+root_tk.geometry(f"{1000}x{700}")
+root_tk.title("map_view_simple_example.py")
+
+# create map widget
+map_widget = tkintermapview.TkinterMapView(root_tk, width=1000, height=700, corner_radius=0)
+map_widget.pack(fill="both", expand=True)
+
+# set other tile server (standard is OpenStreetMap)
+# map_widget.set_tile_server("https://mt0.google.com/vt/lyrs=m&hl=en&x={x}&y={y}&z={z}&s=Ga", max_zoom=22)  # google normal
+# map_widget.set_tile_server("https://mt0.google.com/vt/lyrs=s&hl=en&x={x}&y={y}&z={z}&s=Ga", max_zoom=22)  # google satellite
+
+
+
+
+def show_image(coords):
+    image_window = tkinter.Toplevel()
+    image_window.title("Simulated Sunset")
+    image_window.config(width=256,height=256)
+    #image = skydome.renderFromCamera(coords)
+    #canvas = tkinter.Canvas(image_window,width=256,height=256)
+    #canvas.pack()
+    #img = ImageTk.PhotoImage(width=256,height=256,image=Image.fromarray(image,mode="RGB"))
+    img = ImageTk.PhotoImage(Image.open("highpress_camera.png"))
+
+    panel = tkinter.Label(image_window,image=img)
+    panel.pack(side="bottom", fill = "both", expand = "yes")
+    image_window.mainloop()
+    #canvas.create_image(20,20, anchor="nw", image=img)
+
+def marker_click(marker):
+    print(f"marker clicked - text: {marker.text}  position: {marker.position}")
+
+
+def add_marker_event(coords):
+    print("marker added to {}".format(coords))
+    map_widget.delete_all_marker()
+    new_marker = map_widget.set_marker(coords[0], coords[1], text="sunset observer")
+
+
+def get_altitude(coords):
+    print("getting altitude for {}\n".format(coords))
+    params = get_api_params(coords)
+    r = requests.get(URL,params)
+    data = r.json()
+    print(data['results'][0]['elevation'])
+
+map_widget.add_right_click_menu_command(label="Set as observer location",
+                                        command=add_marker_event,
+                                        pass_coords=True)
+
+map_widget.add_right_click_menu_command(label="get altitude",
+                                        command=get_altitude,
+                                        pass_coords=True)
+
+map_widget.add_right_click_menu_command(label="open_window",
+                                        command=show_image,
+                                        pass_coords=True)
+root_tk.mainloop()

astropy_test.py

@@ -0,0 +1,21 @@
+from astropy.time import Time
+from astropy.coordinates import solar_system_ephemeris, EarthLocation , AltAz
+from astropy.coordinates import get_body_barycentric, get_body
+
+def get_location(long,lat):
+    return EarthLocation.from_geodetic(long,lat)
+
+loc = get_location(45.508,-73.561)
+
+t = Time("2024-10-07 18:16",location=loc)
+
+with solar_system_ephemeris.set('builtin'):
+    sun = get_body('sun', t, loc)
+    moon = get_body('moon', t, loc)
+    saturn = get_body('saturn', t, loc)
+
+new_sun = sun.transform_to(AltAz(obstime=t,location=loc))
+new_moon = moon.transform_to(AltAz(obstime=t,location=loc))
+new_saturn = saturn.transform_to(AltAz(obstime=t,location=loc))
+
+print(saturn)

requirements.txt

@@ -0,0 +1,2 @@
+tkintermapview~=1.29
+astropy~=6.1.4

setup.py

@@ -0,0 +1,14 @@
+from setuptools import setup 
+from Cython.Build import cythonize
+
+
+setup(
+    ext_modules = cythonize("skydome.py", compiler_directives={"language_level": "3"})
+)
+
+#run following
+#python setup.py build_ext --inplace
+
+#and then
+#python
+#import skydome

skydome.py

@@ -0,0 +1,443 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# https://journals.ametsoc.org/view/journals/bams/79/5/1520-0477_1998_079_0831_opoaac_2_0_co_2.xml?tab_body=pdf - for aerosols and clouds
+# 
+# https://iopscience.iop.org/article/10.1088/1757-899X/1269/1/012011/pdf - for Rayleigh extinction
+# 
+# https://sci-hub.se/10.1364/josa.48.1018_1 - temperature dependance Rayleigh scattering
+# 
+# https://sci-hub.se/https://doi.org/10.1364/JOSA.47.000176 - Rayleigh coefficients (index of refraction) - Penndorf
+# 
+# We are going to use the real rayleigh scattering coefficient: 
+# 
+# $\beta_R(h, \lambda) = \frac{8\pi(n^2 - 1)^2}{2N\lambda^4} e^{-\frac{h}{H_R}}$
+# 
+# n = index of refraction, these are calculated in the 4th paper, temperature dependance
+# 
+# N = molecular density, dependent on humidity
+# 
+# $\lambda$ = wavelength, but we only care for RGB (440, 550, 680 nm)
+# 
+# From Penndorf, the dispersion equation is:
+# 
+# $(n_s - 1) \times 10^{-8} = 6432.8 + \frac{2949810}{146 - v^2} + \frac{25540}{41 - v^2}$
+# 
+# It is important to note that $v=1/f$ which is measured in um
+# 
+# Once we have calculated our dispersion index of refraction, we can see what our new value for n will be (it is dependant on temperature and pressure)
+# 
+# $$
+# (n - 1) = (n_s - 1) \left( \frac{(1 + \alpha t_s)}{(1 + \alpha t)} \right) \frac{p}{p_s}
+# $$
+# 
+# Here, we are using $n_s$ from the dispersion equation, $t_s$ is 15 degrees Celsius, and $p_s$ is 760 mmHg. $t$ represents the actual temperature, and $p$ represents the actual pressure. $\alpha$ is the coefficient of linear expansion of air = 0.00367 in units of inverse celcius. We see that the influence of air pressure really is critical
+# 
+# 
+# We also need to find N, the molecular density (depends on temperature, humidity, pressure): We are going to use the value at the altitude where our person is standing and then we scale it using the scale factor. Penndorf uses N_s = 2.54743e19, but i am not sure if it is applicable to our scenario. 
+# 
+# 
+# 
+# 
+
+# Ok, so things work... but they are not super realistic. I think we should be updating the dispersion equation because Penndorf's paper is super old. We want the value of the index of refraction at the location of the observer - this is going to depend on the users inputs. 
+# 
+# https://refractiveindex.info/?shelf=other&book=air&page=Ciddor (1996)
+# 
+# New calculation for refractive index (where lambda is in microns)
+# 
+# This applies for : Standard air: dry air at 15 °C, 101.325 kPa and with 450 ppm CO2 content.
+# 
+# $n - 1 = \frac{0.05792105}{238.0185 - \lambda^{-2}} + \frac{0.00167917}{57.362 - \lambda^{-2}}$
+# 
+# 
+# 
+
+# Something really cool here is that we can use the color of the sunset to predict the weather: "Red at night, sailors delight. Red in morning, sailors take warning"
+# When there is a bright red sunset at night, it means there is likely a high pressure cell off to the west. With west to east weather patterns, this means a high pressure cell is about to move in, high pressure = pleasant weather, "sailors delight". When there is a red sunrise in the east, it means the high pressure cell has already passed and a low pressure cell is likely to replace it. Low pressure = poor weather, "sailors take warning".
+# 
+# Low pressure: 720mmHg
+# High pressure: 780 mmHg
+
+import numpy as np
+import math
+import imageio
+import time
+import random
+import cython
+
+
+"""
+This class allows us to manipulate vectors
+Not all of the functions are used yet - but they might be later when we start adding things
+"""
+
+
+@cython.cclass
+class Vec3:
+    v = cython.declare(cython.float[3],visibility="public")
+    def __cinit__(self, x:cython.float =0, y:cython.float=0, z:cython.float=0):
+        self.v:cython.float[3] = [x, y, z]
+
+    def __init__(self, x:cython.float =0, y:cython.float=0, z:cython.float=0):
+        self.v:cython.float[3] = [x, y, z]
+
+    def __add__(self, other ):
+        return Vec3(self.v[0] + other.v[0],self.v[1] + other.v[1],self.v[2] + other.v[2])
+
+    def __sub__(self, other):
+        return Vec3(self.v[0] - other.v[0],self.v[1] - other.v[1],self.v[2] - other.v[2])
+
+    def __mul__(self, other):
+        if isinstance(other, Vec3):
+            return Vec3(self.v[0] * other.v[0],self.v[1] * other.v[1],self.v[2] * other.v[2])
+        #elif isinstance(other, (int, float)):
+        return Vec3(self.v[0] * other,self.v[1] * other,self.v[2] * other)
+        #raise TypeError("Can only multiply by a Vec3 or scalar (int or float).")
+
+    def __truediv__(self, scalar):
+        return Vec3(self.v[0] / scalar,self.v[1] / scalar,self.v[2] / scalar)
+
+    def length(self)->cython.float:
+        return math.sqrt(self.v[0]**2 + self.v[1]**2 + self.v[2]**2)
+
+    def normalize(self):
+        length = self.length()
+        return self / length if length > 0 else Vec3()
+
+    def to_tuple(self):
+        return (self.v[0],self.v[1],self.v[2])
+
+    def __repr__(self):
+        return f"Vec3({self.v[0]}, {self.v[1]}, {self.v[2]})"
+
+    @cython.cfunc
+    def dot(self, other):
+        return self.v[0]* other.v[0]+self.v[1]* other.v[1]+self.v[2]* other.v[2]
+
+
+"""
+Here we will be calculating what the value for n is based on T, P, but we will
+also be taking dispersion into account - remember, we are considering 440 nm, 550 nm, 680 nm
+
+These equations match the tabulated results made by Penndorf
+"""
+
+wavelengths = Vec3(0.68, 0.55, 0.44)  # These are the wavelengths in um
+
+# The n_s values gives us the different refractive indices for the different wavelenghts
+# The n values give us the new refractive indices depending on the temperature and pressure
+
+@cython.cfunc
+def refraction_calculator(T:cython.int, P:cython.int) -> tuple[Vec3, Vec3]:
+    alpha:cython.float = 0.00367  # in inverse Celsius
+    t_s:cython.int = 15  # degrees Celsius
+    p_s:cython.int = 760  # mmHg
+    
+    n_s_values = []  # for debugging / comparing to the internet values
+    n_values = []
+    n_s:cython.double
+    n:cython.double
+    i:cython.int
+    for i in [wavelengths.v[0], wavelengths.v[1], wavelengths.v[2]]:
+        n_s = 1 + (0.05792105 / (238.0185 - i**-2) + 0.00167917 / (57.362 - i**-2))
+        
+        n = ((n_s - 1) * ((1 + alpha * t_s) / (1 + alpha * T)) * (P / p_s)) + 1
+        
+        n_s_values.append(n_s)
+        n_values.append(n)
+
+    return Vec3(*n_s_values), Vec3(*n_values)
+
+@cython.cfunc
+def beta(n_values:Vec3):
+    
+    N_s:cython.double = 2.54743e7  # um^-3
+    
+    beta_values = []
+
+    for i, n in enumerate(n_values.v):
+        wavelength = wavelengths.v[i]  # Corresponding wavelength value
+        beta_1 = ((8 * math.pi * (n**2 - 1)**2) / ((2 * N_s) * (wavelength)**4))*1e6  # The whole thing should be in m^-1 to match with the rest of the code
+       
+        beta = (beta_1*10)
+        beta_values.append(beta)
+
+    return Vec3(beta_values[0],beta_values[1],beta_values[2])
+    
+
+# Example usage, normal is 25, 750
+T:cython.int = 5  # Temperature in degrees Celsius
+P:cython.int = 760  # Pressure in mmHg
+n_s_values:Vec3
+n_values:cython.double
+n_s_values, n_values = refraction_calculator(T, P)
+beta_values = beta(n_values)
+
+#print("n_s values:", refraction_calculator_result.v)
+#print("n values:", n_values.v)
+#print("Beta values:", beta_values.v)
+
+
+# We will have a scattering enhancement factor for our scattering values
+
+kInfinity = float('inf')
+angle_degrees:cython.int = 90  # Sun direction in degrees (0 to 90, noon to sunset)
+angle_radians:cython.double = math.radians(angle_degrees)  # Convert to radians
+sunDir:Vec3 = Vec3(0, math.cos(angle_radians), -math.sin(angle_radians))
+
+"""
+The values for betaR depend on many things
+For one, they scale proportionally to barometric pressure
+It also varies according to temperature
+
+"""
+
+betaR = Vec3(*beta_values.v)  # Keeps it as a Vec3
+#betaR = Vec3(3.8e-6, 13.5e-6, 33.1e-6)
+
+# R, G, B (the smaller the value, the less it will scatter during the day) in m-1
+# 3.8 e-6, 13.5 e-6, 33.1 e-6 - default
+# The order should be 680 nm, 550 nm, 440 nm
+
+
+"""
+There are many "categories" of aerosols and each has their own extinction coefficient (I removed the factor of 1.1)
+
+Continental Clean - 26e-6
+Continental Average - 75e-6
+Continental Polluted - 175e-6
+
+Urban - 353e-6
+
+Desert - 145e-6
+
+Maritime Clean - 90e-6
+Maritime Polluted - 115e-6
+Maritime Tropic - 43e-6
+
+Arctic - 23e-6
+
+Antarctic - 11e-6
+
+
+"""
+betaM = Vec3(21e-6, 21e-6, 21e-6) # We need to have the Mie scattering be the same in all the directions
+# The greater the value of beta, the smaller the Mie scattering point (responsable for the halo around the sun)
+# If there is more pollution, we get a larger halo and the colors of the sunset become desaturated (more hazy)
+# Default is 21e-6
+
+
+earthRadius:cython.int = 6360e3 #m
+atmosphereRadius:cython.int = 6420e3 #m (60 km higher than earth)
+Hr:cython.int = 7994
+Hm:cython.int = 1200
+
+# The direction will change for each pixel
+@cython.cfunc
+def computeIncidentLight(direction:Vec3) -> tuple[float,float,float]:
+    tmin:cython.float=0
+    tmax:cython.float=kInfinity
+
+    # We can change the origin position if we want, but for now it is 1 meter above the surface
+    orig:Vec3 = Vec3(0, earthRadius + 1, 0)
+    t0:cython.float
+    t1:cython.float
+    t0, t1 = raySphereIntersect(orig, direction, atmosphereRadius)
+    if t1 < 0:
+        return Vec3(0, 0, 0)
+    if t0 > tmin and t0 > 0:
+        tmin = t0
+    if t1 < tmax:
+        tmax = t1
+
+    numSamples:int = 16
+    numSamplesLight:int = 8
+    segmentLength = (tmax - tmin) / numSamples
+    tCurrent = tmin
+    sumR = Vec3()
+    sumM = Vec3()
+    opticalDepthR = 0
+    opticalDepthM = 0
+    mu:cython.float = direction.dot(sunDir) # This is the cosine of the angle between the direction vector (V) and the sun Direction
+
+    phaseR = (3 * (1 + (mu * mu))) / (16 * math.pi)
+    g = 0.76
+    phaseM = 3 / (8 * math.pi) * ((1 - g * g) * (1 + mu * mu)) / ((2 + g * g) * ((1 + g * g - 2 * g * mu) ** 1.5))
+
+    i:cython.int
+    for i in range(numSamples):
+        samplePosition = orig + direction * (tCurrent + segmentLength * 0.5)
+        height = samplePosition.length() - earthRadius
+        hr = math.exp(-height / Hr) * segmentLength
+        hm = math.exp(-height / Hm) * segmentLength
+        opticalDepthR += hr
+        opticalDepthM += hm
+
+        # Sample position is the start, but it should be in the direction of the sun
+        t0Light, t1Light = raySphereIntersect(samplePosition, sunDir, atmosphereRadius)
+
+        if t1Light < 0:
+            continue
+
+        segmentLengthLight = (t1Light - t0Light) / numSamplesLight
+        tCurrentLight = 0
+        opticalDepthLightR = 0
+        opticalDepthLightM = 0
+
+        j:cython.int
+        for j in range(numSamplesLight):
+            samplePositionLight = samplePosition + (sunDir * (tCurrentLight + segmentLengthLight * 0.5))
+            heightLight = samplePositionLight.length() - earthRadius
+            if heightLight < 0:
+                break
+            opticalDepthLightR += (math.exp(-heightLight / Hr) * segmentLengthLight)
+            opticalDepthLightM += (math.exp(-heightLight / Hm) * segmentLengthLight)
+            tCurrentLight += segmentLengthLight
+
+        if j == numSamplesLight - 1:
+            tau = (betaR * (opticalDepthR + opticalDepthLightR)) + (betaM * (opticalDepthM + opticalDepthLightM))
+            attenuation = Vec3(math.exp(-tau.v[0]), math.exp(-tau.v[1]), math.exp(-tau.v[2]))
+            sumR += (attenuation * hr)
+            sumM += (attenuation * hm)
+        tCurrent += (segmentLength)
+
+    # Changing the * 20 number just changes the intensity os the light, it does not much change the colors themselves
+    # Well the colors kinda change, but more they just get super pale or super not pale. 20 should be fine I guess
+    final_color = (sumR * betaR * phaseR + sumM * betaM * phaseM) * 20 
+    return final_color.to_tuple()
+
+@cython.cfunc
+def raySphereIntersect(orig:Vec3, direction:Vec3, radius:cython.double) -> tuple[cython.float,cython.float]:
+    A:cython.float = direction.v[0]**2 + direction.v[1]**2 +direction.v[2]**2
+    B:cython.float = 2 * (direction.v[0]*orig.v[0] +direction.v[1]*orig.v[1] +direction.v[2]*orig.v[2] )
+    C:cython.float = orig.v[0]**2 + orig.v[1]**2 +orig.v[2]**2 - radius ** 2
+
+    t0: cython.float
+    t1: cython.float
+    t0, t1 = solveQuadratic(A, B, C)
+    if t1 < 0:
+        return (kInfinity, kInfinity)
+    #if t0 > t1:
+    #    t0, t1 = t1, t0
+    return (t0, t1)
+
+@cython.cfunc
+def solveQuadratic(a:cython.float, b:cython.float, c:cython.float) -> tuple[cython.float,cython.float]:
+    if b == 0:
+        if a == 0:
+            return (0, 0)
+        x1 = 0
+        x2 = math.sqrt(-c / a)
+        return (x1, x2)
+    discr = b ** 2 - 4 * a * c
+    if discr < 0:
+        return (kInfinity, kInfinity)
+
+    discr:cython.float
+    q:cython.float
+    discrs = math.sqrt(discr)
+    q = (-0.5 * (b - discrs)) if b < 0 else (-0.5 * (b + discrs))
+    return (q/a, c/q)
+
+
+
+@cython.cfunc
+def renderSkydome(filename):
+    width:cython.int
+    height:cython.int
+    width, height = 256,256
+    image = np.zeros((height, width, 3), dtype=float)
+
+    start_time = time.time()  # Start timing
+
+    j:cython.int
+    i:cython.int
+    x:cython.float
+    y:cython.float
+    z2:cython.float
+    phi:float
+    theta:float
+    direction:Vec3
+    for j in range(height):
+        for i in range(width):
+            x = 2 * (i + 0.5) / (width - 1) - 1
+            y = 2 * (j + 0.5) / (height - 1) - 1
+            z2 = x * x + y * y
+            if z2 <= 1:
+                phi = math.atan2(y, x)
+                theta = math.acos(1 - z2)
+                # This changes for each pixel
+                direction = Vec3(math.sin(theta) * math.cos(phi), math.cos(theta), math.sin(theta) * math.sin(phi))
+                color = computeIncidentLight(direction)
+                #color = computeIncidentLight(direction)
+
+
+                # Assign the clipped color directly to the image array
+                image[j][i][0] = np.clip(color, 0, 1)[0]
+                image[j][i][1] = np.clip(color, 0, 1)[1]
+                image[j][i][2] = np.clip(color, 0, 1)[2]
+
+        # Print elapsed time after each row
+        elapsed_time = time.time() - start_time
+        print(f"Rendering row {j + 1}/{height}, elapsed time: {elapsed_time:.2f} seconds")
+        #print(f"Rendering row {j + 1}/{height}")
+
+    # Save result to a PNG image
+    image = np.clip(image, 0, 1) * 255 # change 255
+    imageio.imwrite(filename, image.astype(np.uint8))
+
+
+
+def renderFromCamera(filename):
+
+    width:cython.int
+    height:cython.int
+    width, height = 100, 100
+
+    x:cython.int
+    y:cython.int
+    rayx:cython.float
+    rayy:cython.float
+    image = np.zeros((height, width, 3), dtype=np.float32)
+    
+    aspectRatio = width / float(height)
+    fov = 65 # field of view
+    angle = math.tan(fov * np.pi / 180 * 0.5) # How much the camera can see
+
+    start_time = time.time()  # Start timing
+
+    numPixelSamples:cython.int = 4
+    for y in range(height):
+        for x in range(width):
+#            for m in range(numPixelSamples):
+#                for n in range(numPixelSamples):
+#                    rayx = (2 * (x + (m + random.uniform(0, 1)) / numPixelSamples) / float(width) - 1) * aspectRatio * angle
+#                    rayy = (1 - (y + (n + random.uniform(0, 1)) / numPixelSamples) / float(height) * 2) * angle
+
+            rayx = (2 * x  / float(width) - 1) * aspectRatio * angle
+            rayy = (1 - y / float(height) * 2) * angle
+
+            #This changes for each pixel, it is the direction we are looking in
+            direction = Vec3(rayx, rayy, -1).normalize()
+            color = computeIncidentLight(direction)
+            #image[y, x] = np.array(color)
+            image[y][x] = np.clip(color, 0, 1)
+
+            #image[y, x] /= numPixelSamples * numPixelSamples  # Average color
+
+        elapsed_time = time.time() - start_time
+        print("what")
+        print(f"Renderinggg row {y + 1}/{height}, elapsed time: {elapsed_time:.2f} seconds")  
+
+    image = np.clip(image, 0, 1) * 255
+    return image
+    #imageio.imwrite(filename, image.astype(np.uint8))
+
+    #renderFromCamera("camera_render.png")
+
+
+#renderSkydome("highpress_15C.png")
+#renderFromCamera("highpress_camera.png")
+
+