牧妖的牧场

import numpy as np
from astropy.time import Time
from astropy.coordinates import EarthLocation, AltAz, ICRS, SkyCoord, CartesianRepresentation
import astropy.units as u

# ====== 1. 基本设置 ======
obs_lat = 39.92    # 纬度
obs_lon = 116.38   # 经度
obs_time = Time('2025-06-15T17:00:00')  # UTC

location = EarthLocation(lat=obs_lat*u.deg, lon=obs_lon*u.deg)
altaz_frame = AltAz(obstime=obs_time, location=location)

# ====== 2. 求银道面法线与地平面法线 ======
# 所有运算在赤道坐标系（ICRS）下进行
north_galactic_pole = SkyCoord(l=0*u.deg, b=90*u.deg, frame='galactic')
north_gal_icrs = north_galactic_pole.transform_to(ICRS())
n_gal = np.array(north_gal_icrs.cartesian.xyz.value)

print(f"银道面法线（赤道坐标系）:")
print(f"  赤经 RA  = {north_gal_icrs.ra.deg:.2f}°")
print(f"  赤纬 Dec = {north_gal_icrs.dec.deg:.2f}°")

zenith = SkyCoord(alt=90*u.deg, az=0*u.deg, frame=altaz_frame)
zenith_icrs = zenith.transform_to(ICRS())
n_h = np.array(zenith_icrs.cartesian.xyz.value)

print(f"地平面法线（赤道坐标系）:")
print(f"  赤经 RA  = {zenith_icrs.ra.deg:.2f}°")
print(f"  赤纬 Dec = {zenith_icrs.dec.deg:.2f}°")

# ====== 3. 求银道面与地平面的交线方向 ======
# 这是银道面和地平面的交线方向（沿银道弧在地平线上升下落方向）
# 两平面法向量的叉积就是交线方向
vec_galactic_horizon_cross = np.cross(n_gal, n_h)
vec_galactic_horizon_cross /= np.linalg.norm(vec_galactic_horizon_cross)

# ====== 4. 求银拱在地平坐标中最高点方向 ======
# 这是银拱在地平坐标中最高点所指的方向
# 这个方向位于银道面内，并垂直于上面那条交线
vec_galactic_arch_top = np.cross(vec_galactic_horizon_cross, n_gal)
vec_galactic_arch_top /= np.linalg.norm(vec_galactic_arch_top)

# ====== 5.转换回地平坐标系并输出 ======
sc_cross = SkyCoord(CartesianRepresentation(*vec_galactic_horizon_cross),
                    frame=ICRS())
h_cross = sc_cross.transform_to(altaz_frame)
az_cross = h_cross.az.deg % 360
alt_cross = h_cross.alt.deg
print(f"银道与地平交线方向（升落点）:")
print(f"  方位角 Az  = {az_cross:.2f}°")
print(f"  方位角 Az  = {(az_cross+180)%360:.2f}°")

sc_top = SkyCoord(CartesianRepresentation(*vec_galactic_arch_top),
                    frame=ICRS())
h_top = sc_top.transform_to(altaz_frame)
az_top = h_top.az.deg % 360
alt_top = h_top.alt.deg
print(f"银道弧最高点方向:")
print(f"  方位角 Az  = {az_top:.2f}°")
print(f"  仰角   Alt = {alt_top:.2f}°")

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
from astropy.time import Time
from astropy.coordinates import EarthLocation, AltAz, ICRS, SkyCoord, CartesianRepresentation
import astropy.units as u
import pandas as pd

# ====== 1. 观测点与时间设置 ======
obs_lat = 40    # 纬度
obs_lon = 116.38   # 经度
location = EarthLocation(lat=obs_lat*u.deg, lon=obs_lon*u.deg)

start_time = Time('2025-06-15T00:00:00')  # UTC 时间起点
delta_minutes = 10
num_points = int(24*60 / delta_minutes)
times = start_time + np.linspace(0, 24, num_points) * u.hour

# ====== 2. 银道面法线（固定） ======
north_gal_icrs = SkyCoord(l=0*u.deg, b=90*u.deg, frame='galactic').transform_to(ICRS())
n_gal = np.array(north_gal_icrs.cartesian.xyz.value)

# ====== 3. 逐时间计算最高点方向在地平坐标系下的位置 ======
az_top_list = []
az_bottom_list = []
alt_list = []

for t in times:
    altaz_frame = AltAz(obstime=t, location=location)

    # 地平法线（天顶方向）
    zenith = SkyCoord(alt=90*u.deg, az=0*u.deg, frame=altaz_frame).transform_to(ICRS())
    n_h = np.array(zenith.cartesian.xyz.value)

    # 银道与地平交线方向
    vec_cross = np.cross(n_gal, n_h)
    vec_cross /= np.linalg.norm(vec_cross)

    # 银道弧最高点方向
    vec_top = np.cross(vec_cross, n_gal)
    vec_top /= np.linalg.norm(vec_top)

    # 转换到地平坐标系
    sc_top = SkyCoord(CartesianRepresentation(*vec_top), frame=ICRS()).transform_to(altaz_frame)

    az_top_list.append(sc_top.az.deg % 360)
    az_bottom_list.append((sc_top.az.deg+180) % 360)
    alt_list.append(sc_top.alt.deg)

