source: trunk/src/calc_field.m @ 257

Last change on this file since 257 was 247, checked in by sommeria, 14 years ago

thin plate shell corrected for spatial derivatives
File size: 15.6 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)
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    % new civ data
67    if isfield(DataIn,'X_SubRange')
68        XMax=max(max(DataIn.X_SubRange));
69        YMax=max(max(DataIn.Y_SubRange));
70        XMin=min(min(DataIn.X_SubRange));
71        YMin=min(min(DataIn.Y_SubRange));
72        Mesh=sqrt((YMax-YMin)*(XMax-XMin)/numel(DataIn.X_tps));%2D
73        xI=XMin:Mesh:XMax;
74        yI=YMin:Mesh:YMax;
75        [XI,YI]=meshgrid(xI,yI);
76        XI=reshape(XI,[],1);
77        YI=reshape(YI,[],1);
78        DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
79        for ilist=1:numel(DataOut.ListGlobalAttribute)
80            eval(['DataOut.' DataOut.ListGlobalAttribute{ilist} '=DataIn.' DataIn.ListGlobalAttribute{ilist} ';'])
81        end
82        DataOut.ListVarName={'coord_y','coord_x','FF'};
83        DataOut.VarDimName{1}='coord_y';
84        DataOut.VarDimName{2}='coord_x';
85        DataOut.coord_y=[yI(1) yI(end)];
86        DataOut.coord_x=[xI(1) xI(end)];
87        DataOut.U=zeros(size(XI));
88        DataOut.V=zeros(size(XI));
89        DataOut.vort=zeros(size(XI));
90        DataOut.div=zeros(size(XI));
91        nbval=zeros(size(XI));
92        [NbSubDomain,xx]=size(DataIn.X_SubRange);
93        for isub=1:NbSubDomain
94            nbvec_sub=DataIn.NbSites(isub);         
95                ind_sel=find(XI>=DataIn.X_SubRange(isub,1) & XI<=DataIn.X_SubRange(isub,2) & YI>=DataIn.Y_SubRange(isub,1) & YI<=DataIn.Y_SubRange(isub,2));
96                %rho smoothing parameter
97                epoints = [XI(ind_sel) YI(ind_sel)];% coordinates of interpolation sites
98                ctrs=[DataIn.X_tps(1:nbvec_sub,isub) DataIn.Y_tps(1:nbvec_sub,isub)];%(=initial points) ctrs
99                nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
100%                 PM = [ones(size(epoints,1),1) epoints];
101                switch FieldName{1}
102                    case {'velocity','u','v'}
103                       % DM_eval = DistanceMatrix(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
104                       % EM = tps(1,DM_eval);%values of thin plate
105                        EM = tps_eval(epoints,ctrs);
106                    case{'vort','div'}
107                        [EMDX,EMDY] = tps_eval_dxy(epoints,ctrs);%2D matrix of distances between extrapolation points epoints and spline centres (=site points) ctrs
108     
109                end
110                switch FieldName{1}
111                    case 'velocity'
112                        ListFields={'U', 'V'};
113                        VarAttributes{1}.Role='vector_x';
114                        VarAttributes{2}.Role='vector_y';
115                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
116                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
117                    case 'u'
118                        ListFields={'U'};
119                        VarAttributes{1}.Role='scalar';
120                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
121                    case 'v'
122                        ListFields={'V'};
123                        VarAttributes{1}.Role='scalar';
124                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
125                    case 'vort'
126                        ListFields={'vort'};
127                        VarAttributes{1}.Role='scalar';
128%                         EMX = [DMXY(:,:,1) PM];
129%                         EMY = [DMXY(:,:,2) PM];
130%                         DataIn.U_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
131%                         DataIn.V_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
132%                         DataIn.V_tps(nbvec_sub+2,isub)=0;% X coefficient suppressed for x wise derivatives
133%                         DataIn.U_tps(nbvec_sub+3,isub)=0;% Y coefficient suppressed for x wise derivatives
134                        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);
135                    case 'div'
136                        ListFields={'div'};
137                        VarAttributes{1}.Role='scalar';
138%                         EMX = [DMXY(:,:,1) PM];
139%                         EMY = [DMXY(:,:,2) PM];
140%                         DataIn.U_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
141%                         DataIn.V_tps(nbvec_sub+1,isub)=0;%constant value suppressed by spatial derivatives
142%                         DataIn.V_tps(nbvec_sub+2,isub)=0;% X coefficient suppressed for x wise derivatives
143%                         DataIn.U_tps(nbvec_sub+3,isub)=0;% Y coefficient suppressed for x wise derivatives
144                        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);
145                end
146        end
147        DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang 
148       
149        DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI));
150        nbval(nbval==0)=1;
151        DataOut.U=DataOut.U ./nbval;
152        DataOut.V=DataOut.V ./nbval;
153        DataOut.U=reshape(DataOut.U,numel(yI),numel(xI));
154        DataOut.V=reshape(DataOut.V,numel(yI),numel(xI));
155        DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI));
156        DataOut.div=reshape(DataOut.div,numel(yI),numel(xI));
157        DataOut.ListVarName=[DataOut.ListVarName ListFields];
158        for ilist=3:numel(DataOut.ListVarName)
159            DataOut.VarDimName{ilist}={'coord_y','coord_x'};
160        end
161        DataOut.VarAttribute={[],[]};
162        DataOut.VarAttribute{3}.Role='errorflag';
163        DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];     
164    else   
165        DataOut=DataIn;
166        for ilist=1:length(FieldName)
167            if ~isempty(FieldName{ilist})
168                [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}
169                ListVarName=[ListVarName VarName];
170                ValueList=[ValueList Value];
171                RoleList=[RoleList Role];
172                units_cell=[units_cell units];
173            end
174        end
175        %erase previous data (except coordinates)
176        for ivar=nbcoord+1:length(DataOut.ListVarName)
177            VarName=DataOut.ListVarName{ivar};
178            DataOut=rmfield(DataOut,VarName);
179        end
180        DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
181        if isfield(DataOut,'VarDimName')
182            DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
183        else
184            errormsg='element .VarDimName missing in input data';
185            return
186        end
187        DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
188        %append new data
189        DataOut.ListVarName=[DataOut.ListVarName ListVarName];
190        for ivar=1:length(ListVarName)
191            DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
192            DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
193            DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
194            eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])
195        end
196    end
197end
198
199
200%%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
201function [VarName,ValCell,Role,units_cell]=velocity(DataIn)
202VarName={};
203ValCell={};
204Role={};
205units_cell={};
206if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
207    units=[DataIn.CoordUnit '/' DataIn.TimeUnit];
208else
209    units='pixel';
210end
211if isfield(DataIn,'U')
212    VarName=[VarName {'U'}];
213    ValCell=[ValCell {DataIn.U}];
214    Role=[Role {'vector_x'}];
215    units_cell=[units_cell {units}];
216end
217if isfield(DataIn,'V')
218    VarName=[VarName {'V'}];
219    ValCell=[ValCell {DataIn.V}];
220    Role=[Role {'vector_y'}];
221    units_cell=[units_cell {units}];
222end
223if isfield(DataIn,'W')
224    VarName=[VarName {'W'}];
225    ValCell=[ValCell {DataIn.W}];
226    Role=[Role {'vector_z'}];
227    units_cell=[units_cell {units}];
228end
229if isfield(DataIn,'F')
230    VarName=[VarName {'F'}];
231    ValCell=[ValCell {DataIn.F}];
232    Role=[Role {'warnflag'}];
233    units_cell=[units_cell {[]}];
234end
235if isfield(DataIn,'FF')
236    VarName=[VarName,{'FF'}];
237    ValCell=[ValCell {DataIn.FF}];
238    Role=[Role {'errorflag'}];
239    units_cell=[units_cell {[]}];
240end
241
242%%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%
243function [VarName,ValCell,Role,units]=ima_cor(DataIn)
244VarName={};
245ValCell={};
246Role={};
247units={};
248if isfield(DataIn,'C')
249    VarName{1}='C';
250    ValCell{1}=DataIn.C;
251    Role={'ancillary'};
252    units={[]};
253end
254
255%%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%
256function [VarName,ValCell,Role,units]=norm_vel(DataIn)
257VarName={};
258ValCell={};
259Role={};
260units={};
261if isfield(DataIn,'U') && isfield(DataIn,'V')
262    VarName{1}='norm_vel';
263    ValCell{1}=DataIn.U.*DataIn.U+ DataIn.V.*DataIn.V;
264    if isfield(DataIn,'W') && isequal(size(DataIn.W),size(DataIn.U))
265        ValCell{1}=ValCell{1}+DataIn.W.*DataIn.W;
266    end
267    ValCell{1}=sqrt(ValCell{1});
268    Role{1}='scalar';
269    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
270        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
271    else
272        units={'pixel'};
273    end
274end
275
276
277
278%%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%
279function [VarName,ValCell,Role,units]=vort(DataIn)
280VarName={};
281ValCell={};
282Role={};
283units={};
284if isfield(DataIn,'DjUi')
285    VarName{1}='vort';
286    ValCell{1}=DataIn.DjUi(:,1,2)-DataIn.DjUi(:,2,1);  %vorticity
287    siz=size(ValCell{1});
288    ValCell{1}=reshape(ValCell{1},siz(1),1);
289    Role{1}='scalar';
290    if isfield(DataIn,'TimeUnit')
291        units={[DataIn.TimeUnit '-1']};
292    else
293        units={[]};
294    end
295end
296
297%%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%
298function [VarName,ValCell,Role,units]=div(DataIn)
299VarName={};
300ValCell={};
301Role={};
302units={};
303if isfield(DataIn,'DjUi')
304    VarName{1}='div';
305    ValCell{1}=DataIn.DjUi(:,1,1)+DataIn.DjUi(:,2,2); %DUDX+DVDY
306    siz=size(ValCell{1});
307    ValCell{1}=reshape(ValCell{1},siz(1),1);
308    Role{1}='scalar';
309    if isfield(DataIn,'TimeUnit')
310        units={[DataIn.TimeUnit '-1']};
311    else
312        units={[]};
313    end
314end
315
316%%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
317function [VarName,ValCell,Role,units]=strain(DataIn)
318VarName={};
319ValCell={};
320Role={};
321units={};
322if isfield(DataIn,'DjUi')
323    VarName{1}='strain';
324    ValCell{1}=DataIn.DjUi(:,1,2)+DataIn.DjUi(:,2,1);%DVDX+DUDY
325    siz=size(ValCell{1});
326    ValCell{1}=reshape(ValCell{1},siz(1),1);
327    Role{1}='scalar';
328    if isfield(DataIn,'TimeUnit')
329        units={[DataIn.TimeUnit '-1']};
330    else
331        units={[]};
332    end
333end
334
335%%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
336function [VarName,ValCell,Role,units]=u(DataIn)
337VarName={};
338ValCell={};
339Role={};
340units={};
341if isfield(DataIn,'U')
342    VarName{1}='U';
343    ValCell{1}=DataIn.U;
344    Role{1}='scalar';
345    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
346        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
347    else
348        units={'pixel'};
349    end
350end
351
352%%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
353function [VarName,ValCell,Role,units]=v(DataIn)
354VarName={};
355ValCell={};
356Role={};
357units={};
358if isfield(DataIn,'V')
359    VarName{1}='V';
360    ValCell{1}=DataIn.V;
361    Role{1}='scalar';
362    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
363        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
364    else
365        units={'pixel'};
366    end
367end
368
369%%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
370function [VarName,ValCell,Role,units]=w(DataIn)
371VarName={};
372ValCell={};
373Role={};
374units={};
375if isfield(DataIn,'W')
376    VarName{1}='W';
377    ValCell{1}=DataIn.W;
378    Role{1}='scalar';%will remain unchanged by projection
379    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
380        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
381    else
382        units={'pixel'};
383    end
384end
385
386%%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
387function [VarName,ValCell,Role,units]=w_normal(DataIn)
388VarName={};
389ValCell={};
390Role={};
391units={};
392if isfield(DataIn,'W')
393    VarName{1}='W';
394    ValCell{1}=DataIn.W;
395    Role{1}='vector_z';%will behave like a vector component  by projection
396    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
397        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
398    else
399        units={'pixel'};
400    end
401end
402
403%%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
404function [VarName,ValCell,Role,units]=error(DataIn)
405VarName={};
406ValCell={};
407Role={};
408units={};
409if isfield(DataIn,'E')
410    VarName{1}='E';
411    ValCell{1}=DataIn.E;
412    Role{1}='ancillary'; %TODO CHECK units in actual fields
413    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
414        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
415    else
416        units={'pixel'};
417    end
418end
419
Note: See TracBrowser for help on using the repository browser.