[8] | 1 | %'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them
|
---|
| 2 | %---------------------------------------------------------------------
|
---|
[123] | 3 | % [DataOut,errormsg]=calc_field(FieldName,DataIn)
|
---|
[8] | 4 | %
|
---|
| 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
|
---|
| 11 | % 'der': needs spatial derivatives
|
---|
| 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:
|
---|
[123] | 18 | % FieldName: string representing the name of the field
|
---|
[8] | 19 | % DataIn: structure representing the field, as defined in check_field_srtructure.m
|
---|
[89] | 20 | %
|
---|
[8] | 21 | % FUNCTION related
|
---|
| 22 | % varname_generator.m: determines the field names to read in the netcdf file, depending on the scalar
|
---|
| 23 |
|
---|
| 24 | function [DataOut,errormsg]=calc_field(FieldName,DataIn)
|
---|
| 25 |
|
---|
| 26 | %list of defined scalars to display in menus (in addition to 'ima_cor').
|
---|
| 27 | % a type is associated to each scalar:
|
---|
| 28 | % 'discrete': related to the individual velocity vectors, not interpolated by patch
|
---|
| 29 | % 'vel': calculated from velocity components, continuous field (interpolated with velocity)
|
---|
| 30 | % 'der': needs spatial derivatives
|
---|
| 31 | % 'var': the scalar name corresponds to a field name in the netcdf files
|
---|
| 32 | % a specific variable name for civ1 and civ2 fields are also associated, if
|
---|
[123] | 33 | % the scalar is calculated from other fields, as explicited below
|
---|
| 34 |
|
---|
[8] | 35 | %list_scal={title, type, civ1 var name,civ2 var name}
|
---|
| 36 | list_field={'velocity';...%image correlation corresponding to a vel vector
|
---|
[123] | 37 | 'ima_cor';...%image correlation corresponding to a vel vector
|
---|
[8] | 38 | 'norm_vel';...%norm of the velocity
|
---|
| 39 | 'vort';...%vorticity
|
---|
| 40 | 'div';...%divergence
|
---|
| 41 | 'strain';...%rate of strain
|
---|
| 42 | 'u';... %u velocity component
|
---|
| 43 | 'v';... %v velocity component
|
---|
| 44 | 'w';... %w velocity component
|
---|
| 45 | 'w_normal';... %w velocity component normal to the plane
|
---|
| 46 | 'error'}; %error associated to a vector (for stereo or patch)
|
---|
| 47 | errormsg=[]; %default error message
|
---|
| 48 | if ~exist('FieldName','var')
|
---|
| 49 | DataOut=list_field;% gives the list of possible fields in the absence of input
|
---|
| 50 | else
|
---|
| 51 | if ~exist('DataIn','var')
|
---|
| 52 | DataIn=[];
|
---|
| 53 | end
|
---|
| 54 | if ischar(FieldName)
|
---|
| 55 | FieldName={FieldName};
|
---|
| 56 | end
|
---|
| 57 | if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
|
---|
| 58 | nbcoord=3;
|
---|
| 59 | else
|
---|
| 60 | nbcoord=2;
|
---|
| 61 | end
|
---|
[123] | 62 | DataOut=DataIn; %reproduce global attribute
|
---|
[8] | 63 | ListVarName={};
|
---|
| 64 | ValueList={};
|
---|
| 65 | RoleList={};
|
---|
| 66 | units_cell={};
|
---|
| 67 | for ilist=1:length(FieldName)
|
---|
| 68 | [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}
|
---|
| 69 | ListVarName=[ListVarName VarName];
|
---|
| 70 | ValueList=[ValueList Value];
|
---|
| 71 | RoleList=[RoleList Role];
|
---|
| 72 | units_cell=[units_cell units];
|
---|
[123] | 73 | end
|
---|
[8] | 74 | %erase previous data (except coordinates)
|
---|
| 75 | for ivar=nbcoord+1:length(DataOut.ListVarName)
|
---|
| 76 | VarName=DataOut.ListVarName{ivar};
|
---|
| 77 | DataOut=rmfield(DataOut,VarName);
|
---|
| 78 | end
|
---|
| 79 | DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
|
---|
| 80 | if isfield(DataOut,'VarDimName')
|
---|
| 81 | DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
|
---|
| 82 | else
|
---|
| 83 | errormsg='element .VarDimName missing in input data';
|
---|
| 84 | return
|
---|
| 85 | end
|
---|
| 86 | DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
|
---|
| 87 | %append new data
|
---|
| 88 | DataOut.ListVarName=[DataOut.ListVarName ListVarName];
|
---|
| 89 | for ivar=1:length(ListVarName)
|
---|
| 90 | DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
|
---|
| 91 | DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
|
---|
| 92 | DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
|
---|
| 93 | eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])
|
---|
| 94 | end
|
---|
| 95 | end
|
---|
| 96 |
|
---|
| 97 | %%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
|
---|
| 98 | function [VarName,ValCell,Role,units_cell]=velocity(DataIn)
|
---|
| 99 | VarName={};
|
---|
| 100 | ValCell={};
|
---|
| 101 | Role={};
|
---|
| 102 | units_cell={};
|
---|
| 103 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
| 104 | units=[DataIn.CoordUnit '/' DataIn.TimeUnit];
|
---|
| 105 | else
|
---|
| 106 | units='pixel';
|
---|
| 107 | end
|
---|
| 108 | if isfield(DataIn,'U')
|
---|
| 109 | VarName=[VarName {'U'}];
|
---|
| 110 | ValCell=[ValCell {DataIn.U}];
|
---|
| 111 | Role=[Role {'vector_x'}];
|
---|
| 112 | units_cell=[units_cell {units}];
|
---|
| 113 | end
|
---|
| 114 | if isfield(DataIn,'V')
|
---|
| 115 | VarName=[VarName {'V'}];
|
---|
| 116 | ValCell=[ValCell {DataIn.V}];
|
---|
| 117 | Role=[Role {'vector_y'}];
|
---|
| 118 | units_cell=[units_cell {units}];
|
---|
| 119 | end
|
---|
| 120 | if isfield(DataIn,'W')
|
---|
| 121 | VarName=[VarName {'W'}];
|
---|
| 122 | ValCell=[ValCell {DataIn.W}];
|
---|
| 123 | Role=[Role {'vector_z'}];
|
---|
| 124 | units_cell=[units_cell {units}];
|
---|
| 125 | end
|
---|
| 126 | if isfield(DataIn,'F')
|
---|
| 127 | VarName=[VarName {'F'}];
|
---|
| 128 | ValCell=[ValCell {DataIn.F}];
|
---|
| 129 | Role=[Role {'warnflag'}];
|
---|
| 130 | units_cell=[units_cell {[]}];
|
---|
| 131 | end
|
---|
| 132 | if isfield(DataIn,'FF')
|
---|
| 133 | VarName=[VarName,{'FF'}];
|
---|
| 134 | ValCell=[ValCell {DataIn.FF}];
|
---|
| 135 | Role=[Role {'errorflag'}];
|
---|
| 136 | units_cell=[units_cell {[]}];
|
---|
| 137 | end
|
---|
| 138 |
|
---|
| 139 | %%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%
|
---|
| 140 | function [VarName,ValCell,Role,units]=ima_cor(DataIn)
|
---|
| 141 | VarName={};
|
---|
| 142 | ValCell={};
|
---|
| 143 | Role={};
|
---|
| 144 | units={};
|
---|
| 145 | if isfield(DataIn,'C')
|
---|
| 146 | VarName{1}='C';
|
---|
| 147 | ValCell{1}=DataIn.C;
|
---|
| 148 | Role={'ancillary'};
|
---|
| 149 | units={[]};
|
---|
| 150 | end
|
---|
| 151 |
|
---|
| 152 | %%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%
|
---|
| 153 | function [VarName,ValCell,Role,units]=norm_vel(DataIn)
|
---|
| 154 | VarName={};
|
---|
| 155 | ValCell={};
|
---|
| 156 | Role={};
|
---|
| 157 | units={};
|
---|
| 158 | if isfield(DataIn,'U') && isfield(DataIn,'V')
|
---|
| 159 | VarName{1}='norm_vel';
|
---|
| 160 | ValCell{1}=DataIn.U.*DataIn.U+ DataIn.V.*DataIn.V;
|
---|
| 161 | if isfield(DataIn,'W') && isequal(size(DataIn.W),size(DataIn.U))
|
---|
| 162 | ValCell{1}=ValCell{1}+DataIn.W.*DataIn.W;
|
---|
| 163 | end
|
---|
| 164 | ValCell{1}=sqrt(ValCell{1});
|
---|
| 165 | Role{1}='scalar';
|
---|
| 166 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
| 167 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
| 168 | else
|
---|
| 169 | units={'pixel'};
|
---|
| 170 | end
|
---|
| 171 | end
|
---|
| 172 |
|
---|
| 173 |
|
---|
| 174 |
|
---|
| 175 | %%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%
|
---|
| 176 | function [VarName,ValCell,Role,units]=vort(DataIn)
|
---|
| 177 | VarName={};
|
---|
| 178 | ValCell={};
|
---|
| 179 | Role={};
|
---|
| 180 | units={};
|
---|
| 181 | if isfield(DataIn,'DjUi')
|
---|
| 182 | VarName{1}='vort';
|
---|
| 183 | ValCell{1}=DataIn.DjUi(:,1,2)-DataIn.DjUi(:,2,1); %vorticity
|
---|
| 184 | siz=size(ValCell{1});
|
---|
| 185 | ValCell{1}=reshape(ValCell{1},siz(1),1);
|
---|
| 186 | Role{1}='scalar';
|
---|
| 187 | if isfield(DataIn,'TimeUnit')
|
---|
| 188 | units={[DataIn.TimeUnit '-1']};
|
---|
| 189 | else
|
---|
| 190 | units={[]};
|
---|
| 191 | end
|
---|
| 192 | end
|
---|
| 193 |
|
---|
| 194 | %%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%
|
---|
| 195 | function [VarName,ValCell,Role,units]=div(DataIn)
|
---|
| 196 | VarName={};
|
---|
| 197 | ValCell={};
|
---|
| 198 | Role={};
|
---|
| 199 | units={};
|
---|
| 200 | if isfield(DataIn,'DjUi')
|
---|
| 201 | VarName{1}='div';
|
---|
| 202 | ValCell{1}=DataIn.DjUi(:,1,1)+DataIn.DjUi(:,2,2); %DUDX+DVDY
|
---|
| 203 | siz=size(ValCell{1});
|
---|
| 204 | ValCell{1}=reshape(ValCell{1},siz(1),1);
|
---|
| 205 | Role{1}='scalar';
|
---|
| 206 | if isfield(DataIn,'TimeUnit')
|
---|
| 207 | units={[DataIn.TimeUnit '-1']};
|
---|
| 208 | else
|
---|
| 209 | units={[]};
|
---|
| 210 | end
|
---|
| 211 | end
|
---|
| 212 |
|
---|
| 213 | %%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
|
---|
| 214 | function [VarName,ValCell,Role,units]=strain(DataIn)
|
---|
| 215 | VarName={};
|
---|
| 216 | ValCell={};
|
---|
| 217 | Role={};
|
---|
| 218 | units={};
|
---|
| 219 | if isfield(DataIn,'DjUi')
|
---|
| 220 | VarName{1}='strain';
|
---|
| 221 | ValCell{1}=DataIn.DjUi(:,1,2)+DataIn.DjUi(:,2,1);%DVDX+DUDY
|
---|
| 222 | siz=size(ValCell{1});
|
---|
| 223 | ValCell{1}=reshape(ValCell{1},siz(1),1);
|
---|
| 224 | if isfield(DataIn,'TimeUnit')
|
---|
| 225 | units={[DataIn.TimeUnit '-1']};
|
---|
| 226 | else
|
---|
| 227 | units={[]};
|
---|
| 228 | end
|
---|
| 229 | end
|
---|
| 230 |
|
---|
| 231 | %%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
|
---|
| 232 | function [VarName,ValCell,Role,units]=u(DataIn)
|
---|
| 233 | VarName={};
|
---|
| 234 | ValCell={};
|
---|
| 235 | Role={};
|
---|
| 236 | units={};
|
---|
| 237 | if isfield(DataIn,'U')
|
---|
| 238 | VarName{1}='U';
|
---|
| 239 | ValCell{1}=DataIn.U;
|
---|
| 240 | Role{1}='scalar';
|
---|
| 241 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
| 242 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
| 243 | else
|
---|
| 244 | units={'pixel'};
|
---|
| 245 | end
|
---|
| 246 | end
|
---|
| 247 |
|
---|
| 248 | %%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
|
---|
| 249 | function [VarName,ValCell,Role,units]=v(DataIn)
|
---|
| 250 | VarName={};
|
---|
| 251 | ValCell={};
|
---|
| 252 | Role={};
|
---|
| 253 | units={};
|
---|
| 254 | if isfield(DataIn,'V')
|
---|
| 255 | VarName{1}='V';
|
---|
| 256 | ValCell{1}=DataIn.V;
|
---|
| 257 | Role{1}='scalar';
|
---|
| 258 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
| 259 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
| 260 | else
|
---|
| 261 | units={'pixel'};
|
---|
| 262 | end
|
---|
| 263 | end
|
---|
| 264 |
|
---|
| 265 | %%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
|
---|
| 266 | function [VarName,ValCell,Role,units]=w(DataIn)
|
---|
| 267 | VarName={};
|
---|
| 268 | ValCell={};
|
---|
| 269 | Role={};
|
---|
| 270 | units={};
|
---|
| 271 | if isfield(DataIn,'W')
|
---|
| 272 | VarName{1}='W';
|
---|
| 273 | ValCell{1}=DataIn.W;
|
---|
| 274 | Role{1}='scalar';%will remain unchanged by projection
|
---|
| 275 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
| 276 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
| 277 | else
|
---|
| 278 | units={'pixel'};
|
---|
| 279 | end
|
---|
| 280 | end
|
---|
| 281 |
|
---|
| 282 | %%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
|
---|
| 283 | function [VarName,ValCell,Role,units]=w_normal(DataIn)
|
---|
| 284 | VarName={};
|
---|
| 285 | ValCell={};
|
---|
| 286 | Role={};
|
---|
| 287 | units={};
|
---|
| 288 | if isfield(DataIn,'W')
|
---|
| 289 | VarName{1}='W';
|
---|
| 290 | ValCell{1}=DataIn.W;
|
---|
| 291 | Role{1}='vector_z';%will behave like a vector component by projection
|
---|
| 292 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
| 293 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
| 294 | else
|
---|
| 295 | units={'pixel'};
|
---|
| 296 | end
|
---|
| 297 | end
|
---|
| 298 |
|
---|
| 299 | %%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
|
---|
| 300 | function [VarName,ValCell,Role,units]=error(DataIn)
|
---|
| 301 | VarName={};
|
---|
| 302 | ValCell={};
|
---|
| 303 | Role={};
|
---|
| 304 | units={};
|
---|
| 305 | if isfield(DataIn,'E')
|
---|
| 306 | VarName{1}='E';
|
---|
| 307 | ValCell{1}=DataIn.E;
|
---|
| 308 | Role{1}='ancillary'; %TODO CHECK units in actual fields
|
---|
| 309 | if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
|
---|
| 310 | units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
|
---|
| 311 | else
|
---|
| 312 | units={'pixel'};
|
---|
| 313 | end
|
---|
| 314 | end
|
---|
| 315 |
|
---|