import bpy
import math
from mathutils import Vector

# Blender 4.0.2 diagnostic scene: textured disk + latitude/longitude grid only.
# No stars, labels, camera, lights, or animation; one static Sun is added below.
IMAGE_PATH = ""
DISK_RADIUS = 20000.0
GRID_Z = 60.0
LINE_RADIUS = 18.0

# Clear current scene
bpy.ops.object.select_all(action='SELECT')
bpy.ops.object.delete(use_global=False)
for datablocks in (bpy.data.meshes, bpy.data.curves, bpy.data.materials, bpy.data.cameras, bpy.data.lights):
    for datablock in list(datablocks):
        if datablock.users == 0:
            datablocks.remove(datablock)

# Create circular disk mesh
segments = 256
verts = [(0.0, 0.0, 0.0)]
faces = []
for i in range(segments):
    a = 2.0 * math.pi * i / segments
    verts.append((DISK_RADIUS * math.cos(a), DISK_RADIUS * math.sin(a), 0.0))
for i in range(segments):
    faces.append((0, i + 1, ((i + 1) % segments) + 1))
mesh = bpy.data.meshes.new("Earth disk mesh")
mesh.from_pydata(verts, [], faces)
mesh.update()
disk = bpy.data.objects.new("Earth disk | diameter 40000", mesh)
bpy.context.collection.objects.link(disk)

# UV coordinates for the AE.jpg map
uv_layer = mesh.uv_layers.new(name="MapUV")
for poly in mesh.polygons:
    for loop_index in poly.loop_indices:
        vi = mesh.loops[loop_index].vertex_index
        x, y, _ = mesh.vertices[vi].co
        # Rotate the azimuthal-equidistant map 90 degrees clockwise on the disk.
        # UV mapping around image center: (u, v) -> (v, 1 - u)
        u = 0.5 + x / (2.0 * DISK_RADIUS)
        v = 0.5 + y / (2.0 * DISK_RADIUS)
        uv_layer.data[loop_index].uv = (v, 1.0 - u)

mat = bpy.data.materials.new("Earth map material")
mat.use_nodes = True
nodes = mat.node_tree.nodes
links = mat.node_tree.links
bsdf = nodes.get("Principled BSDF")
if bsdf:
    bsdf.inputs["Roughness"].default_value = 1.0
try:
    image = bpy.data.images.load(IMAGE_PATH, check_existing=True)
    tex = nodes.new("ShaderNodeTexImage")
    tex.image = image
    links.new(tex.outputs["Color"], bsdf.inputs["Base Color"])
except Exception as e:
    print("Could not load map image:", IMAGE_PATH, e)
    mat.diffuse_color = (0.12, 0.35, 0.15, 1.0)
disk.data.materials.append(mat)

# Grid material (viewport and rendered appearance)
grid_mat = bpy.data.materials.new("Latitude longitude grid material")
grid_mat.diffuse_color = (0.02, 0.02, 0.02, 1.0)
grid_mat.use_nodes = True
p = grid_mat.node_tree.nodes.get("Principled BSDF")
if p:
    p.inputs["Base Color"].default_value = (0.015, 0.015, 0.015, 1.0)
    p.inputs["Roughness"].default_value = 1.0

# Curves are deliberately simple and lie just above the disk.
def add_polyline(name, points, cyclic=False):
    curve = bpy.data.curves.new(name + " curve", type='CURVE')
    curve.dimensions = '3D'
    curve.resolution_u = 1
    curve.bevel_depth = LINE_RADIUS
    curve.bevel_resolution = 1
    spline = curve.splines.new('POLY')
    spline.points.add(len(points) - 1)
    for pnt, co in zip(spline.points, points):
        pnt.co = (co[0], co[1], co[2], 1.0)
    spline.use_cyclic_u = cyclic
    obj = bpy.data.objects.new(name, curve)
    bpy.context.collection.objects.link(obj)
    obj.data.materials.append(grid_mat)
    return obj

# Latitude circles for an azimuthal-equidistant map centered on the North Pole.
for lat in [80, 60, 40, 20, 0, -20, -40, -60, -80]:
    r = DISK_RADIUS * (90.0 - lat) / 180.0
    pts = []
    for i in range(361):
        a = 2.0 * math.pi * i / 360.0
        pts.append((r * math.cos(a), r * math.sin(a), GRID_Z))
    add_polyline("Latitude %s deg" % lat, pts, cyclic=False)

