由基本天体物理量测量引发的惨案

一个天体物理量测量有关的实验，强行提高了我对Linux基本操作的熟悉程度以及python使用的能力。所以说，学东西还是要在使用中学。

Linux有关操作

安装和删除各种软件使用

 cd ..  %返回上层目录

 ls    %查看本层有什么文件

 cd 文件名 %进入该文件夹

 ssh 用户名@地址 %用于连接服务器

logout    %断开连接

vim a.py %打开.py文件（此后输入'i'可进入编辑模式，按Esc退出编辑模式，输入'：wq' 保存并退出

%也可安装文本编辑器vscode，管理代码更加方便。

python a.py %编译运行a.py

gunzip *.gz   %星号代表任意，此句解压所有文件

pip install 第三方库名字 %用于为python安装库

sudo apt-get install 软件名 %安装软件（删除用remove）

%可以下载gftp连接服务器，可视化传文件很爽

%有一个可以在服务器安装vscode在本地打开的，好像是ssh -X 

将虚拟机的文件移动到主机

在VMworkstation 17 player虚拟机中安装linux后可以设置一些与windows公用的文件夹。
打开虚拟机后点击左上角的player–管理–虚拟机设置，新窗口中点击上方 选项，点击左边共享文件夹，启用后选择主机中共享文件夹中位置。

在虚拟机中，把文件夹移动到路径 /mnt/hgfs/共享文件夹 中，可以共享。

如何到这个路径？linux点击 files，一路backspace，就能找mnt。

python有关内容

conda

conda可以创建许多环境，每个环境都可以有自己的python版本，自己的库，不互相干扰。

当进行不同目的的计算时，用到的库版本可能不一样，还有可能相互影响。库所要求的python版本也不同，这时它的便利性就体现了出来。

我安装了一个miniconda就能用了。下面是网址，安装方法里面有。
https://docs.conda.io/projects/miniconda/en/latest/

一些基本的操作：

conda create -n 环境名字 python=3.9   %也可在后面跟一些第三方库，具体见官方文档或上网查

conda info -e    %电脑上有哪些环境

conda activate 环境名字    %激活某环境

conda deactivate     %退出环境

没了，就用过这几个

如果下载了vscode，打开.py文件后点击右下角的版本，可以选择相应的环境运行代码。ipynb 是交互式的，好用，在vscode里面需要下载一些拓展才能使用。运行时会提示选择环境。

文件读取

当要处理很多数据在不同文件中时，可以将文件名进行一个排序，写入到某个文件中，保证后续文件处理顺序不会乱。使用下面代码可以获得按字典序排的文件名。即1.gz < 12.gz < 2.gz

import os
with open('SpecFiles.txt', 'w') as file:  #此句不用close，运行完自动关。
    f_path  = '/home/simonz/J034813/RV/mid'
    
    f = sorted(os.listdir(f_path))
    for filename in f:
        if filename.endswith('.gz'):
            file.write(filename+'\n')

读取文件中的东西用下面的句子：

files = open('SpecFiles.txt','r')
RVfile = open('RV.txt','w')
Tfile = open('SpecTime.txt','w')
Lines = files.readlines()
count = 0
t0 = 0
for line in Lines:
    specfile=line.strip()  #括号里是行内分隔符，一行一个东西时不用填。
files.close()
RVfiles.close()
Tfile.close()  #用完记得关闭文件释放内存

使用pandas库进行csv读取也很不错

import pandas as pd
data = pd.read_csv('a.csv')
time = data['time']
flux = data['flux']
flux_err = data['flux_err']

剩下的都是一些奇奇怪怪的用法：

曲线拟合

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize as op

def rv_func(x2,a,b,c):
#第一个是变量，后面的是参数，设好函数形式
    return a*np.sin(4*np.pi*(x2+b/(2*np.pi)))+c 
popt, pcov = op.curve_fit(rv_func,phase,rv_data)  #函数，自变量，因变量；返回值为参数和协方差（也可能不是协方差，反正是用来看拟合好坏的）

画图

三维的：

import matplotlib.pyplot as plt
ax = plt.subplot(projection = '3d')
ax.set_title('luoxiban')
ax.scatter(xx,yy,zz,c='b')
plt.show()

二维的：

plt.plot(x,y,'*')

解方程

import sympy as sp
y0 = sp.Symbol('y0')
m1 = (0.98/(0.98+1.35))#不可见
m2 = (1.35/(0.98+1.35))#可见
q = m2/m1
a = sp.solve(-1/((y0+m2)**2)+q/((y0-m1)**2)+(1+q)*y0, y0)
print(a)

有的解会出现一个极小的虚部，而这个虚数与python自带的虚数是不兼容的。这时可以用sp.re()
和sp.im()进行虚部大小的判断。以及可以用这个转化为python自带的复数。

洛希瓣

