计算不同像素点处对应的动量分辨

博士后出站报告附录

# -*- 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 # 将二维列表转换为一维列表
from scipy import interpolate

maindir = r'D:\CSNS\Summary\Data Treatment\python\Momentum resolution\009_Q_resolution_XYZ_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

# incident beam parameters
L_i = L_1 = 90 # unit m
lambda_i_array = np.array([2,3,4,5,6.267,7,8,9,10])*1e-10 #unit m
theta_i_deg = 90 # unit deg
theta_i_rad = np.radians(theta_i_deg)
phi_i_deg = 0 # unit deg
phi_i_rad = np.radians(phi_i_deg)

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

# Primary Spec
dtheta_i_deg = 4.9
dtheta_i_rad = np.radians(dtheta_i_deg) 

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

theta_Dm_deg = np.linspace(-10,-160,16)
theta_Dp_deg = np.linspace(10,160,16)
theta_D_deg = np.concatenate([theta_Dm_deg,theta_Dp_deg]) # horizontal cover angle of the detectors

# Detector parameters
h_0_min = 0.129
h_0_max = 0.277
h_0s = np.linspace(h_0_min,h_0_max,round(100*(h_0_max-h_0_min))+1)
h_0_inters = []
for i0 in range(len(h_0s)-1): 
    h_0_inter = (h_0s[i0+1]+h_0s[i0])/2
    h_0_inters.append(h_0_inter)
# 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,31) # 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)

dtheta_i_deg = 4.9
dtheta_i_rad = np.radians(dtheta_i_deg)
twotheta_degs = np.linspace(0,180,181)
twotheta_rads = np.radians(twotheta_degs)
cos_2thetas = np.cos(twotheta_rads)
dcos_2thetas = abs(np.sin(twotheta_rads))*dtheta_i_rad
'''
name = ['2theta','dcos_2theta']
unit = ['deg','none']
comment = ['PriSpec','PriSpec']
data = np.c_[twotheta_degs,dcos_2thetas]
filename = '000-dcos-vs-twotheta.dat'
sf.dataoutput(name, unit, comment, data, datapath, filename)
'''

