Changeset 1157 for trunk/src/series/sliding_average.m

Timestamp:

Jul 11, 2024, 4:13:03 PM (4 months ago)

Author:

sommeria

Message:

filter_time added

File:

• trunk/src/series/sliding_average.m (modified) (6 diffs)

trunk/src/series/sliding_average.m

@@ -1 @@
-%'sliding_average_natalya'
+%'sliding_average': calculate the time correlation function at each point
 %------------------------------------------------------------------------
-% function ParamOut=sliding_average(Param)
+% function ParamOut=turb_correlation_time(Param)
 %
 %%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%
@@ -39 +39 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %=======================================================================
-% Copyright 2008-2024, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
+% Copyright 2008-2022, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
 %   http://www.legi.grenoble-inp.fr
-%   Joel.Sommeria - Joel.Sommeria (A) univ-grenoble-alpes.fr
+%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
 %
 %     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)
     ParamOut.FieldName='one';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
-    ParamOut.FieldTransform = 'on';%can use a transform function
+    ParamOut.FieldTransform = 'off';%can use a transform function
     ParamOut.ProjObject='off';%can use projection object39(option 'off'/'on',
     ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
     ParamOut.OutputDirExt='.tfilter';%set the output dir extension
-    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
-%     filecell=get_file_series(Param);%check existence of the first input file
-%     if ~exist(filecell{1,1},'file')
-%         msgbox_uvmat('WARNING','the first input file does not exist')
-%     end
+    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
+   
+    %% check the first input file in the series
+     first_j=[];% 
+    if isfield(Param.IndexRange,'first_j'); first_j=Param.IndexRange.first_j; end
+    PairString='';
+    if isfield(Param.IndexRange,'PairString'); PairString=Param.IndexRange.PairString; end
+    [i1,i2,j1,j2] = get_file_index(Param.IndexRange.first_i,first_j,PairString);
+    FirstFileName=fullfile_uvmat(Param.InputTable{1,1},Param.InputTable{1,2},Param.InputTable{1,3},...
+        Param.InputTable{1,5},Param.InputTable{1,4},i1,i2,j1,j2); 
+    if exist(FirstFileName,'file')
+        FileInfo=get_file_info(FirstFileName);
+        if ~(isfield(FileInfo,'FileType')&& strcmp(FileInfo.FileType,'netcdf'))
+            msgbox_uvmat('ERROR','this fct works only on netcdf files (fields projected on a grid, notraw civ data)')
+        return
+        end
+    else
+        msgbox_uvmat('ERROR',['the input file ' FirstFileName ' does not exist'])
+        return
+    end
+
+    %% setting of intput parameters 
+%     answer = msgbox_uvmat('INPUT_MENU','select the variables to filter',ListVarName');
+%     ParamOut.ListVarName=ListVarName(answer);
+    answer = msgbox_uvmat('INPUT_TXT','set the filering window length','3');
+    ParamOut.WindowLength=str2double(answer);
+    
+%         [ParamOut.ActionInput,errormsg] = set_param_input(ListParam,DefaultValue,ParamIn);
+%         if ~isempty(errormsg)
+%             msgbox_uvmat('ERROR',errormsg)
+%         end
+%     [ParamOut,errormsg] = set_param_input(ListParam,DefaultValue,ParamIn);
     return
 end
@@ -184 +211 @@
 nbmissing=0;
 
-%% initialisation manip Coriolis
-char_index=regexp(SubDir{1},'waves_L1_');
-switch(SubDir{1}(char_index+9))
-    case '1'
-        amplitude=2.5 %oscillation amplitude
-        T=24.46;
-        t0=3 ;% dt=0.5 s, torus at its max x at the beginning of motion, i0=7
-    case '2'
-        amplitude=5 %oscillation amplitude
-        T=24.47;
-        t0=8.5; % dt=1/3 s -> image index of starting motion = 26, % torus at its max x at the beginning of motion
-    case '3'
-        amplitude=10 %oscillation amplitude
-        T=24.45;
-        t0=6.5-T/2;% dt=0.25, torus at its minimum x at the beginning of motion
-    case '4'  
-        amplitude=15 %oscillation amplitude
-        T=24.48;
-        t0=3.4;     %dt=0.2 -> i0=18 image index of starting motion, % torus at its max x at the beginning of motion
-end
-%NbPeriod=2; %number of periods for the sliding average
-NbPeriod=4; %number of periods for the sliding average
-omega=2*2*pi/T;%harmonic
-Lscale=15;%diameter of the torus, length scale for normalisation
-Uscale=amplitude*omega;
-
-DataOut.ListGlobalAttribute= {'Conventions','Time'};
+%% initialisation output data
+
+DataOut.ListGlobalAttribute= {'Conventions'};
 DataOut.Conventions='uvmat';
-DataOut.ListVarName={'coord_y','coord_x','Umean','Vmean','Ucos','Vcos','Usin','Vsin','DUDXsin','DUDXcos','DUDYsin','DVDXsin','DVDXcos'...
-    ,'DVDYsin','Ustokes','Vstokes'};
-DataOut.VarDimName={'coord_y','coord_x',{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},...
-    {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},...
-    {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};
+DataOut.ListVarName={'Time','coord_y','coord_x','Ufilter','Vfilter'};
+DataOut.VarDimName={'Time','coord_y','coord_x',{'Time','coord_y','coord_x'},{'Time','coord_y','coord_x'}};
 
 %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
@@ -226 +226 @@
     return
 end
-[Data,tild,errormsg]=nc2struct(filecell{1,end});    
-Time_end=Data.Time;
-dt=(Time_end-Time_1)/(NbField-1); %time interval 
-NpTime=round(NbPeriod*T/dt+1);
+[Data,tild,errormsg]=nc2struct(filecell{1,2});    
+Time_2=Data.Time;
+dt=(Time_2-Time_1); %time interval 
+NpTime=21; %Nbre de champ pour la moyenne glissante (choisir impair)
 
 OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device));
@@ -246 +246 @@
     %%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
     if index==1 %first field
-        DataOut.coord_x=Field.coord_x/Lscale;
-        DataOut.coord_y=Field.coord_y/Lscale;
+        DataOut.coord_x=Field.coord_x;
+        DataOut.coord_y=Field.coord_y;
         npy=numel(DataOut.coord_y);
         npx=numel(DataOut.coord_x);
-        Umean=zeros(NpTime,npy,npx);
-        Vmean=zeros(NpTime,npy,npx);
-        Ucos=zeros(NpTime,npy,npx);
-        Vcos=zeros(NpTime,npy,npx);
-        Usin=zeros(NpTime,npy,npx);
-        Vsin=zeros(NpTime,npy,npx);
-        DUDXcos=zeros(NpTime,npy,npx);
-        DUDXsin=zeros(NpTime,npy,npx);
-        DUDYsin=zeros(NpTime,npy,npx);
-        DVDXcos=zeros(NpTime,npy,npx);
-        DVDXsin=zeros(NpTime,npy,npx);
-        DVDYsin=zeros(NpTime,npy,npx);
-    end
-    Time(index)=Field.Time-t0;%time from the start of the motion
-    Umean=circshift(Umean,[-1 0 0]); %shift U by ishift along the first index
-    Vmean=circshift(Vmean,[-1 0 0]); %shift U by ishift along the first index
-    Ucos=circshift(Ucos,[-1 0 0]); %shift U by ishift along the first index
-    Vcos=circshift(Vcos,[-1 0 0]); %shift U by ishift along the first index
-    Usin=circshift(Usin,[-1 0 0]); %shift U by ishift along the first index
-    Vsin=circshift(Vsin,[-1 0 0]); %shift U by ishift along the first index
-    DUDXcos=circshift(DUDXcos,[-1 0 0]);
-    DUDXsin=circshift(DUDXsin,[-1 0 0]);
-    DUDYsin=circshift(DUDYsin,[-1 0 0]);        
-    DVDXcos=circshift(DVDXcos,[-1 0 0]);
-    DVDXsin=circshift(DVDXsin,[-1 0 0]);
-    DVDYsin=circshift(DVDYsin,[-1 0 0]);       
-    Umean(end,:,:)=Field.U;
-    Vmean(end,:,:)=Field.V;
-    Ucos(end,:,:)=Field.U*cos(omega*Time(index));
-    Vcos(end,:,:)=Field.V*cos(omega*Time(index));
-    Usin(end,:,:)=Field.U*sin(omega*Time(index));
-    Vsin(end,:,:)=Field.V*sin(omega*Time(index));
-    DUDXcos(end,:,:)=Field.DUDX*cos(omega*Time(index));
-    DUDXsin(end,:,:)=Field.DUDX*sin(omega*Time(index));
-    DUDYsin(end,:,:)=Field.DUDY*sin(omega*Time(index));% ParamOut=[];%default output
-
-    DVDXcos(end,:,:)=Field.DVDX*cos(omega*Time(index));
-    DVDXsin(end,:,:)=Field.DVDX*sin(omega*Time(index));
-    DVDYsin(end,:,:)=Field.DVDY*sin(omega*Time(index));
-    DataOut.Time=(Time(index)-(NpTime-1)*dt/2)/T;%time inperiods from the beginning of the oscillation (torus at max abscissa)
-    DataOut.Umean=(1/Uscale)*squeeze(nanmean(Umean,1));
-    DataOut.Vmean=(1/Uscale)*squeeze(nanmean(Vmean,1));
-    DataOut.Ucos=2*(1/Uscale)*squeeze(nanmean(Ucos,1));
-    DataOut.Vcos=2*(1/Uscale)*squeeze(nanmean(Vcos,1));
-    DataOut.Usin=2*(1/Uscale)*squeeze(nanmean(Usin,1));
-    DataOut.Vsin=2*(1/Uscale)*squeeze(nanmean(Vsin,1));
-    DataOut.DUDXcos=2*squeeze(nanmean(DUDXcos,1));
-    DataOut.DUDXsin=2*squeeze(nanmean(DUDXsin,1));
-    DataOut.DUDYsin=2*squeeze(nanmean(DUDYsin,1));
-    DataOut.DVDXcos=2*squeeze(nanmean(DVDXcos,1));
-    DataOut.DVDXsin=2*squeeze(nanmean(DVDXsin,1));
-    DataOut.DVDYsin=2*squeeze(nanmean(DVDYsin,1));
-    DataOut.Ustokes=(1/omega)*(1/Uscale)*(DataOut.Ucos.*DataOut.DUDXsin+DataOut.Vcos.*DataOut.DUDYsin);
-    DataOut.Vstokes=(1/omega)*(1/Uscale)*(DataOut.Ucos.*DataOut.DVDXsin+DataOut.Vcos.*DataOut.DVDYsin);
+        Ufilter=zeros(NpTime,npy,npx);
+        Vfilter=zeros(NpTime,npy,npx);
+    end
+    Time(index)=Field.Time;%time 
+    Ufilter=circshift(Ufilter,[-1 0 0]); %shift U by ishift along the first index
+    Vfilter=circshift(Vfilter,[-1 0 0]); %shift U by ishift along the first index    
+    Ufilter(end,:,:)=Field.U;
+    Vfilter(end,:,:)=Field.V;
+    DataOut.Time=(Time(index)-(NpTime-1)*dt/2);%mid time 
+    DataOut.Ufilter=squeeze(nanmean(Ufilter,1));
+    DataOut.Vfilter=squeeze(nanmean(Vfilter,1));
+    
 
     % writing the result file as netcdf file
-    i1=i1_series{1}(index);
+    i1=i1_series{1}(index)-ceil(NpTime/2);
     OutputFile=fullfile_uvmat(OutputPath,OutputDir,RootFileOut,'.nc',NomTypeOut,i1);
     errormsg=struct2nc(OutputFile, DataOut);