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


[405]  19 


 20  %% global attributes


[8]  21  test_attr=0;


 22  if isfield(Field,'ListGlobalAttribute')


 23  SubData.ListGlobalAttribute=Field.ListGlobalAttribute;


 24  for ilist=1:numel(Field.ListGlobalAttribute)


 25  AttrName=Field.ListGlobalAttribute{ilist};


[405]  26  SubData.(AttrName)=Field.(AttrName);


[8]  27  end


 28  test_attr=1;


 29  end


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


 52  end


[405]  53 


 54  %% variables


[8]  55  SubData.ListVarName=Field.ListVarName;


 56  SubData.VarDimName=Field.VarDimName;


 57  if isfield(Field,'VarAttribute')


 58  SubData.VarAttribute=Field.VarAttribute;


 59  end


[405]  60  %reproduce the first field Field by default


[8]  61  for ivar=1:numel(Field.ListVarName)


 62  VarName=Field.ListVarName{ivar};


[405]  63  SubData.(VarName)=Field.(VarName);


[8]  64  end


 65 


[405]  66  %% check the two input fields


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


 68  if ~isempty(errormsg)


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


 70  return


 71  end


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


 73  if ~isempty(errormsg)


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


 75  return


 76  end


 77  iselect=find(NbDim==2);


[402]  78  for icell=iselect


 79  if ~isempty(VarTypeCell{icell}.coord_tps)


 80  NbDim(icell)=0;


 81  end


 82  end


 83  iselect=find(NbDim==2);


[8]  84  if ~isequal(numel(iselect),1)


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


 86  return


 87  end


 88  iselect_1=find(NbDim_1==2);


[402]  89  for icell=iselect_1


 90  if ~isempty(VarTypeCell_1{icell}.coord_tps)


 91  NbDim_1(icell)=0;


 92  end


 93  end


 94  iselect_1=find(NbDim_1==2);


[8]  95  if ~isequal(numel(iselect_1),1)


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


 97  return


 98  end


 99  VarType=VarTypeCell{iselect};


 100  VarType_1=VarTypeCell_1{iselect_1};


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


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


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


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


 105  testfalse_1=~isempty(VarType_1.errorflag);


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


 107  if numel(ivar_C)>1


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


 109  return


 110  end


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


 112  if numel(ivar_C_1)>1


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


 114  return


 115  end


 116 


[405]  117  %% substract two vector fields or two scalars


[8]  118  if (testU && testU_1)  (~testU && ~testU_1)


 119  %check coincidence in positions


[179]  120  %unstructured coordinates for the first field


 121  if testX


[8]  122  XName=Field.ListVarName{VarType.coord_x};


 123  YName=Field.ListVarName{VarType.coord_y};


[405]  124  vec_X=Field.(XName);


 125  vec_Y=Field.(YName);


[8]  126  nbpoints=numel(vec_X);


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


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


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


[109]  134 


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


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


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


 147  end


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


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


[8]  150  if testfalse_1


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


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


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


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


 155  vec_X_1=vec_X_1(indsel);


 156  vec_Y_1=vec_Y_1(indsel);


 157  end


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


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


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


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


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


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


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


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


 199  else


 200  AName=Field.ListVarName{ivar_C};


[127]  201  % size(Field.vort)


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


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


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


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


 243  end


 244 


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


 246  if testU && ~testU_1


 247  AName_1=Field_1.ListVarName{ivar_C_1};


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


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


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


[8]  260  else


 261  XAttr=[];


 262  end


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


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


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


[428]  302  if ~testX_1% if the second field has structured coordinates


 303  DimCell={XName_1_1,YName_1_1, {YName_1_1,XName_1_1}};


 304  else


 305  DimCell=[DimCell Field_1.VarDimName(ivar_C_1)];


[8]  306  end


 307  SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {AName_1_1}];


[428]  308  % check that a different dimension name is used for the two fields


 309  if testX_1% if the second field has unstructured coordinates


