source: trunk/src/calc_field.m @ 190

Last change on this file since 190 was 180, checked in by sommeria, 13 years ago

rationalisation of uvmat, introduction of the new function read_field, links with get_field, several bug repairs

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