# Longitude rays every 15 degrees, from the North Pole to the disk edge.
for lon in range(0, 360, 15):
    a = math.radians(lon)
    pts = []
    for i in range(201):
        r = DISK_RADIUS * i / 200.0
        pts.append((r * math.cos(a), r * math.sin(a), GRID_Z))
    add_polyline("Longitude %s deg" % lon, pts)


# Test addition: a single static Sun sphere plus three static star spheres. No camera, lights, or animation.
SUN_RADIUS = 300.0
TROPIC_OF_CANCER_LATITUDE = 23.4366
SUN_ORBIT_RADIUS = DISK_RADIUS * (90.0 - TROPIC_OF_CANCER_LATITUDE) / 180.0
SUN_LOCATION = (SUN_ORBIT_RADIUS, 0.0, 5000.0)
bpy.ops.mesh.primitive_uv_sphere_add(segments=32, ring_count=16, radius=SUN_RADIUS, location=SUN_LOCATION)
sun = bpy.context.object
sun.name = "Sun | diagnostic"
sun_mat = bpy.data.materials.new("Sun emission material")
sun_mat.use_nodes = True
nodes = sun_mat.node_tree.nodes
links = sun_mat.node_tree.links
bsdf = nodes.get("Principled BSDF")
if bsdf:
    bsdf.inputs["Base Color"].default_value = (1.0, 0.45, 0.03, 1.0)
    if "Emission Color" in bsdf.inputs:
        bsdf.inputs["Emission Color"].default_value = (1.0, 0.25, 0.01, 1.0)
        bsdf.inputs["Emission Strength"].default_value = 3.0
    elif "Emission" in bsdf.inputs:
        bsdf.inputs["Emission"].default_value = (1.0, 0.25, 0.01, 1.0)
        if "Emission Strength" in bsdf.inputs:
            bsdf.inputs["Emission Strength"].default_value = 3.0
sun.data.materials.append(sun_mat)

# Fixed Polaris directly above the North Pole at 5000 km model height.
POLARIS_HEIGHT = 5000.0
POLARIS_RADIUS = 180.0
bpy.ops.mesh.primitive_uv_sphere_add(
    segments=24, ring_count=12, radius=POLARIS_RADIUS,
    location=(0.0, 0.0, POLARIS_HEIGHT)
)
polaris = bpy.context.object
polaris.name = "Polaris | fixed at 5000 km above North Pole"
polaris_mat = bpy.data.materials.new("Polaris emission material")
polaris_mat.use_nodes = True
polaris_bsdf = polaris_mat.node_tree.nodes.get("Principled BSDF")
if polaris_bsdf:
    polaris_bsdf.inputs["Base Color"].default_value = (0.85, 0.9, 1.0, 1.0)
    if "Emission Color" in polaris_bsdf.inputs:
        polaris_bsdf.inputs["Emission Color"].default_value = (0.75, 0.85, 1.0, 1.0)
        polaris_bsdf.inputs["Emission Strength"].default_value = 2.0
polaris.data.materials.append(polaris_mat)

# Test addition: three static star spheres, placed at different map latitudes.
STAR_RADIUS = 180.0
STAR_HEIGHT = 5000.0
STAR_SPECS = [
    ("Star | 70N", 70.0, 180.0),
    ("Star | 50N", 50.0, 90.0),
    ("Star | 30N", 30.0, 270.0),
    ("Star | 10N", 10.0, 45.0),
]
star_mat = bpy.data.materials.new("Star emission material")
star_mat.use_nodes = True
star_nodes = star_mat.node_tree.nodes
star_links = star_mat.node_tree.links
star_bsdf = star_nodes.get("Principled BSDF")
if star_bsdf:
    star_bsdf.inputs["Base Color"].default_value = (0.8, 0.9, 1.0, 1.0)
    if "Emission Color" in star_bsdf.inputs:
        star_bsdf.inputs["Emission Color"].default_value = (0.65, 0.8, 1.0, 1.0)
        star_bsdf.inputs["Emission Strength"].default_value = 2.0
    elif "Emission" in star_bsdf.inputs:
        star_bsdf.inputs["Emission"].default_value = (0.65, 0.8, 1.0, 1.0)
        if "Emission Strength" in star_bsdf.inputs:
            star_bsdf.inputs["Emission Strength"].default_value = 2.0
for star_name, latitude, longitude in STAR_SPECS:
    radial_distance = DISK_RADIUS * (90.0 - latitude) / 180.0
    angle = math.radians(longitude)
    star_location = (radial_distance * math.cos(angle), radial_distance * math.sin(angle), STAR_HEIGHT)
    bpy.ops.mesh.primitive_uv_sphere_add(segments=24, ring_count=12, radius=STAR_RADIUS, location=star_location)
    star = bpy.context.object
    star.name = star_name
    star.data.materials.append(star_mat)


