calc_field

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

0001 %'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them
0002 %---------------------------------------------------------------------
0003 
0004 %
0005 % OUTPUT:
0006 % Scal: matlab vector representing the scalar values (length nbvec defined by var_read)
0007             % if no input , Scal=list of programmed scalar names (to put in menus)
0008             % if only the sclar name is put as input, vec_A=type of scalar, which can be:
0009 %                   'vel': scalar calculated solely from velocity components
0010 %                   'der': needs spatial derivatives
0011 %                   'var': the scalar name directly corresponds to a field name in the netcdf files
0012 % error: error flag
0013             % error = 0; OK
0014             % error = 1; the prescribed scalar cannot be read or calculated from available fields
0015 % INPUT:
0016 % ScalName: string representing the name of the scalar
0017 % DataIn: structure representing the field, as defined in check_field_srtructure.m
0018 
0019 % FUNCTION related
0020 % varname_generator.m: determines the field names to read in the netcdf file, depending on the scalar
0021 
0022 function [DataOut,errormsg]=calc_field(FieldName,DataIn)
0023 
0024 %list of defined scalars to display in menus (in addition to 'ima_cor').
0025 % a type is associated to each scalar:
0026 %              'discrete': related to the individual velocity vectors, not interpolated by patch
0027 %              'vel': calculated from velocity components, continuous field (interpolated with velocity)
0028 %              'der': needs spatial derivatives
0029 %              'var': the scalar name corresponds to a field name in the netcdf files
0030 % a specific variable name for civ1 and civ2 fields are also associated, if
0031 % '' the scalar is calculated from other fields, as explicited below
0032 %list_scal={title, type, civ1 var name,civ2 var name}
0033 list_field={'velocity';...%image correlation corresponding to a vel vector
0034             'ima_cor';...%image correlation corresponding to a vel vector
0035            'norm_vel';...%norm of the velocity
0036            'vort';...%vorticity
0037            'div';...%divergence
0038            'strain';...%rate of strain
0039            'u';... %u velocity component
0040            'v';... %v velocity component
0041            'w';... %w velocity component
0042            'w_normal';... %w velocity component normal to the plane
0043           'error'}; %error associated to a vector (for stereo or patch)
0044 errormsg=[]; %default error message
0045 if ~exist('FieldName','var') 
0046     DataOut=list_field;% gives the list of possible fields in the absence of input
0047 else
0048     if ~exist('DataIn','var')
0049         DataIn=[];
0050     end
0051     if ischar(FieldName)
0052         FieldName={FieldName};
0053     end
0054     if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
0055         nbcoord=3;
0056     else
0057         nbcoord=2;
0058     end    
0059     DataOut=DataIn; %reporduce global attribute
0060     ListVarName={};
0061     ValueList={};
0062     RoleList={};
0063     units_cell={};
0064     for ilist=1:length(FieldName)  
0065         [VarName,Value,Role,units]=feval(FieldName{ilist},DataIn);%calculate field with appropriate function named FieldName{ilist}
0066         ListVarName=[ListVarName VarName];
0067         ValueList=[ValueList Value];
0068         RoleList=[RoleList Role];
0069         units_cell=[units_cell units];
0070     end   
0071        %erase previous data (except coordinates)
0072     for ivar=nbcoord+1:length(DataOut.ListVarName)
0073         VarName=DataOut.ListVarName{ivar};
0074         DataOut=rmfield(DataOut,VarName);
0075     end  
0076     DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
0077     if isfield(DataOut,'VarDimName')
0078         DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
0079     else
0080         errormsg='element .VarDimName missing in input data';
0081         return
0082     end
0083     DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
0084     %append new data
0085     DataOut.ListVarName=[DataOut.ListVarName ListVarName];
0086     for ivar=1:length(ListVarName)
0087         DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};      
0088         DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};  
0089         DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
0090         eval(['DataOut.' ListVarName{ivar} '=ValueList{ivar};'])
0091     end
0092 end
0093 
0094 %%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
0095 function [VarName,ValCell,Role,units_cell]=velocity(DataIn)
0096 VarName={};
0097 ValCell={};
0098 Role={};
0099 units_cell={};
0100 if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
0101     units=[DataIn.CoordUnit '/' DataIn.TimeUnit];
0102 else
0103     units='pixel';
0104 end
0105 if isfield(DataIn,'U')
0106     VarName=[VarName {'U'}];
0107     ValCell=[ValCell {DataIn.U}];
0108     Role=[Role {'vector_x'}];
0109     units_cell=[units_cell {units}];
0110 end
0111 if isfield(DataIn,'V')
0112     VarName=[VarName {'V'}];
0113     ValCell=[ValCell {DataIn.V}];
0114     Role=[Role {'vector_y'}];  
0115     units_cell=[units_cell {units}];
0116 end
0117 if isfield(DataIn,'W')
0118     VarName=[VarName {'W'}];
0119     ValCell=[ValCell {DataIn.W}];
0120     Role=[Role {'vector_z'}];
0121     units_cell=[units_cell {units}];
0122 end
0123 if isfield(DataIn,'F')
0124     VarName=[VarName {'F'}];
0125     ValCell=[ValCell {DataIn.F}];
0126     Role=[Role {'warnflag'}];
0127     units_cell=[units_cell {[]}];
0128 end
0129 if isfield(DataIn,'FF')
0130     VarName=[VarName,{'FF'}];
0131     ValCell=[ValCell {DataIn.FF}];
0132     Role=[Role {'errorflag'}];
0133     units_cell=[units_cell {[]}];
0134 end
0135 
0136 %%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%
0137 function [VarName,ValCell,Role,units]=ima_cor(DataIn)
0138 VarName={};
0139 ValCell={};
0140 Role={};
0141 units={};
0142 if isfield(DataIn,'C')
0143     VarName{1}='C';
0144     ValCell{1}=DataIn.C;
0145     Role={'ancillary'};
0146     units={[]};
0147 end
0148 
0149 %%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%
0150 function [VarName,ValCell,Role,units]=norm_vel(DataIn)
0151 VarName={};
0152 ValCell={};
0153 Role={};
0154 units={};
0155 if isfield(DataIn,'U') && isfield(DataIn,'V')
0156     VarName{1}='norm_vel';
0157      ValCell{1}=DataIn.U.*DataIn.U+ DataIn.V.*DataIn.V;
0158      if isfield(DataIn,'W') && isequal(size(DataIn.W),size(DataIn.U))
0159          ValCell{1}=ValCell{1}+DataIn.W.*DataIn.W;
0160      end
0161      ValCell{1}=sqrt(ValCell{1});
0162      Role{1}='scalar';
0163      if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
0164         units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
0165      else
0166         units={'pixel'};
0167      end
0168 end  
0169 
0170 
0171 
0172 %%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%
0173 function [VarName,ValCell,Role,units]=vort(DataIn)
0174 VarName={};
0175 ValCell={};
0176 Role={};
0177 units={};
0178 if isfield(DataIn,'DjUi')
0179     VarName{1}='vort';
0180     ValCell{1}=DataIn.DjUi(:,1,2)-DataIn.DjUi(:,2,1);  %vorticity
0181     siz=size(ValCell{1});
0182     ValCell{1}=reshape(ValCell{1},siz(1),1);
0183     Role{1}='scalar';
0184     if isfield(DataIn,'TimeUnit')
0185         units={[DataIn.TimeUnit '-1']};
0186     else
0187         units={[]};
0188     end
0189 end  
0190 
0191 %%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%
0192 function [VarName,ValCell,Role,units]=div(DataIn)
0193 VarName={};
0194 ValCell={};
0195 Role={};
0196 units={};
0197 if isfield(DataIn,'DjUi')
0198     VarName{1}='div';
0199     ValCell{1}=DataIn.DjUi(:,1,1)+DataIn.DjUi(:,2,2); %DUDX+DVDY
0200     siz=size(ValCell{1});
0201     ValCell{1}=reshape(ValCell{1},siz(1),1);
0202     Role{1}='scalar';
0203     if isfield(DataIn,'TimeUnit')
0204         units={[DataIn.TimeUnit '-1']};
0205     else
0206         units={[]};
0207     end
0208 end  
0209 
0210 %%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
0211 function [VarName,ValCell,Role,units]=strain(DataIn)
0212 VarName={};
0213 ValCell={};
0214 Role={};
0215 units={};
0216 if isfield(DataIn,'DjUi')
0217    VarName{1}='strain'; 
0218    ValCell{1}=DataIn.DjUi(:,1,2)+DataIn.DjUi(:,2,1);%DVDX+DUDY
0219    siz=size(ValCell{1});    
0220    ValCell{1}=reshape(ValCell{1},siz(1),1); 
0221    if isfield(DataIn,'TimeUnit')
0222         units={[DataIn.TimeUnit '-1']};
0223    else
0224         units={[]};
0225    end
0226 end  
0227 
0228 %%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
0229 function [VarName,ValCell,Role,units]=u(DataIn)
0230 VarName={};
0231 ValCell={};
0232 Role={};
0233 units={};
0234 if isfield(DataIn,'U')
0235     VarName{1}='U';
0236     ValCell{1}=DataIn.U;
0237     Role{1}='scalar';
0238     if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
0239         units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
0240     else
0241         units={'pixel'};
0242     end
0243 end
0244 
0245 %%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
0246 function [VarName,ValCell,Role,units]=v(DataIn)
0247 VarName={};
0248 ValCell={};
0249 Role={};
0250 units={};
0251 if isfield(DataIn,'V')
0252     VarName{1}='V';
0253     ValCell{1}=DataIn.V;
0254     Role{1}='scalar';
0255     if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
0256         units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
0257     else
0258         units={'pixel'};
0259     end
0260 end
0261 
0262 %%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
0263 function [VarName,ValCell,Role,units]=w(DataIn)
0264 VarName={};
0265 ValCell={};
0266 Role={};
0267 units={};
0268 if isfield(DataIn,'W')
0269     VarName{1}='W';
0270     ValCell{1}=DataIn.W;
0271     Role{1}='scalar';%will remain unchanged by projection
0272     if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
0273         units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
0274     else
0275         units={'pixel'};
0276     end
0277 end
0278 
0279 %%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
0280 function [VarName,ValCell,Role,units]=w_normal(DataIn)
0281 VarName={};
0282 ValCell={};
0283 Role={};
0284 units={};
0285 if isfield(DataIn,'W')
0286     VarName{1}='W';
0287     ValCell{1}=DataIn.W;
0288     Role{1}='vector_z';%will behave like a vector component  by projection
0289     if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
0290         units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
0291     else
0292         units={'pixel'};
0293     end
0294 end
0295 
0296 %%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
0297 function [VarName,ValCell,Role,units]=error(DataIn)
0298 VarName={};
0299 ValCell={};
0300 Role={};
0301 units={};
0302 if isfield(DataIn,'E')
0303     VarName{1}='E';
0304     ValCell{1}=DataIn.E;
0305     Role{1}='ancillary'; %TODO CHECK units in actual fields
0306     if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
0307         units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
0308     else
0309         units={'pixel'};
0310     end
0311 end
0312