forked from sgarrettroe/data_analysis
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathchooseLineshapeFunction.m
144 lines (121 loc) · 7.63 KB
/
chooseLineshapeFunction.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
function [out] = chooseLineshapeFunction(damping,s)
global wavenumbersToInvPs
c2 = [];
g = [];
switch damping,
case 'voigt'
% c2 = @(t,s) (t==0)/T2 + (Delta1)^2;
% g = @(t,s) t./T2 + (Delta1)^2/2.*t.^2;
c2 = @(t,s) (t==0)/s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2;
g = @(t,s) t./s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/2.*t.^2;
case '1fast'
c2 = @(t,s) (t==0)/s.T2;
g = @(t,s) t./s.T2;
case {'overdamped', '1exp'}
%overdamped exp(-t/tau)
c2 = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau1).*t);
g = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau1))^2.*(exp(-(1/s.tau1).*t)-1+(1/s.tau1).*t);
case {'1exp1fast'}
%overdamped exp(-t/tau)
c2 = @(t,s) (t==0)/s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau1).*t);
g = @(t,s) t./s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau1))^2.*(exp(-(1/s.tau1).*t)-1+(1/s.tau1).*t);
case {'1exp1slow'}
c2 = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau1).*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2;
g = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau1))^2.*(exp(-(1/s.tau1).*t)-1+(1/s.tau1).*t) ...
+ ((s.Delta2_cm*wavenumbersToInvPs*2*pi)^2).*t.^2/2;
case {'1exp1fast1slow'}
c2 = @(t,s) (t==0)/s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau1).*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2;
g = @(t,s) t./s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau1))^2.*(exp(-(1/s.tau1).*t)-1+(1/s.tau1).*t) ...
+ ((s.Delta2_cm*wavenumbersToInvPs*2*pi)^2).*t.^2/2;
case {'2exp1fast'}
% Developed for use with Zhe's SCN- data in ILs, where he has two
% inhomogeneous timescales, one longer, and one shorter.
% Expanded form of '1exp1fast' above --Tom
c2 = @(t,s) (t==0)/s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau1).*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau2).*t);
g = @(t,s) t./s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau1))^2.*(exp(-(1/s.tau1).*t)-1+(1/s.tau1).*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau2))^2.*(exp(-(1/s.tau2).*t)-1+(1/s.tau2).*t);
case {'2exp1fast1slow'}
c2 = @(t,s) (t==0)/s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau1).*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau2).*t) + (s.Delta3_cm*wavenumbersToInvPs*2*pi)^2;
g = @(t,s) t./s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau1))^2.*(exp(-(1/s.tau1).*t)-1+(1/s.tau1).*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau2))^2.*(exp(-(1/s.tau2).*t)-1+(1/s.tau2).*t) + ((s.Delta3_cm*wavenumbersToInvPs*2*pi)^2).*t.^2/2;
case {'2exp1slow'}
c2 = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau1).*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau2).*t) + (s.Delta3_cm*wavenumbersToInvPs*2*pi)^2;
g = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau1))^2.*(exp(-(1/s.tau1).*t)-1+(1/s.tau1).*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau2))^2.*(exp(-(1/s.tau2).*t)-1+(1/s.tau2).*t) + ((s.Delta3_cm*wavenumbersToInvPs*2*pi)^2).*t.^2/2;
case {'3exp1fast'}
c2 = @(t,s) (t==0)/s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau1).*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau2).*t) + (s.Delta3_cm*wavenumbersToInvPs*2*pi)^2.*exp(-(1/s.tau3).*t);
g = @(t,s) t./s.T2 + (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau1))^2.*(exp(-(1/s.tau1).*t)-1+(1/s.tau1).*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau2))^2.*(exp(-(1/s.tau2).*t)-1+(1/s.tau2).*t) ...
+ (s.Delta3_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau3))^2.*(exp(-(1/s.tau3).*t)-1+(1/s.tau3).*t);
case 'critical'
%critically damped (1+2t/tau)exp(-2t/tau)
c2 = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2.*(1+2*(1/s.tau1).*t).*exp(-2*(1/s.tau1).*t);
g = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/4/((1/s.tau1))^2.*exp(-2.*(1/s.tau1).*t) ...
.*(3 + 2*(1/s.tau1)*t + exp(2.*(1/s.tau1).*t).*(4*(1/s.tau1).*t - 3));
case '2expcrit'
g = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/4/((1/s.tau1))^2 ...
.*(3.*exp(-2.*(1/s.tau1).*t) + 2*(1/s.tau1)*t.*exp(-2.*(1/s.tau1).*t) + (4*(1/s.tau1).*t - 3)) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2/4/((1/s.tau2))^2 ...
.*(3.*exp(-2.*(1/s.tau2).*t) + 2*(1/s.tau2)*t.*exp(-2.*(1/s.tau2).*t) + (4*(1/s.tau2).*t - 3));
case '3expcrit'
g = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/4/((1/s.tau1))^2 ...
.*(3.*exp(-2.*(1/s.tau1).*t) + 2*(1/s.tau1)*t.*exp(-2.*(1/s.tau1).*t) + (4*(1/s.tau1).*t - 3)) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2/4/((1/s.tau2))^2 ...
.*(3.*exp(-2.*(1/s.tau2).*t) + 2*(1/s.tau2)*t.*exp(-2.*(1/s.tau2).*t) + (4*(1/s.tau2).*t - 3)) ...
+ (s.Delta3_cm*wavenumbersToInvPs*2*pi)^2/4/((1/s.tau3))^2 ...
.*(3.*exp(-2.*(1/s.tau3).*t) + 2*(1/s.tau3)*t.*exp(-2.*(1/s.tau3).*t) + (4*(1/s.tau3).*t - 3));
case '3exp1offcrit'
g = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/4/((1/s.tau1))^2 ...
.*(3.*exp(-2.*(1/s.tau1).*t) + 2*(1/s.tau1)*t.*exp(-2.*(1/s.tau1).*t) + (4*(1/s.tau1).*t - 3)) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2/4/((1/s.tau2))^2 ...
.*(3.*exp(-2.*(1/s.tau2).*t) + 2*(1/s.tau2)*t.*exp(-2.*(1/s.tau2).*t) + (4*(1/s.tau2).*t - 3)) ...
+ (s.Delta3_cm*wavenumbersToInvPs*2*pi)^2/4/((1/s.tau3))^2 ...
.*(3.*exp(-2.*(1/s.tau3).*t) + 2*(1/s.tau3)*t.*exp(-2.*(1/s.tau3).*t) + (4*(1/s.tau3).*t - 3)) ...
+ (s.Delta4_cm*wavenumbersToInvPs*2*pi)^2.*t.^2/2;
case '2exp'
g = @(t,s) (s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau1))^2 ...
.*(exp(-(1/s.tau1).*t) - 1 + (1/s.tau1)*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau2))^2 ...
.*(exp(-(1/s.tau2).*t) - 1 + (1/s.tau2)*t);
case '3exp'
g = @(t,s)(s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau1))^2 ...
.*(exp(-(1/s.tau1).*t) - 1 + (1/s.tau1)*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau2))^2 ...
.*(exp(-(1/s.tau2).*t) - 1 + (1/s.tau2)*t) ...
+ (s.Delta3_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau3))^2 ...
.*(exp(-(1/s.tau3).*t) - 1 + (1/s.tau3)*t);
case '3exp1off'
g = @(t,s)(s.Delta1_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau1))^2 ...
.*(exp(-(1/s.tau1).*t) - 1 + (1/s.tau1)*t) ...
+ (s.Delta2_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau2))^2 ...
.*(exp(-(1/s.tau2).*t) - 1 + (1/s.tau2)*t) ...
+ (s.Delta3_cm*wavenumbersToInvPs*2*pi)^2/((1/s.tau3))^2 ...
.*(exp(-(1/s.tau3).*t) - 1 + (1/s.tau3)*t) ...
+ (s.Delta4_cm*wavenumbersToInvPs*2*pi)^2.*t.^2/2;
case 'hynesform'
%fit c2 to hynes type function, double integrate with Mathematica
%to find g(t)
a1 = 0.3232;
k11 = 30.01; %ps-1
k12 = 17.41; %ps-1
a2 = 0.3378;
k2 = 8.270; %ps-1
a3 = 0.3455;
k3 = 1.897; %ps-1
%yuck
g = @(t,s) Delta^2*exp(-t.*(k12+k2+k3))/(k2^2*k3^2*(k11^2+k12^2)^2) ...
.*(a1*k2^2*k3^2.*exp((k2+k3).*t).*(cos(k11.*t).*(k12^2-k11^2)-2*k11*k12*sin(k11*t) ...
+ exp(k12.*t).*(k11^2*(k12.*t+1)+k12^2*(k12.*t-1))) ...
+ (k11^2+k12^2)^2.*exp(k12.*t).*(a3*k2^2*exp(k2.*t).*(exp(k3.*t).*(k3.*t-1)+1) ...
+a2*k3^2*exp(k3.*t).*(exp(k2.*t).*(k2*t-1)+1)));
otherwise
error('damping value is unknown');
end
out.g = g;
out.c2 = c2;