source: trunk/src/calc_field.m @ 89

Last change on this file since 89 was 89, checked in by sommeria, 14 years ago

many bug corrections and cleaning. Activation of the BW option in uvmat. Improvement of the interaction of get_field with uvmat.

File size: 9.6 KB
Line 
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
22function [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}
33list_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)
44errormsg=[]; %default error message
45if ~exist('FieldName','var')
46    DataOut=list_field;% gives the list of possible fields in the absence of input
47else
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
92end
93
94%%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
95function [VarName,ValCell,Role,units_cell]=velocity(DataIn)
96VarName={};
97ValCell={};
98Role={};
99units_cell={};
100if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
101    units=[DataIn.CoordUnit '/' DataIn.TimeUnit];
102else
103    units='pixel';
104end
105if isfield(DataIn,'U')
106    VarName=[VarName {'U'}];
107    ValCell=[ValCell {DataIn.U}];
108    Role=[Role {'vector_x'}];
109    units_cell=[units_cell {units}];
110end
111if isfield(DataIn,'V')
112    VarName=[VarName {'V'}];
113    ValCell=[ValCell {DataIn.V}];
114    Role=[Role {'vector_y'}]; 
115    units_cell=[units_cell {units}];
116end
117if isfield(DataIn,'W')
118    VarName=[VarName {'W'}];
119    ValCell=[ValCell {DataIn.W}];
120    Role=[Role {'vector_z'}];
121    units_cell=[units_cell {units}];
122end
123if isfield(DataIn,'F')
124    VarName=[VarName {'F'}];
125    ValCell=[ValCell {DataIn.F}];
126    Role=[Role {'warnflag'}];
127    units_cell=[units_cell {[]}];
128end
129if isfield(DataIn,'FF')
130    VarName=[VarName,{'FF'}];
131    ValCell=[ValCell {DataIn.FF}];
132    Role=[Role {'errorflag'}];
133    units_cell=[units_cell {[]}];
134end
135
136%%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%
137function [VarName,ValCell,Role,units]=ima_cor(DataIn)
138VarName={};
139ValCell={};
140Role={};
141units={};
142if isfield(DataIn,'C')
143    VarName{1}='C';
144    ValCell{1}=DataIn.C;
145    Role={'ancillary'};
146    units={[]};
147end
148
149%%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%
150function [VarName,ValCell,Role,units]=norm_vel(DataIn)
151VarName={};
152ValCell={};
153Role={};
154units={};
155if 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
168end 
169
170
171
172%%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%
173function [VarName,ValCell,Role,units]=vort(DataIn)
174VarName={};
175ValCell={};
176Role={};
177units={};
178if 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
189end 
190
191%%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%
192function [VarName,ValCell,Role,units]=div(DataIn)
193VarName={};
194ValCell={};
195Role={};
196units={};
197if 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
208end 
209
210%%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
211function [VarName,ValCell,Role,units]=strain(DataIn)
212VarName={};
213ValCell={};
214Role={};
215units={};
216if 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
226end 
227
228%%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
229function [VarName,ValCell,Role,units]=u(DataIn)
230VarName={};
231ValCell={};
232Role={};
233units={};
234if 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
243end
244
245%%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
246function [VarName,ValCell,Role,units]=v(DataIn)
247VarName={};
248ValCell={};
249Role={};
250units={};
251if 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
260end
261
262%%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
263function [VarName,ValCell,Role,units]=w(DataIn)
264VarName={};
265ValCell={};
266Role={};
267units={};
268if 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
277end
278
279%%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
280function [VarName,ValCell,Role,units]=w_normal(DataIn)
281VarName={};
282ValCell={};
283Role={};
284units={};
285if 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
294end
295
296%%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
297function [VarName,ValCell,Role,units]=error(DataIn)
298VarName={};
299ValCell={};
300Role={};
301units={};
302if 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
311end
312
Note: See TracBrowser for help on using the repository browser.