source: trunk/src/calc_field.m @ 425

Last change on this file since 425 was 412, checked in by sommeria, 12 years ago

bug corrected in file indexing in civ
bugs corrected in calc_field for new matlab PIV

File size: 16.7 KB
Line 
1
2
3%'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them
4%---------------------------------------------------------------------
5% [DataOut,errormsg]=calc_field(FieldList,DataIn,Coord_interp)
6%
7% OUTPUT:
8% Scal: matlab vector representing the scalar values (length nbvec defined by var_read)
9%      if no input, Scal=list of programmed scalar names (to put in menus)
10%      if only the field name is put as input, vec_A=type of scalar, which can be:
11%                   'discrete': related to the individual velocity vectors, not interpolated by patch
12%                   'vel': scalar calculated solely from velocity components
13%                   'der': needs spatial derivatives
14%                   'var': the scalar name directly corresponds to a field name in the netcdf files
15% error: error flag
16%      error = 0; OK
17%      error = 1; the prescribed scalar cannot be read or calculated from available fields
18%
19% INPUT:
20% FieldList: cell array of strings representing the name(s) of the field(s) to calculate
21% DataIn: structure representing the field, as defined in check_field_srtructure.m
22% Coord_interp(:,nb_coord) optional set of coordinates to interpolate the field (use with thin plate shell)
23%
24% FUNCTION related
25% varname_generator.m: determines the field names to read in the netcdf
26% file, depending on the scalar
27
28function [DataOut,errormsg]=calc_field(FieldList,DataIn,Coord_interp)
29
30%list of defined scalars to display in menus (in addition to 'ima_cor').
31% a type is associated to each scalar:
32%              'discrete': related to the individual velocity vectors, not interpolated by patch
33%              'vel': calculated from velocity components, continuous field (interpolated with velocity)
34%              'der': needs spatial derivatives
35%              'var': the scalar name corresponds to a field name in the netcdf files
36% a specific variable name for civ1 and civ2 fields are also associated, if
37% the scalar is calculated from other fields, as explicited below
38
39%% list of field options implemented
40FieldOptions={'velocity';...%image correlation corresponding to a vel vector
41    'ima_cor';...%image correlation corresponding to a vel vector
42    'norm_vel';...%norm of the velocity
43    'vort';...%vorticity
44    'div';...%divergence
45    'strain';...%rate of strain
46    'u';... %u velocity component
47    'v';... %v velocity component
48    'w';... %w velocity component
49    'w_normal';... %w velocity component normal to the plane
50    'error'}; %error associated to a vector (for stereo or patch)
51errormsg=[]; %default error message
52if ~exist('FieldList','var')
53    DataOut=FieldOptions;% gives the list of possible field inputs in the absence of input
54    return
55end
56
57%% check input
58if ~exist('DataIn','var')
59    DataIn=[];
60end
61if ischar(FieldList)
62    FieldList={FieldList};%convert a string input to a cell with one string element
63end
64check_grid=0;
65check_der=0;
66check_calc=ones(size(FieldList));
67for ilist=1:length(FieldList)
68    switch FieldList{ilist}
69        case {'u','v'}
70            check_grid=1;% needs a regular grid
71        case{'vort','div','strain'}% needs spatial derivatives spatial derivatives
72            check_der=1;
73        case {'velocity','norm_vel','ima_cor'};
74        otherwise
75            check_calc(ilist)=0;
76    end
77end
78FieldList=FieldList(check_calc==1);
79% if isempty(FieldList)
80%     DataOut=DataIn;
81%     return
82% end
83if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
84    nbcoord=3;
85else
86    nbcoord=2;
87end
88ListVarName={};
89ValueList={};
90RoleList={};
91units_cell={};
92
93%% interpolation with new civ data
94if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var') || check_grid ||check_der)
95    %reproduce global attributes
96    DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute;
97    for ilist=1:numel(DataOut.ListGlobalAttribute)
98        DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
99    end
100
101    %create a default grid if needed
102    if  ~exist('Coord_interp','var')
103            XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates
104            YMax=max(max(DataIn.SubRange(2,:,:)));
105            XMin=min(min(DataIn.SubRange(1,:,:)));
106            YMin=min(min(DataIn.SubRange(2,:,:)));
107        if ~isfield(DataIn,'Mesh')
108            DataIn.Mesh=sqrt(2*(XMax-XMin)*(YMax-YMin)/numel(DataIn.Coord_tps));
109            % adjust the mesh to a value 1, 2 , 5 *10^n
110            ord=10^(floor(log10(DataIn.Mesh)));%order of magnitude
111            if DataIn.Mesh/ord>=5
112                DataIn.Mesh=5*ord;
113            elseif DataIn.Mesh/ord>=2
114                DataIn.Mesh=2*ord;
115            else
116                DataIn.Mesh=ord;
117            end
118        end
119        coord_x=XMin:DataIn.Mesh:XMax;
120        coord_y=YMin:DataIn.Mesh:YMax;
121        npx=length(coord_x);
122        npy=length(coord_y);
123        DataOut.coord_x=[XMin XMax];
124        DataOut.coord_y=[YMin YMax];
125        [XI,YI]=meshgrid(coord_x,coord_y);
126        XI=reshape(XI,[],1);
127        YI=reshape(YI,[],1);
128        Coord_interp=[XI YI];
129    end
130   
131    %initialise output
132    nb_sites=size(Coord_interp,1);
133    nb_coord=size(Coord_interp,2);
134    nbval=zeros(nb_sites,1);
135    NbSubDomain=size(DataIn.SubRange,3);
136    DataOut.ListVarName={'coord_y','coord_x'};
137    DataOut.VarDimName={{'coord_y'},{'coord_x'}};
138    DataOut.VarAttribute{1}=[];
139    DataOut.VarAttribute{2}=[];
140    for ilist=1:length(FieldList)
141        switch FieldList{ilist}
142            case 'velocity'
143                DataOut.U=zeros(nb_sites,1);
144                DataOut.V=zeros(nb_sites,1);
145            otherwise
146                DataOut.(FieldList{ilist})=zeros(nb_sites,1);
147        end
148    end
149   
150    % interpolate data in each subdomain
151    for isub=1:NbSubDomain
152        nbvec_sub=DataIn.NbSites(isub);
153        check_range=(Coord_interp >=ones(nb_sites,1)*DataIn.SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*DataIn.SubRange(:,2,isub)');
154        ind_sel=find(sum(check_range,2)==nb_coord);
155        %rho smoothing parameter
156        %                 epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites
157        %                 ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
158        nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
159        if check_grid
160            EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
161        end
162        if check_der
163            [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
164        end
165        ListVar={};
166        for ilist=1:length(FieldList)
167            var_count=numel(ListVar);
168            switch FieldList{ilist}
169                case 'velocity'
170                    ListVar=[ListVar {'U', 'V'}];
171                    VarAttribute{var_count+1}.Role='vector_x';
172                    VarAttribute{var_count+2}.Role='vector_y';
173                    DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
174                    DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
175                case 'u'
176                    ListVar=[ListVar {'u'}];
177                    VarAttribute{var_count+1}.Role='scalar';
178                    DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
179                case 'v'
180                    ListVar=[ListVar {'v'}];
181                    VarAttribute{var_count+1}.Role='scalar';
182                    DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
183                case 'norm_vel'
184                    ListVar=[ListVar {'norm_vel'}];
185                    VarAttribute{var_count+1}.Role='scalar';
186                    U=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
187                    V=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
188                    DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V);
189                case 'vort'
190                    ListVar=[ListVar {'vort'}];
191                    VarAttribute{var_count+1}.Role='scalar';
192                    DataOut.vort(ind_sel)=DataOut.vort(ind_sel)-EMDY *DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
193                case 'div'
194                    ListVar=[ListVar {'div'}];
195                    VarAttribute{var_count+1}.Role='scalar';
196                    DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDY *DataIn.V_tps(1:nbvec_sub+3,isub);
197                case 'strain'
198                    ListVar=[ListVar {'strain'}];
199                    VarAttribute{var_count+1}.Role='scalar';
200                    DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
201            end
202        end
203    end
204    DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
205    nbval(nbval==0)=1;% to avoid division by zero for averaging
206    if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed
207        DataOut.ListVarName=[DataOut.ListVarName {'FF'}];
208        DataOut.VarDimName=[DataOut.VarDimName {{'coord_y','coord_x'}}];
209        DataOut.VarAttribute{length(DataOut.ListVarName)}.Role='errorflag';
210    end
211    DataOut.ListVarName=[DataOut.ListVarName ListVar];
212    for ifield=1:numel(ListVar)
213        VarDimName{ifield}={'coord_y','coord_x'};
214        DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval;
215        DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx);
216    end
217    DataOut.FF=reshape(DataOut.FF,npy,npx);
218    DataOut.VarDimName=[DataOut.VarDimName VarDimName];     
219    DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];
220else
221
222    %% civx data
223    DataOut=DataIn;
224    for ilist=1:length(FieldList)
225        if ~isempty(FieldList{ilist})
226            [VarName,Value,Role,units]=feval(FieldList{ilist},DataIn);%calculate field with appropriate function named FieldList{ilist}
227            ListVarName=[ListVarName VarName];
228            ValueList=[ValueList Value];
229            RoleList=[RoleList Role];
230            units_cell=[units_cell units];
231        end
232    end
233    %erase previous data (except coordinates)
234    for ivar=nbcoord+1:length(DataOut.ListVarName)
235        VarName=DataOut.ListVarName{ivar};
236        DataOut=rmfield(DataOut,VarName);
237    end
238    DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
239    if isfield(DataOut,'VarDimName')
240        DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
241    else
242        errormsg='element .VarDimName missing in input data';
243        return
244    end
245    DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
246    %append new data
247    DataOut.ListVarName=[DataOut.ListVarName ListVarName];
248    for ivar=1:length(ListVarName)
249        DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
250        DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
251        DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
252        DataOut.(ListVarName{ivar})=ValueList{ivar};
253    end
254end
255
256
257
258%%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
259function [VarName,ValCell,Role,units_cell]=velocity(DataIn)
260VarName={};
261ValCell={};
262Role={};
263units_cell={};
264if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
265    units=[DataIn.CoordUnit '/' DataIn.TimeUnit];
266else
267    units='pixel';
268end
269if isfield(DataIn,'U')
270    VarName=[VarName {'U'}];
271    ValCell=[ValCell {DataIn.U}];
272    Role=[Role {'vector_x'}];
273    units_cell=[units_cell {units}];
274end
275if isfield(DataIn,'V')
276    VarName=[VarName {'V'}];
277    ValCell=[ValCell {DataIn.V}];
278    Role=[Role {'vector_y'}];
279    units_cell=[units_cell {units}];
280end
281if isfield(DataIn,'W')
282    VarName=[VarName {'W'}];
283    ValCell=[ValCell {DataIn.W}];
284    Role=[Role {'vector_z'}];
285    units_cell=[units_cell {units}];
286end
287if isfield(DataIn,'F')
288    VarName=[VarName {'F'}];
289    ValCell=[ValCell {DataIn.F}];
290    Role=[Role {'warnflag'}];
291    units_cell=[units_cell {[]}];
292end
293if isfield(DataIn,'FF')
294    VarName=[VarName,{'FF'}];
295    ValCell=[ValCell {DataIn.FF}];
296    Role=[Role {'errorflag'}];
297    units_cell=[units_cell {[]}];
298end
299
300%%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%
301function [VarName,ValCell,Role,units]=ima_cor(DataIn)
302VarName={};
303ValCell={};
304Role={};
305units={};
306if isfield(DataIn,'C')
307    VarName{1}='C';
308    ValCell{1}=DataIn.C;
309    Role={'ancillary'};
310    units={[]};
311end
312
313%%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%
314function [VarName,ValCell,Role,units]=norm_vel(DataIn)
315VarName={};
316ValCell={};
317Role={};
318units={};
319if isfield(DataIn,'U') && isfield(DataIn,'V')
320    VarName{1}='norm_vel';
321    ValCell{1}=DataIn.U.*DataIn.U+ DataIn.V.*DataIn.V;
322    if isfield(DataIn,'W') && isequal(size(DataIn.W),size(DataIn.U))
323        ValCell{1}=ValCell{1}+DataIn.W.*DataIn.W;
324    end
325    ValCell{1}=sqrt(ValCell{1});
326    Role{1}='scalar';
327    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
328        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
329    else
330        units={'pixel'};
331    end
332end
333
334
335
336%%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%
337function [VarName,ValCell,Role,units]=vort(DataIn)
338VarName={};
339ValCell={};
340Role={};
341units={};
342if isfield(DataIn,'DjUi')
343    VarName{1}='vort';
344    ValCell{1}=DataIn.DjUi(:,1,2)-DataIn.DjUi(:,2,1);  %vorticity
345    siz=size(ValCell{1});
346    ValCell{1}=reshape(ValCell{1},siz(1),1);
347    Role{1}='scalar';
348    if isfield(DataIn,'TimeUnit')
349        units={[DataIn.TimeUnit '-1']};
350    else
351        units={[]};
352    end
353end
354
355%%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%
356function [VarName,ValCell,Role,units]=div(DataIn)
357VarName={};
358ValCell={};
359Role={};
360units={};
361if isfield(DataIn,'DjUi')
362    VarName{1}='div';
363    ValCell{1}=DataIn.DjUi(:,1,1)+DataIn.DjUi(:,2,2); %DUDX+DVDY
364    siz=size(ValCell{1});
365    ValCell{1}=reshape(ValCell{1},siz(1),1);
366    Role{1}='scalar';
367    if isfield(DataIn,'TimeUnit')
368        units={[DataIn.TimeUnit '-1']};
369    else
370        units={[]};
371    end
372end
373
374%%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
375function [VarName,ValCell,Role,units]=strain(DataIn)
376VarName={};
377ValCell={};
378Role={};
379units={};
380if isfield(DataIn,'DjUi')
381    VarName{1}='strain';
382    ValCell{1}=DataIn.DjUi(:,1,2)+DataIn.DjUi(:,2,1);%DVDX+DUDY
383    siz=size(ValCell{1});
384    ValCell{1}=reshape(ValCell{1},siz(1),1);
385    Role{1}='scalar';
386    if isfield(DataIn,'TimeUnit')
387        units={[DataIn.TimeUnit '-1']};
388    else
389        units={[]};
390    end
391end
392
393%%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
394function [VarName,ValCell,Role,units]=u(DataIn)
395VarName={};
396ValCell={};
397Role={};
398units={};
399if isfield(DataIn,'U')
400    VarName{1}='U';
401    ValCell{1}=DataIn.U;
402    Role{1}='scalar';
403    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
404        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
405    else
406        units={'pixel'};
407    end
408end
409
410%%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
411function [VarName,ValCell,Role,units]=v(DataIn)
412VarName={};
413ValCell={};
414Role={};
415units={};
416if isfield(DataIn,'V')
417    VarName{1}='V';
418    ValCell{1}=DataIn.V;
419    Role{1}='scalar';
420    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
421        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
422    else
423        units={'pixel'};
424    end
425end
426
427%%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
428function [VarName,ValCell,Role,units]=w(DataIn)
429VarName={};
430ValCell={};
431Role={};
432units={};
433if isfield(DataIn,'W')
434    VarName{1}='W';
435    ValCell{1}=DataIn.W;
436    Role{1}='scalar';%will remain unchanged by projection
437    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
438        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
439    else
440        units={'pixel'};
441    end
442end
443
444%%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
445function [VarName,ValCell,Role,units]=w_normal(DataIn)
446VarName={};
447ValCell={};
448Role={};
449units={};
450if isfield(DataIn,'W')
451    VarName{1}='W';
452    ValCell{1}=DataIn.W;
453    Role{1}='vector_z';%will behave like a vector component  by projection
454    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
455        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
456    else
457        units={'pixel'};
458    end
459end
460
461%%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
462function [VarName,ValCell,Role,units]=error(DataIn)
463VarName={};
464ValCell={};
465Role={};
466units={};
467if isfield(DataIn,'E')
468    VarName{1}='E';
469    ValCell{1}=DataIn.E;
470    Role{1}='ancillary'; %TODO CHECK units in actual fields
471    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
472        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
473    else
474        units={'pixel'};
475    end
476end
477
Note: See TracBrowser for help on using the repository browser.