source: trunk/src/sub_field.m @ 404

Last change on this file since 404 was 402, checked in by sommeria, 12 years ago

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