[179]  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]  8  %


 9  % function SubData=sub_field(Field,Field_1)


 10  %


 11  % OUPUT:


 12  % SubData: structure representing the resulting field


 13  %


 14  % INPUT:


[179]  15  % Field: matlab structure representing the first field


 16  % Field_1:matlab structure representing the second field


[8]  17 


 18  function [SubData,errormsg]=sub_field(Field,Field_1)


 19  test_attr=0;


 20  if 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;


 27  end


 28  if 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


 50  end


 51  SubData.ListVarName=Field.ListVarName;


 52  SubData.VarDimName=Field.VarDimName;


 53  if isfield(Field,'VarAttribute')


 54  SubData.VarAttribute=Field.VarAttribute;


 55  end


 56  %reproduce Field by default


 57  for ivar=1:numel(Field.ListVarName)


 58  VarName=Field.ListVarName{ivar};


 59  eval(['SubData.' VarName '=Field.' VarName ';'])


 60  end


 61 


 62  %fields


 63  [CellVarIndex,NbDim,VarTypeCell,errormsg]=find_field_indices(Field);


 64  if ~isempty(errormsg)


 65  errormsg=['invalid first input to sub_field:' errormsg];


 66  return


 67  end


 68  [CellVarIndex_1,NbDim_1,VarTypeCell_1,errormsg]=find_field_indices(Field_1);


 69  if ~isempty(errormsg)


 70  errormsg=['invalid second input to sub_field:' errormsg];


 71  return


 72  end


 73  iselect=find(NbDim==2);


 74  if ~isequal(numel(iselect),1)


 75  errormsg='invalid first input to sub_field: it must contain a single 2D field cell';


 76  return


 77  end


 78  iselect_1=find(NbDim_1==2);


 79  if ~isequal(numel(iselect_1),1)


 80  errormsg='invalid second input to sub_field: it must contain a single 2D field cell';


 81  return


 82  end


 83  VarType=VarTypeCell{iselect};


 84  VarType_1=VarTypeCell_1{iselect_1};


 85  testX=~isempty(VarType.coord_x)&& ~isempty(VarType.coord_y);%unstructured coordiantes


 86  testX_1=~isempty(VarType_1.coord_x)&& ~isempty(VarType_1.coord_y);%unstructured coordiantes


 87  testU=~isempty(VarType.vector_x)&& ~isempty(VarType.vector_y);%vector field


 88  testU_1=~isempty(VarType_1.vector_x)&& ~isempty(VarType_1.vector_y);%vector field


 89  testfalse_1=~isempty(VarType_1.errorflag);


 90  ivar_C=[VarType.scalar VarType.image VarType.color VarType.ancillary]; %defines index (indices) for the scalar or ancillary fields


 91  if numel(ivar_C)>1


 92  errormsg='too many scalar fields in the first input of sub_field.m';


 93  return


 94  end


 95  ivar_C_1=[VarType_1.scalar VarType_1.image VarType_1.color VarType_1.ancillary]; %defines index (indices) for the scalar or ancillary fields


 96  if numel(ivar_C_1)>1


 97  errormsg='too many scalar fields in the second input of sub_field.m';


 98  return


 99  end


 100 


 101  %substract two vector fields or two scalars


 102  if (testU && testU_1)  (~testU && ~testU_1)


 103  %check coincidence in positions


[179]  104  %unstructured coordinates for the first field


 105  if testX


[8]  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);


[109]  111  vec_X=reshape(vec_X,nbpoints,1);


 112  vec_Y=reshape(vec_Y,nbpoints,1);


[8]  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 ';'])


[109]  118 


[8]  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 ';'])


[109]  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


[8]  130  [vec_X_1,vec_Y_1]=meshgrid(x_1,y_1);


 131  end


[179]  132  vec_X_1=reshape(vec_X_1,[],1);


 133  vec_Y_1=reshape(vec_Y_1,[],1);


[8]  134  if testfalse_1


 135  FFName_1=Field_1.ListVarName{VarType_1.errorflag};


 136  eval(['vec_FF_1=Field_1.' FFName_1 ';'])


[179]  137  vec_FF_1=reshape(vec_FF_1,[],1);


[8]  138  indsel=find(~vec_FF_1);


 139  vec_X_1=vec_X_1(indsel);


 140  vec_Y_1=vec_Y_1(indsel);


 141  end


[179]  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 


[8]  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


[109]  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


[8]  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 ';'])


[109]  179  vec_U=reshape(vec_U,numel(vec_U),1);


 180  vec_V=reshape(vec_V,numel(vec_V),1);


[8]  181  eval(['SubData.' UName '=vec_Uvec_U_1;'])


 182  eval(['SubData.' VName '=vec_Vvec_V_1;'])


 183  else


 184  AName=Field.ListVarName{ivar_C};


[127]  185  % size(Field.vort)


[8]  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


[89]  214  eval(['SubData.' UName '=Field.' UName 'U_1;'])


 215  eval(['SubData.' VName '=Field.' VName 'V_1;'])


[8]  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


 227  end


 228 


 229  % merge a vector field and a scalar as second input


 230  if testU && ~testU_1


 231  AName_1=Field_1.ListVarName{ivar_C_1};


[179]  232  if isfield(Field_1,'VarAttribute') & numel(Field_1.VarAttribute)>=ivar_C_1


[8]  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


[179]  243  XAttr=Field_1.VarAttribute{VarType_1.coord_x};


[8]  244  else


 245  XAttr=[];


 246  end


 247  if numel(Field_1.VarAttribute)>=VarType_1.coord_y


[179]  248  YAttr=Field_1.VarAttribute{VarType_1.coord_y};


[8]  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)


[180]  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


[8]  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 ';'])


 305  end


 306 


 307  %merge a scalar as the first input and a vector field as second input


 308  if ~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];


[165]  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


[8]  363  end


[165]  364  SubData.VarAttribute=[SubData.VarAttribute {XAttr} {YAttr} {UAttr} {VAttr}];


[8]  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 ';'])


 370  end


 371 

