[404]  1 


 2 


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


[8]  4  %


[399]  5  % [DataOut,errormsg]=calc_field(FieldList,DataIn,Coord_interp)


[8]  6  %


[124]  7  % OUTPUT:


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


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


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


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


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


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


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


 15  % error: error flag


[89]  16  % error = 0; OK


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


 18  %


[8]  19  % INPUT:


[399]  20  % FieldList: cell array of strings representing the name(s) of the field(s) to calculate


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


[399]  22  % Coord_interp(:,nb_coord) optional set of coordinates to interpolate the field (use with thin plate shell)


[89]  23  %


[8]  24  % FUNCTION related


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


 26  % file, depending on the scalar


[8]  27 


[399]  28  function [DataOut,errormsg]=calc_field(FieldList,DataIn,Coord_interp)


[8]  29 


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


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


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


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


 34  % 'der': needs spatial derivatives


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


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


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


 38 


[404]  39  %% list of field options implemented


 40  FieldOptions={'velocity';...%image correlation corresponding to a vel vector


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


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


 43  'vort';...%vorticity


 44  'div';...%divergence


 45  'strain';...%rate of strain


 46  'u';... %u velocity component


 47  'v';... %v velocity component


 48  'w';... %w velocity component


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


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


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


[399]  52  if ~exist('FieldList','var')


[404]  53  DataOut=FieldOptions;% gives the list of possible field inputs in the absence of input


 54  return


 55  end


 56 


 57  %% check input


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


 59  DataIn=[];


 60  end


 61  if ischar(FieldList)


 62  FieldList={FieldList};%convert a string input to a cell with one string element


 63  end


 64  check_grid=0;


 65  check_der=0;


 66  check_calc=ones(size(FieldList));


 67  for ilist=1:length(FieldList)


 68  switch FieldList{ilist}


 69  case {'u','v'}


 70  check_grid=1;% needs a regular grid


 71  case{'vort','div','strain'}% needs spatial derivatives spatial derivatives


 72  check_der=1;


[408]  73  case {'velocity','norm_vel','ima_cor'};


[404]  74  otherwise


 75  check_calc(ilist)=0;


 76  end


 77  end


 78  FieldList=FieldList(check_calc==1);


[405]  79  % if isempty(FieldList)


 80  % DataOut=DataIn;


 81  % return


 82  % end


[404]  83  if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))


 84  nbcoord=3;


[8]  85  else


[404]  86  nbcoord=2;


 87  end


 88  ListVarName={};


 89  ValueList={};


 90  RoleList={};


 91  units_cell={};


 92 


 93  %% interpolation with new civ data


 94  if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var')  check_grid check_der)


[406]  95  %reproduce global attributes


 96  DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute;


[404]  97  for ilist=1:numel(DataOut.ListGlobalAttribute)


 98  DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});


[8]  99  end


[406]  100 


[405]  101  %create a default grid if needed


 102  if ~exist('Coord_interp','var')


[406]  103  XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates


 104  YMax=max(max(DataIn.SubRange(2,:,:)));


 105  XMin=min(min(DataIn.SubRange(1,:,:)));


 106  YMin=min(min(DataIn.SubRange(2,:,:)));


[405]  107  if ~isfield(DataIn,'Mesh')


 108  DataIn.Mesh=sqrt(2*(XMaxXMin)*(YMaxYMin)/numel(DataIn.Coord_tps));


 109  % adjust the mesh to a value 1, 2 , 5 *10^n


 110  ord=10^(floor(log10(DataIn.Mesh)));%order of magnitude


 111  if DataIn.Mesh/ord>=5


 112  DataIn.Mesh=5*ord;


 113  elseif DataIn.Mesh/ord>=2


 114  DataIn.Mesh=2*ord;


 115  else


 116  DataIn.Mesh=ord;


 117  end


 118  end


 119  coord_x=XMin:DataIn.Mesh:XMax;


 120  coord_y=YMin:DataIn.Mesh:YMax;


[406]  121  npx=length(coord_x);


 122  npy=length(coord_y);


 123  DataOut.coord_x=[XMin XMax];


 124  DataOut.coord_y=[YMin YMax];


[405]  125  [XI,YI]=meshgrid(coord_x,coord_y);


 126  XI=reshape(XI,[],1);


 127  YI=reshape(YI,[],1);


 128  Coord_interp=[XI YI];


 129  end


[406]  130 


 131  %initialise output


[404]  132  nb_sites=size(Coord_interp,1);


 133  nb_coord=size(Coord_interp,2);


