[575]  1  %'calc_field': defines fields (velocity, vort, div...) from civx data (old conventions) 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);


[575]  123  if ~isfield(DataIn,'CoordMesh')


 124  DataIn.Coord


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


[405]  126  % adjust the mesh to a value 1, 2 , 5 *10^n


[575]  127  ord=10^(floor(log10(DataIn.CoordMesh)));%order of magnitude


 128  if DataIn.CoordMesh/ord>=5


 129  DataIn.CoordMesh=5*ord;


 130  elseif DataIn.CoordMesh/ord>=2


 131  DataIn.CoordMesh=2*ord;


[405]  132  else


[575]  133  DataIn.CoordMesh=ord;


[405]  134  end


 135  end


[575]  136  coord_x=XMin:DataIn.CoordMesh:XMax;% increase the recommanded mesh to visualisation


 137  coord_y=YMin:DataIn.CoordMesh:YMax;


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


 139  DataOut.coord_y=[YMin YMax];


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


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


[405]  142  end


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


 144  npy=size(Coord_interp,1);


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


[515]  146 


[406]  147  %initialise output


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


 149  nb_coord=size(Coord_interp,2);


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


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


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


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


 154  DataOut.VarAttribute{1}=[];


 155  DataOut.VarAttribute{2}=[];


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


 157  switch FieldList{ilist}


 158  case 'velocity'


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


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


 161  otherwise


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


[389]  163  end


[404]  164  end


[406]  165 


 166  % interpolate data in each subdomain


[515]  167  for icell=icell_tps


 168  for isub=1:NbSubDomain


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


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


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


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


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


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


 175  if check_grid


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


 177  end


 178  if check_der


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


 180  end


 181  ListVar={};


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


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


 184  for ilist=1:length(FieldList)


 185  var_count=numel(ListVar);


 186  switch FieldList{ilist}


 187  case 'velocity'


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


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


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


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


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


 193  case 'u'


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


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


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


 197  case 'v'


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


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


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


 201  case 'norm_vel'


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


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


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


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


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


 207  case 'vort'


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


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


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


 211  case 'div'


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


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


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


 215  case 'strain'


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


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


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


 219  end


 220  end


[246]  221  end


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


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


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


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


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


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


[404]  228  end


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


 230  for ifield=1:numel(ListVar)


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


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


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


[246]  234  end


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


 236  DataOut.VarDimName=[DataOut.VarDimName VarDimName];


 237  DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];


[404]  238  end


 239  else


 240 


 241  %% civx data


 242  DataOut=DataIn;


 243  for ilist=1:length(FieldList)


 244  if ~isempty(FieldList{ilist})


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


 246  ListVarName=[ListVarName VarName];


 247  ValueList=[ValueList Value];


 248  RoleList=[RoleList Role];


 249  units_cell=[units_cell units];


[246]  250  end


[124]  251  end


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


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


 254  VarName=DataOut.ListVarName{ivar};


 255  DataOut=rmfield(DataOut,VarName);


 256  end


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


 258  if isfield(DataOut,'VarDimName')


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


 260  else


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


 262  return


 263  end


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


 265  %append new data


 266  DataOut.ListVarName=[DataOut.ListVarName ListVarName];


 267  for ivar=1:length(ListVarName)


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


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


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


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


 272  end


[8]  273  end


 274 


[246]  275 


[404]  276 


[8]  277  %%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%


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


 279  VarName={};


 280  ValCell={};


 281  Role={};


 282  units_cell={};


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


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


 285  else


 286  units='pixel';


 287  end


 288  if isfield(DataIn,'U')


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


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


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


 292  units_cell=[units_cell {units}];


 293  end


 294  if isfield(DataIn,'V')


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


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


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


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


 299  end


 300  if isfield(DataIn,'W')


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


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


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


 304  units_cell=[units_cell {units}];


 305  end


 306  if isfield(DataIn,'F')


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


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


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


 310  units_cell=[units_cell {[]}];


 311  end


 312  if isfield(DataIn,'FF')


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


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


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


 316  units_cell=[units_cell {[]}];


 317  end


 318 


 319  %%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%


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


 321  VarName={};


 322  ValCell={};


 323  Role={};


 324  units={};


 325  if isfield(DataIn,'C')


 326  VarName{1}='C';


 327  ValCell{1}=DataIn.C;


 328  Role={'ancillary'};


 329  units={[]};


 330  end


 331 


 332  %%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%


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


 334  VarName={};


 335  ValCell={};


 336  Role={};


 337  units={};


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


 339  VarName{1}='norm_vel';


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


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


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


 343  end


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


 345  Role{1}='scalar';


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


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


[124]  348  else


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


[124]  350  end


 351  end


[8]  352 


 353 


 354 


 355  %%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%


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


 357  VarName={};


 358  ValCell={};


 359  Role={};


 360  units={};


 361  if isfield(DataIn,'DjUi')


 362  VarName{1}='vort';


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


 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  %%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%


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


 376  VarName={};


 377  ValCell={};


 378  Role={};


 379  units={};


 380  if isfield(DataIn,'DjUi')


 381  VarName{1}='div';


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


 383  siz=size(ValCell{1});


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


 385  Role{1}='scalar';


 386  if isfield(DataIn,'TimeUnit')


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


 388  else


 389  units={[]};


 390  end


[124]  391  end


[8]  392 


 393  %%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%


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


 395  VarName={};


 396  ValCell={};


 397  Role={};


 398  units={};


 399  if isfield(DataIn,'DjUi')


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


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


 402  siz=size(ValCell{1});


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


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


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


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


[124]  407  else


[8]  408  units={[]};


[124]  409  end


 410  end


[8]  411 


 412  %%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%


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


 414  VarName={};


 415  ValCell={};


 416  Role={};


 417  units={};


 418  if isfield(DataIn,'U')


 419  VarName{1}='U';


 420  ValCell{1}=DataIn.U;


 421  Role{1}='scalar';


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


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


 424  else


 425  units={'pixel'};


 426  end


 427  end


 428 


 429  %%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%


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


 431  VarName={};


 432  ValCell={};


 433  Role={};


 434  units={};


 435  if isfield(DataIn,'V')


 436  VarName{1}='V';


 437  ValCell{1}=DataIn.V;


 438  Role{1}='scalar';


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


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


 441  else


 442  units={'pixel'};


 443  end


 444  end


 445 


 446  %%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%


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


 448  VarName={};


 449  ValCell={};


 450  Role={};


 451  units={};


 452  if isfield(DataIn,'W')


 453  VarName{1}='W';


 454  ValCell{1}=DataIn.W;


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


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


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


 458  else


 459  units={'pixel'};


 460  end


 461  end


 462 


 463  %%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%


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


 465  VarName={};


 466  ValCell={};


 467  Role={};


 468  units={};


 469  if isfield(DataIn,'W')


 470  VarName{1}='W';


 471  ValCell{1}=DataIn.W;


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


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


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


 475  else


 476  units={'pixel'};


 477  end


 478  end


 479 


 480  %%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%


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


 482  VarName={};


 483  ValCell={};


 484  Role={};


 485  units={};


 486  if isfield(DataIn,'E')


 487  VarName{1}='E';


 488  ValCell{1}=DataIn.E;


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


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


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


 492  else


 493  units={'pixel'};


 494  end


 495  end


 496 

