[124]  1  %'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them


[8]  2  %


[123]  3  % [DataOut,errormsg]=calc_field(FieldName,DataIn)


[8]  4  %


[124]  5  % OUTPUT:


 6  % Scal: matlab vector representing the scalar values (length nbvec defined by var_read)


[123]  7  % if no input, Scal=list of programmed scalar names (to put in menus)


 8  % if only the field name is put as input, vec_A=type of scalar, which can be:


 9  % 'discrete': related to the individual velocity vectors, not interpolated by patch


[8]  10  % 'vel': scalar calculated solely from velocity components


[124]  11  % 'der': needs spatial derivatives


[8]  12  % 'var': the scalar name directly corresponds to a field name in the netcdf files


 13  % error: error flag


[89]  14  % error = 0; OK


 15  % error = 1; the prescribed scalar cannot be read or calculated from available fields


 16  %


[8]  17  % INPUT:


[123]  18  % FieldName: string representing the name of the field


[8]  19  % DataIn: structure representing the field, as defined in check_field_srtructure.m


[89]  20  %


[8]  21  % FUNCTION related


[124]  22  % varname_generator.m: determines the field names to read in the netcdf


 23  % file, depending on the scalar


[8]  24 


 25  function [DataOut,errormsg]=calc_field(FieldName,DataIn)


 26 


 27  %list of defined scalars to display in menus (in addition to 'ima_cor').


[124]  28  % a type is associated to each scalar:


[8]  29  % 'discrete': related to the individual velocity vectors, not interpolated by patch


[124]  30  % 'vel': calculated from velocity components, continuous field (interpolated with velocity)


 31  % 'der': needs spatial derivatives


[8]  32  % 'var': the scalar name corresponds to a field name in the netcdf files


 33  % a specific variable name for civ1 and civ2 fields are also associated, if


[123]  34  % the scalar is calculated from other fields, as explicited below


 35 


[8]  36  %list_scal={title, type, civ1 var name,civ2 var name}


 37  list_field={'velocity';...%image correlation corresponding to a vel vector


[124]  38  'ima_cor';...%image correlation corresponding to a vel vector


 39  'norm_vel';...%norm of the velocity


 40  'vort';...%vorticity


 41  'div';...%divergence


 42  'strain';...%rate of strain


 43  'u';... %u velocity component


 44  'v';... %v velocity component


 45  'w';... %w velocity component


 46  'w_normal';... %w velocity component normal to the plane


 47  'error'}; %error associated to a vector (for stereo or patch)


[8]  48  errormsg=[]; %default error message


[124]  49  if ~exist('FieldName','var')


[8]  50  DataOut=list_field;% gives the list of possible fields in the absence of input


 51  else


 52  if ~exist('DataIn','var')


 53  DataIn=[];


 54  end


 55  if ischar(FieldName)


 56  FieldName={FieldName};


 57  end


 58  if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))


 59  nbcoord=3;


 60  else


 61  nbcoord=2;


[124]  62  end


[123]  63  DataOut=DataIn; %reproduce global attribute


[8]  64  ListVarName={};


 65  ValueList={};


 66  RoleList={};


 67  units_cell={};


[124]  68  for ilist=1:length(FieldName)


[8]  69  [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}


 70  ListVarName=[ListVarName VarName];


 71  ValueList=[ValueList Value];


 72  RoleList=[RoleList Role];


 73  units_cell=[units_cell units];


[124]  74  end


 75  %erase previous data (except coordinates)


[8]  76  for ivar=nbcoord+1:length(DataOut.ListVarName)


 77  VarName=DataOut.ListVarName{ivar};


 78  DataOut=rmfield(DataOut,VarName);


[124]  79  end


[8]  80  DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);


 81  if isfield(DataOut,'VarDimName')


 82  DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);


 83  else


 84  errormsg='element .VarDimName missing in input data';


 85  return


 86  end


 87  DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);


 88  %append new data


 89  DataOut.ListVarName=[DataOut.ListVarName ListVarName];


 90  for ivar=1:length(ListVarName)


