Changeset 959


Ignore:
Timestamp:
Jun 24, 2016, 8:49:08 PM (8 years ago)
Author:
sommeria
Message:

various

Location:
trunk/src
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/proj_field.m

    r955 r959  
    16531653                    Mrot=rodrigues(PlaneAngle);
    16541654                    newsummit=zeros(3,8);% initialize the rotated summit coordinates
    1655                     ObjectData.Coord=[ObjectData.Coord(1,1); ObjectData.Coord(1,2); 0];%TODO: set in set_object
     1655                    ObjectData.Coord=ObjectData.Coord';
     1656                    if size(ObjectData.Coord,2)<3
     1657                    ObjectData.Coord=[ObjectData.Coord; 0];%add z origin at z=0 by default
     1658                    end
    16561659                    for isummit=1:8% TODO: introduce a function for rotation of n points (to use also for scattered data)
    16571660                        newsummit(:,isummit)=Mrot*(summit(:,isummit)-(ObjectData.Coord));
    16581661                    end
    1659                     coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:));
     1662                    coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:));% set of coordinqtes in the projection plane
    16601663                    coord_y_proj=min(newsummit(2,:)):InterpMesh: max(newsummit(2,:));
    1661                     coord_z_proj=-width:InterpMesh:width;
    1662                     Mrot_inverse=rodrigues(-PlaneAngle);
     1664                    coord_z_proj=-width:width;
     1665                    Mrot_inverse=rodrigues(-PlaneAngle);% inverse rotation matrix
    16631666                    Origin=Mrot_inverse*[coord_x_proj(1);coord_y_proj(1);coord_z_proj(1)]+ObjectData.Coord;
    1664                     ix=Mrot_inverse*[coord_x_proj(end)-coord_x_proj(1);0;0];% unit vector along the new x coordinates
    1665                     iy=Mrot_inverse*[0;coord_y_proj(end)-coord_y_proj(1);0];% unit vector y coordinates
    1666                     iz=Mrot_inverse*[0;0;coord_z_proj(end)-coord_z_proj(1)];% x unit vector z coordinates
    1667                     [Grid_x,Grid_y,Grid_z]=meshgrid(1:numel(coord_x_proj),1:numel(coord_y_proj),1:numel(coord_z_proj));
    1668                     if ismatrix(Grid_x)
     1667                    npx=numel(coord_x_proj);
     1668                    npy=numel(coord_y_proj);
     1669                    npz=numel(coord_z_proj);
     1670                    ix=Mrot_inverse*[coord_x_proj(end)-coord_x_proj(1);0;0]/(npx-1);% unit vector along the new x coordinates
     1671                    iy=Mrot_inverse*[0;coord_y_proj(end)-coord_y_proj(1);0]/(npy-1);% unit vector y coordinates
     1672                    iz=Mrot_inverse*[0;0;coord_z_proj(end)-coord_z_proj(1)]/(npz-1);% x unit vector z coordinates
     1673                    [Grid_x,Grid_y,Grid_z]=meshgrid(0:npx-1,0:npy-1,0:npz-1);
     1674                    if ismatrix(Grid_x)% add a singleton in case of a single z value
    16691675                        Grid_x=shiftdim(Grid_x,-1);
    16701676                        Grid_y=shiftdim(Grid_y,-1);
    … …  
    16811687                            VarName=FieldData.ListVarName{ivar};
    16821688                            ListVarName=[ListVarName VarName];
     1689                            VarDimName=[VarDimName {{'coord_y','coord_x'}}];
    16831690                            VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
    16841691                            FieldData.(VarName)=permute(FieldData.(VarName),[2 3 1]);
    1685                             ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI);
     1692                            ProjData.coord_x=coord_x_proj;
     1693                            ProjData.coord_y=coord_y_proj;
     1694                            ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI,'*linear');
     1695                            ProjData.(VarName)=nanmean(ProjData.(VarName),3);
     1696                            ProjData.(VarName)=squeeze(ProjData.(VarName));
    16861697                    end
    16871698                end

  • trunk/src/script_readlvm.m

    r950 r959  
    11%% get the input file
    2 fileinput=uigetfile_uvmat('pick an input file','/fsnet/project/coriolis/2016/16MILESTONE/Data');
     2project='/fsnet/project/coriolis/2015/15MINI_MEDDY/PROBES';
     3if ~exist(project,'dir')
     4    project='U:\project\coriolis\2015\15MINI_MEDDY\PROBES';%windows
     5end
     6fileinput=uigetfile_uvmat('pick an input file',project);
    37[Path,Name,Ext]=fileparts(fileinput);
    48
    … …  
    711if strcmp(Ext,'.lvm')
    812    disp(['reading ' fileinput '...'])
    9 Data=read_lvm(fileinput)
     13    Data=read_lvm(fileinput)
    1014elseif strcmp(Ext,'.nc')
    1115    disp(['reading ' fileinput '...'])
    … …  
    1519end
    1620
    17 %% save as netcdf file if it has been opened as .lvm
     21%% save netcdf file as .lvm
    1822if strcmp(Ext,'.lvm')
    1923    OutputFile=fullfile(Path,[Name '.nc']);
    … …  
    2731end
    2832
    29 %% use calibration data stored in an xml file
    30 ProbeDoc=[];
    31 XmlFile=fullfile(Path,[Name '.xml']);
     33%% use calibration stored in a specified xml file; always use the same file
     34ProbeDoc=[];
     35%XmlFile=fullfile(Path,[Name '.xml']);
     36XmlFile=fullfile(Path,'calib.xml'); 
    3237if exist(XmlFile,'file')
    33     ProbeDoc=xml2struct(XmlFile);
    34 end
     38    ProbeDoc=xml2struct(XmlFile)
     39else
     40    disp('no calibration file .xml detected')
     41end
     42
    3543% if isfield(Data,'Position'), C2,C4,C6
    36 %     Min=1;
    37 %     Data.Position=Data.Position+PositionMin;
     44%     Min=1; Data.Position=Data.Position+PositionMin;
    3845% end
    39 % a5=-2.134088,b5=1010.1611,
    40 % Data.C2=Data.C2;
    41 % Data.C5=a5*Data.C5+b5;
    42 % Data.C6=Data.C2;
     46% a5=-2.134088,b5=1010.1611, Data.C2=Data.C2; Data.C5=a5*Data.C5+b5; Data.C6=Data.C2;
    4347% ylabelstring='density drho (kg/m3)';
    4448
    45 %% treqnsform conductivity probe signals: calibration + filter at 10 Hz
    46 ylabelstring='conductivity signal (volts)';
    47 clist=0;% counter of conductivity probes
     49
     50%% transform temperature probe signals
     51if isfield(ProbeDoc,'T5')&& ~isempty(ProbeDoc.T5)  % if temperature calibration exists; see calibT.m
     52   
     53    Data.T5=exp((Data.T5 - ProbeDoc.T5.a)./ProbeDoc.T5.b);% transform volt signal into temperature
     54    Data.T5=filter(ones(1,60)/60,1,Data.T5); % filter the signal to 4 Hz
     55    figure(6); set(6,'name','temperature'); plot(Data.Time,Data.T5);
     56    %plot(Data.Time,Data.T5,Data.Time,20+Data.Position/100)
     57    title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', temperature'])
     58    xlabel('Time(s)'); ylabel('temperature (degree C)'); %ylim([20 37])
     59    grid on
     60end
     61
     62%% check camera signal
     63ind_start=find(Data.Trig_Cam>3.5,1,'first')
     64disp(['camera starts at time ' num2str(Data.Time(ind_start))])
     65%% transform and filter conductivity probe signals into [temperature-corrected] density
     66ylabelstring='conductivity signal (volts)'; clist=0;% counter of conductivity probes
    4867for ilist=1:numel(Data.ListVarName)
    4968    if isequal(Data.ListVarName{ilist}(1),'C');% if the var name begins by 'C'
    … …  
    5170        CName{clist}=Data.ListVarName{ilist};
    5271        if isfield(ProbeDoc,CName{clist})&& ~isempty(ProbeDoc.(CName{clist}))
    53             a=ProbeDoc.(CName{clist}).a;
    54             b=ProbeDoc.(CName{clist}).b;
    55             Data.(CName{clist})=a*Data.(CName{clist})+b;% transform volt signal into density
     72       
     73            a=ProbeDoc.(CName{clist}).a; b=ProbeDoc.(CName{clist}).b;
     74            % Data.(CName{clist})=a*Data.(CName{clist})+b;% volts STRAIGHT into density (if const T)
     75           
     76            %% BUT now need to modify conductivity, density due to temperature effect
     77            %% first put into conductivity - using just C5 conductivity calibration for now (22/7/15)
     78            ac1 = 0.5766e4; bc1 = 2.9129e4; %% this was found later, w/ different gain. Use the one below instead
     79            %% ac1 = 0.4845e4; bc1 = 2.064e4; %% this was w/ the gain when we did the experiments; but it doesn't work?
     80            Data.(CName{clist})=ac1*Data.(CName{clist})+bc1; %% voltage translated into conductivity via calibration
     81             
     82            %% read in temperature...
     83            refT = 23; T = Data.T5;
     84            Hewitt_fit = [2.2794885e-11 -6.2634979e-9 1.5439826e-7 7.8601061e-5 2.1179818e-2];
     85            bfit = reshape(polyval(Hewitt_fit,T(:)),size(T)); bref = reshape(polyval(Hewitt_fit,refT(:)),size(refT));
     86            sigTsig18 = (1+bfit.*(T-18)); sigREFsig18 = (1+bref.*(refT-18));
     87            modfactor = sigTsig18./sigREFsig18;
     88            modfactor = 1./modfactor;  %% now in sigREF/sigT, which (* sigT) below to get sigREF, which can convert to rho
     89            Data.(CName{clist})=Data.(CName{clist}).*modfactor %% = temp corrected conductivity (= as if at reference temp)
     90           
     91            %% now we need to put the modified conductivity back into a (modified) voltage, then estimate density directly.
     92            ac2 = 1.7343e-4; bc2 = -5.0048;
     93            Data.(CName{clist})=ac2*Data.(CName{clist})+bc2; % temp-corrected voltage
     94            Data.(CName{clist})=a*Data.(CName{clist})+b; % now finally voltage into density
    5695            Data.(CName{clist})=filter(ones(1,20)/20,1,Data.(CName{clist})); % filter the signal to 10 Hz
    5796            ylabelstring='density drho (g/cm3)';
     97                 
    5898        end
    5999    end
    60100end
    61101
    62 %% plot conductivity probe signals
    63 figure(1)
    64 set(1,'name','conductivity')
     102%% plot conductivity signals
     103figure(1); set(1,'name','conductivity')
     104bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz
    65105plot_string='plot(';
    66106for clist=1:numel(CName)
     107    Data.(CName{clist})=filter(ones(1,bandwidth2)/bandwidth2,1,Data.(CName{clist}));%low pass filter
    67108    plot_string=[plot_string 'Data.Time,Data.' CName{clist} ','];
    68109end
    … …  
    74115xlabel('Time(s)')
    75116ylabel(ylabelstring)
     117        ylim([1 1.02])
    76118grid on
    77119
    78120if isfield(Data,'Position')
    79     %% plot  motor position
     121    %% plot motor position
    80122    figure(2)
    81123    set(2,'name','position')
    … …  
    88130    hold on
    89131    plot(Data.Time,Data.Speed)
    90    
    91     %%plot first profile
    92     figure(4)
    93     index1=find(Data.Speed<0,1); %start of descent
    94     Speed=Data.Speed(index1:end);
    95     index2=index1-1+find(Speed>0,1); %end of descent
    96     plot(Data.Position(index1:index2),Data.C1(index1:index2))
    97    
     132       
    98133    %% plot conductivity probe profiles (limited to downward motion, Data.Speed<0)
    99134    if ~isempty(ProbeDoc)
    … …  
    119154        xlabel('Z(cm)')
    120155        ylabel(ylabelstring)
     156        ylim([1 1.02])
    121157        grid on
    122158    end
    123     %     plot(Z,Data.C2(Data.Speed<0),Z,Data.C3(Data.Speed<0),Z,Data.C4(Data.Speed<0),Z,Data.C5(Data.Speed<0))
    124     %     legend({'C2';'C3';'C4';'C5'})
    125     %     title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', profiles conductivity probes'])
    126     %     xlabel('Z(cm)')
    127     %     %ylabel('signal (Volt)')
    128     %     ylabel(ylabelstring)
    129     %     grid on
    130    
    131    
    132     %     set(3,'name','profiles')
    133     %     PositionMin=1; %position of probes at the bottom (Z=1 cm)
    134     %     Data.Position=Data.Position-min(Data.Position)+PositionMin;
    135     %     Z=Data.Position(Data.Speed<0);
    136     %     plot(Z,Data.C2(Data.Speed<0),Z,Data.C5(Data.Speed<0),Z,Data.C6(Data.Speed<0))
    137     %     legend({'C2';'C5';'C6'})
    138     %     title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', profiles conductivity probes'])
    139     %     xlabel('Z(cm)')
    140     %     %ylabel('signal (Volt)')
    141     %     ylabel(ylabelstring)
    142     %     grid on
    143    
    144 end
    145 
    146 
    147 
    148 
    149 %% plot temperature signals
    150 % figure(4)
    151 % set(4,'name','temperature')
    152 % plot(Data.Time,Data.T2,Data.Time,Data.T5,Data.Time,Data.T6)
    153 % legend({'T2';'T5';'T6'})
    154 % title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', signal conductivity probes'])
    155 % xlabel('Time(s)')
    156 % ylabel('signal (Volt)')
    157 % grid on
     159       
     160end
     161
     162%%%%
     163figure(4)
     164bandwidth1=480; % corresponds to 2 cm (1 cm/s with 240 pts/s)
     165bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz
     166C5_filter_low=filter(ones(1,bandwidth2)/bandwidth2,1,Data.C5);%low pass filter
     167C5_filter=filter(ones(1,bandwidth1)/bandwidth1,1,C5_filter_low);%low pass filter
     168C5_filter=C5_filter_low-C5_filter;% high pass filter
     169C5_filter(Data.Speed>-0.1)=NaN;
     170plot(Data.Time,C5_filter)
     171ylim([-0.001 0.001])
     172hold on
     173plot(Data.Time,Data.Position/100000)
     174grid on
     175hold off
     176title('C5 filtered')
     177
     178%%%%
     179figure(5)
     180bandwidth1=480; % corresponds to 2 cm (1 cm/s with 240 pts/s)
     181bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz
     182C3_filter_low=filter(ones(1,bandwidth2)/bandwidth2,1,Data.C3);%low pass filter
     183C3_filter=filter(ones(1,bandwidth1)/bandwidth1,1,C3_filter_low);%low pass filter
     184C3_filter=C3_filter_low-C3_filter;% high pass filter
     185C3_filter(Data.Speed>-0.1)=NaN;
     186plot(Data.Time,C3_filter)
     187ylim([-0.001 0.001])
     188hold on
     189plot(Data.Time,Data.Position/100000)
     190grid on
     191hold off
     192title('C3 filtered')
    158193
    159194%% plot velocity (ADV) signals
    … …  
    178213
    179214
     215

  • trunk/src/series/aver_stat.m

    r924 r959  
    111111        end
    112112    end
     113    % determine volume scan mode
     114        prompt = {'volume scan mode (Yes/No)'};
     115    dlg_title = 'determine volume scan';
     116    num_lines= 1;
     117    def     = { 'No'};
     118     answer=msgbox_uvmat('INPUT_Y-N','volume scan mode (OK/No)?');
     119%     answer = inputdlg(prompt,dlg_title,num_lines,def);
     120    if isempty(answer)
     121        return
     122    end
     123    %check input consistency
     124    if strcmp(answer,'Yes')
     125            ParamOut.NbSlice=1;% set NbSlice to 1 ( for i index)
     126            ParamOut.ActionInput.CheckVolume=1;
     127    end
    113128    return
    114129end
    … …  
    248263    end
    249264end
    250 %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
    251 for index=1:NbField
    252     update_waitbar(WaitbarHandle,index/NbField)
    253     if ~isempty(RUNHandle)&& ~strcmp(get(RUNHandle,'BusyAction'),'queue')
    254         disp('program stopped by user')
    255         break
    256     end
     265
     266%% set volume scan
     267NbSlice_j=1;
     268index_series=1:size(filecell,2);
     269first_j_out=first_j;
     270last_j_out=last_j;
     271if isfield(Param,'ActionInput') && isfield(Param.ActionInput,'CheckVolume') ...
     272        && Param.ActionInput.CheckVolume
     273   index_j=Param.IndexRange.first_j:Param.IndexRange.incr_j:Param.IndexRange.last_j;
     274   NbSlice_j=numel(index_j);
     275   index_i=index_series(1:NbSlice_j:end);
     276end
     277
     278%%% loop on slices (volume scan)
     279for islice=index_j
     280    if NbSlice_j>1
     281    first_j_out=islice;
     282    last_j_out=islice;
    257283   
    258     %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
    259     for iview=1:NbView
    260         % reading input file(s)
    261         [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index));
    262         if ~isempty(errormsg)
    263             errormsg=['error of input reading: ' errormsg];
     284    end
     285    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     286    %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
     287    for index=index_i+index_j(islice)
     288        update_waitbar(WaitbarHandle,index/NbField)
     289        if ~isempty(RUNHandle)&& ~strcmp(get(RUNHandle,'BusyAction'),'queue')
     290            disp('program stopped by user')
    264291            break
    265292        end
    266         if ~isempty(NbSlice_calib)
    267             Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
    268         end
    269     end
    270     %%%%%%%%%%%%%%%% end loop on views (input lines) %%%%%%%%%%%%%%%%
    271     %%%%%%%%%%%% END STANDARD PART  %%%%%%%%%%%%
    272     % EDIT FROM HERE
    273    
    274     if isempty(errormsg)
    275         Field=Data{1}; % default input field structure
    276         nbfiles=nbfiles+1; %increment the file counter
    277         %% coordinate transform (or other user defined transform)
    278         if ~isempty(transform_fct)
    279             switch nargin(transform_fct)
    280                 case 4
    281                     if length(Data)==2
    282                         Field=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});
     293       
     294        %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%%
     295        for iview=1:NbView
     296            % reading input file(s)
     297            %InputFile=fullfile_uvmat(RootPath{iview},SubDir{iview},RootFile{iview},FileExt{iview},NomType{iview},first_i,last_i,first_j_out,last_j_out);
     298            [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index));
     299            if ~isempty(errormsg)
     300                errormsg=['error of input reading: ' errormsg];
     301                break
     302            end
     303            if ~isempty(NbSlice_calib)
     304                Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform
     305            end
     306        end
     307        %%%%%%%%%%%%%%%% end loop on views (input lines) %%%%%%%%%%%%%%%%
     308       
     309        if isempty(errormsg)
     310            Field=Data{1}; % default input field structure
     311            nbfiles=nbfiles+1; %increment the file counter
     312            %% coordinate transform (or other user defined transform)
     313            if ~isempty(transform_fct)
     314                switch nargin(transform_fct)
     315                    case 4
     316                        if length(Data)==2
     317                            Field=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2});
     318                        else
     319                            Field=transform_fct(Data{1},XmlData{1});
     320                        end
     321                    case 3
     322                        if length(Data)==2
     323                            Field=transform_fct(Data{1},XmlData{1},Data{2});
     324                        else
     325                            Field=transform_fct(Data{1},XmlData{1});
     326                        end
     327                    case 2
     328                        Field=transform_fct(Data{1},XmlData{1});
     329                    case 1
     330                        Field=transform_fct(Data{1});
     331                end
     332            end
     333           
     334            %% field projection on an object
     335            if Param.CheckObject
     336                if strcmp(Param.ProjObject.ProjMode,'interp_tps')
     337                    Field=tps_coeff_field(Field,check_proj_tps);% calculate tps coefficients if needed
     338                end
     339                [Field,errormsg]=proj_field(Field,Param.ProjObject,VarMesh);
     340                if ~isempty(errormsg)
     341                    disp_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg],checkrun)
     342                    return
     343                end
     344            end
     345           
     346            %%%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
     347            if nbfiles==1 %first field
     348                time_1=[];
     349                if isfield(Field,'Time')
     350                    time_1=Field.Time(1);
     351                end
     352                DataOut=Field;%outcome reproduces the first (projected) field by default
     353                DataOut.Conventions='uvmat'; %suppress Conventions='uvmat/civdata' for civ input files
     354                if isfield(Param,'ProjObject')&& ismember(Param.ProjObject.ProjMode,{'inside','outside'})%case of histograms
     355                    for ivar=1:numel(Field.ListVarName)% list of variable names before projection (histogram)
     356                        VarName=Field.ListVarName{ivar};
     357                        if isfield(Data{1},VarName)
     358                            DataOut.(VarName)=Field.(VarName);
     359                            DataOut.([VarName 'Histo'])=zeros(size(DataOut.(VarName)));
     360                            VarMesh=DataOut.(VarName)(2)-DataOut.(VarName)(1);
     361                        end
     362                    end
     363                    disp(['mesh for histogram = ' num2str(VarMesh)])
     364                else
     365                    errorvar=zeros(numel(Field.ListVarName));%index of errorflag associated to each variable
     366                    if isfield(Field,'VarAttribute')
     367                        for ivar=1:numel(Field.ListVarName)
     368                            VarName=Field.ListVarName{ivar};
     369                            DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero
     370                            NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero
     371                           
     372                            for iivar=1:length(Field.VarAttribute)
     373                                if isequal(Field.VarDimName{iivar},Field.VarDimName{ivar})&& isfield(Field.VarAttribute{iivar},'Role')...
     374                                        && strcmp(Field.VarAttribute{iivar}.Role,'errorflag')
     375                                    errorvar(ivar)=iivar; % index of the errorflag variable corresponding to ivar
     376                                end
     377                            end
     378                        end
     379                        DataOut.ListVarName(errorvar(errorvar~=0))=[]; %remove errorflag from result
     380                        DataOut.VarDimName(errorvar(errorvar~=0))=[]; %remove errorflag from result
     381                        DataOut.VarAttribute(errorvar(errorvar~=0))=[]; %remove errorflag from result
    283382                    else
    284                         Field=transform_fct(Data{1},XmlData{1});
     383                        for ivar=1:numel(Field.ListVarName)
     384                            VarName=Field.ListVarName{ivar};
     385                            DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero
     386                            NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero
     387                        end
    285388                    end
    286                 case 3
    287                     if length(Data)==2
    288                         Field=transform_fct(Data{1},XmlData{1},Data{2});
     389                   
     390                end
     391            end   %current field
     392            for ivar=1:length(DataOut.ListVarName)
     393                VarName=DataOut.ListVarName{ivar};
     394                sizmean=size(DataOut.(VarName));
     395                siz=size(Field.(VarName));
     396                if isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'})
     397                    if isfield(Data{1},VarName)
     398                        MaxValue=max(DataOut.(VarName));% current max of histogram absissa
     399                        MinValue=min(DataOut.(VarName));% current min of histogram absissa
     400                        %                     VarMesh=Field.VarAttribute{ivar}.Mesh;
     401                        MaxIndex=round(MaxValue/VarMesh);
     402                        MinIndex=round(MinValue/VarMesh);
     403                        MaxIndex_new=round(max(Field.(VarName)/VarMesh));% max of the current field
     404                        MinIndex_new=round(min(Field.(VarName)/VarMesh));
     405                        if MaxIndex_new>MaxIndex% the variable max for the current field exceeds the previous one
     406                            DataOut.(VarName)=[DataOut.(VarName) VarMesh*(MaxIndex+1:MaxIndex_new)];% append the new variable values
     407                            DataOut.([VarName 'Histo'])=[DataOut.([VarName 'Histo']) zeros(1,MaxIndex_new-MaxIndex)]; % append the new histo values
     408                        end
     409                        if MinIndex_new <= MinIndex-1
     410                            DataOut.(VarName)=[VarMesh*(MinIndex_new:MinIndex-1) DataOut.(VarName)];% insert the new variable values
     411                            DataOut.([VarName 'Histo'])=[zeros(1,MinIndex-MinIndex_new) DataOut.([VarName 'Histo'])];% insert the new histo values
     412                            ind_start=1;
     413                        else
     414                            ind_start=MinIndex_new-MinIndex+1;
     415                        end
     416                        DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)=...
     417                            DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)+Field.([VarName 'Histo']);
     418                    end
     419                elseif ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean)
     420                    disp_uvmat('ERROR',['unequal size of input field ' VarName ', need to project  on a grid'],checkrun)
     421                    return
     422                else
     423                    if errorvar(ivar)==0
     424                        check_bad=isnan(Field.(VarName));%=0 for NaN data values, 1 else
    289425                    else
    290                         Field=transform_fct(Data{1},XmlData{1});
     426                        check_bad=isnan(Field.(VarName)) | Field.(Field.ListVarName{errorvar(ivar)})~=0;%=0 for NaN or error flagged data values, 1 else
    291427                    end
    292                 case 2
    293                     Field=transform_fct(Data{1},XmlData{1});
    294                 case 1
    295                     Field=transform_fct(Data{1});
    296             end
    297         end
    298        
    299         %% field projection on an object
    300         if Param.CheckObject
    301             if strcmp(Param.ProjObject.ProjMode,'interp_tps')
    302                 Field=tps_coeff_field(Field,check_proj_tps);% calculate tps coefficients if needed
    303             end
    304             [Field,errormsg]=proj_field(Field,Param.ProjObject,VarMesh);
    305             if ~isempty(errormsg)
    306                 disp_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg],checkrun)
    307                 return
    308             end
    309         end
    310              
    311         %%%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
    312         if nbfiles==1 %first field
    313             time_1=[];
    314             if isfield(Field,'Time')
    315                 time_1=Field.Time(1);
    316             end
    317             DataOut=Field;%outcome reproduces the first (projected) field by default
    318             DataOut.Conventions='uvmat'; %suppress Conventions='uvmat/civdata' for civ input files         
    319             if isfield(Param,'ProjObject')&& ismember(Param.ProjObject.ProjMode,{'inside','outside'})%case of histograms
    320                 for ivar=1:numel(Field.ListVarName)% list of variable names before projection (histogram)
    321                     VarName=Field.ListVarName{ivar};
    322                     if isfield(Data{1},VarName)
    323                         DataOut.(VarName)=Field.(VarName);
    324                         DataOut.([VarName 'Histo'])=zeros(size(DataOut.(VarName)));
    325                         VarMesh=DataOut.(VarName)(2)-DataOut.(VarName)(1);
    326                     end
    327                 end
    328                 disp(['mesh for histogram = ' num2str(VarMesh)])
    329             else
    330                 errorvar=zeros(numel(Field.ListVarName));%index of errorflag associated to each variable
    331                 if isfield(Field,'VarAttribute')
    332                     for ivar=1:numel(Field.ListVarName)
    333                         VarName=Field.ListVarName{ivar};
    334                         DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero
    335                         NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero
    336                        
    337                         for iivar=1:length(Field.VarAttribute)
    338                             if isequal(Field.VarDimName{iivar},Field.VarDimName{ivar})&& isfield(Field.VarAttribute{iivar},'Role')...
    339                                     && strcmp(Field.VarAttribute{iivar}.Role,'errorflag')
    340                                 errorvar(ivar)=iivar; % index of the errorflag variable corresponding to ivar
    341                             end
    342                         end
    343                     end
    344                     DataOut.ListVarName(errorvar(errorvar~=0))=[]; %remove errorflag from result
    345                     DataOut.VarDimName(errorvar(errorvar~=0))=[]; %remove errorflag from result
    346                     DataOut.VarAttribute(errorvar(errorvar~=0))=[]; %remove errorflag from result
    347                 else
    348                                    for ivar=1:numel(Field.ListVarName)
    349                                        VarName=Field.ListVarName{ivar};
    350                         DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero
    351                         NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero
    352                                    end
    353                 end
    354                
    355             end
    356         end   %current field
    357         for ivar=1:length(DataOut.ListVarName)
    358             VarName=DataOut.ListVarName{ivar};
    359             sizmean=size(DataOut.(VarName));
    360             siz=size(Field.(VarName));
    361             if isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'})
    362                 if isfield(Data{1},VarName)
    363                     MaxValue=max(DataOut.(VarName));% current max of histogram absissa
    364                     MinValue=min(DataOut.(VarName));% current min of histogram absissa
    365 %                     VarMesh=Field.VarAttribute{ivar}.Mesh;
    366                     MaxIndex=round(MaxValue/VarMesh);
    367                     MinIndex=round(MinValue/VarMesh);
    368                     MaxIndex_new=round(max(Field.(VarName)/VarMesh));% max of the current field
    369                     MinIndex_new=round(min(Field.(VarName)/VarMesh));
    370                     if MaxIndex_new>MaxIndex% the variable max for the current field exceeds the previous one
    371                         DataOut.(VarName)=[DataOut.(VarName) VarMesh*(MaxIndex+1:MaxIndex_new)];% append the new variable values
    372                         DataOut.([VarName 'Histo'])=[DataOut.([VarName 'Histo']) zeros(1,MaxIndex_new-MaxIndex)]; % append the new histo values                   
    373                     end
    374                     if MinIndex_new <= MinIndex-1
    375                         DataOut.(VarName)=[VarMesh*(MinIndex_new:MinIndex-1) DataOut.(VarName)];% insert the new variable values
    376                         DataOut.([VarName 'Histo'])=[zeros(1,MinIndex-MinIndex_new) DataOut.([VarName 'Histo'])];% insert the new histo values
    377                         ind_start=1;
    378                     else
    379                         ind_start=MinIndex_new-MinIndex+1;
    380                     end
    381                     DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)=...
    382                         DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)+Field.([VarName 'Histo']);   
    383                 end
    384             elseif ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean)
    385                 disp_uvmat('ERROR',['unequal size of input field ' VarName ', need to project  on a grid'],checkrun)
    386                 return
    387             else
    388                 if errorvar(ivar)==0
    389                     check_bad=isnan(Field.(VarName));%=0 for NaN data values, 1 else
    390                 else
    391                     check_bad=isnan(Field.(VarName)) | Field.(Field.ListVarName{errorvar(ivar)})~=0;%=0 for NaN or error flagged data values, 1 else
    392                 end
    393                 Field.(VarName)(check_bad)=0; %set to zero NaN or data marked by error flag
    394                 DataOut.(VarName)=DataOut.(VarName)+ double(Field.(VarName)); % update the sum
    395                 NbData.(VarName)=NbData.(VarName)+ ~check_bad;% records the number of data for each point
    396             end
    397         end
    398         %%%%%%%%%%%%   END MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
    399     else
    400         disp(errormsg)
    401     end
    402 end
    403 %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%%
     428                    Field.(VarName)(check_bad)=0; %set to zero NaN or data marked by error flag
     429                    DataOut.(VarName)=DataOut.(VarName)+ double(Field.(VarName)); % update the sum
     430                    NbData.(VarName)=NbData.(VarName)+ ~check_bad;% records the number of data for each point
     431                end
     432            end
     433            %%%%%%%%%%%%   END MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
     434        else
     435            disp(errormsg)
     436        end
     437    end
     438    %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%%
     439
    404440if ~(isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'}))
    405441    for ivar=1:length(Field.ListVarName)
    … …  
    426462
    427463%% writing the result file
    428 OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j,last_j);
     464OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j_out,last_j_out);
    429465if strcmp(FileExtOut,'.png') %case of images
    430466    if isequal(FileInfo{1}.BitDepth,16)||(numel(FileInfo)==2 &&isequal(FileInfo{2}.BitDepth,16))
    … …  
    444480    end
    445481end  % end averaging  loop
     482end
    446483
    447484%% open the result file with uvmat (in RUN mode)

  • trunk/src/series/merge_proj.m

    r956 r959  
    103103    WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series
    104104end
    105 tild=phys([],[]);% test to provoke the inclusion of the function phys at compilation
     105
    106106%% define the directory for result file (with path=RootPath{1})
    107107OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files
Note: See TracChangeset for help on using the changeset viewer.