# Tokyo observation camera, aimed at the fixed Polaris above the North Pole.
TOKYO_LATITUDE = 35.7
TOKYO_LONGITUDE = 139.77
TOKYO_RADIUS = DISK_RADIUS * (90.0 - TOKYO_LATITUDE) / 180.0
TOKYO_CAMERA_HEIGHT = 100.0

cam_data = bpy.data.cameras.new("Tokyo Observation Camera")
cam = bpy.data.objects.new("Tokyo Observation Camera", cam_data)
bpy.context.collection.objects.link(cam)
tokyo_angle = math.radians(TOKYO_LONGITUDE)
cam.location = (
    TOKYO_RADIUS * math.cos(tokyo_angle),
    TOKYO_RADIUS * math.sin(tokyo_angle),
    TOKYO_CAMERA_HEIGHT
)
cam_data.type = 'PERSP'
cam_data.lens = 50.0
cam_data.clip_start = 0.1
cam_data.clip_end = 100000.0
bpy.context.scene.camera = cam

# Point the Tokyo camera at the fixed Polaris, now that both objects exist.
direction_to_polaris = Vector(polaris.location) - cam.location
cam.rotation_euler = direction_to_polaris.to_track_quat('-Z', 'Y').to_euler()


# Diagnostic addition: animate the Sun and three stars on circular XY orbits.
# One revolution takes 5 seconds at 24 fps (frames 1 through 121).
scene = bpy.context.scene
scene.render.fps = 24
scene.frame_start = 1
scene.frame_end = 121
scene.render.resolution_x = 1280
scene.render.resolution_y = 720
scene.render.resolution_percentage = 100

def animate_circular_orbit(obj, radius, height, start_angle_degrees=0.0):
    # Insert a keyframe at EVERY frame. The earlier version used only four
    # quarter-turn points, so linear interpolation connected them as a square.
    # Dense points calculated from sine/cosine make the trajectory circular.
    for frame in range(scene.frame_start, scene.frame_end + 1):
        fraction = (frame - scene.frame_start) / (scene.frame_end - scene.frame_start)
        angle = math.radians(start_angle_degrees) - 2.0 * math.pi * fraction
        obj.location = (radius * math.cos(angle), radius * math.sin(angle), height)
        obj.keyframe_insert(data_path="location", frame=frame)
    # Linear interpolation between adjacent frame samples avoids Bézier overshoot.
    if obj.animation_data and obj.animation_data.action:
        for fcurve in obj.animation_data.action.fcurves:
            for key in fcurve.keyframe_points:
                key.interpolation = 'LINEAR'

# Sun follows a circle centered exactly on the North Pole (map origin).
# Its radius corresponds to the Tropic of Cancer in this azimuthal-equidistant map.
animate_circular_orbit(sun, radius=SUN_ORBIT_RADIUS, height=5000.0, start_angle_degrees=0.0)


# Diagnostic reference ring: exact North-Pole-centered Tropic of Cancer orbit.
# This ring is independent of the Sun animation and makes visual alignment testable.
ORBIT_RING_Z = 80.0
ORBIT_RING_RADIUS = SUN_ORBIT_RADIUS
ring_mat = bpy.data.materials.new("Tropic of Cancer orbit reference | red")
ring_mat.diffuse_color = (1.0, 0.0, 0.0, 1.0)
ring_mat.use_nodes = True
ring_bsdf = ring_mat.node_tree.nodes.get("Principled BSDF")
if ring_bsdf:
    ring_bsdf.inputs["Base Color"].default_value = (1.0, 0.0, 0.0, 1.0)
    ring_bsdf.inputs["Roughness"].default_value = 1.0
ring_curve = bpy.data.curves.new("Tropic of Cancer reference orbit curve", type='CURVE')
ring_curve.dimensions = '3D'
ring_curve.resolution_u = 1
ring_curve.bevel_depth = 35.0
ring_curve.bevel_resolution = 1
ring_spline = ring_curve.splines.new('POLY')
ring_spline.points.add(360)
for i, pnt in enumerate(ring_spline.points):
    a = 2.0 * math.pi * i / 360.0
    pnt.co = (ORBIT_RING_RADIUS * math.cos(a), ORBIT_RING_RADIUS * math.sin(a), ORBIT_RING_Z, 1.0)
ring_spline.use_cyclic_u = True
ring_obj = bpy.data.objects.new("REFERENCE ORBIT | Tropic of Cancer | centered on North Pole", ring_curve)
bpy.context.collection.objects.link(ring_obj)
ring_obj.data.materials.append(ring_mat)