[180]  310  for icell=1:numel(DimCell)


 311  if isequal(DimCell{icell}{1},SubData.VarDimName{1}{1})


 312  DimCell{icell}{1}=[DimCell{icell}{1} '_1'];


 313  end


 314  end


[428]  315  end


[8]  316  SubData.VarDimName=[SubData.VarDimName DimCell];


 317  if isfield(Field_1,'VarAttribute')


 318  SubData.VarAttribute=[SubData.VarAttribute {AAttr}];


 319  end


 320  eval(['SubData.' XName_1_1 '=Field_1.' XName_1 ';'])


 321  eval(['SubData.' YName_1_1 '=Field_1.' YName_1 ';'])


 322  eval(['SubData.' AName_1_1 '=Field_1.' AName_1 ';'])


 323  end


 324 


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


 326  if ~testU && testU_1


 327  UName_1=Field_1.ListVarName{VarType_1.vector_x};


 328  VName_1=Field_1.ListVarName{VarType_1.vector_y};


 329  UAttr=Field_1.VarAttribute{VarType_1.vector_x};


 330  VAttr=Field_1.VarAttribute{VarType_1.vector_y};


 331  if testX_1 %unstructured coordinate for the second field


 332  XName_1=Field_1.ListVarName{VarType_1.coord_x};


 333  YName_1=Field_1.ListVarName{VarType_1.coord_y};


 334 


 335  XAttr=Field_1.VarAttribute{VarType_1.coord_x};


 336  YAttr=Field_1.VarAttribute{VarType_1.coord_y};


 337  % SubData.ListVarName=[SubData.ListVarName {XName_1} {YName_1}];


 338  DimCell=Field_1.VarDimName([VarType_1.coord_x VarType_1.coord_y ]);


 339  else


 340  XName_1=Field_1.ListVarName{VarType_1.coord(2)};


 341  YName_1=Field_1.ListVarName{VarType_1.coord(1)};


 342  if numel(Field_1.VarAttribute)>=VarType_1.coord(2)


 343  XAttr=Field_1.VarAttribute{VarType_1.coord(2)};


 344  else


 345  XAttr=[];


 346  end


 347  if numel(Field_1.VarAttribute)>=VarType_1.coord(1)


 348  YAttr=Field_1.VarAttribute{VarType_1.coord(1)};


 349  else


 350  YAttr=[];


 351  end


 352  end


 353  %check for the existence of the same variable name


 354  XName_1_1=XName_1; %default


 355  YName_1_1=YName_1; %default


 356  UName_1_1=UName_1; %default


 357  VName_1_1=VName_1; %default


 358  for iprev=1:numel(SubData.ListVarName)


 359  switch SubData.ListVarName{iprev}


 360  case XName_1


 361  XName_1_1=[XName_1 '_1'];


 362  case YName_1


 363  YName_1_1=[YName_1 '_1'];


 364  case UName_1


 365  UName_1_1=[UName_1 '_1'];


 366  case VName_1


 367  VName_1_1=[VName_1 '_1'];


 368  end


[405]  369  end


 370 


[8]  371  if ~testX_1


 372  DimCell=[{XName_1_1} {YName_1_1}];


 373  end


 374  SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {UName_1_1} {VName_1_1}];


 375  DimCell=[DimCell Field_1.VarDimName([VarType_1.vector_x VarType_1.vector_y ])];


 376  SubData.VarDimName=[SubData.VarDimName DimCell];


[165]  377  if isfield(SubData,'VarAttribute')


 378  if ~(numel(SubData.VarAttribute)==numel(SubData.ListVarName))


 379  for ivar=numel(SubData.VarAttribute)+1:numel(SubData.ListVarName)4


 380  SubData.VarAttribute{ivar}=[];


 381  end


[8]  382  end


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


[8]  384  end


[405]  385  SubData.(XName_1_1)=Field_1.(XName_1);


 386  SubData.(YName_1_1)=Field_1.(YName_1);


 387  SubData.(UName_1_1)=Field_1.(UName_1);


 388  SubData.(VName_1_1)=Field_1.(VName_1);


[8]  389  end


 390 

