import bpy
import math
from mathutils import Vector

# Blender 4.0.2: Flat-Earth / azimuthal-equidistant disk with Orion stars.
# Distances are treated as kilometres (1 Blender unit = 1 km).
IMAGE_PATH = ""
DISK_RADIUS = 20000.0
GRID_Z = 60.0
LINE_RADIUS = 18.0
STAR_HEIGHT = 5000.0
STAR_RADIUS = 325.0

# Clear the 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 km", mesh)
bpy.context.collection.objects.link(disk)

# UV coordinates for AE.jpg; rotate the map 90 degrees counter-clockwise.
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
        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)

# Latitude/longitude grid.
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
grid_bsdf = grid_mat.node_tree.nodes.get("Principled BSDF")
if grid_bsdf:
    grid_bsdf.inputs["Base Color"].default_value = (0.015, 0.015, 0.015, 1.0)
    grid_bsdf.inputs["Roughness"].default_value = 1.0

def add_polyline(name, points, cyclic=False, material=grid_mat):
    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(material)
    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=True)

# Longitude rays every 15 degrees.
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)

# Two equal-sized stars on the equator, both at altitude 5000 km.
# Star 1: same longitude as Tokyo, representing the due-south meridian position.
# Star 2: 90 degrees west of Tokyo, the equatorial subsolar longitude when
# the Sun is setting due west in Tokyo at an equinox.
TOKYO_LAT = 35.6762
TOKYO_LON = 139.6503
STAR_HEIGHT = 5000.0
STAR_RADIUS = 325.0


def ae_xy(latitude, longitude):
    # Match the map/grid convention in the supplied script.
    r = DISK_RADIUS * (90.0 - latitude) / 180.0
    a = math.radians(longitude)
    return (r * math.cos(a), r * math.sin(a))


def make_emission_material(name, color, strength=0.7):
    material = bpy.data.materials.new(name)
    material.diffuse_color = (*color, 1.0)
    material.use_nodes = True
    bsdf = material.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 = strength
        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 = strength
    return material


def add_split_star(name, latitude, longitude):
    x, y = ae_xy(latitude, longitude)
    bpy.ops.mesh.primitive_uv_sphere_add(
        segments=48, ring_count=24, radius=STAR_RADIUS,
        location=(x, y, STAR_HEIGHT))
    star = bpy.context.object
    star.name = name
    bright_orange = make_emission_material(name + " | lower bright orange", (1.0, 0.32, 0.025))
    dark_orange = make_emission_material(name + " | upper dark orange", (0.42, 0.075, 0.008))
    star.data.materials.append(bright_orange)
    star.data.materials.append(dark_orange)
    # Face assignment is based on each polygon's local Z coordinate:
    # lower half is bright orange; upper half is dark orange.
    for poly in star.data.polygons:
        poly.material_index = 0 if poly.center.z < 0 else 1
        poly.use_smooth = True
    return star


# Remove the former Orion stars by virtue of the scene reset above; create only these two.
south_star = add_split_star(
    "Equator Star 1 | Tokyo south culmination | altitude 5000 km",
    0.0, TOKYO_LON)
sunset_longitude = ((TOKYO_LON - 90.0 + 180.0) % 360.0) - 180.0
sunset_star = add_split_star(
    "Equator Star 2 | Tokyo equinox sunset position | altitude 5000 km",
    0.0, sunset_longitude)


def add_camera(name, location, target, lens=50.0):
    cam_data = bpy.data.cameras.new(name)
    cam = bpy.data.objects.new(name, cam_data)
    bpy.context.collection.objects.link(cam)
    cam.location = location
    cam_data.type = 'PERSP'
    cam_data.lens = lens
    cam_data.clip_start = 0.1
    cam_data.clip_end = 100000.0
    direction = Vector(target) - cam.location
    cam.rotation_euler = direction.to_track_quat('-Z', 'Y').to_euler()
    return cam


tokyo_x, tokyo_y = ae_xy(TOKYO_LAT, TOKYO_LON)
south_camera = add_camera(
    "Camera Tokyo | looking at south-culminating equator star",
    (tokyo_x, tokyo_y, 100.0), south_star.location, lens=50.0)
sunset_camera = add_camera(
    "Camera Tokyo | looking at equinox sunset-position star",
    (tokyo_x, tokyo_y, 100.0), sunset_star.location, lens=50.0)

# Third camera is directly below the sunset-position star's nadir and looks up at it.
sunset_x, sunset_y = ae_xy(0.0, sunset_longitude)
nadir_camera = add_camera(
    "Camera directly below sunset star | looking upward",
    (sunset_x, sunset_y, 100.0), sunset_star.location, lens=50.0)

scene = bpy.context.scene
scene.camera = south_camera
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

# Predictable viewport.
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.ops.object.select_all(action='DESELECT')
disk.select_set(True)
bpy.context.view_layer.objects.active = disk

print("Created AE disk/grid, two equatorial stars at altitude 5000 km, and three cameras.")
print("Active camera: Tokyo south-culmination star view.")
print("Choose the other cameras in the Outliner and assign them to scene.camera for alternate views.")
print("Tokyo longitude:", TOKYO_LON, "equinox sunset-star longitude:", sunset_longitude)
print("Map image path:", IMAGE_PATH)