-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathClosedCalc.m
More file actions
69 lines (69 loc) · 2 KB
/
Copy pathClosedCalc.m
File metadata and controls
69 lines (69 loc) · 2 KB
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
%% Finction to run a simulation with some set of parameters and return the fit
% accuracy
function [raw,t,freqAxis] = ClosedCalc(base)
import HypWright.*
import HypWright.Models.*
%% Initilaize variable
if isempty(base)
base = struct();
end
Default = struct('gamma', 67.262e6, 'readBandwidth', 4096, 'rfBandwidth', 5000,...
'nPoints', 2048, 'endTime', 100, 'T1a', 56, 'T2a', 0.02, 'T1b', 30,...
'T2b', 0.02, 'ppma', -7e-6,'ppmb', 7e-6, 'Kab', 0.1, 'flipAngle', 20,...
'TR', 2,'A', TwoSiteExchange(), 'verbose', false,'B0',3.0);
tmpNames = fieldnames(Default);
for i = 1:numel(tmpNames)
if ~isfield(base,tmpNames{i})
base.(tmpNames{i}) = Default.(tmpNames{i});
end
end
base.flipAngle = base.flipAngle*pi/180;
gamma = base.gamma;
readBandwidth = base.readBandwidth;
rfBandwidth = base.rfBandwidth;
nPoints = base.nPoints;
endTime = base.endTime;
T1a = base.T1a;
T2a = base.T2a;
T1b = base.T1b;
ppma = base.ppma;
ppmb = base.ppmb;
Kab = base.Kab;
flipAngle = base.flipAngle;
TR = base.TR;
B0 = base.B0;
verbose = base.verbose;
%TODO add input validation;
%% Init World
world = HypWright.World.getWorld;
world.initWorld()
world.setB0([0;0;B0])
Spin = TwoSiteExchangeGroup([0;0;1;0;0;0],[0;0;0;0;0;0],...
T1a,T2a,ppma,T1b,T2a,ppmb,gamma,1,Kab,[]);
V = Voxel([0;0;0],Spin);
world.addVoxel(V);
%% Build Pulse Sequence
PS = PulseSequence;
t = 0:TR:endTime;
ADC = zeros(length(t),nPoints);
for i = 1:length(t)
Pulse = SincPulse(t(i)+0.001,rfBandwidth,flipAngle/(gamma),gamma*B0,[],...
sprintf('Excitation%d',1));
PS.addPulse(Pulse)
ADC(i,:) = Pulse.endTime:1/readBandwidth:Pulse.endTime+(nPoints-1)/readBandwidth;
end
world.setPulseSequence(PS)
%% Calculate
world.calculate(t(end)+10);
FID = zeros(length(t),nPoints);
for i = 1:length(t)
[FID(i,:), freqAxis] = world.evaluate(ADC(i,:),-gamma*B0);
end
raw = FID.';
t = linspace(ADC(1,floor(end/2)),ADC(end,floor(end/2)),size(ADC,1));
if (verbose)
figure
surf(freqAxis,t,abs(fftshift(fft(FID,[],2),2)));
drawnow
end
end