[406]  134  nbval=zeros(nb_sites,1);


 135  NbSubDomain=size(DataIn.SubRange,3);


 136  DataOut.ListVarName={'coord_y','coord_x'};


 137  DataOut.VarDimName={{'coord_y'},{'coord_x'}};


 138  DataOut.VarAttribute{1}=[];


 139  DataOut.VarAttribute{2}=[];


[404]  140  for ilist=1:length(FieldList)


 141  switch FieldList{ilist}


 142  case 'velocity'


 143  DataOut.U=zeros(nb_sites,1);


 144  DataOut.V=zeros(nb_sites,1);


 145  otherwise


 146  DataOut.(FieldList{ilist})=zeros(nb_sites,1);


[389]  147  end


[404]  148  end


[406]  149 


 150  % interpolate data in each subdomain


[404]  151  for isub=1:NbSubDomain


 152  nbvec_sub=DataIn.NbSites(isub);


 153  check_range=(Coord_interp >=ones(nb_sites,1)*DataIn.SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*DataIn.SubRange(:,2,isub)');


 154  ind_sel=find(sum(check_range,2)==nb_coord);


 155  %rho smoothing parameter


 156  % epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites


 157  % ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs


 158  nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)


[405]  159  if check_grid


[404]  160  EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'


[246]  161  end


[404]  162  if check_der


 163  [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'


 164  end


[412]  165  ListVar={};


[399]  166  for ilist=1:length(FieldList)


[412]  167  var_count=numel(ListVar);


[404]  168  switch FieldList{ilist}


 169  case 'velocity'


[412]  170  ListVar=[ListVar {'U', 'V'}];


[411]  171  VarAttribute{var_count+1}.Role='vector_x';


 172  VarAttribute{var_count+2}.Role='vector_y';


[404]  173  DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);


 174  DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);


 175  case 'u'


[412]  176  ListVar=[ListVar {'u'}];


[411]  177  VarAttribute{var_count+1}.Role='scalar';


[406]  178  DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);


[404]  179  case 'v'


[412]  180  ListVar=[ListVar {'v'}];


[411]  181  VarAttribute{var_count+1}.Role='scalar';


[406]  182  DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);


[405]  183  case 'norm_vel'


[412]  184  ListVar=[ListVar {'norm_vel'}];


[411]  185  VarAttribute{var_count+1}.Role='scalar';


 186  U=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);


[404]  187  V=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);


 188  DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V);


 189  case 'vort'


[412]  190  ListVar=[ListVar {'vort'}];


[411]  191  VarAttribute{var_count+1}.Role='scalar';


 192  DataOut.vort(ind_sel)=DataOut.vort(ind_sel)EMDY *DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);


[404]  193  case 'div'


[412]  194  ListVar=[ListVar {'div'}];


[411]  195  VarAttribute{var_count+1}.Role='scalar';


[404]  196  DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDY *DataIn.V_tps(1:nbvec_sub+3,isub);


 197  case 'strain'


[412]  198  ListVar=[ListVar {'strain'}];


[411]  199  VarAttribute{var_count+1}.Role='scalar';


[404]  200  DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);


[246]  201  end


 202  end


[404]  203  end


[411]  204  DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang


[412]  205  nbval(nbval==0)=1;% to avoid division by zero for averaging


[411]  206  if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed


 207  DataOut.ListVarName=[DataOut.ListVarName {'FF'}];


 208  DataOut.VarDimName=[DataOut.VarDimName {{'coord_y','coord_x'}}];


 209  DataOut.VarAttribute{length(DataOut.ListVarName)}.Role='errorflag';


 210  end


[412]  211  DataOut.ListVarName=[DataOut.ListVarName ListVar];


 212  for ifield=1:numel(ListVar)


[411]  213  VarDimName{ifield}={'coord_y','coord_x'};


[412]  214  DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval;


 215  DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx);


[411]  216  end


 217  DataOut.FF=reshape(DataOut.FF,npy,npx);


 218  DataOut.VarDimName=[DataOut.VarDimName VarDimName];


 219  DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];


[404]  220  else


 221 


 222  %% civx data


 223  DataOut=DataIn;


 224  for ilist=1:length(FieldList)


 225  if ~isempty(FieldList{ilist})


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


 227  ListVarName=[ListVarName VarName];


 228  ValueList=[ValueList Value];


 229  RoleList=[RoleList Role];


 230  units_cell=[units_cell units];


[246]  231  end


[124]  232  end


[404]  233  %erase previous data (except coordinates)


 234  for ivar=nbcoord+1:length(DataOut.ListVarName)


 235  VarName=DataOut.ListVarName{ivar};


 236  DataOut=rmfield(DataOut,VarName);


 237  end


 238  DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);


 239  if isfield(DataOut,'VarDimName')


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


 241  else


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


 243  return


 244  end


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


 246  %append new data


 247  DataOut.ListVarName=[DataOut.ListVarName ListVarName];


 248  for ivar=1:length(ListVarName)


 249  DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};


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


 251  DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};


 252  DataOut.(ListVarName{ivar})=ValueList{ivar};


 253  end


