source: trunk/src/sub_field.m @ 323

Last change on this file since 323 was 180, checked in by sommeria, 14 years ago

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