星の位置を四角い紙に投影するためのコード

star.py

#%%
#データ読み込み
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
from scipy.spatial.transform import Rotation as R

#データのもってきかた
#https://heasarc.gsfc.nasa.gov/cgi-bin/W3Browse/w3query.pl?&tablehead=name%3Dheasarc_hipparcos%26description%3DHipparcos+Main+Catalog%26url%3Dhttp%3A%2F%2Fheasarc.gsfc.nasa.gov%2FW3Browse%2Fstar-catalog%2Fhipparcos.html%26archive%3D%26radius%3D1%26mission%3DSTAR%2BCATALOG%26priority%3D3&mission=STAR+CATALOG&Action=More+Options&Action=Parameter+Search&ConeAdd=1
#vmag,ra_deg,dec_regを出力項目として選ぶ
#下の方で
#Limit Result to:をNO Limitに設定
#Output formatをexcelに設定して「Search」を実行

#参考：
# https://coelostat.hatenablog.com/entry/2016/02/16/223337
# http://www.kusastro.kyoto-u.ac.jp/~iwamuro/LECTURE/OBS/coord.html

df=pd.read_excel("browse_results.xls",sheet_name=1)
pos=df[["ra_deg","dec_deg"]].to_numpy()
vmag=df["vmag"].to_numpy()
# %%
#投影条件の設定

#光源からLmm離れた所に幅W,高さH (mm)の紙をおいてそこに穴を空ける事を考える
#光源からの距離
L=75
#遮光版の大きさ
W=150
H=150
#紙のセンター位置
Cwh=np.array((W/2,H/2))

# %%
#デバグ用に3D天球儀表示
if False:
    %matplotlib qt
    Crd=np.array((0,60))
    fig = plt.figure()
    ax = fig.add_subplot(projection='3d')

    for vv,pp in zip(vmag,pos):
        #暗いのは切る
        if 3<vv:continue
        p=np.copy(pp)
        p[0]-=Crd[0]
        #x:右 y:奥行き z:上　
        r=(0,1,0)
        r=R.from_euler("ZX",p,degrees=True).apply(r)
        r=R.from_euler("x",-Crd[1],degrees=True).apply(r)

        dc=cm.hsv( (pp[1]+90)/180)
        rc=cm.hsv( pp[0]/360)
        ax.scatter(r[0],r[1],r[2],s=100*2.5**(-vv),color=dc )
    plt.show()

#%%

#見ている角度
Crd=np.array((80,0)) #オリオン座    
Crd=np.array((100,0)) #冬の大三角    
Crd=np.array((-70,20)) # 夏の大三角   
Crd=np.array((180,60)) #北斗七星

fp=open("star.nc","w")
fp.write(
"""
G90
M3 S0
F 1000
""")
def vmag_to_size(vmag,rate=20):
    #h_max等星⇒0
    #h_min等星⇒1
    h_max=4
    h_min=1.5

    v=(h_max-vmag)/(h_max-h_min)
    v=np.max((0,v))
    v=np.min((1,v))
    return np.sqrt((rate*rate)**v)


def plot(fp,_x,_y,vmag):
    _size=vmag_to_size(vmag)*0.05
    size='{:0.2f}'.format(_size)
    x='{:0.3f}'.format(_x)
    y='{:0.3f}'.format(_y)
    xp='{:0.3f}'.format(_x+_size)
    yp='{:0.3f}'.format(_y+_size)
    fp.write(
f"""
G0 X{x} Y{y}
G4 P0.2
S500
G1 X{xp} Y{y}
G1 X{xp} Y{yp} 
G1 X{x} Y{yp}
G1 X{x} Y{y}
S0
""")


plt.figure(figsize=(5,5))
plt.ylim([0,H])
plt.xlim([0,W])
vpos=[]
count=0
for vv,pp in zip(vmag,pos):
    #暗いのは切る
    if 6<vv:continue
    p=np.copy(pp)
    p[0]-=Crd[0]
    #x:右 y:奥行き z:上　
    r=(0,1,0)   
    r=R.from_euler("ZX",p,degrees=True).apply(r)
    r=R.from_euler("x",-Crd[1],degrees=True).apply(r)
    #裏側なので描画しない
    if r[1]<=0:continue
    
    wh=r[[0,2]]*(L/r[1]) + Cwh
    if np.isnan(wh).any() or np.min(wh)<0 or W<wh[0] or H<wh[1]:continue

    dc=cm.hsv( (pp[1]+90)/180)
    rc=cm.hsv( pp[0]/360)
    vpos.append((wh[0],wh[1],vv))
    plt.scatter(wh[0],wh[1],s=vmag_to_size(vv)*2,color=dc )
    count+=1
print(f"count:{count}")


vpos=np.array(vpos)
vidx=[0]
ridx=list(range(1,vpos.shape[0]))
while(len(ridx)):
    d=vpos[ridx,:2]-vpos[vidx[-1],:2]
    i=np.argmin(np.linalg.norm(d,axis=1))
    vidx.append(ridx.pop(i))
plt.plot(vpos[vidx,0],vpos[vidx,1],alpha=0.5)    
plt.show()

for x,y,v in vpos[vidx,:]:
    plot(fp,x,y,v)




fp.write(
"""
M5
""")
fp.close()

# %%