source: trunk/src/calc_field.m @ 162

Last change on this file since 162 was 156, checked in by sommeria, 13 years ago

many bug repairs and corrections for mouse action
create_grid: option for black marjkers for grid detection

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