利用二分法解方程求解次级谱仪

博士后出站报告附录

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

maindir = r'D:\001_CSNS\Summary\Data Treatment\TOF_DATA_Treatment\Python\20220406_TOF_PSD_data_2-10 AA_5 detectors\Second_Spec_data'
datapath = sf.mkdir(maindir+r'\Norm Data Output')
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 = np.linspace(2e-10,10e-10,5) #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
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;
y_off = np.linspace(-0.2,+0.2,201) # the distance from the pixel to the center of detector;
#################################

R = 2.5 # unit m, radius of analyzer sphere
r_0 = 0.27 #unit m, radius of detector cylinder


# Detector pixel information
theta_Dp_deg = np.linspace(+40,+140,100)
theta_Dp_rad = np.radians(theta_Dp_deg)
x_off = r_0*theta_Dp_rad # arc positon of detector pixel
h_0 = y_off +h_d # vertical distance from pixel to origin

# Array assignment
phi_deg = np.linspace(0,180,181)
phi_rad = np.radians(phi_deg)

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

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

name_t_f = []
unit_t_f = []
comment_t_f = []
data_t_f = []

for k in range(len(h_array)):
    h = h_array[k]
    psi_rad = sf.arctan((h_0-h)/r_0) 
    psi_deg = np.degrees(psi_rad)
    rho = np.sqrt(r_0**2+(h_0-h)**2) # the distance from pixel to center of analyzer
    phi_rad_01 = []
    phi_rad_02 = []
    phi_deg_01 = []
    phi_deg_02 = []
    for i in range(len(y_off)): # y is y_off
        y_i = []
        def Delta_y(phi):  # define a multi variable function, search for the zero point of the function.
            L_2i = np.sqrt(R**2+h**2-2*R*h*np.cos(phi))
            L_3i = np.sqrt(R**2+rho[i]**2-2*R*rho[i]*np.sin(phi-psi_rad[i]))
            y = L_3i*(R**2+L_2i**2-h**2)-L_2i*(R**2+L_3i**2-rho[i]**2)
            return y
        for j in range(len(phi_rad)):  # x is phi    
            y_ij = Delta_y(phi_rad[j])
            y_i.append(y_ij)  
        a = y_i.index(min(y_i)) 
        phi_rad_01i = sf.biSection(phi_rad[0], phi_rad[a], 0.001, Delta_y)
        phi_rad_02i = sf.biSection(phi_rad[a], phi_rad[-1], 0.001, Delta_y)
        phi_rad_01.append(phi_rad_01i)
        phi_rad_02.append(phi_rad_02i)
        phi_deg_01i = np.degrees(phi_rad_01i)
        phi_deg_02i = np.degrees(phi_rad_02i)
        phi_deg_01.append(phi_deg_01i)
        phi_deg_02.append(phi_deg_02i)

    phi_deg_01 = np.array(phi_deg_01)
    phi_deg_02 = np.array(phi_deg_02)       
    name = _flatten(['y_off','phi_01','phi_02'])
    unit = _flatten(['m','deg','deg'])
    comment = _flatten([(str('%.2f' % (100*(0.087-h)))+'cm',)*len(name)])
    data = np.c_[y_off,phi_deg_01,phi_deg_02]
    filename = str(k+1).rjust(3,'1')+'_phi-solution_'+str('%.2f' % (100*(0.087-h)))+'cm.dat'
    sf.dataoutput(name,unit,comment,data,datapath,filename)

    alpha_rad = sf.arctan(h*np.sin(phi_rad_01)/(R-h*np.cos(phi_rad_01)))
    alpha_deg = np.degrees(alpha_rad)
    theta_B_deg = 90 - alpha_deg
    theta_B_rad = np.radians(theta_B_deg)
    gamma_deg = phi_deg_01+alpha_deg-90
    beta_deg = 180-phi_deg_01+alpha_deg
    y_ana = h - R*np.cos(phi_rad_01)
    z_ana = R*np.sin(phi_rad_01)

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

    lambda_f = 2*d*np.sin(theta_B_rad) 
    v_f = h_P/(lambda_f*m_n)
    t_f = L_f/v_f

    # alldata output
    name = _flatten(['gamma','alpha','theta_B','phi','beta','h_0','y_off','y_ana','z_ana','L_f','L_2','t_f'])
    unit = _flatten([('deg',)*5,('m',)*6,'s'])
    comment = _flatten([(str('%.2f' % (100*(0.087-h)))+'cm',)*len(name)])
    data = np.c_[gamma_deg,alpha_deg,theta_B_deg,phi_deg_01,beta_deg,h_0,y_off,y_ana,z_ana,L_f,L_2,t_f]
    filename = str(k+1).rjust(3,'2')+'_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]

    # theta_B output
    name_t_f = name_t_f+['h_0','t_f',]
    unit_t_f = unit_t_f+['m','s',]
    comment_t_f = comment_t_f+[str('%.2f' % (100*(0.087-h)))+'cm',str('%.2f' % (100*(0.087-h)))+'cm',]
    data_t_f = data_t_f+[h_0,t_f]
    print ('Loop '+str(k)+'-finished!')

data_theta = np.transpose(data_theta)
filename = str(len(h_array)+1).rjust(3,'3')+'_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,'4')+'_L_f.dat'
sf.dataoutput(name_L_f,unit_L_f,comment_L_f,data_L_f,datapath,filename)

comment_t_f = _flatten(comment_t_f)
data_t_f = np.transpose(data_t_f)
filename = str(len(h_array)+1).rjust(3,'4')+'_t_f.dat'
sf.dataoutput(name_t_f,unit_t_f,comment_t_f,data_t_f,datapath,filename)

print ('Run Successfully!')