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 counter-clockwise on the disk.
        # UV rotation around image center: (u, v) -> (1 - v, 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)


# Static celestial trajectories only: no celestial spheres and no animation.
TROPIC_OF_CANCER_LATITUDE = 23.4366
SUN_ORBIT_RADIUS = DISK_RADIUS * (90.0 - TROPIC_OF_CANCER_LATITUDE) / 180.0
SUN_HEIGHT = 5000.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),
]

# Tokyo observation camera, positioned on the map at Tokyo's latitude/longitude.
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

def add_orbit_curve(name, radius, height, color, bevel=28.0):
    curve = bpy.data.curves.new(name + " curve", type='CURVE')
    curve.dimensions = '3D'
    curve.resolution_u = 1
    curve.bevel_depth = bevel
    curve.bevel_resolution = 2
    spline = curve.splines.new('POLY')
    steps = 720
    spline.points.add(steps - 1)
    for i, pnt in enumerate(spline.points):
        angle = 2.0 * math.pi * i / steps
        pnt.co = (
            radius * math.cos(angle),
            radius * math.sin(angle),
            height,
            1.0
        )
    spline.use_cyclic_u = True
    obj = bpy.data.objects.new(name, curve)
    bpy.context.collection.objects.link(obj)
    mat = bpy.data.materials.new(name + " material")
    mat.diffuse_color = (*color, 1.0)
    mat.use_nodes = True
    bsdf = mat.node_tree.nodes.get("Principled BSDF")
    if bsdf:
        bsdf.inputs["Base Color"].default_value = (*color, 1.0)
        if "Emission Color" in bsdf.inputs:
            bsdf.inputs["Emission Color"].default_value = (*color, 1.0)
            bsdf.inputs["Emission Strength"].default_value = 1.5
        elif "Emission" in bsdf.inputs:
            bsdf.inputs["Emission"].default_value = (*color, 1.0)
            if "Emission Strength" in bsdf.inputs:
                bsdf.inputs["Emission Strength"].default_value = 1.5
    obj.data.materials.append(mat)
    return obj

# Polaris is retained as a visible sphere, fixed above the North Pole.
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 sphere | height 5000"
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.75, 0.85, 1.0, 1.0)
    if "Emission Color" in polaris_bsdf.inputs:
        polaris_bsdf.inputs["Emission Color"].default_value = (0.65, 0.8, 1.0, 1.0)
        polaris_bsdf.inputs["Emission Strength"].default_value = 2.0
polaris.data.materials.append(polaris_mat)

# Draw the Sun's path and the four stars' paths above the map.
# The paths are static circles centered on the North Pole.
add_orbit_curve(
    "Trajectory | Sun | Tropic of Cancer",
    SUN_ORBIT_RADIUS, SUN_HEIGHT, (1.0, 0.35, 0.02), bevel=35.0
)

orbit_colors = [
    (0.9, 0.1, 0.9),   # 70N
    (0.1, 0.35, 1.0),  # 50N
    (0.0, 0.8, 0.9),   # 30N
    (0.2, 0.8, 0.25),  # 10N
]
for (star_name, latitude, longitude), color in zip(STAR_SPECS, orbit_colors):
    orbit_radius = DISK_RADIUS * (90.0 - latitude) / 180.0
    add_orbit_curve(
        "Trajectory | " + star_name,
        orbit_radius, STAR_HEIGHT, color, bevel=28.0
    )

# Aim the Tokyo camera at Polaris' fixed geometric position (not a visible sphere).
polaris_target = Vector((0.0, 0.0, 5000.0))
direction_to_polaris = polaris_target - cam.location
cam.rotation_euler = direction_to_polaris.to_track_quat('-Z', 'Y').to_euler()

scene = bpy.context.scene
scene.frame_start = 1
scene.frame_end = 1
scene.frame_set(1)
scene.render.resolution_x = 1280
scene.render.resolution_y = 720
scene.render.resolution_percentage = 100

# 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 map, grid, Tokyo camera, static trajectory curves, and fixed Polaris sphere (no other celestial spheres or animation).")
print("Image path:", IMAGE_PATH)
print("Disk vertices:", len(mesh.vertices), "polygons:", len(mesh.polygons))
print("Grid objects:", 9 + 24)
print("Sun trajectory radius:", SUN_ORBIT_RADIUS)
print("Tokyo camera location:", tuple(cam.location))
print("Polaris location:", tuple(polaris.location))
print("Trajectory curves:", 1 + 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.")