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


[8]  2  %


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


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


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


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


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


[89]  21  %


[8]  22  % FUNCTION related


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


 24  % file, depending on the scalar


[517]  25  function [FieldOptions,ColorList]=calc_field


 26  %function [DataOut,errormsg]=calc_field(FieldList,DataIn,Coord_interp)


[8]  27 


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


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


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


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


 32  % 'der': needs spatial derivatives


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


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


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


 36 


[404]  37  %% list of field options implemented


[515]  38  FieldOptions={'vec(U,V)';...%image correlation corresponding to a vel vector


 39  'C';...%image correlation corresponding to a vel vector


 40  'norm(U,V)';...%norm of the velocity


 41  'curl(U,V)';...%vorticity


 42  'div(U,V)';...%divergence


 43  'strain(U,V)';...%rate of strain


 44  'U';... %u velocity component


 45  'V';... %v velocity component


[124]  46  'w';... %w velocity component


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


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


[517]  49  ColorList={'C';...%image correlation corresponding to a vel vector


 50  'norm(U,V)';...%norm of the velocity


 51  'U';... %u velocity component


 52  'V';... %v velocity component


 53  }


 54  % errormsg=[]; %default error message


 55  % if ~exist('FieldList','var')


 56  % DataOut=FieldOptions;% gives the list of possible field inputs in the absence of input


 57  % return


 58  % end


 59  return


[404]  60  %% check input


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


 62  DataIn=[];


 63  end


 64  if ischar(FieldList)


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


 66  end


 67  check_grid=0;


 68  check_der=0;


 69  check_calc=ones(size(FieldList));


 70  for ilist=1:length(FieldList)


 71  switch FieldList{ilist}


[491]  72  case {'u','v','velocity','norm_vel','ima_cor'}


[404]  73  check_grid=1;% needs a regular grid


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


 75  check_der=1;


[491]  76  % case {'velocity','norm_vel','ima_cor'};


[404]  77  otherwise


 78  check_calc(ilist)=0;


 79  end


 80  end


 81  FieldList=FieldList(check_calc==1);


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


 83  nbcoord=3;


[8]  84  else


[404]  85  nbcoord=2;


 86  end


 87  ListVarName={};


 88  ValueList={};


 89  RoleList={};


 90  units_cell={};


 91 


[515]  92  %% reproduce global attributes


 93  DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute;


 94  for ilist=1:numel(DataOut.ListGlobalAttribute)


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


 96  end


 97 


[404]  98  %% interpolation with new civ data


[515]  99  [CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(DataIn);


 100  icell_tps=[];


 101  for icell=1:numel(CellVarIndex)


 102  VarType=VarTypeCell{icell};


 103  if NbDimVec(icell)>=2 && ~isempty(VarType.subrange_tps)


 104  icell_tps=[icell_tps icell];


[8]  105  end


[515]  106  end


[406]  107 


[515]  108  %if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var')  check_grid check_der)


 109  if ~isempty(icell_tps)


[405]  110  %create a default grid if needed


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


[515]  112  for ind=1:numel(icell_tps)


 113  SubRange=DataIn.(DataIn.ListVarName{VarType{icell_tps(ind)}.subrange_tps});


 114  XMax(ind)=max(max(SubRange(1,:,:)));% extrema of the coordinates


 115  YMax(ind)=max(max(SubRange(2,:,:)));


 116  XMin(ind)=min(min(SubRange(1,:,:)));


 117  YMin(ind)=min(min(SubRange(2,:,:)));


 118  end


 119  XMax=max(XMax);


 120  YMax=max(YMax);


 121  XMin=min(XMin);


 122  YMin=min(YMin);


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


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


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


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


 127  if DataIn.Mesh/ord>=5


 128  DataIn.Mesh=5*ord;


 129  elseif DataIn.Mesh/ord>=2


 130  DataIn.Mesh=2*ord;


 131  else


 132  DataIn.Mesh=ord;


 133  end


 134  end


[501]  135  coord_x=XMin:DataIn.Mesh:XMax;% increase the recommanded mesh to visualisation


[405]  136  coord_y=YMin:DataIn.Mesh:YMax;


[406]  137  DataOut.coord_x=[XMin XMax];


 138  DataOut.coord_y=[YMin YMax];


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


[491]  140  Coord_interp=cat(3,XI,YI);%[XI YI];


[405]  141  end


[491]  142  npx=size(Coord_interp,2);


 143  npy=size(Coord_interp,1);


 144  Coord_interp=reshape(Coord_interp,npx*npy,size(Coord_interp,3));


[515]  145 


[406]  146  %initialise output


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


 148  nb_coord=size(Coord_interp,2);


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


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


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


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


 153  DataOut.VarAttribute{1}=[];


 154  DataOut.VarAttribute{2}=[];


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


 156  switch FieldList{ilist}


 157  case 'velocity'


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


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


 160  otherwise


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


[389]  162  end


[404]  163  end


[406]  164 


 165  % interpolate data in each subdomain


[515]  166  for icell=icell_tps


 167  for isub=1:NbSubDomain


 168  nbvec_sub=DataIn.(DataIn.ListVarName{VarType{icell}.nbsites_tps});


 169  SubRange=DataIn.(DataIn.ListVarName{VarType{icell}.subrange_tps});


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


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


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


 173  Coord_tps=DataIn.(DataIn.ListVarName{VarType{icell}.coord_tps});


 174  if check_grid


 175  EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'


 176  end


 177  if check_der


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


 179  end


 180  ListVar={};


 181  U_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_x_tps});


 182  V_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_y_tps});


 183  for ilist=1:length(FieldList)


 184  var_count=numel(ListVar);


 185  switch FieldList{ilist}


 186  case 'velocity'


 187  ListVar=[ListVar {'U', 'V'}];


 188  VarAttribute{var_count+1}.Role='vector_x';


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


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


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


 192  case 'u'


 193  ListVar=[ListVar {'u'}];


 194  VarAttribute{var_count+1}.Role='scalar';


 195  DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);


 196  case 'v'


 197  ListVar=[ListVar {'v'}];


 198  VarAttribute{var_count+1}.Role='scalar';


 199  DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);


 200  case 'norm_vel'


 201  ListVar=[ListVar {'norm_vel'}];


 202  VarAttribute{var_count+1}.Role='scalar';


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


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


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


 206  case 'vort'


 207  ListVar=[ListVar {'vort'}];


 208  VarAttribute{var_count+1}.Role='scalar';


 209  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);


 210  case 'div'


 211  ListVar=[ListVar {'div'}];


 212  VarAttribute{var_count+1}.Role='scalar';


 213  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);


 214  case 'strain'


 215  ListVar=[ListVar {'strain'}];


 216  VarAttribute{var_count+1}.Role='scalar';


 217  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);


 218  end


 219  end


