1 | function [alpha1 alpha2 alpha3] = physics_mcf(varargin) |
---|
2 | % PHYSICS_MCF - Computes momentum compaction up to 3rd order |
---|
3 | % |
---|
4 | % INPUTS |
---|
5 | % 1. Option Flag Display, NoDisplay |
---|
6 | % |
---|
7 | % OUPUTS |
---|
8 | % 2. pvalue - momentum compaction factor : alpha1, alpha2, alpha3 |
---|
9 | % |
---|
10 | % See Also mcf2, getmcf |
---|
11 | |
---|
12 | % |
---|
13 | %% Written by Laurent S. Nadolski |
---|
14 | |
---|
15 | DisplayFlag = 1; |
---|
16 | |
---|
17 | global THERING |
---|
18 | |
---|
19 | for i = length(varargin):-1:1 |
---|
20 | if strcmpi(varargin{i},'Display') |
---|
21 | DisplayFlag = 1; |
---|
22 | varargin{i} = []; |
---|
23 | elseif strcmpi(varargin{i},'NoDisplay') |
---|
24 | DisplayFlag = 0; |
---|
25 | varargin{i} = []; |
---|
26 | end |
---|
27 | end |
---|
28 | |
---|
29 | % energy variation for fitting mcf |
---|
30 | deltamin = -1e-4; |
---|
31 | deltamax = -deltamin; |
---|
32 | %delta = deltamin:5e-5:deltamax; |
---|
33 | delta = linspace(deltamin, deltamax, 7); |
---|
34 | delta4eval = linspace(deltamin, deltamax, 21); |
---|
35 | alpha = zeros(size(delta)); |
---|
36 | |
---|
37 | for k =1:length(delta) |
---|
38 | alpha(k) = mcf2(THERING,delta(k)); |
---|
39 | end |
---|
40 | |
---|
41 | % fit curve upto 3rd order |
---|
42 | porder =2; |
---|
43 | pvalue = polyfit(delta, alpha, porder); |
---|
44 | palpha = polyval(pvalue, delta4eval); |
---|
45 | |
---|
46 | |
---|
47 | % warning, alpha value are effective mcf, local slope for offmomentum |
---|
48 | % particle, so alpha = alpha1 + 2*alpha2*dp * 2*alpha3*dp^2 |
---|
49 | % This is the derivative of DT/T = alpha1*dp + alpha2*dp^2 + alpha3*dp^3; |
---|
50 | alpha1 = pvalue(end); |
---|
51 | alpha2 = pvalue(end-1)/2; |
---|
52 | alpha3 = pvalue(end-2)/3; |
---|
53 | |
---|
54 | if (0) % for debug |
---|
55 | fprintf('Momentum compaction factor \n mcf = %3.2e + 2x %3.2e dp + 3x %3.2e dp**2\n', ... |
---|
56 | alpha1, alpha2, alpha3) |
---|
57 | fprintf('Momentum compaction factor \n mcf = %3.2e (%3.2e) + %3.2e (%3.2e) dp + %3.2e (%3.2e) dp**2\n', ... |
---|
58 | [polymodel.Coefficients(end:-1:end-porder); polymodel.ParameterStd(end:-1:end-porder)]) |
---|
59 | end |
---|
60 | |
---|
61 | if DisplayFlag |
---|
62 | figure |
---|
63 | plot(delta*100,alpha,'k.') |
---|
64 | grid on; hold on; |
---|
65 | plot(delta4eval*100,palpha,'r') |
---|
66 | title(sprintf('Momentum compaction factor \n alpha1= %3.2e alpha2= %3.2e alpha3= %3.2e', ... |
---|
67 | alpha1, alpha2, alpha3)) |
---|
68 | xlabel('Energy (%)') |
---|
69 | ylabel('Momemtum compaction factor'); |
---|
70 | legend('Data', 'polyfit') |
---|
71 | end |
---|