Ignore:
Timestamp:
Jul 11, 2024, 4:13:03 PM (3 months ago)
Author:
sommeria
Message:

filter_time added

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
    22%------------------------------------------------------------------------
    3 % function ParamOut=sliding_average(Param)
     3% function ParamOut=turb_correlation_time(Param)
    44%
    55%%%%%%%%%%% GENERAL TO ALL SERIES ACTION FCTS %%%%%%%%%%%%%%%%%%%%%%%%%%%
     
    3939%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    4040%=======================================================================
    41 % Copyright 2008-2024, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
     41% Copyright 2008-2022, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
    4242%   http://www.legi.grenoble-inp.fr
    43 %   Joel.Sommeria - Joel.Sommeria (A) univ-grenoble-alpes.fr
     43%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
    4444%
    4545%     This file is part of the toolbox UVMAT.
     
    6565    ParamOut.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
    6666    ParamOut.FieldName='one';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
    67     ParamOut.FieldTransform = 'on';%can use a transform function
     67    ParamOut.FieldTransform = 'off';%can use a transform function
    6868    ParamOut.ProjObject='off';%can use projection object39(option 'off'/'on',
    6969    ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
    7070    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);
    76103    return
    77104end
     
    184211nbmissing=0;
    185212
    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
     215DataOut.ListGlobalAttribute= {'Conventions'};
    213216DataOut.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'}};
     217DataOut.ListVarName={'Time','coord_y','coord_x','Ufilter','Vfilter'};
     218DataOut.VarDimName={'Time','coord_y','coord_x',{'Time','coord_y','coord_x'},{'Time','coord_y','coord_x'}};
    219219
    220220%%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
     
    226226    return
    227227end
    228 [Data,tild,errormsg]=nc2struct(filecell{1,end});   
    229 Time_end=Data.Time;
    230 dt=(Time_end-Time_1)/(NbField-1); %time interval
    231 NpTime=round(NbPeriod*T/dt+1);
     228[Data,tild,errormsg]=nc2struct(filecell{1,2});   
     229Time_2=Data.Time;
     230dt=(Time_2-Time_1); %time interval
     231NpTime=21; %Nbre de champ pour la moyenne glissante (choisir impair)
    232232
    233233OutputPath=fullfile(Param.OutputPath,num2str(Param.Experiment),num2str(Param.Device));
     
    246246    %%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
    247247    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;
    250250        npy=numel(DataOut.coord_y);
    251251        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   
    306264
    307265    % writing the result file as netcdf file
    308     i1=i1_series{1}(index);
     266    i1=i1_series{1}(index)-ceil(NpTime/2);
    309267    OutputFile=fullfile_uvmat(OutputPath,OutputDir,RootFileOut,'.nc',NomTypeOut,i1);
    310268    errormsg=struct2nc(OutputFile, DataOut);
Note: See TracChangeset for help on using the changeset viewer.