[515]  1 


[575]  2  %'calc_field_tps': defines fields (velocity, vort, div...) from civ data and calculate them with tps interpolation


[515]  3  %


[581]  4  % [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord_tps,NbSites,SubRange,FieldVar,FieldName,Coord_interp)


[515]  5  %


 6  % OUTPUT:


 7  % DataOut: structure representing the output fields


 8  %


 9  % INPUT:


 10  % Coord_tps:


 11  % NbSites


 12  % SubRange


 13  % FieldVar


[581]  14  % FieldName: cell array representing the list of operations (eg div, rot..)


[515]  15  % Coord_interp: coordiantes of sites on which the fields need to be calculated


 16 


[581]  17  function [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord_tps,NbSites,SubRange,FieldVar,FieldName,Coord_interp)


[515]  18 


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


 20  % a type is associated to each scalar:


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


 22  % 'vel': calculated from velocity components, continuous field (interpolated with velocity)


 23  % 'der': needs spatial derivatives


 24  % 'var': the scalar name corresponds to a field name in the netcdf files


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


 26  % the scalar is calculated from other fields, as explicited below


 27  errormsg='';


 28 


 29  %% nbre of subdomains


 30  if ndims(Coord_interp)==3


 31  nb_coord=size(Coord_interp,3);


 32  npx=size(Coord_interp,2);


 33  npy=size(Coord_interp,1);


 34  nb_sites=npx*npy;


 35  Coord_interp=reshape(Coord_interp,nb_sites,nb_coord);


 36  else


 37  nb_coord=size(Coord_interp,2);


 38  nb_sites=size(Coord_interp,1);


 39  end


 40  NbSubDomain=size(Coord_tps,3);


 41  nbval=zeros(nb_sites,1);


 42 


 43  %% list of operations


 44  check_grid=0;


 45  check_der=0;


[581]  46  for ilist=1:length(FieldName)


 47  FieldNameType=regexprep(FieldName{ilist},'(.+','');


 48  switch FieldNameType


[515]  49  case 'vec'


[521]  50  check_grid=1;


[515]  51  DataOut.U=zeros(nb_sites,1);


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


 53  VarAttribute{1}.Role='vector_x';


 54  VarAttribute{2}.Role='vector_y';


 55  case {'U','V','norm'}


 56  check_grid=1;


[581]  57  DataOut.(FieldNameType)=zeros(nb_sites,1);


[515]  58  VarAttribute{1}.Role='scalar';


 59  case {'curl','div','strain'}


 60  check_der=1;


[581]  61  DataOut.(FieldNameType)=zeros(nb_sites,1);


[515]  62  VarAttribute{1}.Role='scalar';


 63  end


 64  end


 65  Attr_FF.Role='errorflag';


 66  VarAttribute=[VarAttribute {Attr_FF}];


 67 


 68  %% loop on subdomains


 69  for isub=1:NbSubDomain


 70  nbvec_sub=NbSites(isub);


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


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


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


 74  if check_grid


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


 76  end


 77  if check_der


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


 79  end


 80  ListVar={};


[581]  81  for ilist=1:length(FieldName)


[515]  82  var_count=numel(ListVar);


[581]  83  switch FieldName{ilist}


[515]  84  case 'vec(U,V)'


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


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


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


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


 89  DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);


 90  case 'U'


 91  ListVar=[ListVar {'U'}];


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


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


 94  case 'V'


 95  ListVar=[ListVar {'V'}];


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


 97  DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);


 98  case 'norm(U,V)'


 99  ListVar=[ListVar {'norm'}];


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


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


 102  V=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);


 103  DataOut.norm(ind_sel)=sqrt(U.*U+V.*V);


 104  case 'curl(U,V)'


 105  ListVar=[ListVar {'curl'}];


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


 107  DataOut.curl(ind_sel)=DataOut.curl(ind_sel)EMDY *FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);


 108  case 'div(U,V)'


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


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


 111  DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*FieldVar(1:nbvec_sub+3,isub,1)+EMDY *FieldVar(1:nbvec_sub+3,isub,2);


 112  case 'strain(U,V)'


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


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


 115  DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);


 116  end


 117  end


 118  end


 119  DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang


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


 121  ListFieldOut=fieldnames(DataOut);


 122  for ifield=1:numel(ListFieldOut)


 123  DataOut.(ListFieldOut{ifield})=DataOut.(ListFieldOut{ifield})./nbval;


 124  DataOut.(ListFieldOut{ifield})=reshape(DataOut.(ListFieldOut{ifield}),npy,npx);


 125  end


 126 


 127 


 128 


 129 

