diff --git a/Shaders/std/sky.glsl b/Shaders/std/sky.glsl
index 921eb7e8..6b85b4cc 100644
--- a/Shaders/std/sky.glsl
+++ b/Shaders/std/sky.glsl
@@ -29,6 +29,8 @@
 #define nishita_mie_dir 0.76 // Aerosols anisotropy ("direction")
 #define nishita_mie_dir_sq 0.5776 // Squared aerosols anisotropy
 
+#define sun_limb_darkening_col vec3(0.397, 0.503, 0.652)
+
 /* ray-sphere intersection that assumes
  * the sphere is centered at the origin.
  * No intersection when result.x > result.y */
@@ -135,4 +137,19 @@ vec3 nishita_atmosphere(const vec3 r, const vec3 r0, const vec3 pSun, const floa
 	return nishita_sun_intensity * (pRlh * nishita_rayleigh_coeff * totalRlh + pMie * nishita_mie_coeff * totalMie);
 }
 
+vec3 sun_disk(const vec3 n, const vec3 light_dir, const float disk_size, const float intensity) {
+	// Normalized SDF
+	float dist = distance(n, light_dir) / disk_size;
+
+	// Darken the edges of the sun
+	// Reference: https://media.contentapi.ea.com/content/dam/eacom/frostbite/files/s2016-pbs-frostbite-sky-clouds-new.pdf
+	// (Physically Based Sky, Atmosphere and Cloud Rendering in Frostbite by Sebastien Hillaire)
+	// Page 28, Page 60 (Code from [Nec96])
+	float invDist = 1.0 - dist;
+	float mu = sqrt(invDist * invDist);
+	vec3 limb_darkening = 1.0 - (1.0 - pow(vec3(mu), sun_limb_darkening_col));
+
+	return 1 + (1.0 - step(1.0, dist)) * nishita_sun_intensity * intensity * limb_darkening;
+}
+
 #endif
diff --git a/blender/arm/material/cycles_nodes/nodes_texture.py b/blender/arm/material/cycles_nodes/nodes_texture.py
index 4aef494c..b095785d 100644
--- a/blender/arm/material/cycles_nodes/nodes_texture.py
+++ b/blender/arm/material/cycles_nodes/nodes_texture.py
@@ -1,3 +1,4 @@
+import math
 import os
 from typing import Union
 
@@ -382,7 +383,25 @@ def parse_sky_nishita(node: bpy.types.ShaderNodeTexSky, state: ParserState) -> v
     # d_ozone = node.ozone_density
     density = c.to_vec2((d_air, d_dust))
 
-    return f'nishita_atmosphere(n, vec3(0, 0, {ray_origin_z}), sunDir, {planet_radius}, {density})'
+    sun = ''
+    if node.sun_disc:
+        # The sun size is calculated relative in terms of the distance
+        # between the sun position and the sky dome normal at every
+        # pixel (see sun_disk() in sky.glsl).
+        #
+        # An isosceles triangle is created with the camera at the
+        # opposite side of the base with node.sun_size being the vertex
+        # angle from which the base angle theta is calculated. Iron's
+        # skydome geometry roughly resembles a unit sphere, so the leg
+        # size is set to 1. The base size is the doubled normal-relative
+        # target size.
+
+        # sun_size is already in radians despite being degrees in the UI
+        theta = 0.5 * (math.pi - node.sun_size)
+        size = math.cos(theta)
+        sun = f'* sun_disk(n, sunDir, {size}, {node.sun_intensity})'
+
+    return f'nishita_atmosphere(n, vec3(0, 0, {ray_origin_z}), sunDir, {planet_radius}, {density}){sun}'
 
 
 def parse_tex_environment(node: bpy.types.ShaderNodeTexEnvironment, out_socket: bpy.types.NodeSocket, state: ParserState) -> vec3str: