1 | %'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them
|
---|
2 | %---------------------------------------------------------------------
|
---|
3 | %
|
---|
4 | % OUTPUT:
|
---|
5 | % Scal: matlab vector representing the scalar values (length nbvec defined by var_read)
|
---|
6 | % if no input , Scal=list of programmed scalar names (to put in menus)
|
---|
7 | % if only the sclar name is put as input, vec_A=type of scalar, which can be:
|
---|
8 | % 'vel': scalar calculated solely from velocity components
|
---|
9 | % 'der': needs spatial derivatives
|
---|
10 | % 'var': the scalar name directly corresponds to a field name in the netcdf files
|
---|
11 | % error: error flag
|
---|
12 | % error = 0; OK
|
---|
13 | % error = 1; the prescribed scalar cannot be read or calculated from available fields
|
---|
14 | %
|
---|
15 | % INPUT:
|
---|
16 | % ScalName: string representing the name of the scalar
|
---|
17 | % DataIn: structure representing the field, as defined in check_field_srtructure.m
|
---|
18 | %
|
---|
19 | % FUNCTION related
|
---|
20 | % varname_generator.m: determines the field names to read in the netcdf file, depending on the scalar
|
---|
21 |
|
---|
22 | function [DataOut,errormsg]=calc_field(FieldName,DataIn)
|
---|
23 |
|
---|
24 | %list of defined scalars to display in menus (in addition to 'ima_cor').
|
---|
25 | % a type is associated to each scalar:
|
---|
26 | % 'discrete': related to the individual velocity vectors, not interpolated by patch
|
---|
27 | % 'vel': calculated from velocity components, continuous field (interpolated with velocity)
|
---|
28 | % 'der': needs spatial derivatives
|
---|
29 | % 'var': the scalar name corresponds to a field name in the netcdf files
|
---|
30 | % a specific variable name for civ1 and civ2 fields are also associated, if
|
---|
31 | % '' the scalar is calculated from other fields, as explicited below
|
---|
32 | %list_scal={title, type, civ1 var name,civ2 var name}
|
---|
33 | list_field={'velocity';...%image correlation corresponding to a vel vector
|
---|
34 | 'ima_cor';...%image correlation corresponding to a vel vector
|
---|
35 | 'norm_vel';...%norm of the velocity
|
---|
36 | 'vort';...%vorticity
|
---|
37 | 'div';...%divergence
|
---|
38 | 'strain';...%rate of strain
|
---|
39 | 'u';... %u velocity component
|
---|
40 | 'v';... %v velocity component
|
---|
41 | 'w';... %w velocity component
|
---|
42 | 'w_normal';... %w velocity component normal to the plane
|
---|
43 | 'error'}; %error associated to a vector (for stereo or patch)
|
---|
44 | errormsg=[]; %default error message
|
---|
45 | if ~exist('FieldName','var')
|
---|
46 | DataOut=list_field;% gives the list of possible fields in the absence of input
|
---|
47 | else
|
---|
48 | if ~exist('DataIn','var')
|
---|
49 | DataIn=[];
|
---|
50 | end
|
---|
51 | if ischar(FieldName)
|
---|
52 | FieldName={FieldName};
|
---|
53 | end
|
---|
54 | if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
|
---|
55 | nbcoord=3;
|
---|
56 | else
|
---|
57 | nbcoord=2;
|
---|
58 | end
|
---|
59 | DataOut=DataIn; %reporduce global attribute
|
---|
60 | ListVarName={};
|
---|
61 | ValueList={};
|
---|
62 | RoleList={};
|
---|
63 | units_cell={};
|
---|
64 | for ilist=1:length(FieldName)
|
---|
65 | [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}
|
---|
66 | ListVarName=[ListVarName VarName];
|
---|
67 | ValueList=[ValueList Value];
|
---|
68 | RoleList=[RoleList Role];
|
---|
69 | units_cell=[units_cell units];
|
---|
70 | end
|
---|
71 | %erase previous data (except coordinates)
|
---|
72 | for ivar=nbcoord+1:length(DataOut.ListVarName)
|
---|
73 | VarName=DataOut.ListVarName{ivar};
|
---|
74 | DataOut=rmfield(DataOut,VarName);
|
---|
75 | end
|
---|
76 | DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
|
---|
77 | if isfield(DataOut,'VarDimName')
|
---|
78 | DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
|
---|
79 | else
|
---|
80 | errormsg='element .VarDimName missing in input data';
|
---|
81 | return
|
---|
82 | end
|
---|
83 | DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
|
---|
84 | %append new data
|
---|
85 | DataOut.ListVarName=[DataOut.ListVarName ListVarName];
|
---|
86 | for ivar=1:length(ListVarName)
|
---|
87 | DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
|
---|
88 | DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
|
---|
89 | DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
|
---|
90 | eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])
|
---|
91 | end
|
---|
92 | end
|
---|
93 |
|
---|
94 | %%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
|
---|
95 | function [VarName,ValCell,Role,units_cell]=velocity(DataIn)
|
---|
96 | VarName={};
|
---|
97 | ValCell={};
|
---|
98 | Role={};
|
---|
99 | units_cell={};
|
---|
100 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
101 | units=[DataIn.CoordUnit '/' DataIn.TimeUnit];
|
---|
102 | else
|
---|
103 | units='pixel';
|
---|
104 | end
|
---|
105 | if isfield(DataIn,'U')
|
---|
106 | VarName=[VarName {'U'}];
|
---|
107 | ValCell=[ValCell {DataIn.U}];
|
---|
108 | Role=[Role {'vector_x'}];
|
---|
109 | units_cell=[units_cell {units}];
|
---|
110 | end
|
---|
111 | if isfield(DataIn,'V')
|
---|
112 | VarName=[VarName {'V'}];
|
---|
113 | ValCell=[ValCell {DataIn.V}];
|
---|
114 | Role=[Role {'vector_y'}];
|
---|
115 | units_cell=[units_cell {units}];
|
---|
116 | end
|
---|
117 | if isfield(DataIn,'W')
|
---|
118 | VarName=[VarName {'W'}];
|
---|
119 | ValCell=[ValCell {DataIn.W}];
|
---|
120 | Role=[Role {'vector_z'}];
|
---|
121 | units_cell=[units_cell {units}];
|
---|
122 | end
|
---|
123 | if isfield(DataIn,'F')
|
---|
124 | VarName=[VarName {'F'}];
|
---|
125 | ValCell=[ValCell {DataIn.F}];
|
---|
126 | Role=[Role {'warnflag'}];
|
---|
127 | units_cell=[units_cell {[]}];
|
---|
128 | end
|
---|
129 | if isfield(DataIn,'FF')
|
---|
130 | VarName=[VarName,{'FF'}];
|
---|
131 | ValCell=[ValCell {DataIn.FF}];
|
---|
132 | Role=[Role {'errorflag'}];
|
---|
133 | units_cell=[units_cell {[]}];
|
---|
134 | end
|
---|
135 |
|
---|
136 | %%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%
|
---|
137 | function [VarName,ValCell,Role,units]=ima_cor(DataIn)
|
---|
138 | VarName={};
|
---|
139 | ValCell={};
|
---|
140 | Role={};
|
---|
141 | units={};
|
---|
142 | if isfield(DataIn,'C')
|
---|
143 | VarName{1}='C';
|
---|
144 | ValCell{1}=DataIn.C;
|
---|
145 | Role={'ancillary'};
|
---|
146 | units={[]};
|
---|
147 | end
|
---|
148 |
|
---|
149 | %%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%
|
---|
150 | function [VarName,ValCell,Role,units]=norm_vel(DataIn)
|
---|
151 | VarName={};
|
---|
152 | ValCell={};
|
---|
153 | Role={};
|
---|
154 | units={};
|
---|
155 | if isfield(DataIn,'U') && isfield(DataIn,'V')
|
---|
156 | VarName{1}='norm_vel';
|
---|
157 | ValCell{1}=DataIn.U.*DataIn.U+ DataIn.V.*DataIn.V;
|
---|
158 | if isfield(DataIn,'W') && isequal(size(DataIn.W),size(DataIn.U))
|
---|
159 | ValCell{1}=ValCell{1}+DataIn.W.*DataIn.W;
|
---|
160 | end
|
---|
161 | ValCell{1}=sqrt(ValCell{1});
|
---|
162 | Role{1}='scalar';
|
---|
163 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
164 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
165 | else
|
---|
166 | units={'pixel'};
|
---|
167 | end
|
---|
168 | end
|
---|
169 |
|
---|
170 |
|
---|
171 |
|
---|
172 | %%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%
|
---|
173 | function [VarName,ValCell,Role,units]=vort(DataIn)
|
---|
174 | VarName={};
|
---|
175 | ValCell={};
|
---|
176 | Role={};
|
---|
177 | units={};
|
---|
178 | if isfield(DataIn,'DjUi')
|
---|
179 | VarName{1}='vort';
|
---|
180 | ValCell{1}=DataIn.DjUi(:,1,2)-DataIn.DjUi(:,2,1); %vorticity
|
---|
181 | siz=size(ValCell{1});
|
---|
182 | ValCell{1}=reshape(ValCell{1},siz(1),1);
|
---|
183 | Role{1}='scalar';
|
---|
184 | if isfield(DataIn,'TimeUnit')
|
---|
185 | units={[DataIn.TimeUnit '-1']};
|
---|
186 | else
|
---|
187 | units={[]};
|
---|
188 | end
|
---|
189 | end
|
---|
190 |
|
---|
191 | %%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%
|
---|
192 | function [VarName,ValCell,Role,units]=div(DataIn)
|
---|
193 | VarName={};
|
---|
194 | ValCell={};
|
---|
195 | Role={};
|
---|
196 | units={};
|
---|
197 | if isfield(DataIn,'DjUi')
|
---|
198 | VarName{1}='div';
|
---|
199 | ValCell{1}=DataIn.DjUi(:,1,1)+DataIn.DjUi(:,2,2); %DUDX+DVDY
|
---|
200 | siz=size(ValCell{1});
|
---|
201 | ValCell{1}=reshape(ValCell{1},siz(1),1);
|
---|
202 | Role{1}='scalar';
|
---|
203 | if isfield(DataIn,'TimeUnit')
|
---|
204 | units={[DataIn.TimeUnit '-1']};
|
---|
205 | else
|
---|
206 | units={[]};
|
---|
207 | end
|
---|
208 | end
|
---|
209 |
|
---|
210 | %%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
|
---|
211 | function [VarName,ValCell,Role,units]=strain(DataIn)
|
---|
212 | VarName={};
|
---|
213 | ValCell={};
|
---|
214 | Role={};
|
---|
215 | units={};
|
---|
216 | if isfield(DataIn,'DjUi')
|
---|
217 | VarName{1}='strain';
|
---|
218 | ValCell{1}=DataIn.DjUi(:,1,2)+DataIn.DjUi(:,2,1);%DVDX+DUDY
|
---|
219 | siz=size(ValCell{1});
|
---|
220 | ValCell{1}=reshape(ValCell{1},siz(1),1);
|
---|
221 | if isfield(DataIn,'TimeUnit')
|
---|
222 | units={[DataIn.TimeUnit '-1']};
|
---|
223 | else
|
---|
224 | units={[]};
|
---|
225 | end
|
---|
226 | end
|
---|
227 |
|
---|
228 | %%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
|
---|
229 | function [VarName,ValCell,Role,units]=u(DataIn)
|
---|
230 | VarName={};
|
---|
231 | ValCell={};
|
---|
232 | Role={};
|
---|
233 | units={};
|
---|
234 | if isfield(DataIn,'U')
|
---|
235 | VarName{1}='U';
|
---|
236 | ValCell{1}=DataIn.U;
|
---|
237 | Role{1}='scalar';
|
---|
238 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
239 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
240 | else
|
---|
241 | units={'pixel'};
|
---|
242 | end
|
---|
243 | end
|
---|
244 |
|
---|
245 | %%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
|
---|
246 | function [VarName,ValCell,Role,units]=v(DataIn)
|
---|
247 | VarName={};
|
---|
248 | ValCell={};
|
---|
249 | Role={};
|
---|
250 | units={};
|
---|
251 | if isfield(DataIn,'V')
|
---|
252 | VarName{1}='V';
|
---|
253 | ValCell{1}=DataIn.V;
|
---|
254 | Role{1}='scalar';
|
---|
255 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
256 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
257 | else
|
---|
258 | units={'pixel'};
|
---|
259 | end
|
---|
260 | end
|
---|
261 |
|
---|
262 | %%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
|
---|
263 | function [VarName,ValCell,Role,units]=w(DataIn)
|
---|
264 | VarName={};
|
---|
265 | ValCell={};
|
---|
266 | Role={};
|
---|
267 | units={};
|
---|
268 | if isfield(DataIn,'W')
|
---|
269 | VarName{1}='W';
|
---|
270 | ValCell{1}=DataIn.W;
|
---|
271 | Role{1}='scalar';%will remain unchanged by projection
|
---|
272 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
273 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
274 | else
|
---|
275 | units={'pixel'};
|
---|
276 | end
|
---|
277 | end
|
---|
278 |
|
---|
279 | %%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
|
---|
280 | function [VarName,ValCell,Role,units]=w_normal(DataIn)
|
---|
281 | VarName={};
|
---|
282 | ValCell={};
|
---|
283 | Role={};
|
---|
284 | units={};
|
---|
285 | if isfield(DataIn,'W')
|
---|
286 | VarName{1}='W';
|
---|
287 | ValCell{1}=DataIn.W;
|
---|
288 | Role{1}='vector_z';%will behave like a vector component by projection
|
---|
289 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
290 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
291 | else
|
---|
292 | units={'pixel'};
|
---|
293 | end
|
---|
294 | end
|
---|
295 |
|
---|
296 | %%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
|
---|
297 | function [VarName,ValCell,Role,units]=error(DataIn)
|
---|
298 | VarName={};
|
---|
299 | ValCell={};
|
---|
300 | Role={};
|
---|
301 | units={};
|
---|
302 | if isfield(DataIn,'E')
|
---|
303 | VarName{1}='E';
|
---|
304 | ValCell{1}=DataIn.E;
|
---|
305 | Role{1}='ancillary'; %TODO CHECK units in actual fields
|
---|
306 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
307 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
308 | else
|
---|
309 | units={'pixel'};
|
---|
310 | end
|
---|
311 | end
|
---|
312 |
|
---|