[8]  254  end


 255 


[246]  256 


[404]  257 


[8]  258  %%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%


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


 260  VarName={};


 261  ValCell={};


 262  Role={};


 263  units_cell={};


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


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


 266  else


 267  units='pixel';


 268  end


 269  if isfield(DataIn,'U')


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


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


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


 273  units_cell=[units_cell {units}];


 274  end


 275  if isfield(DataIn,'V')


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


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


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


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


 280  end


 281  if isfield(DataIn,'W')


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


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


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


 285  units_cell=[units_cell {units}];


 286  end


 287  if isfield(DataIn,'F')


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


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


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


 291  units_cell=[units_cell {[]}];


 292  end


 293  if isfield(DataIn,'FF')


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


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


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


 297  units_cell=[units_cell {[]}];


 298  end


 299 


 300  %%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%


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


 302  VarName={};


 303  ValCell={};


 304  Role={};


 305  units={};


 306  if isfield(DataIn,'C')


 307  VarName{1}='C';


 308  ValCell{1}=DataIn.C;


 309  Role={'ancillary'};


 310  units={[]};


 311  end


 312 


 313  %%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%


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


 315  VarName={};


 316  ValCell={};


 317  Role={};


 318  units={};


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


 320  VarName{1}='norm_vel';


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


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


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


 324  end


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


 326  Role{1}='scalar';


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


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


[124]  329  else


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


[124]  331  end


 332  end


[8]  333 


 334 


 335 


 336  %%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%


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


 338  VarName={};


 339  ValCell={};


 340  Role={};


 341  units={};


 342  if isfield(DataIn,'DjUi')


 343  VarName{1}='vort';


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


 345  siz=size(ValCell{1});


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


 347  Role{1}='scalar';


 348  if isfield(DataIn,'TimeUnit')


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


 350  else


 351  units={[]};


 352  end


[124]  353  end


[8]  354 


 355  %%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%


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


 357  VarName={};


 358  ValCell={};


 359  Role={};


 360  units={};


 361  if isfield(DataIn,'DjUi')


 362  VarName{1}='div';


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


 364  siz=size(ValCell{1});


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


 366  Role{1}='scalar';


 367  if isfield(DataIn,'TimeUnit')


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


 369  else


 370  units={[]};


 371  end


[124]  372  end


[8]  373 


 374  %%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%


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


 376  VarName={};


 377  ValCell={};


 378  Role={};


 379  units={};


 380  if isfield(DataIn,'DjUi')


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


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


 383  siz=size(ValCell{1});


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


[156]  385  Role{1}='scalar';


[124]  386  if isfield(DataIn,'TimeUnit')


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


[124]  388  else


[8]  389  units={[]};


[124]  390  end


 391  end


[8]  392 


 393  %%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%


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


 395  VarName={};


 396  ValCell={};


 397  Role={};


 398  units={};


 399  if isfield(DataIn,'U')


 400  VarName{1}='U';


 401  ValCell{1}=DataIn.U;


 402  Role{1}='scalar';


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


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


 405  else


 406  units={'pixel'};


 407  end


 408  end


 409 


 410  %%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%


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


 412  VarName={};


 413  ValCell={};


 414  Role={};


 415  units={};


 416  if isfield(DataIn,'V')


 417  VarName{1}='V';


 418  ValCell{1}=DataIn.V;


 419  Role{1}='scalar';


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


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


 422  else


 423  units={'pixel'};


 424  end


 425  end


 426 


 427  %%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%


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


 429  VarName={};


 430  ValCell={};


 431  Role={};


 432  units={};


 433  if isfield(DataIn,'W')


 434  VarName{1}='W';


 435  ValCell{1}=DataIn.W;


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


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


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


 439  else


 440  units={'pixel'};


 441  end


 442  end


 443 


 444  %%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%


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


 446  VarName={};


 447  ValCell={};


 448  Role={};


 449  units={};


 450  if isfield(DataIn,'W')


 451  VarName{1}='W';


 452  ValCell{1}=DataIn.W;


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


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


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


 456  else


 457  units={'pixel'};


 458  end


 459  end


 460 


 461  %%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%


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


 463  VarName={};


 464  ValCell={};


 465  Role={};


 466  units={};


 467  if isfield(DataIn,'E')


 468  VarName{1}='E';


 469  ValCell{1}=DataIn.E;


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


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


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


 473  else


 474  units={'pixel'};


 475  end


 476  end


 477 