# Stars orbit at their existing radial distances, each remaining at height 5000.
for star_name, latitude, longitude in STAR_SPECS:
    star_obj = bpy.data.objects.get(star_name)
    if star_obj:
        star_radius = DISK_RADIUS * (90.0 - latitude) / 180.0
        animate_circular_orbit(star_obj, radius=star_radius, height=STAR_HEIGHT, start_angle_degrees=longitude)

scene.frame_set(1)

# Numeric validation at the current frame: the Sun's XY distance must remain fixed.
scene.frame_set(1)
print("\n--- ORBIT VALIDATION ---")
print("North Pole / map origin: (0, 0)")
print("Tropic of Cancer latitude:", TROPIC_OF_CANCER_LATITUDE)
print("Reference orbit radius:", ORBIT_RING_RADIUS)
for frame in (1, 31, 61, 91, 121):
    scene.frame_set(frame)
    loc = sun.matrix_world.translation
    radial_distance = math.hypot(loc.x, loc.y)
    print(f"Frame {frame}: Sun XY=({loc.x:.3f}, {loc.y:.3f}), radius={radial_distance:.3f}, difference={radial_distance-ORBIT_RING_RADIUS:.6f}")
scene.frame_set(1)

# Configure viewport clipping and a predictable top view.
for screen in bpy.data.screens:
    for area in screen.areas:
        if area.type == 'VIEW_3D':
            space = area.spaces.active
            space.clip_start = 0.1
            space.clip_end = 100000.0
            space.region_3d.view_distance = 50000.0
            space.region_3d.view_location = (0.0, 0.0, 0.0)
            space.region_3d.view_perspective = 'ORTHO'
            space.shading.type = 'MATERIAL'

bpy.context.view_layer.objects.active = disk
bpy.ops.object.select_all(action='DESELECT')
disk.select_set(True)
print("Created disk + map + grid; Sun and stars orbit clockwise; fixed Polaris; Tokyo observation camera aimed at Polaris.")
print("Image path:", IMAGE_PATH)
print("Disk vertices:", len(mesh.vertices), "polygons:", len(mesh.polygons))
print("Grid objects:", 9 + 24)
print("Sun location:", tuple(sun.location), "sphere radius:", SUN_RADIUS)
print("Sun orbit radius from North Pole:", SUN_ORBIT_RADIUS)
print("Tokyo camera location:", tuple(cam.location))
print("Polaris location:", tuple(polaris.location))
print("Star objects:", len(STAR_SPECS))


# --- ON-SCREEN ORBIT DIAGNOSTIC (no console required) ---
# Displays numeric checks as a text object in the scene.
# The Sun and reference ring should have the same XY radius at every sampled frame.
diag_lines = [
    "ORBIT DIAGNOSTIC",
    "Expected Sun orbit radius: %.3f" % ORBIT_RING_RADIUS,
]
for check_frame in (1, 31, 61, 91, 121):
    scene.frame_set(check_frame)
    world_loc = sun.matrix_world.translation
    actual_r = math.hypot(world_loc.x, world_loc.y)
    diag_lines.append(
        "Frame %d: radius %.3f | delta %.6f" %
        (check_frame, actual_r, actual_r - ORBIT_RING_RADIUS)
    )
scene.frame_set(1)

diag_curve = bpy.data.curves.new("Orbit diagnostic text data", type='FONT')
diag_curve.body = "\n".join(diag_lines)
diag_curve.size = 500.0
diag_curve.extrude = 0.0
diag_obj = bpy.data.objects.new("ORBIT DIAGNOSTIC TEXT", diag_curve)
bpy.context.collection.objects.link(diag_obj)
diag_obj.location = (-19000.0, 17000.0, 6000.0)
diag_obj.rotation_euler = (0.0, 0.0, 0.0)
diag_mat = bpy.data.materials.new("Orbit diagnostic text material")
diag_mat.diffuse_color = (1.0, 1.0, 1.0, 1.0)
diag_mat.use_nodes = True
diag_bsdf = diag_mat.node_tree.nodes.get("Principled BSDF")
if diag_bsdf:
    diag_bsdf.inputs["Base Color"].default_value = (1.0, 1.0, 1.0, 1.0)
    if "Emission Color" in diag_bsdf.inputs:
        diag_bsdf.inputs["Emission Color"].default_value = (1.0, 1.0, 1.0, 1.0)
        diag_bsdf.inputs["Emission Strength"].default_value = 2.0
diag_obj.data.materials.append(diag_mat)

print("On-screen diagnostic text created.")