最后贴一个画洛希瓣的代码，m2为亮的一半，填充率为f，质量比m2/m1为q，倾角thet：

import sympy as sp
from decimal import Decimal
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import math as mt
y0 = sp.Symbol('y0')
m1 = (0.98/(0.98+1.35))#不可见
m2 = (1.35/(0.98+1.35))#可见
q = m2/m1
a = sp.solve(-1/((y0+m2)**2)+q/((y0-m1)**2)+(1+q)*y0, y0)
print(a)
for i in range(0,len(a)):
    if sp.re(a[i])<m1 and sp.re(a[i])>-m2:
        L1 = float(sp.re(a[i]))
print(L1)
#fig = plt.figure()
#ax = fig.add_subplot(111, projection='3d')


#phi = np.linspace(0, np.pi, 100)
#theta, phi = np.meshgrid(theta, phi)
f = 0.6
L1 = m1-f*(m1-L1)
ans = -1/mt.sqrt((L1+m2)**2)-q/mt.sqrt((L1-m1)**2)-(1+q)*(L1**2)/2
ans = round(ans, 4)
print(L1)

y1 = sp.Symbol('y1')
b = sp.solve(-1/(y1+m2)-q/(y1-m1)-(1+q)*(y1**2)/2-ans, y1)
for i in range(0,len(b)):
    if sp.re(b[i])>m1:
        L2 = float(sp.re(b[i]))
        break
print(type(L1))
y = np.linspace(L1, L2, 100)
phi = np.linspace(0, 2*np.pi,72)#共有72个，0到71，最后一个比2pi大一点
dphi=2*np.pi/72
z = [0]
#for i in range(0,len(y)):
for i in range(1, len(y)):
    temp = 0
    former = 1000
    nowans = 0
    while 1:
        nowans = -1/mt.sqrt((y[i]+m2)**2 + temp**2)-q/mt.sqrt((y[i]-m1)**2 + temp**2)-(1+q)*(temp**2 + y[i]**2)/2-ans
        if nowans > 0:
            z.append(temp)
            break
        former = nowans
        temp =  temp+0.001
s = [[0 for _ in range(len(phi))] for _ in range(len(z))]#本来想直接用一个dphi*z，但是好像并不是平均的。如果能严格要求保留小数位数或许可以 
dirx = [[0 for _ in range(len(phi))] for _ in range(len(z))]
diry = [[0 for _ in range(len(phi))] for _ in range(len(z))]
dirz = [[0 for _ in range(len(phi))] for _ in range(len(z))]
for i in range(1,len(y)): #面积i为i到i-1段的部分，dirx[i][k]为si对应phik到phik+1部分,phi71=phi0因而dir的k有0到70
    dy = y[i]-y[i-1] 
    dz = z[i]-z[i-1]
    guiy = dy/mt.sqrt(dy**2+dz**2)
    guiz = dz/mt.sqrt(dy**2+dz**2)
    for k in range(0,len(phi)-1):
        dirx[i][k] = guiy*mt.sin((phi[k]+phi[k+1])/2)
        diry[i][k] = -guiz
        dirz[i][k] = guiy*mt.cos((phi[k]+phi[k+1])/2)
        s[i][k] = mt.sqrt(dz**2+dy**2)*(z[i]+z[i-1])*(phi[k+1]-phi[k])/2 

#计算面积
S=[]
for i in range(len(phi)):
    tp = 0
    # 双星旋转，等效于视线绕z轴旋转，倾角为\theta，假设phi在0时最大，phi在90时最小（不一定？）
    #先假设一个thet
    thet = 55/180*np.pi
    nx = mt.sin(thet)*mt.cos(phi[i]) # phii=0时用这个
    ny = mt.sin(thet)*mt.sin(phi[i]) # phii=pi/2时用这个
    nz = mt.cos(thet)
    #对于
    for i in range(1,len(y)):
        for k in range(0,len(phi)-1):
            fangxda = nx*dirx[i][k]+nz*dirz[i][k]+ny*diry[i][k]
            if fangxda>0:
                tp = tp+fangxda*s[i][k]
    S.append(tp)
print(max(S)/min(S))
plt.plot(phi, S)
plt.show()

xx=[]
yy=[]
zz=[]


for i in range(0,len(y)):
    for j in range(0,len(phi)):
        yy.append(y[i])
        zz.append(z[i]*mt.cos(phi[j]))
        xx.append(z[i]*mt.sin(phi[j]))
ax = plt.subplot(projection = '3d')
ax.set_title('luoxiban')
ax.scatter(xx,yy,zz,c='b')
plt.show()
## ans = -1/mt.sqrt((y+m2)**2 + z**2 + x**2)-q/mt.sqrt((y-m1)**2 + z**2 + x**2)-(1+q)*(z**2 + x**2 + y**2)/2
#最y最大在1.1以内，左边用a【0】