source: trunk/src/calc_field.m @ 34

Last change on this file since 34 was 19, checked in by gostiaux, 14 years ago

the private files have been moved down to the root folder

File size: 9.6 KB
Line 
1%'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them 
2%---------------------------------------------------------------------
3
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 sclar name is put as input, vec_A=type of scalar, which can be:
9%                   'vel': scalar calculated solely from velocity components
10%                   'der': needs spatial derivatives     
11%                   'var': the scalar name directly corresponds to a field name in the netcdf files
12% error: error flag
13            % error = 0; OK
14            % error = 1; the prescribed scalar cannot be read or calculated from available fields
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.