source: trunk/src/calc_field.m @ 387

Last change on this file since 387 was 383, checked in by sommeria, 13 years ago

bugs corrected , new Matlab patch introduced for patch2
File size: 14.9 KB
Line 
1%'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them
2%---------------------------------------------------------------------
3% [DataOut,errormsg]=calc_field(FieldName,DataIn)
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 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
10%                   'vel': scalar calculated solely from velocity components
11%                   'der': needs spatial derivatives
12%                   'var': the scalar name directly corresponds to a field name in the netcdf files
13% error: error flag
14%      error = 0; OK
15%      error = 1; the prescribed scalar cannot be read or calculated from available fields
16%
17% INPUT:
18% FieldName: string representing the name of the field
19% DataIn: structure representing the field, as defined in check_field_srtructure.m
20%
21% FUNCTION related
22% varname_generator.m: determines the field names to read in the netcdf
23% file, depending on the scalar
24
25function [DataOut,errormsg]=calc_field(FieldName,DataIn,VelType,XI,YI)
26
27%list of defined scalars to display in menus (in addition to 'ima_cor').
28% a type is associated to each scalar:
29%              'discrete': related to the individual velocity vectors, not interpolated by patch
30%              'vel': calculated from velocity components, continuous field (interpolated with velocity)
31%              'der': needs spatial derivatives
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
34% the scalar is calculated from other fields, as explicited below
35
36list_field={'velocity';...%image correlation corresponding to a vel vector
37    'ima_cor';...%image correlation corresponding to a vel vector
38    'norm_vel';...%norm of the velocity
39    'vort';...%vorticity
40    'div';...%divergence
41    'strain';...%rate of strain
42    'u';... %u velocity component
43    'v';... %v velocity component
44    'w';... %w velocity component
45    'w_normal';... %w velocity component normal to the plane
46    'error'}; %error associated to a vector (for stereo or patch)
47errormsg=[]; %default error message
48if ~exist('FieldName','var')
49    DataOut=list_field;% gives the list of possible fields in the absence of input
50else
51    if ~exist('DataIn','var')
52        DataIn=[];
53    end
54    if ischar(FieldName)
55        FieldName={FieldName};%convert a string input to a cell with one string element
56    end
57    if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
58        nbcoord=3;
59    else
60        nbcoord=2;
61    end
62    ListVarName={};
63    ValueList={};
64    RoleList={};
65    units_cell={};
66   
67    %% interpolation with new civ data
68    %if isfield(DataIn,'X_SubRange')
69    if isfield(DataIn,'SubRange') &&(strcmp(VelType,'filter1')||strcmp(VelType,'filter2'))
70        XMax=max(max(DataIn.SubRange(1,:,:)));
71        YMax=max(max(DataIn.SubRange(2,:,:)));
72        XMin=min(min(DataIn.SubRange(1,:,:)));
73        YMin=min(min(DataIn.SubRange(2,:,:)));
74        if ~(exist('XI','var')&&exist('YI','var'))
75            Mesh=sqrt((YMax-YMin)*(XMax-XMin)/numel(DataIn.Coord_tps(:,1,:)));%2D
76            xI=XMin:Mesh:XMax;
77            yI=YMin:Mesh:YMax;
78            [XI,YI]=meshgrid(xI,yI);
79            XI=reshape(XI,[],1);
80            YI=reshape(YI,[],1);
81        end
82        DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
83        for ilist=1:numel(DataOut.ListGlobalAttribute)
84            eval(['DataOut.' DataOut.ListGlobalAttribute{ilist} '=DataIn.' DataIn.ListGlobalAttribute{ilist} ';'])
85        end
86        DataOut.ListVarName={'coord_y','coord_x','FF'};
87        DataOut.VarDimName{1}='coord_y';
88        DataOut.VarDimName{2}='coord_x';
89        DataOut.coord_y=[yI(1) yI(end)];
90        DataOut.coord_x=[xI(1) xI(end)];
91        DataOut.U=zeros(size(XI));
92        DataOut.V=zeros(size(XI));
93        DataOut.vort=zeros(size(XI));       
94        DataOut.div=zeros(size(XI));
95        nbval=zeros(size(XI));
96        NbSubDomain=size(DataIn.SubRange,3);
97        for isub=1:NbSubDomain
98            if isfield(DataIn,'NbSites')
99                nbvec_sub=DataIn.NbSites(isub);
100            else
101                nbvec_sub=numel(find(DataIn.Indices_tps(:,isub)));
102            end
103                ind_sel=find(XI>=DataIn.SubRange(1,1,isub) & XI<=DataIn.SubRange(1,2,isub) & YI>=DataIn.SubRange(2,1,isub) & YI<=DataIn.SubRange(2,2,isub));
104                %rho smoothing parameter
105                epoints = [XI(ind_sel) YI(ind_sel)];% coordinates of interpolation sites
106                ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
107                nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
108%                 PM = [ones(size(epoints,1),1) epoints];
109                switch FieldName{1}
110                    case {'velocity','u','v'}
111                       % DM_eval = DistanceMatrix(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
112                       % EM = tps(1,DM_eval);%values of thin plate
113                        EM = tps_eval(epoints,ctrs);
114                    case{'vort','div'}
115                        [EMDX,EMDY] = tps_eval_dxy(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
116     
117                end
118                switch FieldName{1}
119                    case 'velocity'
120                        ListFields={'U', 'V'};
121                        VarAttributes{1}.Role='vector_x';
122                        VarAttributes{2}.Role='vector_y';
123                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
124                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
125                    case 'u'
126                        ListFields={'U'};
127                        VarAttributes{1}.Role='scalar';
128                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
129                    case 'v'
130                        ListFields={'V'};
131                        VarAttributes{1}.Role='scalar';
132                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
133                    case 'vort'
134                        ListFields={'vort'};
135                        VarAttributes{1}.Role='scalar';
136                        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);
137                    case 'div'
138                        ListFields={'div'};
139                        VarAttributes{1}.Role='scalar';
140                        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);
141                end
142        end
143        DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang 
144       
145        DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI));
146        nbval(nbval==0)=1;
147        DataOut.U=DataOut.U ./nbval;
148        DataOut.V=DataOut.V ./nbval;
149        DataOut.U=reshape(DataOut.U,numel(yI),numel(xI));
150        DataOut.V=reshape(DataOut.V,numel(yI),numel(xI));
151        DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI));
152        DataOut.div=reshape(DataOut.div,numel(yI),numel(xI));
153        DataOut.ListVarName=[DataOut.ListVarName ListFields];
154        for ilist=3:numel(DataOut.ListVarName)
155            DataOut.VarDimName{ilist}={'coord_y','coord_x'};
156        end
157        DataOut.VarAttribute={[],[]};
158        DataOut.VarAttribute{3}.Role='errorflag';
159        DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];   
160    else   
161       
162    %% civx data   
163        DataOut=DataIn;
164        for ilist=1:length(FieldName)
165            if ~isempty(FieldName{ilist})
166                [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}
167                ListVarName=[ListVarName VarName];
168                ValueList=[ValueList Value];
169                RoleList=[RoleList Role];
170                units_cell=[units_cell units];
171            end
172        end
173        %erase previous data (except coordinates)
174        for ivar=nbcoord+1:length(DataOut.ListVarName)
175            VarName=DataOut.ListVarName{ivar};
176            DataOut=rmfield(DataOut,VarName);
177        end
178        DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
179        if isfield(DataOut,'VarDimName')
180            DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
181        else
182            errormsg='element .VarDimName missing in input data';
183            return
184        end
185        DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
186        %append new data
187        DataOut.ListVarName=[DataOut.ListVarName ListVarName];
188        for ivar=1:length(ListVarName)
189            DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
190            DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
191            DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
192            eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])
193        end
194    end
195end
196
197
198%%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
199function [VarName,ValCell,Role,units_cell]=velocity(DataIn)
200VarName={};
201ValCell={};
202Role={};
203units_cell={};
204if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
205    units=[DataIn.CoordUnit '/' DataIn.TimeUnit];
206else
207    units='pixel';
208end
209if isfield(DataIn,'U')
210    VarName=[VarName {'U'}];
211    ValCell=[ValCell {DataIn.U}];
212    Role=[Role {'vector_x'}];
213    units_cell=[units_cell {units}];
214end
215if isfield(DataIn,'V')
216    VarName=[VarName {'V'}];
217    ValCell=[ValCell {DataIn.V}];
218    Role=[Role {'vector_y'}];
219    units_cell=[units_cell {units}];
220end
221if isfield(DataIn,'W')
222    VarName=[VarName {'W'}];
223    ValCell=[ValCell {DataIn.W}];
224    Role=[Role {'vector_z'}];
225    units_cell=[units_cell {units}];
226end
227if isfield(DataIn,'F')
228    VarName=[VarName {'F'}];
229    ValCell=[ValCell {DataIn.F}];
230    Role=[Role {'warnflag'}];
231    units_cell=[units_cell {[]}];
232end
233if isfield(DataIn,'FF')
234    VarName=[VarName,{'FF'}];
235    ValCell=[ValCell {DataIn.FF}];
236    Role=[Role {'errorflag'}];
237    units_cell=[units_cell {[]}];
238end
239
240%%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%
241function [VarName,ValCell,Role,units]=ima_cor(DataIn)
242VarName={};
243ValCell={};
244Role={};
245units={};
246if isfield(DataIn,'C')
247    VarName{1}='C';
248    ValCell{1}=DataIn.C;
249    Role={'ancillary'};
250    units={[]};
251end
252
253%%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%
254function [VarName,ValCell,Role,units]=norm_vel(DataIn)
255VarName={};
256ValCell={};
257Role={};
258units={};
259if isfield(DataIn,'U') && isfield(DataIn,'V')
260    VarName{1}='norm_vel';
261    ValCell{1}=DataIn.U.*DataIn.U+ DataIn.V.*DataIn.V;
262    if isfield(DataIn,'W') && isequal(size(DataIn.W),size(DataIn.U))
263        ValCell{1}=ValCell{1}+DataIn.W.*DataIn.W;
264    end
265    ValCell{1}=sqrt(ValCell{1});
266    Role{1}='scalar';
267    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
268        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
269    else
270        units={'pixel'};
271    end
272end
273
274
275
276%%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%
277function [VarName,ValCell,Role,units]=vort(DataIn)
278VarName={};
279ValCell={};
280Role={};
281units={};
282if isfield(DataIn,'DjUi')
283    VarName{1}='vort';
284    ValCell{1}=DataIn.DjUi(:,1,2)-DataIn.DjUi(:,2,1);  %vorticity
285    siz=size(ValCell{1});
286    ValCell{1}=reshape(ValCell{1},siz(1),1);
287    Role{1}='scalar';
288    if isfield(DataIn,'TimeUnit')
289        units={[DataIn.TimeUnit '-1']};
290    else
291        units={[]};
292    end
293end
294
295%%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%
296function [VarName,ValCell,Role,units]=div(DataIn)
297VarName={};
298ValCell={};
299Role={};
300units={};
301if isfield(DataIn,'DjUi')
302    VarName{1}='div';
303    ValCell{1}=DataIn.DjUi(:,1,1)+DataIn.DjUi(:,2,2); %DUDX+DVDY
304    siz=size(ValCell{1});
305    ValCell{1}=reshape(ValCell{1},siz(1),1);
306    Role{1}='scalar';
307    if isfield(DataIn,'TimeUnit')
308        units={[DataIn.TimeUnit '-1']};
309    else
310        units={[]};
311    end
312end
313
314%%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
315function [VarName,ValCell,Role,units]=strain(DataIn)
316VarName={};
317ValCell={};
318Role={};
319units={};
320if isfield(DataIn,'DjUi')
321    VarName{1}='strain';
322    ValCell{1}=DataIn.DjUi(:,1,2)+DataIn.DjUi(:,2,1);%DVDX+DUDY
323    siz=size(ValCell{1});
324    ValCell{1}=reshape(ValCell{1},siz(1),1);
325    Role{1}='scalar';
326    if isfield(DataIn,'TimeUnit')
327        units={[DataIn.TimeUnit '-1']};
328    else
329        units={[]};
330    end
331end
332
333%%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
334function [VarName,ValCell,Role,units]=u(DataIn)
335VarName={};
336ValCell={};
337Role={};
338units={};
339if isfield(DataIn,'U')
340    VarName{1}='U';
341    ValCell{1}=DataIn.U;
342    Role{1}='scalar';
343    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
344        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
345    else
346        units={'pixel'};
347    end
348end
349
350%%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
351function [VarName,ValCell,Role,units]=v(DataIn)
352VarName={};
353ValCell={};
354Role={};
355units={};
356if isfield(DataIn,'V')
357    VarName{1}='V';
358    ValCell{1}=DataIn.V;
359    Role{1}='scalar';
360    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
361        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
362    else
363        units={'pixel'};
364    end
365end
366
367%%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
368function [VarName,ValCell,Role,units]=w(DataIn)
369VarName={};
370ValCell={};
371Role={};
372units={};
373if isfield(DataIn,'W')
374    VarName{1}='W';
375    ValCell{1}=DataIn.W;
376    Role{1}='scalar';%will remain unchanged by projection
377    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
378        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
379    else
380        units={'pixel'};
381    end
382end
383
384%%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
385function [VarName,ValCell,Role,units]=w_normal(DataIn)
386VarName={};
387ValCell={};
388Role={};
389units={};
390if isfield(DataIn,'W')
391    VarName{1}='W';
392    ValCell{1}=DataIn.W;
393    Role{1}='vector_z';%will behave like a vector component  by projection
394    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
395        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
396    else
397        units={'pixel'};
398    end
399end
400
401%%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
402function [VarName,ValCell,Role,units]=error(DataIn)
403VarName={};
404ValCell={};
405Role={};
406units={};
407if isfield(DataIn,'E')
408    VarName{1}='E';
409    ValCell{1}=DataIn.E;
410    Role{1}='ancillary'; %TODO CHECK units in actual fields
411    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
412        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
413    else
414        units={'pixel'};
415    end
416end
417
Note: See TracBrowser for help on using the repository browser.