利用几何关系求解次级谱仪

博士后出站报告附录

# -*- coding: utf-8 -*-
"""
Created on Tue Dec 21 10:50:15 2021

@author: zhouhuibin
"""

import matplotlib.pyplot as plt
import numpy as np
import math
import os
from zhbcode import selffunction as sf
from tkinter import _flatten # 将二维列表转换为一维列表

maindir = r'D:\CSNS\Summary\Data Treatment\python\Norm Data Output\direct'
datapath = sf.mkdir(maindir+r'\data')
figpath = sf.mkdir(maindir+r'\figs')

# Physics constant
m_n = 1.675e-27 # unit kg, neutron mass
h_P = 6.626e-34 # unit J s, planck constant

# Instrument parameters
L_i = L_1 = 90 # unit m
lambda_i = 10e-10 #unit m

# Moderator
t_0 = 46.887e-6 #unit s, time delay for 6.267 AA

# Analyzer parameter
d = 3.13545e-10 # unit m, Spaing of Si-111 plane
lambda_off = 2*d # wavelength when theta is 90 deg;

# Parameters of secondary spectrometer
h_sample = 0.03 # the height of sample is 3cm
h_ana_center = 0.087 # unit m, distance from the center of analyzer sphere to the center of sample(i.e. origin)
h_array = np.linspace(h_ana_center-h_sample/2,h_ana_center+h_sample/2,5) # unit m, the distance from scattering point to spheric center of analyzer
R = 2.5 # unit m, radius of analyzer sphere
r_0 = 0.27 #unit m, radius of detector cylinder
h_d = 0.15 # unit m, the vertical offset of the detector center; elevation angle:15 cm -- 15.5 deg; 17.33 cm -- 20 deg;

# array assignment
gamma_deg = np.linspace(0,20,201)
gamma_rad = np.radians(gamma_deg)

name_theta = []
unit_theta = []
comment_theta = []
data_theta = []

name_L_f = []
unit_L_f = []
comment_L_f = []
data_L_f = []

name_E = []
unit_E = []
comment_E = []
data_E = []

for i in range(len(h_array)):

    h = h_array[i]
    alpha_rad = sf.arcsin(h*np.cos(gamma_rad)/R)
    alpha_deg = np.degrees(alpha_rad)
    theta_B_deg = 90-alpha_deg
    theta_B_rad = np.radians(theta_B_deg)
    phi_deg = gamma_deg - alpha_deg + 90
    phi_rad = np.radians(phi_deg)
    beta_deg = 180 - phi_deg + alpha_deg
    beta_rad = np.radians(beta_deg)

    y_ana = h - R*np.cos(phi_rad)
    z_ana = R*np.sin(phi_rad)
    h_0 = y_ana - (z_ana - r_0)*sf.cot(beta_rad)
    y_off = h_0 - h_d
    rho = np.sqrt(r_0**2+(h_0-h)**2)

    psi_rad = sf.arctan((h_0-h)/r_0) 
    psi_deg = np.degrees(psi_rad)

    L_2 = np.sqrt(R**2+h**2-2*R*h*np.cos(phi_rad))
    L_3 = np.sqrt(R**2+rho**2-2*R*rho*np.sin(phi_rad-psi_rad))
    L_f = L_2+L_3

    # alldata output
    name = _flatten(['gamma','alpha','theta_B','phi','beta','h_0','y_off','y_ana','z_ana'])
    unit = _flatten([('deg',)*5,('m',)*4])
    comment = _flatten([(str('%.2f' % (100*(0.087-h)))+'cm',)*len(name)])
    data = np.c_[gamma_deg,alpha_deg,theta_B_deg,phi_deg,beta_deg,h_0,y_off,y_ana,z_ana]
    filename = str(i+1).rjust(3,'0')+'_alldata_'+str('%.2f' % (100*(0.087-h)))+'cm.dat'
    sf.dataoutput(name,unit,comment,data,datapath,filename)

    # theta_B output
    name_theta = name_theta+['h_0','theta_B',]
    unit_theta = unit_theta+['m','deg',]
    comment_theta = comment_theta+[str('%.2f' % (100*(0.087-h)))+'cm',str('%.2f' % (100*(0.087-h)))+'cm',]
    data_theta = data_theta+[h_0,theta_B_deg]

    # L_f output
    name_L_f = name_L_f+['h_0','L_2','L_3','L_f',]
    unit_L_f = unit_L_f+['m','m','m','m',]
    comment_L_f = comment_L_f+[(str('%.2f' % (100*(0.087-h)))+'cm',)*4]
    data_L_f = data_L_f+[h_0,L_2,L_3,L_f]

    # Energy transfer
    E_i = 0.5*h_P**2/(m_n*lambda_i**2) # unit J
    lambda_f = lambda_off*np.sin(theta_B_rad) 
    E_f = 0.5*h_P**2/(m_n*lambda_f**2) # unit J
    E = E_i - E_f # unit J
    # energy output
    name_E = name_E+['h_0','E',]
    unit_E = unit_E+['m','meV',]
    comment_E = comment_E+[(str('%.2f' % (100*(0.087-h)))+'cm',)*2]
    data_E = data_E+[h_0,E*1e3/(1.6e-19)]  
# Energy output
comment_E = _flatten(comment_E)
data_E = np.transpose(data_E)
filename = str(i+1).rjust(3,'3')+'_Energy transfer_E_i_'+str('%.3f' % (1e10*lambda_i)+'AA.dat')
sf.dataoutput(name_E,unit_E,comment_E,data_E,datapath,filename)

# theta_B output    
data_theta = np.transpose(data_theta)
filename = str(len(h_array)+1).rjust(3,'1')+'_theta_B.dat'
sf.dataoutput(name_theta,unit_theta,comment_theta,data_theta,datapath,filename)

# L_f output
comment_L_f = _flatten(comment_L_f)
data_L_f = np.transpose(data_L_f)
filename = str(len(h_array)+1).rjust(3,'2')+'_L_f.dat'
sf.dataoutput(name_L_f,unit_L_f,comment_L_f,data_L_f,datapath,filename)

print ('Run Successfully!')