Changeset 1157 for trunk/src/series/sliding_average.m
- Timestamp:
- Jul 11, 2024, 4:13:03 PM (3 months ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/series/sliding_average.m
r1127 r1157 1 %'sliding_average _natalya'1 %'sliding_average': calculate the time correlation function at each point 2 2 %------------------------------------------------------------------------ 3 % function ParamOut= sliding_average(Param)3 % function ParamOut=turb_correlation_time(Param) 4 4 % 5 5 %%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%% … … 39 39 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 40 40 %======================================================================= 41 % Copyright 2008-202 4, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France41 % Copyright 2008-2022, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France 42 42 % http://www.legi.grenoble-inp.fr 43 % Joel.Sommeria - Joel.Sommeria (A) univ-grenoble-alpes.fr43 % Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr 44 44 % 45 45 % This file is part of the toolbox UVMAT. … … 65 65 ParamOut.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two', 'off' by default) 66 66 ParamOut.FieldName='one';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default) 67 ParamOut.FieldTransform = 'o n';%can use a transform function67 ParamOut.FieldTransform = 'off';%can use a transform function 68 68 ParamOut.ProjObject='off';%can use projection object39(option 'off'/'on', 69 69 ParamOut.Mask='off';%can use mask option (option 'off'/'on', 'off' by default) 70 70 ParamOut.OutputDirExt='.tfilter';%set the output dir extension 71 ParamOut.OutputFileMode='NbSlice';% '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice 72 % filecell=get_file_series(Param);%check existence of the first input file 73 % if ~exist(filecell{1,1},'file') 74 % msgbox_uvmat('WARNING','the first input file does not exist') 75 % end 71 ParamOut.OutputFileMode='NbSlice'; % '=NbInput': 1 output file per input file index, '=NbInput_i': 1 file per input file index i, '=NbSlice': 1 file per slice 72 73 %% check the first input file in the series 74 first_j=[];% 75 if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end 76 PairString=''; 77 if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end 78 [i1,i2,j1,j2] = get_file_index(Param.IndexRange.first_i,first_j,PairString); 79 FirstFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},... 80 Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2); 81 if exist(FirstFileName,'file') 82 FileInfo=get_file_info(FirstFileName); 83 if ~(isfield(FileInfo,'FileType')&& strcmp(FileInfo.FileType,'netcdf')) 84 msgbox_uvmat('ERROR','this fct works only on netcdf files (fields projected on a grid, notraw civ data)') 85 return 86 end 87 else 88 msgbox_uvmat('ERROR',['the input file ' FirstFileName ' does not exist']) 89 return 90 end 91 92 %% setting of intput parameters 93 % answer = msgbox_uvmat('INPUT_MENU','select the variables to filter',ListVarName'); 94 % ParamOut.ListVarName=ListVarName(answer); 95 answer = msgbox_uvmat('INPUT_TXT','set the filering window length','3'); 96 ParamOut.WindowLength=str2double(answer); 97 98 % [ParamOut.ActionInput,errormsg] = set_param_input(ListParam,DefaultValue,ParamIn); 99 % if ~isempty(errormsg) 100 % msgbox_uvmat('ERROR',errormsg) 101 % end 102 % [ParamOut,errormsg] = set_param_input(ListParam,DefaultValue,ParamIn); 76 103 return 77 104 end … … 184 211 nbmissing=0; 185 212 186 %% initialisation manip Coriolis 187 char_index=regexp(SubDir{1},'waves_L1_'); 188 switch(SubDir{1}(char_index+9)) 189 case '1' 190 amplitude=2.5 %oscillation amplitude 191 T=24.46; 192 t0=3 ;% dt=0.5 s, torus at its max x at the beginning of motion, i0=7 193 case '2' 194 amplitude=5 %oscillation amplitude 195 T=24.47; 196 t0=8.5; % dt=1/3 s -> image index of starting motion = 26, % torus at its max x at the beginning of motion 197 case '3' 198 amplitude=10 %oscillation amplitude 199 T=24.45; 200 t0=6.5-T/2;% dt=0.25, torus at its minimum x at the beginning of motion 201 case '4' 202 amplitude=15 %oscillation amplitude 203 T=24.48; 204 t0=3.4; %dt=0.2 -> i0=18 image index of starting motion, % torus at its max x at the beginning of motion 205 end 206 %NbPeriod=2; %number of periods for the sliding average 207 NbPeriod=4; %number of periods for the sliding average 208 omega=2*2*pi/T;%harmonic 209 Lscale=15;%diameter of the torus, length scale for normalisation 210 Uscale=amplitude*omega; 211 212 DataOut.ListGlobalAttribute= {'Conventions','Time'}; 213 %% initialisation output data 214 215 DataOut.ListGlobalAttribute= {'Conventions'}; 213 216 DataOut.Conventions='uvmat'; 214 DataOut.ListVarName={'coord_y','coord_x','Umean','Vmean','Ucos','Vcos','Usin','Vsin','DUDXsin','DUDXcos','DUDYsin','DVDXsin','DVDXcos'... 215 ,'DVDYsin','Ustokes','Vstokes'}; 216 DataOut.VarDimName={'coord_y','coord_x',{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},... 217 {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},... 218 {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}}; 217 DataOut.ListVarName={'Time','coord_y','coord_x','Ufilter','Vfilter'}; 218 DataOut.VarDimName={'Time','coord_y','coord_x',{'Time','coord_y','coord_x'},{'Time','coord_y','coord_x'}}; 219 219 220 220 %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%% … … 226 226 return 227 227 end 228 [Data,tild,errormsg]=nc2struct(filecell{1, end});229 Time_ end=Data.Time;230 dt=(Time_ end-Time_1)/(NbField-1); %time interval231 NpTime= round(NbPeriod*T/dt+1);228 [Data,tild,errormsg]=nc2struct(filecell{1,2}); 229 Time_2=Data.Time; 230 dt=(Time_2-Time_1); %time interval 231 NpTime=21; %Nbre de champ pour la moyenne glissante (choisir impair) 232 232 233 233 OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device)); … … 246 246 %%%%%%%%%%% MAIN RUNNING OPERATIONS %%%%%%%%%%%% 247 247 if index==1 %first field 248 DataOut.coord_x=Field.coord_x /Lscale;249 DataOut.coord_y=Field.coord_y /Lscale;248 DataOut.coord_x=Field.coord_x; 249 DataOut.coord_y=Field.coord_y; 250 250 npy=numel(DataOut.coord_y); 251 251 npx=numel(DataOut.coord_x); 252 Umean=zeros(NpTime,npy,npx); 253 Vmean=zeros(NpTime,npy,npx); 254 Ucos=zeros(NpTime,npy,npx); 255 Vcos=zeros(NpTime,npy,npx); 256 Usin=zeros(NpTime,npy,npx); 257 Vsin=zeros(NpTime,npy,npx); 258 DUDXcos=zeros(NpTime,npy,npx); 259 DUDXsin=zeros(NpTime,npy,npx); 260 DUDYsin=zeros(NpTime,npy,npx); 261 DVDXcos=zeros(NpTime,npy,npx); 262 DVDXsin=zeros(NpTime,npy,npx); 263 DVDYsin=zeros(NpTime,npy,npx); 264 end 265 Time(index)=Field.Time-t0;%time from the start of the motion 266 Umean=circshift(Umean,[-1 0 0]); %shift U by ishift along the first index 267 Vmean=circshift(Vmean,[-1 0 0]); %shift U by ishift along the first index 268 Ucos=circshift(Ucos,[-1 0 0]); %shift U by ishift along the first index 269 Vcos=circshift(Vcos,[-1 0 0]); %shift U by ishift along the first index 270 Usin=circshift(Usin,[-1 0 0]); %shift U by ishift along the first index 271 Vsin=circshift(Vsin,[-1 0 0]); %shift U by ishift along the first index 272 DUDXcos=circshift(DUDXcos,[-1 0 0]); 273 DUDXsin=circshift(DUDXsin,[-1 0 0]); 274 DUDYsin=circshift(DUDYsin,[-1 0 0]); 275 DVDXcos=circshift(DVDXcos,[-1 0 0]); 276 DVDXsin=circshift(DVDXsin,[-1 0 0]); 277 DVDYsin=circshift(DVDYsin,[-1 0 0]); 278 Umean(end,:,:)=Field.U; 279 Vmean(end,:,:)=Field.V; 280 Ucos(end,:,:)=Field.U*cos(omega*Time(index)); 281 Vcos(end,:,:)=Field.V*cos(omega*Time(index)); 282 Usin(end,:,:)=Field.U*sin(omega*Time(index)); 283 Vsin(end,:,:)=Field.V*sin(omega*Time(index)); 284 DUDXcos(end,:,:)=Field.DUDX*cos(omega*Time(index)); 285 DUDXsin(end,:,:)=Field.DUDX*sin(omega*Time(index)); 286 DUDYsin(end,:,:)=Field.DUDY*sin(omega*Time(index));% ParamOut=[];%default output 287 288 DVDXcos(end,:,:)=Field.DVDX*cos(omega*Time(index)); 289 DVDXsin(end,:,:)=Field.DVDX*sin(omega*Time(index)); 290 DVDYsin(end,:,:)=Field.DVDY*sin(omega*Time(index)); 291 DataOut.Time=(Time(index)-(NpTime-1)*dt/2)/T;%time inperiods from the beginning of the oscillation (torus at max abscissa) 292 DataOut.Umean=(1/Uscale)*squeeze(nanmean(Umean,1)); 293 DataOut.Vmean=(1/Uscale)*squeeze(nanmean(Vmean,1)); 294 DataOut.Ucos=2*(1/Uscale)*squeeze(nanmean(Ucos,1)); 295 DataOut.Vcos=2*(1/Uscale)*squeeze(nanmean(Vcos,1)); 296 DataOut.Usin=2*(1/Uscale)*squeeze(nanmean(Usin,1)); 297 DataOut.Vsin=2*(1/Uscale)*squeeze(nanmean(Vsin,1)); 298 DataOut.DUDXcos=2*squeeze(nanmean(DUDXcos,1)); 299 DataOut.DUDXsin=2*squeeze(nanmean(DUDXsin,1)); 300 DataOut.DUDYsin=2*squeeze(nanmean(DUDYsin,1)); 301 DataOut.DVDXcos=2*squeeze(nanmean(DVDXcos,1)); 302 DataOut.DVDXsin=2*squeeze(nanmean(DVDXsin,1)); 303 DataOut.DVDYsin=2*squeeze(nanmean(DVDYsin,1)); 304 DataOut.Ustokes=(1/omega)*(1/Uscale)*(DataOut.Ucos.*DataOut.DUDXsin+DataOut.Vcos.*DataOut.DUDYsin); 305 DataOut.Vstokes=(1/omega)*(1/Uscale)*(DataOut.Ucos.*DataOut.DVDXsin+DataOut.Vcos.*DataOut.DVDYsin); 252 Ufilter=zeros(NpTime,npy,npx); 253 Vfilter=zeros(NpTime,npy,npx); 254 end 255 Time(index)=Field.Time;%time 256 Ufilter=circshift(Ufilter,[-1 0 0]); %shift U by ishift along the first index 257 Vfilter=circshift(Vfilter,[-1 0 0]); %shift U by ishift along the first index 258 Ufilter(end,:,:)=Field.U; 259 Vfilter(end,:,:)=Field.V; 260 DataOut.Time=(Time(index)-(NpTime-1)*dt/2);%mid time 261 DataOut.Ufilter=squeeze(nanmean(Ufilter,1)); 262 DataOut.Vfilter=squeeze(nanmean(Vfilter,1)); 263 306 264 307 265 % writing the result file as netcdf file 308 i1=i1_series{1}(index) ;266 i1=i1_series{1}(index)-ceil(NpTime/2); 309 267 OutputFile=fullfile_uvmat(OutputPath,OutputDir,RootFileOut,'.nc',NomTypeOut,i1); 310 268 errormsg=struct2nc(OutputFile, DataOut);
Note: See TracChangeset
for help on using the changeset viewer.