[515]  1 


 2  %'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them


 3  %


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


 5  %


 6  % OUTPUT:


 7  % DataOut: structure representing the output fields


 8  %


 9  % INPUT:


 10  % Coord_tps:


 11  % NbSites


 12  % SubRange


 13  % FieldVar


 14  % Operation: cell array representing the list of operations (eg div, rot..)


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


 16 


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


 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;


 46  for ilist=1:length(Operation)


 47  OperationType=regexprep(Operation{ilist},'(.+','');


 48  switch OperationType


 49  case 'vec'


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


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


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


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


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


 55  check_grid=1;


 56  DataOut.(OperationType)=zeros(nb_sites,1);


 57  VarAttribute{1}.Role='scalar';


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


 59  check_der=1;


 60  DataOut.(OperationType)=zeros(nb_sites,1);


 61  VarAttribute{1}.Role='scalar';


 62  end


 63  end


 64  Attr_FF.Role='errorflag';


 65  VarAttribute=[VarAttribute {Attr_FF}];


 66 


 67  %% loop on subdomains


 68  for isub=1:NbSubDomain


 69  nbvec_sub=NbSites(isub);


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


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


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


 73  if check_grid


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


 75  end


 76  if check_der


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


 78  end


 79  ListVar={};


 80  for ilist=1:length(Operation)


 81  var_count=numel(ListVar);


 82  switch Operation{ilist}


 83  case 'vec(U,V)'


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


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


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


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


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


 89  case 'U'


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


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


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


 93  case 'V'


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


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


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


 97  case 'norm(U,V)'


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


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


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


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


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


 103  case 'curl(U,V)'


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


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


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


 107  case 'div(U,V)'


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


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


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


 111  case 'strain(U,V)'


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


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


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


 115  end


 116  end


 117  end


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


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


 120  ListFieldOut=fieldnames(DataOut);


 121  for ifield=1:numel(ListFieldOut)


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


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


 124  end


 125 


 126 


 127 


 128 