[246]  220  end


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


 222  nbval(nbval==0)=1;% to avoid division by zero for averaging


 223  if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed


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


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


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


[404]  227  end


[515]  228  DataOut.ListVarName=[DataOut.ListVarName ListVar];


 229  for ifield=1:numel(ListVar)


 230  VarDimName{ifield}={'coord_y','coord_x'};


 231  DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval;


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


[246]  233  end


[515]  234  DataOut.FF=reshape(DataOut.FF,npy,npx);


 235  DataOut.VarDimName=[DataOut.VarDimName VarDimName];


 236  DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];


[404]  237  end


 238  else


 239 


 240  %% civx data


 241  DataOut=DataIn;


 242  for ilist=1:length(FieldList)


 243  if ~isempty(FieldList{ilist})


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


 245  ListVarName=[ListVarName VarName];


 246  ValueList=[ValueList Value];


 247  RoleList=[RoleList Role];


 248  units_cell=[units_cell units];


[246]  249  end


[124]  250  end


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


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


 253  VarName=DataOut.ListVarName{ivar};


 254  DataOut=rmfield(DataOut,VarName);


 255  end


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


 257  if isfield(DataOut,'VarDimName')


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


 259  else


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


 261  return


 262  end


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


 264  %append new data


 265  DataOut.ListVarName=[DataOut.ListVarName ListVarName];


 266  for ivar=1:length(ListVarName)


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


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


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


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


 271  end


[8]  272  end


 273 


[246]  274 


[404]  275 


