Subversion Repositories lagranto.ecmwf

Rev

Details | Last modification | View Log | RSS feed

Rev Author Line No. Line
3 michaesp 1
% ----------------------------------------------------------------------
2
% Read a trajectory file
3
% ----------------------------------------------------------------------
4
 
5
function out = read_trajectory(filename);
6
 
7
% Open the trajectory file
8
fid = fopen(filename,'r')
9
 
10
% Read the info line (tra.info)
11
line     = textscan(fid,'%s',1,'delimiter','\n');
12
tra.info = line{1};
13
 
14
% Read the line with field description (skip empty lines)
15
textscan(fid,'%s',1,'delimiter','\n');
16
line     = textscan(fid,'%s',1,'delimiter','\n');
17
textscan(fid,'%s',1,'delimiter','\n');
18
line     = char(line{1});
19
 
20
% Get the names of the fields (tra.field, tra.nfield)
21
indl       = 1;
22
indr       = 1;
23
tra.nfield = 0;
24
while ( indl <= length(line) )
25
    while ( line(indl) == ' '  )
26
       indl=indl+1;
27
    end
28
    indr=indl;
29
    while ( (line(indr) ~= ' ') & (indr < length(line)) )
30
       indr=indr+1;
31
    end
32
    tra.nfield            = tra.nfield + 1;
33
    tra.field(tra.nfield) = cellstr(line(indl:indr));
34
    indl = indr + 1;
35
end
36
 
37
% Get the number of times for a single trajectory (tra.ntime)
38
fpos = ftell(fid);
39
textscan(fid,'%s',1,'delimiter','\n');
40
line = textscan(fid,'%s',1,'delimiter','\n');
41
tra.ntime = 0;
42
while ( ~ strcmp(cellstr(line{1}),'') )
43
    line = textscan(fid,'%s',1,'delimiter','\n');   
44
    tra.ntime = tra.ntime + 1;
45
end    
46
fseek(fid,fpos,'bof');
47
 
48
% Read the trajectories
49
format = '';
50
for i=1:tra.nfield
51
   format = [ format '%n' ];
52
end
53
lines = textscan(fid,format);
54
 
55
% Set the total number of trajectories
56
tra.ntra = length(lines{1})/tra.ntime;
57
 
58
% Attribute the field names to the columns
59
for i=1:tra.nfield
60
    name = char(tra.field(i));
61
    tra.(name) = lines{i};
62
end
63
 
64
% Attribute a unique label to each trajectory
65
tra.label = 1+floor( (0:tra.ntra*tra.ntime-1) / tra.ntime )'; 
66
 
67
% Close trajectory file
68
fclose(fid);
69
 
70
% Return the output structure
71
out = tra;
72
 
73