[124]  91  DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};


 92  DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};


[8]  93  DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};


 94  eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])


 95  end


 96  end


 97 


 98  %%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%


 99  function [VarName,ValCell,Role,units_cell]=velocity(DataIn)


 100  VarName={};


 101  ValCell={};


 102  Role={};


 103  units_cell={};


 104  if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')


 105  units=[DataIn.CoordUnit '/' DataIn.TimeUnit];


 106  else


 107  units='pixel';


 108  end


 109  if isfield(DataIn,'U')


 110  VarName=[VarName {'U'}];


 111  ValCell=[ValCell {DataIn.U}];


 112  Role=[Role {'vector_x'}];


 113  units_cell=[units_cell {units}];


 114  end


 115  if isfield(DataIn,'V')


 116  VarName=[VarName {'V'}];


 117  ValCell=[ValCell {DataIn.V}];


[124]  118  Role=[Role {'vector_y'}];


[8]  119  units_cell=[units_cell {units}];


 120  end


 121  if isfield(DataIn,'W')


 122  VarName=[VarName {'W'}];


 123  ValCell=[ValCell {DataIn.W}];


 124  Role=[Role {'vector_z'}];


 125  units_cell=[units_cell {units}];


 126  end


 127  if isfield(DataIn,'F')


 128  VarName=[VarName {'F'}];


 129  ValCell=[ValCell {DataIn.F}];


 130  Role=[Role {'warnflag'}];


 131  units_cell=[units_cell {[]}];


 132  end


 133  if isfield(DataIn,'FF')


 134  VarName=[VarName,{'FF'}];


 135  ValCell=[ValCell {DataIn.FF}];


 136  Role=[Role {'errorflag'}];


 137  units_cell=[units_cell {[]}];


 138  end


 139 


 140  %%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%


 141  function [VarName,ValCell,Role,units]=ima_cor(DataIn)


 142  VarName={};


 143  ValCell={};


 144  Role={};


 145  units={};


 146  if isfield(DataIn,'C')


 147  VarName{1}='C';


 148  ValCell{1}=DataIn.C;


 149  Role={'ancillary'};


 150  units={[]};


 151  end


 152 


 153  %%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%


 154  function [VarName,ValCell,Role,units]=norm_vel(DataIn)


 155  VarName={};


 156  ValCell={};


 157  Role={};


 158  units={};


 159  if isfield(DataIn,'U') && isfield(DataIn,'V')


 160  VarName{1}='norm_vel';


[124]  161  ValCell{1}=DataIn.U.*DataIn.U+ DataIn.V.*DataIn.V;


 162  if isfield(DataIn,'W') && isequal(size(DataIn.W),size(DataIn.U))


 163  ValCell{1}=ValCell{1}+DataIn.W.*DataIn.W;


 164  end


 165  ValCell{1}=sqrt(ValCell{1});


 166  Role{1}='scalar';


 167  if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')


[8]  168  units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};


[124]  169  else


[8]  170  units={'pixel'};


[124]  171  end


 172  end


[8]  173 


 174 


 175 


 176  %%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%


 177  function [VarName,ValCell,Role,units]=vort(DataIn)


 178  VarName={};


 179  ValCell={};


 180  Role={};


 181  units={};


 182  if isfield(DataIn,'DjUi')


 183  VarName{1}='vort';


 184  ValCell{1}=DataIn.DjUi(:,1,2)DataIn.DjUi(:,2,1); %vorticity


 185  siz=size(ValCell{1});


 186  ValCell{1}=reshape(ValCell{1},siz(1),1);


 187  Role{1}='scalar';


 188  if isfield(DataIn,'TimeUnit')


 189  units={[DataIn.TimeUnit '1']};


 190  else


 191  units={[]};


 192  end


