-
Notifications
You must be signed in to change notification settings - Fork 7
/
magia_read_input.m
75 lines (64 loc) · 1.61 KB
/
magia_read_input.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
function input = magia_read_input(filename)
if(isempty(filename))
error('Could not read input data.');
end
if(iscell(filename))
filename = filename{1};
end
if(~exist(filename,'file'))
error('%s: Could not read file %s because it does not exist.',subject,filename);
end
[~,~,ext] = fileparts(filename);
if(strcmpi(ext,'dft'))
input = read_input_dft(filename);
else
fid = fopen(filename,'r');
% Get rid of comments (#) in the beginning
l = fgetl(fid);
while(strcmp(l(1),'#'))
l = fgetl(fid);
if(isempty(l))
l = fgetl(fid);
break;
end
end
% Parse the input data
i = 0;
while(l~=-1)
if(ischar(l))
i = i + 1;
h = strsplit(l,' ');
if(length(h)==1)
h = strsplit(l,' ');
end
h(cellfun(@isempty,h)) = [];
input(i,1) = str2double(h{1}); % time
input(i,2) = str2double(h{2}); % radioactivity
end
l = fgetl(fid);
end
fclose(fid);
end
% Remove measurements before time zero
k = input(:,1) >= 0;
input = input(k,:);
k = ~isnan(input(:,2));
input = input(k,:);
maxT = max(input(:,1));
if(maxT > 300)
input(:,1) = input(:,1)/60;
end
if(input(1,1)>0)
input = [0 0;input];
end
zero_time_idx = input(:,1) == 0;
if(sum(zero_time_idx)>1)
last_zero_time_idx = find(zero_time_idx,1,'last');
input = input(last_zero_time_idx:end,:);
end
% Throw an error if there are duplicate time points
T = tabulate(input(:,1));
if(any(T(:,2) > 1))
error('Found duplicate time-points from the input file %s',filename);
end
end