[8]  276  %%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%


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


 278  VarName={};


 279  ValCell={};


 280  Role={};


 281  units_cell={};


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


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


 284  else


 285  units='pixel';


 286  end


 287  if isfield(DataIn,'U')


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


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


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


 291  units_cell=[units_cell {units}];


 292  end


 293  if isfield(DataIn,'V')


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


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


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


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


 298  end


 299  if isfield(DataIn,'W')


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


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


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


 303  units_cell=[units_cell {units}];


 304  end


 305  if isfield(DataIn,'F')


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


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


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


 309  units_cell=[units_cell {[]}];


 310  end


 311  if isfield(DataIn,'FF')


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


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


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


 315  units_cell=[units_cell {[]}];


 316  end


 317 


 318  %%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%


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


 320  VarName={};


 321  ValCell={};


 322  Role={};


 323  units={};


 324  if isfield(DataIn,'C')


 325  VarName{1}='C';


 326  ValCell{1}=DataIn.C;


 327  Role={'ancillary'};


 328  units={[]};


 329  end


 330 


 331  %%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%


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


 333  VarName={};


 334  ValCell={};


 335  Role={};


 336  units={};


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


 338  VarName{1}='norm_vel';


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


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


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


 342  end


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


 344  Role{1}='scalar';


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


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


[124]  347  else


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


[124]  349  end


 350  end


[8]  351 


 352 


 353 


 354  %%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%


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


 356  VarName={};


 357  ValCell={};


 358  Role={};


 359  units={};


 360  if isfield(DataIn,'DjUi')


 361  VarName{1}='vort';


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


 363  siz=size(ValCell{1});


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


 365  Role{1}='scalar';


 366  if isfield(DataIn,'TimeUnit')


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


 368  else


 369  units={[]};


 370  end


[124]  371  end


[8]  372 


 373  %%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%


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


 375  VarName={};


 376  ValCell={};


 377  Role={};


 378  units={};


 379  if isfield(DataIn,'DjUi')


 380  VarName{1}='div';


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


 382  siz=size(ValCell{1});


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


 384  Role{1}='scalar';


 385  if isfield(DataIn,'TimeUnit')


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


 387  else


 388  units={[]};


 389  end


[124]  390  end


[8]  391 


 392  %%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%


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


 394  VarName={};


 395  ValCell={};


 396  Role={};


 397  units={};


 398  if isfield(DataIn,'DjUi')


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


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


 401  siz=size(ValCell{1});


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


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


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


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


[124]  406  else


[8]  407  units={[]};


[124]  408  end


 409  end


[8]  410 


 411  %%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%


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


 413  VarName={};


 414  ValCell={};


 415  Role={};


 416  units={};


 417  if isfield(DataIn,'U')


 418  VarName{1}='U';


 419  ValCell{1}=DataIn.U;


 420  Role{1}='scalar';


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


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


 423  else


 424  units={'pixel'};


 425  end


 426  end


 427 


 428  %%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%


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


 430  VarName={};


 431  ValCell={};


 432  Role={};


 433  units={};


 434  if isfield(DataIn,'V')


 435  VarName{1}='V';


 436  ValCell{1}=DataIn.V;


 437  Role{1}='scalar';


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


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


 440  else


 441  units={'pixel'};


 442  end


 443  end


 444 


 445  %%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%


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


 447  VarName={};


 448  ValCell={};


 449  Role={};


 450  units={};


 451  if isfield(DataIn,'W')


 452  VarName{1}='W';


 453  ValCell{1}=DataIn.W;


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


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


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


 457  else


 458  units={'pixel'};


 459  end


 460  end


 461 


 462  %%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%


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


 464  VarName={};


 465  ValCell={};


 466  Role={};


 467  units={};


 468  if isfield(DataIn,'W')


 469  VarName{1}='W';


 470  ValCell{1}=DataIn.W;


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


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


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


 474  else


 475  units={'pixel'};


 476  end


 477  end


 478 


 479  %%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%


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


 481  VarName={};


 482  ValCell={};


 483  Role={};


 484  units={};


 485  if isfield(DataIn,'E')


 486  VarName{1}='E';


 487  ValCell{1}=DataIn.E;


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


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


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


 491  else


 492  units={'pixel'};


 493  end


 494  end


 495 

