source: trunk/src/sub_field.m @ 511

Last change on this file since 511 was 494, checked in by sommeria, 12 years ago

various bugs corrected after testing in Windows OS. Introduction
of filter tps

File size: 15.7 KB
Line 
1%'sub_field': combines two input fields
2%
3% the two fields are subtstracted when of the same nature (scalar or
4% vector), if the coordinates do not coincide, the second field is
5% interpolated on the cooridintes of the first one
6%
7% when scalar and vectors are combined, the fields are just merged in a single matlab structure for common visualisation
8%-----------------------------------------------------------------------
9% function SubData=sub_field(Field,Field_1)
10%
11% OUPUT:
12% SubData: structure representing the resulting field
13%
14% INPUT:
15% Field: matlab structure representing the first field
16% Field_1:matlab structure representing the second field
17
18function [SubData,errormsg]=sub_field(Field,Field_1)
19
20%% global attributes
21test_attr=0;
22if isfield(Field,'ListGlobalAttribute')
23    SubData.ListGlobalAttribute=Field.ListGlobalAttribute;
24    for ilist=1:numel(Field.ListGlobalAttribute)
25        AttrName=Field.ListGlobalAttribute{ilist};
26        SubData.(AttrName)=Field.(AttrName);
27    end
28    test_attr=1;
29end
30if isfield(Field_1,'ListGlobalAttribute')
31    for ilist=1:numel(Field_1.ListGlobalAttribute)
32        test_1=1;
33        AttrName=Field_1.ListGlobalAttribute{ilist};
34        if test_attr
35            for i_prev=1:numel(Field.ListGlobalAttribute)
36                if isequal(Field.ListGlobalAttribute{i_prev},AttrName)
37                    test_1=0; %attribute already written
38                    eval(['Val=Field.' AttrName ';'])                 
39                    eval(['Val_1=Field_1.' AttrName ';'])
40                    if isequal(Val,Val_1)           
41                        break% data already written
42                    else
43                        eval(['SubData.' AttrName '_1=Field_1.' AttrName ';'])
44                    end
45                end
46            end
47        end
48        if test_1
49            eval(['SubData.' AttrName '=Field_1.' AttrName ';'])
50        end
51    end
52end
53
54%% variables
55SubData.ListVarName=Field.ListVarName;
56SubData.VarDimName=Field.VarDimName;
57if isfield(Field,'VarAttribute')
58    SubData.VarAttribute=Field.VarAttribute;
59end
60%reproduce the first field Field by default
61for ivar=1:numel(Field.ListVarName)
62   VarName=Field.ListVarName{ivar};
63   SubData.(VarName)=Field.(VarName);
64end
65
66%% check the two input fields     
67[CellVarIndex,NbDim,VarTypeCell,errormsg]=find_field_indices(Field);
68if ~isempty(errormsg)
69    errormsg=['invalid  first input to sub_field:' errormsg];
70    return
71end
72[CellVarIndex_1,NbDim_1,VarTypeCell_1,errormsg]=find_field_indices(Field_1);
73if ~isempty(errormsg)
74    errormsg=['invalid second input to sub_field:' errormsg];
75    return
76end
77iselect=find(NbDim==2);
78for icell=iselect
79    if ~isempty(VarTypeCell{icell}.coord_tps)
80        NbDim(icell)=0;
81    end
82end
83iselect=find(NbDim==2);
84if ~isequal(numel(iselect),1)
85    errormsg='invalid  first input to sub_field: it must  contain a single 2D field cell';
86    return
87end
88iselect_1=find(NbDim_1==2);
89for icell=iselect_1
90    if ~isempty(VarTypeCell_1{icell}.coord_tps)
91        NbDim_1(icell)=0;
92    end
93end
94iselect_1=find(NbDim_1==2);
95if ~isequal(numel(iselect_1),1)
96    errormsg='invalid  second input to sub_field: it must  contain a single 2D field cell';
97    return
98end
99VarType=VarTypeCell{iselect};
100VarType_1=VarTypeCell_1{iselect_1};
101testX=~isempty(VarType.coord_x)&& ~isempty(VarType.coord_y);%unstructured coordiantes
102testX_1=~isempty(VarType_1.coord_x)&& ~isempty(VarType_1.coord_y);%unstructured coordiantes
103testU=~isempty(VarType.vector_x)&& ~isempty(VarType.vector_y);%vector field
104testU_1=~isempty(VarType_1.vector_x)&& ~isempty(VarType_1.vector_y);%vector field
105testfalse_1=~isempty(VarType_1.errorflag);
106ivar_C=[VarType.scalar VarType.image VarType.color VarType.ancillary]; %defines index (indices) for the scalar or ancillary fields
107if numel(ivar_C)>1
108    errormsg='too many scalar fields in the first input of sub_field.m';
109    return
110end
111ivar_C_1=[VarType_1.scalar VarType_1.image VarType_1.color VarType_1.ancillary]; %defines index (indices) for the scalar or ancillary fields
112if numel(ivar_C_1)>1
113    errormsg='too many scalar fields in the second input of sub_field.m';
114    return
115end
116
117%% substract two vector fields or two scalars
118if (testU && testU_1) || (~testU && ~testU_1)
119   %check coincidence in positions
120   %unstructured coordinates for the first field
121   if testX 
122       XName=Field.ListVarName{VarType.coord_x};
123       YName=Field.ListVarName{VarType.coord_y};
124       vec_X=Field.(XName);
125       vec_Y=Field.(YName);
126       nbpoints=numel(vec_X);
127       vec_X=reshape(vec_X,nbpoints,1);
128       vec_Y=reshape(vec_Y,nbpoints,1);
129       if testX_1 %unstructured coordinates for the second field
130            X_1_Name=Field_1.ListVarName{VarType_1.coord_x};
131            Y_1_Name=Field_1.ListVarName{VarType_1.coord_y};
132            eval(['vec_X_1=Field_1.' X_1_Name ';'])
133            eval(['vec_Y_1=Field_1.' Y_1_Name ';'])
134
135       else   %structured coordinates for the second field
136           y_1_Name=Field_1.ListVarName{VarType_1.coord(1)};
137           x_1_Name=Field_1.ListVarName{VarType_1.coord(2)};
138           eval(['y_1=Field_1.' y_1_Name ';'])
139           eval(['x_1=Field_1.' x_1_Name ';'])
140           if isequal(numel(x_1),2) 
141               x_1=linspace(x_1(1),x_1(2),nbpoints_x_1);
142           end
143           if isequal(numel(y_1),2) 
144               y_1=linspace(y_1(1),y_1(2),nbpoints_y_1);
145           end
146           [vec_X_1,vec_Y_1]=meshgrid(x_1,y_1);
147       end
148       vec_X_1=reshape(vec_X_1,[],1);
149       vec_Y_1=reshape(vec_Y_1,[],1);
150       if testfalse_1
151           FFName_1=Field_1.ListVarName{VarType_1.errorflag};         
152           eval(['vec_FF_1=Field_1.' FFName_1 ';'])
153           vec_FF_1=reshape(vec_FF_1,[],1);
154           indsel=find(~vec_FF_1);
155           vec_X_1=vec_X_1(indsel);
156           vec_Y_1=vec_Y_1(indsel);
157       end
158       if testU % vector fields
159            U_1_Name=Field_1.ListVarName{VarType_1.vector_x};
160            V_1_Name=Field_1.ListVarName{VarType_1.vector_y};
161            eval(['vec_U_1=Field_1.' U_1_Name ';'])
162            eval(['vec_V_1=Field_1.' V_1_Name ';'])
163            nbpoints_x_1=size(vec_U_1,2);
164            nbpoints_y_1=size(vec_U_1,1);
165            vec_U_1=reshape(vec_U_1,nbpoints_x_1*nbpoints_y_1,1);
166            vec_V_1=reshape(vec_V_1,nbpoints_x_1*nbpoints_y_1,1);
167            if testfalse_1
168                vec_U_1=vec_U_1(indsel);
169                vec_V_1=vec_V_1(indsel);
170            end           
171       else %(~testU && ~testU_1)
172           A_1_Name=Field_1.ListVarName{ivar_C_1};
173           eval(['vec_A_1=Field_1.' A_1_Name ';'])   
174           nbpoints_x_1=size(vec_A_1,2);
175           nbpoints_y_1=size(vec_A_1,1);%TODO: use a faster interpolation method for a regular grid (size(x)=2)
176           vec_A_1=reshape(vec_A_1,nbpoints_x_1*nbpoints_y_1,1);
177           if testfalse_1
178                vec_A_1=vec_A_1(indsel);
179           end
180       end
181
182       if ~isequal(vec_X_1,vec_X) && ~isequal(vec_Y_1,vec_Y) % if the unstructured positions are not the same
183           if testU
184               vec_U_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_U_1,vec_X,vec_Y);  %interpolate vectors in the second field
185               vec_V_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_V_1,vec_X,vec_Y);  %interpolate vectors in the second field   
186           else
187               vec_A_1=griddata_uvmat(vec_X_1,vec_Y_1,double(vec_A_1),vec_X,vec_Y);  %interpolate vectors in the second field
188           end
189       end
190       if testU
191           UName=Field.ListVarName{VarType.vector_x};
192           VName=Field.ListVarName{VarType.vector_y}; 
193           eval(['vec_U=Field.' UName ';'])
194           eval(['vec_V=Field.' VName ';'])       
195           vec_U=reshape(vec_U,numel(vec_U),1);
196           vec_V=reshape(vec_V,numel(vec_V),1);
197           eval(['SubData.' UName '=vec_U-vec_U_1;'])
198           eval(['SubData.' VName '=vec_V-vec_V_1;'])
199       else
200           AName=Field.ListVarName{ivar_C};
201          % size(Field.vort)
202           eval(['SubData.' AName '=Field.' AName '-vec_A_1;'])
203       end
204   else  %structured coordiantes
205       XName=Field.ListVarName{VarType.coord(2)};
206       YName=Field.ListVarName{VarType.coord(1)};
207       eval(['x=Field.' XName ';'])
208       eval(['y=Field.' YName ';'])
209       if testX_1 %unstructured coordinates for the second field
210           errormsg='the second input scalar is not on a regular grid: comparison option not implemented';
211           return
212       else
213           XName_1=Field.ListVarName{VarType_1.coord(2)};
214           YName_1=Field.ListVarName{VarType_1.coord(1)};
215           eval(['x_1=Field_1.' XName_1 ';'])
216           eval(['y_1=Field_1.' YName_1 ';'])
217       end
218       if testU % vector fields
219           UName=Field.ListVarName{VarType.vector_x};
220           VName=Field.ListVarName{VarType.vector_y};
221           U_1_Name=Field_1.ListVarName{VarType_1.vector_x};
222           V_1_Name=Field_1.ListVarName{VarType_1.vector_y};
223           eval(['U_1=Field_1.' U_1_Name ';'])
224           eval(['V_1=Field_1.' V_1_Name ';'])
225           if ~isequal(x_1,x)||~isequal(y_1,y)
226                [X_1,Y_1]=meshgrid(x_1,y_1);
227                U_1 =interp2(X_1,Y_1,U_1,x,y');
228                V_1 =interp2(X_1,Y_1,V_1,x,y');
229           end
230           eval(['SubData.' UName '=Field.' UName '-U_1;'])
231           eval(['SubData.' VName '=Field.' VName '-V_1;'])
232       else
233           AName=Field.ListVarName{ivar_C};
234           A_1_Name=Field_1.ListVarName{ivar_C_1};
235           eval(['A_1=double(Field_1.' A_1_Name ');'])
236           if ~isequal(x_1,x)||~isequal(y_1,y)
237                [X_1,Y_1]=meshgrid(x_1,y_1);
238                A_1 =interp2(X_1,Y_1,A_1,x,y');
239           end
240           eval(['SubData.' AName '=double(Field.' AName ')-A_1;'])
241       end
242   end
243end
244
245% merge a vector field and a scalar as second input
246if testU && ~testU_1
247    AName_1=Field_1.ListVarName{ivar_C_1};
248    if isfield(Field_1,'VarAttribute') & numel(Field_1.VarAttribute)>=ivar_C_1
249        AAttr=Field_1.VarAttribute{ivar_C_1} ;
250    else
251        AAttr=[];
252    end
253    if testX_1 %unstructured coordinate
254       XName_1=Field_1.ListVarName{VarType_1.coord_x};
255       YName_1=Field_1.ListVarName{VarType_1.coord_y};
256       DimCell=Field_1.VarDimName([VarType_1.coord_x VarType_1.coord_y ]);
257       if isfield(Field_1,'VarAttribute')
258           if numel(Field_1.VarAttribute)>=VarType_1.coord_x
259                XAttr=Field_1.VarAttribute{VarType_1.coord_x};
260           else
261                XAttr=[];
262           end
263           if numel(Field_1.VarAttribute)>=VarType_1.coord_y
264               YAttr=Field_1.VarAttribute{VarType_1.coord_y};
265           else
266               YAttr=[];
267           end
268           SubData.VarAttribute=[SubData.VarAttribute {XAttr} {YAttr}];
269       end
270    else
271       XName_1=Field_1.ListVarName{VarType_1.coord(2)};
272       YName_1=Field_1.ListVarName{VarType_1.coord(1)};
273%        DimCell=[{YName_1} {XName_1}];
274       if isfield(Field_1,'VarAttribute')
275           if numel(Field_1.VarAttribute)>=VarType_1.coord(2)
276                XAttr=Field_1.VarAttribute{VarType_1.coord(2)} ;
277           else
278                XAttr=[];
279           end
280           if numel(Field_1.VarAttribute)>=VarType_1.coord(1)
281               YAttr=Field_1.VarAttribute{VarType_1.coord(1)} ;
282           else
283               YAttr=[];
284           end
285           SubData.VarAttribute=[SubData.VarAttribute {YAttr} {XAttr}];
286       end
287    end 
288    %look for previously used variable names
289    XName_1_1=XName_1;%default
290    YName_1_1=YName_1;%default
291    AName_1_1=AName_1;%default
292    for iprev=1:numel(SubData.ListVarName)
293        switch SubData.ListVarName{iprev}
294            case XName_1
295                XName_1_1=[XName_1 '_1'];
296            case YName_1
297                YName_1_1=[YName_1 '_1'];
298            case AName_1
299                AName_1_1=[AName_1 '_1'];
300        end
301    end     
302    if ~testX_1% if the second field has structured coordinates
303        ilist_1=find(strcmp(AName_1,Field_1.ListVarName));
304        DimCell_1=Field_1.VarDimName{ilist_1};
305        if numel(DimCell_1)>2
306             DimCell={XName_1_1,YName_1_1, [{YName_1_1,XName_1_1} DimCell_1(3:end)]};
307        else
308          DimCell={XName_1_1,YName_1_1, {YName_1_1,XName_1_1}};
309        end
310    else
311        DimCell=[DimCell Field_1.VarDimName(ivar_C_1)];
312    end
313    SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {AName_1_1}];
314    % check that a  different dimension name is used for the two fields
315     if testX_1% if the second field has unstructured coordinates
316        for icell=1:numel(DimCell)
317            if isequal(DimCell{icell}{1},SubData.VarDimName{1}{1})
318                DimCell{icell}{1}=[DimCell{icell}{1} '_1'];
319            end
320        end
321     end
322    SubData.VarDimName=[SubData.VarDimName DimCell];
323    if isfield(Field_1,'VarAttribute')
324        SubData.VarAttribute=[SubData.VarAttribute {AAttr}];
325    end
326    eval(['SubData.' XName_1_1 '=Field_1.' XName_1 ';'])
327    eval(['SubData.' YName_1_1 '=Field_1.' YName_1 ';'])
328    eval(['SubData.' AName_1_1 '=Field_1.' AName_1 ';'])
329end
330
331%merge a scalar as the first input and a vector field as second input
332if ~testU && testU_1
333    UName_1=Field_1.ListVarName{VarType_1.vector_x};
334    VName_1=Field_1.ListVarName{VarType_1.vector_y};
335    UAttr=Field_1.VarAttribute{VarType_1.vector_x};
336    VAttr=Field_1.VarAttribute{VarType_1.vector_y};
337    if testX_1 %unstructured coordinate for the second field
338       XName_1=Field_1.ListVarName{VarType_1.coord_x};
339       YName_1=Field_1.ListVarName{VarType_1.coord_y};
340       
341       XAttr=Field_1.VarAttribute{VarType_1.coord_x};
342       YAttr=Field_1.VarAttribute{VarType_1.coord_y};
343%        SubData.ListVarName=[SubData.ListVarName {XName_1} {YName_1}];
344       DimCell=Field_1.VarDimName([VarType_1.coord_x VarType_1.coord_y ]);
345    else
346       XName_1=Field_1.ListVarName{VarType_1.coord(2)};
347       YName_1=Field_1.ListVarName{VarType_1.coord(1)};
348       if numel(Field_1.VarAttribute)>=VarType_1.coord(2)
349           XAttr=Field_1.VarAttribute{VarType_1.coord(2)};
350       else
351           XAttr=[];
352       end
353       if numel(Field_1.VarAttribute)>=VarType_1.coord(1)
354           YAttr=Field_1.VarAttribute{VarType_1.coord(1)};
355       else
356           YAttr=[];
357       end     
358    end 
359    %check for the existence of the same  variable name
360    XName_1_1=XName_1; %default
361    YName_1_1=YName_1; %default
362    UName_1_1=UName_1; %default
363    VName_1_1=VName_1; %default
364    for iprev=1:numel(SubData.ListVarName)
365        switch SubData.ListVarName{iprev}
366            case XName_1
367                XName_1_1=[XName_1 '_1'];
368            case YName_1
369                YName_1_1=[YName_1 '_1'];
370            case UName_1
371                UName_1_1=[UName_1 '_1'];
372            case VName_1
373                VName_1_1=[VName_1 '_1'];
374        end
375    end   
376   
377    if ~testX_1
378          DimCell=[{XName_1_1} {YName_1_1}];
379    end
380    SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {UName_1_1} {VName_1_1}];
381    DimCell=[DimCell Field_1.VarDimName([VarType_1.vector_x VarType_1.vector_y ])];
382    SubData.VarDimName=[SubData.VarDimName DimCell];
383    if isfield(SubData,'VarAttribute')
384        if ~(numel(SubData.VarAttribute)==numel(SubData.ListVarName))
385            for ivar=numel(SubData.VarAttribute)+1:numel(SubData.ListVarName)-4
386                SubData.VarAttribute{ivar}=[];
387            end
388        end
389        SubData.VarAttribute=[SubData.VarAttribute {XAttr} {YAttr} {UAttr} {VAttr}];
390    end
391    SubData.(XName_1_1)=Field_1.(XName_1);
392    SubData.(YName_1_1)=Field_1.(YName_1);
393    SubData.(UName_1_1)=Field_1.(UName_1);
394    SubData.(VName_1_1)=Field_1.(VName_1);
395end
396 
Note: See TracBrowser for help on using the repository browser.