for k in range(len(lambda_i_array)):
    lambda_i = lambda_i_array[k]
    lambda_f = 6.267e-10 # unit AA

    Qs = np.sqrt(4*np.pi**2*(1/lambda_i**2+1/lambda_f**2-2*cos_2thetas/(lambda_i*lambda_f)))
    dQ1s = abs(sf.dydx(cos_2thetas, Qs)*dcos_2thetas)

    name_Q = ['h0']
    unit_Q = ['m']
    comment_Q = [f'{lambda_i*1e10}AA']
    data_Q = [h_0_inters]

    name_Q_reso1 = ['h0']
    unit_Q_reso1 = ['m']
    comment_Q_reso1 = [f'{lambda_i*1e10}AA']
    data_Q_reso1 = [h_0_inters]

    name_Q_reso2 = ['h0']
    unit_Q_reso2 = ['m']
    comment_Q_reso2 = [f'{lambda_i*1e10}AA']
    data_Q_reso2 = [h_0_inters]

    name_Q_reso = ['h0']
    unit_Q_reso = ['m']
    comment_Q_reso = [f'{lambda_i*1e10}AA']
    data_Q_reso = [h_0_inters]

    phi_f_X = []
    h_0_inters_Y = []
    Q_reso1s_Z = []
    Q_reso1_ranges_Z = []
    Q_reso2s_Z = []
    Q_resos_Z = []
    Q_Z = []

    for j in range(len(theta_Dp_deg)):
        phi_f_deg = theta_Dp_deg[j]
        phi_f_rad = np.radians(phi_f_deg)   
        Q_centers = []
        Q_reso1s = []
        Q_reso1_ranges = []
        Q_reso2s = []
        Q_resos = []
        for i0 in range(len(h_0s)-1): # 对于每一个h_0s[i0]，我们都可以得到一个Q_center和Q_reso
            h_0_inter = (h_0s[i0+1]+h_0s[i0])/2
            Qs = []
            Q_reso1s_dis = []
            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)

                lambda_f = lambda_off*np.sin(theta_B_rad)
                theta_f_deg = 90 - gamma_deg
                theta_f_rad = np.radians(theta_f_deg)
                cos_2theta = np.sin(theta_i_rad)*np.sin(theta_f_rad)*np.cos(phi_i_rad-phi_f_rad)+np.cos(theta_i_rad)*np.cos(theta_f_rad) 
                twotheta_rad = sf.arccos(cos_2theta)
                twotheta_deg = np.degrees(twotheta_rad)
                Q = np.sqrt(4*np.pi**2*((1/(lambda_i**2))+(1/(lambda_f**2))-2*cos_2theta/(lambda_i*lambda_f)))
                Q_reso1_array = interpolate.interp1d(twotheta_rads, dQ1s, kind='linear')(twotheta_rad)

                for i1 in range(len(h_0)):
                    if h_0[i1]>=h_0s[i0] and h_0[i1]<=h_0s[i0+1]:
                        Qs.append(Q[i1])
                        Q_reso1s_dis.append(Q_reso1_array[i1])

            Q_reso1_min = min(Q_reso1s_dis)
            Q_reso1_max = max(Q_reso1s_dis)
            Q_reso1 = (Q_reso1_min + Q_reso1_max)/2
            Q_reso1_range = Q_reso1_max - Q_reso1_min
            Q_reso1s.append(Q_reso1)
            Q_reso1_ranges.append(Q_reso1_range)

            Q_min = min(Qs)
            Q_max = max(Qs)
            Q_center = (Q_max + Q_min)/2
            Q_reso2 = Q_max-Q_min
            Q_centers.append(Q_center)
            Q_reso2s.append(Q_reso2)

        Q_resos = np.sqrt(np.array(Q_reso1s)**2+np.array(Q_reso2s)**2)
        Q_resos = np.array(Q_resos)
        Q_reso1s = np.array(Q_reso1s)
        Q_reso1_ranges = np.array(Q_reso1_ranges)
        Q_reso2s = np.array(Q_reso2s)


        phi_f_X = np.concatenate([phi_f_X,np.ones(len(h_0_inters))*phi_f_deg])
        h_0_inters_Y = np.concatenate([h_0_inters_Y,h_0_inters])
        Q_reso1s_Z = np.concatenate([Q_reso1s_Z,Q_reso1s])
        Q_reso1_ranges_Z = np.concatenate([Q_reso1_ranges_Z,Q_reso1_ranges])
        Q_reso2s_Z = np.concatenate([Q_reso2s_Z,Q_reso2s])
        Q_resos_Z = np.concatenate([Q_resos_Z,Q_resos])
        Q_Z = np.concatenate([Q_Z,Q_centers])

        name_Q = name_Q+['Q','dQ_half',]
        unit_Q = unit_Q+['1/A','1/A',]
        comment_Q = comment_Q+[(f'{phi_f_deg}deg',)*2]
        data_Q = data_Q+[np.array(Q_centers)*1e-10,np.array(Q_resos)*1e-10/2]

        name_Q_reso1 = name_Q_reso1+['Q_reso1','Q_reso1_range_half',]
        unit_Q_reso1 = unit_Q_reso1+['1/A','1/A']
        comment_Q_reso1 = comment_Q_reso1+[(f'{phi_f_deg}deg',)*2]
        data_Q_reso1 = data_Q_reso1+[np.array(Q_reso1s)*1e-10,np.array(Q_reso1_ranges)*1e-10/2]

        name_Q_reso2 = name_Q_reso2+['Q_reso2',]
        unit_Q_reso2 = unit_Q_reso2+['1/A',]
        comment_Q_reso2 = comment_Q_reso2+[(f'{phi_f_deg}deg',)*1]
        data_Q_reso2 = data_Q_reso2+[np.array(Q_reso2s)*1e-10]

        name_Q_reso = name_Q_reso+['Q_reso',]
        unit_Q_reso = unit_Q_reso+['1/A',]
        comment_Q_reso = comment_Q_reso+[(f'{phi_f_deg}deg',)*1]
        data_Q_reso = data_Q_reso+[np.array(Q_resos)*1e-10]

    # momentum output
    name_Q = _flatten(name_Q)
    unit_Q = _flatten(unit_Q)
    comment_Q = _flatten(comment_Q)
    data_Q = np.transpose(data_Q)
    filename = str(k+1).rjust(3,'0')+'_momentum transfer_'+f'lambda_i = {lambda_i*1e10}AA.dat'
    sf.dataoutput(name_Q,unit_Q,comment_Q,data_Q,datapath,filename)

    # momentum resolution1
    name_Q_reso1 = _flatten(name_Q_reso1)
    unit_Q_reso1 = _flatten(unit_Q_reso1)
    comment_Q_reso1 = _flatten(comment_Q_reso1)
    data_Q_reso1 = np.transpose(data_Q_reso1)
    filename = str(k+1).rjust(3,'0')+'_momentum resolution1_'+f'lambda_i = {lambda_i*1e10}AA.dat'
    sf.dataoutput(name_Q_reso1,unit_Q_reso1,comment_Q_reso1,data_Q_reso1,datapath,filename)

    # momentum resolution2
    name_Q_reso2 = _flatten(name_Q_reso2)
    unit_Q_reso2 = _flatten(unit_Q_reso2)
    comment_Q_reso2 = _flatten(comment_Q_reso2)
    data_Q_reso2 = np.transpose(data_Q_reso2)
    filename = str(k+1).rjust(3,'0')+'_momentum resolution2_'+f'lambda_i = {lambda_i*1e10}AA.dat'
    sf.dataoutput(name_Q_reso2,unit_Q_reso2,comment_Q_reso2,data_Q_reso2,datapath,filename)

    # momentum resolution
    name_Q_reso = _flatten(name_Q_reso)
    unit_Q_reso = _flatten(unit_Q_reso)
    comment_Q_reso = _flatten(comment_Q_reso)
    data_Q_reso = np.transpose(data_Q_reso)
    filename = str(k+1).rjust(3,'0')+'_momentum resolution_'+f'lambda_i = {lambda_i*1e10}AA.dat'
    sf.dataoutput(name_Q_reso,unit_Q_reso,comment_Q_reso,data_Q_reso,datapath,filename)

    # Q and Q resolution-XYZ
    name_Q_all = ['phi_f','h_0','Q','Q_reso_half','Q_reso1','Q_reso1_range_half','Q_reso2']
    unit_Q_all = ['deg','m','1/A','1/A','1/A','1/A','1/A']
    comment_Q_all = _flatten([(f'{lambda_i*1e10}AA',)*7])
    data_Q_all = np.c_[phi_f_X,h_0_inters_Y,Q_Z*1e-10,Q_resos_Z*1e-10/2,Q_reso1s_Z*1e-10,Q_reso1_ranges_Z*1e-10/2,Q_reso2s_Z*1e-10]
    filename = str(k+1).rjust(3,'0')+'_momentum_XYZ_Data_'+f'lambda_i = {lambda_i*1e10}AA.dat'
    sf.dataoutput(name_Q_all,unit_Q_all,comment_Q_all,data_Q_all,datapath,filename)
    print (f'Loop-{k} finished')
print ('Run Successfully!')