# 去除az_list跳变(优化折线图显示)
def remove_jumps(az_list):
    result = [az_list[0]]
    for i in range(1, len(az_list)):
        if abs(az_list[i] - az_list[i-1]) > 100: # 如果相邻az差值超过100
            result[-1] = None  # 将前一个值置为 None
        result.append(az_list[i])
    return result
az_top_list=remove_jumps(az_top_list)
az_bottom_list=remove_jumps(az_bottom_list)

# ====== 4. 可视化 ======

font_path = "/usr/share/fonts/wenquanyi/wqy-zenhei/wqy-zenhei.ttc"
font_prop = FontProperties(fname=font_path)

time_list = [t.datetime for t in times]
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), sharey=False)

ax1.plot(time_list, alt_list)
ax1.set_title("银道弧最高点仰角随时间变化", fontproperties=font_prop)
ax1.set_xlabel("时间 (UTC)", fontproperties=font_prop)
ax1.set_ylabel("仰角 Alt (°)", fontproperties=font_prop)
ax1.grid(True)

ax2.plot(time_list, az_top_list, label="最高点 Az")
ax2.plot(time_list, az_bottom_list, label="最低点 Az")
ax2.set_title("银道弧最高点方位角随时间变化", fontproperties=font_prop)
ax2.set_xlabel("时间 (UTC)", fontproperties=font_prop)
ax2.set_ylabel("方位角 Az (°)", fontproperties=font_prop)
ax2.grid(True)
ax2.legend(prop=font_prop)

fig.autofmt_xdate()
fig.tight_layout()
plt.show()
plt.savefig("movement.png")

// 可以通过https://editor.p5js.org/直接生成

let angle = 0;
let observerLat = 40;
function setup() {
  createCanvas(600, 600, WEBGL);
  angleMode(RADIANS);
  frameRate(60);
}

function draw() {
  background(255);
  lights();

  // 摄像机设置
  let camRadius = 1000;
  let theta = radians(200);
  let camX = camRadius * sin(theta);
  let camY = camRadius * -cos(theta);
  let camZ = 750;

  camera(camX, camY, camZ, 0, 0, 0, 0, 0, -1);

  // 地平面
  noStroke();
  fill(180);
  box(450, 450, 5);
  rotateY(radians(90-observerLat));
  rotateZ(radians(angle));

  // 球体
  noStroke();
  fill(50, 50, 100, 255);
  sphere(200);

  // 赤道坐标极轴（Z轴）
  stroke(255, 50, 255);  // 极轴
  strokeWeight(2);
  line(0, 0, -300, 0, 0, 300);  // Z轴线

  // 渐变银拱
  noFill();
  strokeWeight(3);
  drawGalacticArc(220);
  // drawXYZAxes(250);
  angle += 2;
  // if (frameCount <= 180) {
  // saveCanvas('frame_' + nf(frameCount, 4), 'png');
  // }
}

function drawXYZAxes(length = 300) {
  strokeWeight(2);

  // X 轴（红色）
  stroke(255, 0, 0);
  line(0, 0, 0, length, 0, 0);

  // Y 轴（绿色）
  stroke(0, 255, 0);
  line(0, 0, 0, 0, length, 0);

  // Z 轴（蓝色）
  stroke(0, 0, 255);
  line(0, 0, 0, 0, 0, length);
}

function raDecToVec(ra_h, dec_d) {
  let ra = radians(ra_h * 15);
  let dec = radians(dec_d);
  let x = cos(dec) * cos(ra);
  let y = cos(dec) * sin(ra);
  let z = sin(dec);
  return createVector(x, y, z);
}

function rotateAroundAxis(vec, axis, angle) {
  let cosA = cos(angle);
  let sinA = sin(angle);
  let my_dot = vec.dot(axis);
  let my_cross = p5.Vector.cross(axis, vec);

  return p5.Vector.add(
    p5.Vector.add(
      vec.copy().mult(cosA),
      my_cross.mult(sinA)
    ),
    axis.copy().mult(my_dot * (1 - cosA))
  ).normalize();
}

function drawGalacticArc(radius) {
  let steps = 200;

  // 北银极方向（圆心）
  let pole = raDecToVec(192.8595 / 15, 27.1283);  // RA:小时, Dec:度

  // 银心方向（起点）
  let start = raDecToVec(266.4168 / 15, -29.0078);

  // 用单位向量绘制
  beginShape();
  for (let i = 0; i <= steps; i++) {
    let t = i / steps;
    let angle = t * TWO_PI;

    // Rodrigues’ rotation: rotate `start` around `pole`
    let v = rotateAroundAxis(start, pole, angle).mult(radius);

    // 渐变色：红到蓝
    let r = lerp(255, 0, t);
    let g = 0;
    let b = lerp(0, 255, t);
    stroke(r, g, b);

    vertex(v.x, v.y, v.z);
  }
  endShape();
}