一个天体物理量测量有关的实验,强行提高了我对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】