source: trunk/src/sub_field.m @ 405

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

many bugs corrected: composition of 2 input fields, arrow colors from different scalar fields...
File size: 15.2 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
303          DimCell=[{XName_1_1} {YName_1_1}];
304    end
305    SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {AName_1_1}];
306    DimCell=[DimCell Field_1.VarDimName(ivar_C_1)]; %(TODO: check for dimension names)
307    if testX_1
308        for icell=1:numel(DimCell)
309            if isequal(DimCell{icell}{1},SubData.VarDimName{1}{1})
310                DimCell{icell}{1}=[DimCell{icell}{1} '_1'];
311            end
312        end
313    end
314    SubData.VarDimName=[SubData.VarDimName DimCell];
315    if isfield(Field_1,'VarAttribute')
316        SubData.VarAttribute=[SubData.VarAttribute {AAttr}];
317    end
318    eval(['SubData.' XName_1_1 '=Field_1.' XName_1 ';'])
319    eval(['SubData.' YName_1_1 '=Field_1.' YName_1 ';'])
320    eval(['SubData.' AName_1_1 '=Field_1.' AName_1 ';'])
321end
322
323%merge a scalar as the first input and a vector field as second input
324if ~testU && testU_1
325    UName_1=Field_1.ListVarName{VarType_1.vector_x};
326    VName_1=Field_1.ListVarName{VarType_1.vector_y};
327    UAttr=Field_1.VarAttribute{VarType_1.vector_x};
328    VAttr=Field_1.VarAttribute{VarType_1.vector_y};
329    if testX_1 %unstructured coordinate for the second field
330       XName_1=Field_1.ListVarName{VarType_1.coord_x};
331       YName_1=Field_1.ListVarName{VarType_1.coord_y};
332       
333       XAttr=Field_1.VarAttribute{VarType_1.coord_x};
334       YAttr=Field_1.VarAttribute{VarType_1.coord_y};
335%        SubData.ListVarName=[SubData.ListVarName {XName_1} {YName_1}];
336       DimCell=Field_1.VarDimName([VarType_1.coord_x VarType_1.coord_y ]);
337    else
338       XName_1=Field_1.ListVarName{VarType_1.coord(2)};
339       YName_1=Field_1.ListVarName{VarType_1.coord(1)};
340       if numel(Field_1.VarAttribute)>=VarType_1.coord(2)
341           XAttr=Field_1.VarAttribute{VarType_1.coord(2)};
342       else
343           XAttr=[];
344       end
345       if numel(Field_1.VarAttribute)>=VarType_1.coord(1)
346           YAttr=Field_1.VarAttribute{VarType_1.coord(1)};
347       else
348           YAttr=[];
349       end     
350    end 
351    %check for the existence of the same  variable name
352    XName_1_1=XName_1; %default
353    YName_1_1=YName_1; %default
354    UName_1_1=UName_1; %default
355    VName_1_1=VName_1; %default
356    for iprev=1:numel(SubData.ListVarName)
357        switch SubData.ListVarName{iprev}
358            case XName_1
359                XName_1_1=[XName_1 '_1'];
360            case YName_1
361                YName_1_1=[YName_1 '_1'];
362            case UName_1
363                UName_1_1=[UName_1 '_1'];
364            case VName_1
365                VName_1_1=[VName_1 '_1'];
366        end
367    end   
368   
369    if ~testX_1
370          DimCell=[{XName_1_1} {YName_1_1}];
371    end
372    SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {UName_1_1} {VName_1_1}];
373    DimCell=[DimCell Field_1.VarDimName([VarType_1.vector_x VarType_1.vector_y ])];
374    SubData.VarDimName=[SubData.VarDimName DimCell];
375    if isfield(SubData,'VarAttribute')
376        if ~(numel(SubData.VarAttribute)==numel(SubData.ListVarName))
377            for ivar=numel(SubData.VarAttribute)+1:numel(SubData.ListVarName)-4
378                SubData.VarAttribute{ivar}=[];
379            end
380        end
381        SubData.VarAttribute=[SubData.VarAttribute {XAttr} {YAttr} {UAttr} {VAttr}];
382    end
383    SubData.(XName_1_1)=Field_1.(XName_1);
384    SubData.(YName_1_1)=Field_1.(YName_1);
385    SubData.(UName_1_1)=Field_1.(UName_1);
386    SubData.(VName_1_1)=Field_1.(VName_1);
387end
388 
Note: See TracBrowser for help on using the repository browser.