Changeset 1100 for trunk/src/series/turb_stat.m
 Timestamp:
 Apr 23, 2021, 4:59:15 PM (3 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/src/series/turb_stat.m
r1095 r1100 202 202 203 203 %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%% 204 for index=1:NbField 205 update_waitbar(WaitbarHandle,index/NbField) 206 if ~isempty(RUNHandle)&& ~strcmp(get(RUNHandle,'BusyAction'),'queue') 207 disp('program stopped by user') 208 break 209 end 210 [Field,tild,errormsg] = read_field(filecell{1,index},FileType{iview},InputFields{iview},frame_index{iview}(index)); 211 212 %%%%%%%%%%%% MAIN RUNNING OPERATIONS %%%%%%%%%%%% 213 if index==1 %initiate the output data structure in the first field 214 [CellInfo,NbDim,errormsg]=find_field_cells(Field); 215 YName='coord_y';%default 216 XName='coord_x';%default 217 for icell=1:numel(NbDim) 218 if NbDim(icell)==2 && strcmp(CellInfo{icell}.CoordType,'grid') 219 YName=CellInfo{icell}.YName; 220 XName=CellInfo{icell}.XName; 221 break 204 % for i_slice=1:Param.IndexRange.NbSlice 205 % i_slice 206 ind_first=Param.IndexRange.first_i 207 for index_i=ind_first:Param.IndexRange.NbSlice:Param.IndexRange.last_i 208 if ~isempty(RUNHandle)&& ~strcmp(get(RUNHandle,'BusyAction'),'queue') 209 disp('program stopped by user') 210 break 211 end 212 for index_j=Param.IndexRange.first_j:Param.IndexRange.last_j 213 InputFile=fullfile_uvmat(RootPath{1},SubDir{1},RootFile{1},FileExt{1},NomType{1},index_i,index_i,index_j,index_j); 214 [Field,tild,errormsg] = read_field(InputFile,FileType{iview},InputFields{iview}); 215 216 %[Field,tild,errormsg] = read_field(filecell{1,index},FileType{iview},InputFields{iview},frame_index{iview}(index)); 217 218 %%%%%%%%%%%% MAIN RUNNING OPERATIONS %%%%%%%%%%%% 219 if index_i==ind_first && index_j==Param.IndexRange.first_j %initiate the output data structure in the first field 220 [CellInfo,NbDim,errormsg]=find_field_cells(Field); 221 YName='coord_y';%default 222 XName='coord_x';%default 223 for icell=1:numel(NbDim) 224 if NbDim(icell)==2 && strcmp(CellInfo{icell}.CoordType,'grid') 225 YName=CellInfo{icell}.YName; 226 XName=CellInfo{icell}.XName; 227 break 228 end 229 end 230 DataOut.ListVarName={YName, XName ,'UMean' , 'VMean','u2Mean','v2Mean','u2Mean_1','v2Mean_1','uvMean','CurlMean','DivMean','Curl2Mean','Div2Mean','Counter'}; 231 DataOut.VarDimName={YName,XName,{YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName},... 232 {YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName}}; 233 DataOut.(YName)=Field.(YName); 234 DataOut.(XName)=Field.(XName); 235 Uprev=Field.U;% store the current field for next iteration 236 Vprev=Field.V; 237 if isfield(Field,'FF') 238 FFprev=Field.FF;% possible flag for false data 239 else 240 %FFprev=true(size(Field.U)); 241 FFprev=isnan(Field.U); 242 end 243 end 244 FF=isnan(Field.U);%Field.U<60Field.U>30;% threshold on U 245 DataOut.Counter=DataOut.Counter+ (~FF);% add 1 to the couter for non NaN point 246 Counter_1=Counter_1+(~FF & ~FFprev); 247 Field.U(FF)=0;% set to 0 the nan values 248 Field.V(FF)=0; 249 DataOut.UMean=DataOut.UMean+Field.U; %increment the sum 250 DataOut.VMean=DataOut.VMean+Field.V; %increment the sum 251 252 U2Mean=U2Mean+(Field.U).*(Field.U); %increment the U squared sum 253 V2Mean=V2Mean+(Field.V).*(Field.V); %increment the V squared sum 254 UVMean=UVMean+(Field.U).*(Field.V); %increment the sum 255 U2Mean_1=U2Mean_1+(Field.U).*Uprev; %increment the U squared sum 256 V2Mean_1=V2Mean_1+(Field.V).*Vprev; %increment the V squared sum 257 Uprev=Field.U; %store for next iteration 258 Vprev=Field.V; 259 FFprev=FF; 260 if isfield(Field,'curl') && isfield(Field,'div') 261 Field.curl(FF)=0;% set to 0 the nan values 262 Field.div(FF)=0; 263 DataOut.CurlMean=DataOut.CurlMean+Field.curl; 264 DataOut.DivMean=DataOut.DivMean+Field.div; 265 DataOut.Curl2Mean=DataOut.Curl2Mean+Field.curl.*Field.curl; 266 DataOut.Div2Mean=DataOut.Div2Mean+Field.div.*Field.div; 222 267 end 223 268 end 224 DataOut.ListVarName={YName, XName ,'UMean' , 'VMean','u2Mean','v2Mean','u2Mean_1','v2Mean_1','uvMean','CurlMean','DivMean','Curl2Mean','Div2Mean','Counter'};225 DataOut.VarDimName={YName,XName,{YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName},...226 {YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName}};227 DataOut.(YName)=Field.(YName);228 DataOut.(XName)=Field.(XName);229 Uprev=Field.U;% store the current field for next iteration230 Vprev=Field.V;231 if isfield(Field,'FF')232 FFprev=Field.FF;% possible flag for false data233 else234 %FFprev=true(size(Field.U));235 FFprev=isnan(Field.U);236 end237 269 end 238 FF=isnan(Field.U);%Field.U<60Field.U>30;% threshold on U 239 DataOut.Counter=DataOut.Counter+ (~FF);% add 1 to the couter for non NaN point 240 Counter_1=Counter_1+(~FF & ~FFprev); 241 Field.U(FF)=0;% set to 0 the nan values 242 Field.V(FF)=0; 243 DataOut.UMean=DataOut.UMean+Field.U; %increment the sum 244 DataOut.VMean=DataOut.VMean+Field.V; %increment the sum 270 %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%% 245 271 246 U2Mean=U2Mean+(Field.U).*(Field.U); %increment the U squared sum 247 V2Mean=V2Mean+(Field.V).*(Field.V); %increment the V squared sum 248 UVMean=UVMean+(Field.U).*(Field.V); %increment the sum 249 U2Mean_1=U2Mean_1+(Field.U).*Uprev; %increment the U squared sum 250 V2Mean_1=V2Mean_1+(Field.V).*Vprev; %increment the V squared sum 251 Uprev=Field.U; %store for next iteration 252 Vprev=Field.V; 253 FFprev=FF; 254 if isfield(Field,'curl') && isfield(Field,'div') 255 Field.curl(FF)=0;% set to 0 the nan values 256 Field.div(FF)=0; 257 DataOut.CurlMean=DataOut.CurlMean+Field.curl; 258 DataOut.DivMean=DataOut.DivMean+Field.div; 259 DataOut.Curl2Mean=DataOut.Curl2Mean+Field.curl.*Field.curl; 260 DataOut.Div2Mean=DataOut.Div2Mean+Field.div.*Field.div; 272 DataOut.Counter(DataOut.Counter==0)=1;% put counter to 1 when it is zero 273 DataOut.UMean=DataOut.UMean./DataOut.Counter; % normalize the mean 274 DataOut.VMean=DataOut.VMean./DataOut.Counter; % normalize the mean 275 U2Mean=U2Mean./DataOut.Counter; % normalize the mean 276 V2Mean=V2Mean./DataOut.Counter; % normalize the mean 277 UVMean=UVMean./DataOut.Counter; % normalize the mean 278 U2Mean_1=U2Mean_1./Counter_1; % normalize the mean 279 V2Mean_1=V2Mean_1./Counter_1; % normalize the mean 280 DataOut.u2Mean=U2MeanDataOut.UMean.*DataOut.UMean; % normalize the mean 281 DataOut.v2Mean=V2MeanDataOut.VMean.*DataOut.VMean; % normalize the mean 282 DataOut.uvMean=UVMeanDataOut.UMean.*DataOut.VMean; % normalize the mean \ 283 DataOut.u2Mean_1=U2Mean_1DataOut.UMean.*DataOut.UMean; % normalize the mean 284 DataOut.v2Mean_1=V2Mean_1DataOut.VMean.*DataOut.VMean; % normalize the mean 285 DataOut.CurlMean=DataOut.CurlMean./DataOut.Counter; 286 DataOut.DivMean=DataOut.DivMean./DataOut.Counter; 287 DataOut.Curl2Mean=DataOut.Curl2Mean./DataOut.CounterDataOut.CurlMean.*DataOut.CurlMean; 288 DataOut.Div2Mean=DataOut.Div2Mean./DataOut.CounterDataOut.DivMean.*DataOut.DivMean; 289 290 %% writing the result file as netcdf file 291 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,ind_first,ind_first,first_j,last_j); 292 %case of netcdf input file , determine global attributes 293 errormsg=struct2nc(OutputFile,DataOut); %save result file 294 if isempty(errormsg) 295 disp([OutputFile ' written']); 296 else 297 disp(['error in writting result file: ' errormsg]) 261 298 end 262 end 263 %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%% 264 265 DataOut.Counter(DataOut.Counter==0)=1;% put counter to 1 when it is zero 266 DataOut.UMean=DataOut.UMean./DataOut.Counter; % normalize the mean 267 DataOut.VMean=DataOut.VMean./DataOut.Counter; % normalize the mean 268 U2Mean=U2Mean./DataOut.Counter; % normalize the mean 269 V2Mean=V2Mean./DataOut.Counter; % normalize the mean 270 UVMean=UVMean./DataOut.Counter; % normalize the mean 271 U2Mean_1=U2Mean_1./Counter_1; % normalize the mean 272 V2Mean_1=V2Mean_1./Counter_1; % normalize the mean 273 DataOut.u2Mean=U2MeanDataOut.UMean.*DataOut.UMean; % normalize the mean 274 DataOut.v2Mean=V2MeanDataOut.VMean.*DataOut.VMean; % normalize the mean 275 DataOut.uvMean=UVMeanDataOut.UMean.*DataOut.VMean; % normalize the mean \ 276 DataOut.u2Mean_1=U2Mean_1DataOut.UMean.*DataOut.UMean; % normalize the mean 277 DataOut.v2Mean_1=V2Mean_1DataOut.VMean.*DataOut.VMean; % normalize the mean 278 DataOut.CurlMean=DataOut.CurlMean./DataOut.Counter; 279 DataOut.DivMean=DataOut.DivMean./DataOut.Counter; 280 DataOut.Curl2Mean=DataOut.Curl2Mean./DataOut.CounterDataOut.CurlMean.*DataOut.CurlMean; 281 DataOut.Div2Mean=DataOut.Div2Mean./DataOut.CounterDataOut.DivMean.*DataOut.DivMean; 282 283 %% calculate the profiles 284 % npx=numel(DataOut.coord_x); 285 % band=ceil(npx/5) :floor(4*npx/5);% keep only the central band 286 % for ivar=3:numel(DataOut.ListVarName)1 287 % VarName=DataOut.ListVarName{ivar};% name of the variable 288 % DataOut.ListVarName=[DataOut.ListVarName {[VarName 'Profile']}];%append the name of the profile variable 289 % DataOut.VarDimName=[DataOut.VarDimName {'(YName)'}]; 290 % DataOut.([VarName 'Profile'])=mean(DataOut.(VarName)(:,band),2); %take the mean profile of U, excluding the edges 299 291 300 % end 292 293 %% writing the result file as netcdf file294 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j,last_j);295 %case of netcdf input file , determine global attributes296 errormsg=struct2nc(OutputFile,DataOut); %save result file297 if isempty(errormsg)298 disp([OutputFile ' written']);299 else300 disp(['error in writting result file: ' errormsg])301 end302 303 304 301 %% open the result file with uvmat (in RUN mode) 305 302 if checkrun && isequal(Param.IndexRange.NbSlice,1)
Note: See TracChangeset
for help on using the changeset viewer.