source: trunk/src/calc_field.m @ 142

Last change on this file since 142 was 124, checked in by gostiaux, 14 years ago

smart indent

File size: 9.6 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)
[8]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];
[124]74    end
75    %erase previous data (except coordinates)
[8]76    for ivar=nbcoord+1:length(DataOut.ListVarName)
77        VarName=DataOut.ListVarName{ivar};
78        DataOut=rmfield(DataOut,VarName);
[124]79    end
[8]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)
[124]91        DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
92        DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
[8]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}];
[124]118    Role=[Role {'vector_y'}];
[8]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';
[124]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')
[8]168        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
[124]169    else
[8]170        units={'pixel'};
[124]171    end
172end
[8]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
[124]193end
[8]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
[124]212end
[8]213
214%%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
215function [VarName,ValCell,Role,units]=strain(DataIn)
216VarName={};
217ValCell={};
218Role={};
219units={};
220if isfield(DataIn,'DjUi')
[124]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    if isfield(DataIn,'TimeUnit')
[8]226        units={[DataIn.TimeUnit '-1']};
[124]227    else
[8]228        units={[]};
[124]229    end
230end
[8]231
232%%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
233function [VarName,ValCell,Role,units]=u(DataIn)
234VarName={};
235ValCell={};
236Role={};
237units={};
238if isfield(DataIn,'U')
239    VarName{1}='U';
240    ValCell{1}=DataIn.U;
241    Role{1}='scalar';
242    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
243        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
244    else
245        units={'pixel'};
246    end
247end
248
249%%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
250function [VarName,ValCell,Role,units]=v(DataIn)
251VarName={};
252ValCell={};
253Role={};
254units={};
255if isfield(DataIn,'V')
256    VarName{1}='V';
257    ValCell{1}=DataIn.V;
258    Role{1}='scalar';
259    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
260        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
261    else
262        units={'pixel'};
263    end
264end
265
266%%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
267function [VarName,ValCell,Role,units]=w(DataIn)
268VarName={};
269ValCell={};
270Role={};
271units={};
272if isfield(DataIn,'W')
273    VarName{1}='W';
274    ValCell{1}=DataIn.W;
275    Role{1}='scalar';%will remain unchanged by projection
276    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
277        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
278    else
279        units={'pixel'};
280    end
281end
282
283%%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
284function [VarName,ValCell,Role,units]=w_normal(DataIn)
285VarName={};
286ValCell={};
287Role={};
288units={};
289if isfield(DataIn,'W')
290    VarName{1}='W';
291    ValCell{1}=DataIn.W;
[124]292    Role{1}='vector_z';%will behave like a vector component  by projection
[8]293    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
294        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
295    else
296        units={'pixel'};
297    end
298end
299
300%%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
301function [VarName,ValCell,Role,units]=error(DataIn)
302VarName={};
303ValCell={};
304Role={};
305units={};
306if isfield(DataIn,'E')
307    VarName{1}='E';
308    ValCell{1}=DataIn.E;
309    Role{1}='ancillary'; %TODO CHECK units in actual fields
310    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
311        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
312    else
313        units={'pixel'};
314    end
315end
316
Note: See TracBrowser for help on using the repository browser.