[124]  193  end


[8]  194 


 195  %%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%


 196  function [VarName,ValCell,Role,units]=div(DataIn)


 197  VarName={};


 198  ValCell={};


 199  Role={};


 200  units={};


 201  if isfield(DataIn,'DjUi')


 202  VarName{1}='div';


 203  ValCell{1}=DataIn.DjUi(:,1,1)+DataIn.DjUi(:,2,2); %DUDX+DVDY


 204  siz=size(ValCell{1});


 205  ValCell{1}=reshape(ValCell{1},siz(1),1);


 206  Role{1}='scalar';


 207  if isfield(DataIn,'TimeUnit')


 208  units={[DataIn.TimeUnit '1']};


 209  else


 210  units={[]};


 211  end


[124]  212  end


[8]  213 


 214  %%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%


 215  function [VarName,ValCell,Role,units]=strain(DataIn)


 216  VarName={};


 217  ValCell={};


 218  Role={};


 219  units={};


 220  if isfield(DataIn,'DjUi')


[124]  221  VarName{1}='strain';


 222  ValCell{1}=DataIn.DjUi(:,1,2)+DataIn.DjUi(:,2,1);%DVDX+DUDY


 223  siz=size(ValCell{1});


 224  ValCell{1}=reshape(ValCell{1},siz(1),1);


 225  if isfield(DataIn,'TimeUnit')


[8]  226  units={[DataIn.TimeUnit '1']};


[124]  227  else


[8]  228  units={[]};


[124]  229  end


 230  end


[8]  231 


 232  %%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%


 233  function [VarName,ValCell,Role,units]=u(DataIn)


 234  VarName={};


 235  ValCell={};


 236  Role={};


 237  units={};


 238  if isfield(DataIn,'U')


 239  VarName{1}='U';


 240  ValCell{1}=DataIn.U;


 241  Role{1}='scalar';


 242  if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')


 243  units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};


 244  else


 245  units={'pixel'};


 246  end


 247  end


 248 


 249  %%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%


 250  function [VarName,ValCell,Role,units]=v(DataIn)


 251  VarName={};


 252  ValCell={};


 253  Role={};


 254  units={};


 255  if isfield(DataIn,'V')


 256  VarName{1}='V';


 257  ValCell{1}=DataIn.V;


 258  Role{1}='scalar';


 259  if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')


 260  units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};


 261  else


 262  units={'pixel'};


 263  end


 264  end


 265 


 266  %%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%


 267  function [VarName,ValCell,Role,units]=w(DataIn)


 268  VarName={};


 269  ValCell={};


 270  Role={};


 271  units={};


 272  if isfield(DataIn,'W')


 273  VarName{1}='W';


 274  ValCell{1}=DataIn.W;


 275  Role{1}='scalar';%will remain unchanged by projection


 276  if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')


 277  units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};


 278  else


 279  units={'pixel'};


 280  end


 281  end


 282 


 283  %%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%


 284  function [VarName,ValCell,Role,units]=w_normal(DataIn)


 285  VarName={};


 286  ValCell={};


 287  Role={};


 288  units={};


 289  if isfield(DataIn,'W')


 290  VarName{1}='W';


 291  ValCell{1}=DataIn.W;


[124]  292  Role{1}='vector_z';%will behave like a vector component by projection


[8]  293  if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')


 294  units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};


 295  else


 296  units={'pixel'};


 297  end


 298  end


 299 


 300  %%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%


 301  function [VarName,ValCell,Role,units]=error(DataIn)


 302  VarName={};


 303  ValCell={};


 304  Role={};


 305  units={};


 306  if isfield(DataIn,'E')


 307  VarName{1}='E';


 308  ValCell{1}=DataIn.E;


 309  Role{1}='ancillary'; %TODO CHECK units in actual fields


 310  if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')


 311  units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};


 312  else


 313  units={'pixel'};


 314  end


 315  end


 316 

