Index: /trunk/src/proj_field.m =================================================================== --- /trunk/src/proj_field.m (revision 958) +++ /trunk/src/proj_field.m (revision 959) @@ -1653,18 +1653,24 @@ Mrot=rodrigues(PlaneAngle); newsummit=zeros(3,8);% initialize the rotated summit coordinates - ObjectData.Coord=[ObjectData.Coord(1,1); ObjectData.Coord(1,2); 0];%TODO: set in set_object + ObjectData.Coord=ObjectData.Coord'; + if size(ObjectData.Coord,2)<3 + ObjectData.Coord=[ObjectData.Coord; 0];%add z origin at z=0 by default + end for isummit=1:8% TODO: introduce a function for rotation of n points (to use also for scattered data) newsummit(:,isummit)=Mrot*(summit(:,isummit)-(ObjectData.Coord)); end - coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:)); + coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:));% set of coordinqtes in the projection plane coord_y_proj=min(newsummit(2,:)):InterpMesh: max(newsummit(2,:)); - coord_z_proj=-width:InterpMesh:width; - Mrot_inverse=rodrigues(-PlaneAngle); + coord_z_proj=-width:width; + Mrot_inverse=rodrigues(-PlaneAngle);% inverse rotation matrix Origin=Mrot_inverse*[coord_x_proj(1);coord_y_proj(1);coord_z_proj(1)]+ObjectData.Coord; - ix=Mrot_inverse*[coord_x_proj(end)-coord_x_proj(1);0;0];% unit vector along the new x coordinates - iy=Mrot_inverse*[0;coord_y_proj(end)-coord_y_proj(1);0];% unit vector y coordinates - iz=Mrot_inverse*[0;0;coord_z_proj(end)-coord_z_proj(1)];% x unit vector z coordinates - [Grid_x,Grid_y,Grid_z]=meshgrid(1:numel(coord_x_proj),1:numel(coord_y_proj),1:numel(coord_z_proj)); - if ismatrix(Grid_x) + npx=numel(coord_x_proj); + npy=numel(coord_y_proj); + npz=numel(coord_z_proj); + ix=Mrot_inverse*[coord_x_proj(end)-coord_x_proj(1);0;0]/(npx-1);% unit vector along the new x coordinates + iy=Mrot_inverse*[0;coord_y_proj(end)-coord_y_proj(1);0]/(npy-1);% unit vector y coordinates + iz=Mrot_inverse*[0;0;coord_z_proj(end)-coord_z_proj(1)]/(npz-1);% x unit vector z coordinates + [Grid_x,Grid_y,Grid_z]=meshgrid(0:npx-1,0:npy-1,0:npz-1); + if ismatrix(Grid_x)% add a singleton in case of a single z value Grid_x=shiftdim(Grid_x,-1); Grid_y=shiftdim(Grid_y,-1); @@ -1681,7 +1687,12 @@ VarName=FieldData.ListVarName{ivar}; ListVarName=[ListVarName VarName]; + VarDimName=[VarDimName {{'coord_y','coord_x'}}]; VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes FieldData.(VarName)=permute(FieldData.(VarName),[2 3 1]); - ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI); + ProjData.coord_x=coord_x_proj; + ProjData.coord_y=coord_y_proj; + ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI,'*linear'); + ProjData.(VarName)=nanmean(ProjData.(VarName),3); + ProjData.(VarName)=squeeze(ProjData.(VarName)); end end Index: /trunk/src/script_readlvm.m =================================================================== --- /trunk/src/script_readlvm.m (revision 958) +++ /trunk/src/script_readlvm.m (revision 959) @@ -1,4 +1,8 @@ %% get the input file -fileinput=uigetfile_uvmat('pick an input file','/fsnet/project/coriolis/2016/16MILESTONE/Data'); +project='/fsnet/project/coriolis/2015/15MINI_MEDDY/PROBES'; +if ~exist(project,'dir') + project='U:\project\coriolis\2015\15MINI_MEDDY\PROBES';%windows +end +fileinput=uigetfile_uvmat('pick an input file',project); [Path,Name,Ext]=fileparts(fileinput); @@ -7,5 +11,5 @@ if strcmp(Ext,'.lvm') disp(['reading ' fileinput '...']) -Data=read_lvm(fileinput) + Data=read_lvm(fileinput) elseif strcmp(Ext,'.nc') disp(['reading ' fileinput '...']) @@ -15,5 +19,5 @@ end -%% save as netcdf file if it has been opened as .lvm +%% save netcdf file as .lvm if strcmp(Ext,'.lvm') OutputFile=fullfile(Path,[Name '.nc']); @@ -27,23 +31,38 @@ end -%% use calibration data stored in an xml file -ProbeDoc=[]; -XmlFile=fullfile(Path,[Name '.xml']); +%% use calibration stored in a specified xml file; always use the same file +ProbeDoc=[]; +%XmlFile=fullfile(Path,[Name '.xml']); +XmlFile=fullfile(Path,'calib.xml'); if exist(XmlFile,'file') - ProbeDoc=xml2struct(XmlFile); -end + ProbeDoc=xml2struct(XmlFile) +else + disp('no calibration file .xml detected') +end + % if isfield(Data,'Position'), C2,C4,C6 -% Min=1; -% Data.Position=Data.Position+PositionMin; +% Min=1; Data.Position=Data.Position+PositionMin; % end -% a5=-2.134088,b5=1010.1611, -% Data.C2=Data.C2; -% Data.C5=a5*Data.C5+b5; -% Data.C6=Data.C2; +% a5=-2.134088,b5=1010.1611, Data.C2=Data.C2; Data.C5=a5*Data.C5+b5; Data.C6=Data.C2; % ylabelstring='density drho (kg/m3)'; -%% treqnsform conductivity probe signals: calibration + filter at 10 Hz -ylabelstring='conductivity signal (volts)'; -clist=0;% counter of conductivity probes + +%% transform temperature probe signals +if isfield(ProbeDoc,'T5')&& ~isempty(ProbeDoc.T5) % if temperature calibration exists; see calibT.m + + Data.T5=exp((Data.T5 - ProbeDoc.T5.a)./ProbeDoc.T5.b);% transform volt signal into temperature + Data.T5=filter(ones(1,60)/60,1,Data.T5); % filter the signal to 4 Hz + figure(6); set(6,'name','temperature'); plot(Data.Time,Data.T5); + %plot(Data.Time,Data.T5,Data.Time,20+Data.Position/100) + title([Data.Experiment ', ' Data.FileName ', Time=' Data.DateTime ', temperature']) + xlabel('Time(s)'); ylabel('temperature (degree C)'); %ylim([20 37]) + grid on +end + +%% check camera signal +ind_start=find(Data.Trig_Cam>3.5,1,'first') +disp(['camera starts at time ' num2str(Data.Time(ind_start))]) +%% transform and filter conductivity probe signals into [temperature-corrected] density +ylabelstring='conductivity signal (volts)'; clist=0;% counter of conductivity probes for ilist=1:numel(Data.ListVarName) if isequal(Data.ListVarName{ilist}(1),'C');% if the var name begins by 'C' @@ -51,18 +70,40 @@ CName{clist}=Data.ListVarName{ilist}; if isfield(ProbeDoc,CName{clist})&& ~isempty(ProbeDoc.(CName{clist})) - a=ProbeDoc.(CName{clist}).a; - b=ProbeDoc.(CName{clist}).b; - Data.(CName{clist})=a*Data.(CName{clist})+b;% transform volt signal into density + + a=ProbeDoc.(CName{clist}).a; b=ProbeDoc.(CName{clist}).b; + % Data.(CName{clist})=a*Data.(CName{clist})+b;% volts STRAIGHT into density (if const T) + + %% BUT now need to modify conductivity, density due to temperature effect + %% first put into conductivity - using just C5 conductivity calibration for now (22/7/15) + ac1 = 0.5766e4; bc1 = 2.9129e4; %% this was found later, w/ different gain. Use the one below instead + %% ac1 = 0.4845e4; bc1 = 2.064e4; %% this was w/ the gain when we did the experiments; but it doesn't work? + Data.(CName{clist})=ac1*Data.(CName{clist})+bc1; %% voltage translated into conductivity via calibration + + %% read in temperature... + refT = 23; T = Data.T5; + Hewitt_fit = [2.2794885e-11 -6.2634979e-9 1.5439826e-7 7.8601061e-5 2.1179818e-2]; + bfit = reshape(polyval(Hewitt_fit,T(:)),size(T)); bref = reshape(polyval(Hewitt_fit,refT(:)),size(refT)); + sigTsig18 = (1+bfit.*(T-18)); sigREFsig18 = (1+bref.*(refT-18)); + modfactor = sigTsig18./sigREFsig18; + modfactor = 1./modfactor; %% now in sigREF/sigT, which (* sigT) below to get sigREF, which can convert to rho + Data.(CName{clist})=Data.(CName{clist}).*modfactor %% = temp corrected conductivity (= as if at reference temp) + + %% now we need to put the modified conductivity back into a (modified) voltage, then estimate density directly. + ac2 = 1.7343e-4; bc2 = -5.0048; + Data.(CName{clist})=ac2*Data.(CName{clist})+bc2; % temp-corrected voltage + Data.(CName{clist})=a*Data.(CName{clist})+b; % now finally voltage into density Data.(CName{clist})=filter(ones(1,20)/20,1,Data.(CName{clist})); % filter the signal to 10 Hz ylabelstring='density drho (g/cm3)'; + end end end -%% plot conductivity probe signals -figure(1) -set(1,'name','conductivity') +%% plot conductivity signals +figure(1); set(1,'name','conductivity') +bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz plot_string='plot('; for clist=1:numel(CName) + Data.(CName{clist})=filter(ones(1,bandwidth2)/bandwidth2,1,Data.(CName{clist}));%low pass filter plot_string=[plot_string 'Data.Time,Data.' CName{clist} ',']; end @@ -74,8 +115,9 @@ xlabel('Time(s)') ylabel(ylabelstring) + ylim([1 1.02]) grid on if isfield(Data,'Position') - %% plot motor position + %% plot motor position figure(2) set(2,'name','position') @@ -88,12 +130,5 @@ hold on plot(Data.Time,Data.Speed) - - %%plot first profile - figure(4) - index1=find(Data.Speed<0,1); %start of descent - Speed=Data.Speed(index1:end); - index2=index1-1+find(Speed>0,1); %end of descent - plot(Data.Position(index1:index2),Data.C1(index1:index2)) - + %% plot conductivity probe profiles (limited to downward motion, Data.Speed<0) if ~isempty(ProbeDoc) @@ -119,41 +154,41 @@ xlabel('Z(cm)') ylabel(ylabelstring) + ylim([1 1.02]) grid on end - % 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)) - % legend({'C2';'C3';'C4';'C5'}) - % title([Data.Experiment ', ' Data.FileName ', Time=' Data.DateTime ', profiles conductivity probes']) - % xlabel('Z(cm)') - % %ylabel('signal (Volt)') - % ylabel(ylabelstring) - % grid on - - - % set(3,'name','profiles') - % PositionMin=1; %position of probes at the bottom (Z=1 cm) - % Data.Position=Data.Position-min(Data.Position)+PositionMin; - % Z=Data.Position(Data.Speed<0); - % plot(Z,Data.C2(Data.Speed<0),Z,Data.C5(Data.Speed<0),Z,Data.C6(Data.Speed<0)) - % legend({'C2';'C5';'C6'}) - % title([Data.Experiment ', ' Data.FileName ', Time=' Data.DateTime ', profiles conductivity probes']) - % xlabel('Z(cm)') - % %ylabel('signal (Volt)') - % ylabel(ylabelstring) - % grid on - -end - - - - -%% plot temperature signals -% figure(4) -% set(4,'name','temperature') -% plot(Data.Time,Data.T2,Data.Time,Data.T5,Data.Time,Data.T6) -% legend({'T2';'T5';'T6'}) -% title([Data.Experiment ', ' Data.FileName ', Time=' Data.DateTime ', signal conductivity probes']) -% xlabel('Time(s)') -% ylabel('signal (Volt)') -% grid on + +end + +%%%% +figure(4) +bandwidth1=480; % corresponds to 2 cm (1 cm/s with 240 pts/s) +bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz +C5_filter_low=filter(ones(1,bandwidth2)/bandwidth2,1,Data.C5);%low pass filter +C5_filter=filter(ones(1,bandwidth1)/bandwidth1,1,C5_filter_low);%low pass filter +C5_filter=C5_filter_low-C5_filter;% high pass filter +C5_filter(Data.Speed>-0.1)=NaN; +plot(Data.Time,C5_filter) +ylim([-0.001 0.001]) +hold on +plot(Data.Time,Data.Position/100000) +grid on +hold off +title('C5 filtered') + +%%%% +figure(5) +bandwidth1=480; % corresponds to 2 cm (1 cm/s with 240 pts/s) +bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz +C3_filter_low=filter(ones(1,bandwidth2)/bandwidth2,1,Data.C3);%low pass filter +C3_filter=filter(ones(1,bandwidth1)/bandwidth1,1,C3_filter_low);%low pass filter +C3_filter=C3_filter_low-C3_filter;% high pass filter +C3_filter(Data.Speed>-0.1)=NaN; +plot(Data.Time,C3_filter) +ylim([-0.001 0.001]) +hold on +plot(Data.Time,Data.Position/100000) +grid on +hold off +title('C3 filtered') %% plot velocity (ADV) signals @@ -178,2 +213,3 @@ + Index: /trunk/src/series/aver_stat.m =================================================================== --- /trunk/src/series/aver_stat.m (revision 958) +++ /trunk/src/series/aver_stat.m (revision 959) @@ -111,4 +111,19 @@ end end + % determine volume scan mode + prompt = {'volume scan mode (Yes/No)'}; + dlg_title = 'determine volume scan'; + num_lines= 1; + def = { 'No'}; + answer=msgbox_uvmat('INPUT_Y-N','volume scan mode (OK/No)?'); +% answer = inputdlg(prompt,dlg_title,num_lines,def); + if isempty(answer) + return + end + %check input consistency + if strcmp(answer,'Yes') + ParamOut.NbSlice=1;% set NbSlice to 1 ( for i index) + ParamOut.ActionInput.CheckVolume=1; + end return end @@ -248,158 +263,179 @@ end end -%%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%% -for index=1:NbField - update_waitbar(WaitbarHandle,index/NbField) - if ~isempty(RUNHandle)&& ~strcmp(get(RUNHandle,'BusyAction'),'queue') - disp('program stopped by user') - break - end + +%% set volume scan +NbSlice_j=1; +index_series=1:size(filecell,2); +first_j_out=first_j; +last_j_out=last_j; +if isfield(Param,'ActionInput') && isfield(Param.ActionInput,'CheckVolume') ... + && Param.ActionInput.CheckVolume + index_j=Param.IndexRange.first_j:Param.IndexRange.incr_j:Param.IndexRange.last_j; + NbSlice_j=numel(index_j); + index_i=index_series(1:NbSlice_j:end); +end + +%%% loop on slices (volume scan) +for islice=index_j + if NbSlice_j>1 + first_j_out=islice; + last_j_out=islice; - %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%% - for iview=1:NbView - % reading input file(s) - [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index)); - if ~isempty(errormsg) - errormsg=['error of input reading: ' errormsg]; + end + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%% + for index=index_i+index_j(islice) + update_waitbar(WaitbarHandle,index/NbField) + if ~isempty(RUNHandle)&& ~strcmp(get(RUNHandle,'BusyAction'),'queue') + disp('program stopped by user') break end - if ~isempty(NbSlice_calib) - Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform - end - end - %%%%%%%%%%%%%%%% end loop on views (input lines) %%%%%%%%%%%%%%%% - %%%%%%%%%%%% END STANDARD PART %%%%%%%%%%%% - % EDIT FROM HERE - - if isempty(errormsg) - Field=Data{1}; % default input field structure - nbfiles=nbfiles+1; %increment the file counter - %% coordinate transform (or other user defined transform) - if ~isempty(transform_fct) - switch nargin(transform_fct) - case 4 - if length(Data)==2 - Field=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2}); + + %%%%%%%%%%%%%%%% loop on views (input lines) %%%%%%%%%%%%%%%% + for iview=1:NbView + % reading input file(s) + %InputFile=fullfile_uvmat(RootPath{iview},SubDir{iview},RootFile{iview},FileExt{iview},NomType{iview},first_i,last_i,first_j_out,last_j_out); + [Data{iview},tild,errormsg] = read_field(filecell{iview,index},FileType{iview},InputFields{iview},frame_index{iview}(index)); + if ~isempty(errormsg) + errormsg=['error of input reading: ' errormsg]; + break + end + if ~isempty(NbSlice_calib) + Data{iview}.ZIndex=mod(i1_series{iview}(index)-1,NbSlice_calib{iview})+1;%Zindex for phys transform + end + end + %%%%%%%%%%%%%%%% end loop on views (input lines) %%%%%%%%%%%%%%%% + + if isempty(errormsg) + Field=Data{1}; % default input field structure + nbfiles=nbfiles+1; %increment the file counter + %% coordinate transform (or other user defined transform) + if ~isempty(transform_fct) + switch nargin(transform_fct) + case 4 + if length(Data)==2 + Field=transform_fct(Data{1},XmlData{1},Data{2},XmlData{2}); + else + Field=transform_fct(Data{1},XmlData{1}); + end + case 3 + if length(Data)==2 + Field=transform_fct(Data{1},XmlData{1},Data{2}); + else + Field=transform_fct(Data{1},XmlData{1}); + end + case 2 + Field=transform_fct(Data{1},XmlData{1}); + case 1 + Field=transform_fct(Data{1}); + end + end + + %% field projection on an object + if Param.CheckObject + if strcmp(Param.ProjObject.ProjMode,'interp_tps') + Field=tps_coeff_field(Field,check_proj_tps);% calculate tps coefficients if needed + end + [Field,errormsg]=proj_field(Field,Param.ProjObject,VarMesh); + if ~isempty(errormsg) + disp_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg],checkrun) + return + end + end + + %%%%%%%%%%%% MAIN RUNNING OPERATIONS %%%%%%%%%%%% + if nbfiles==1 %first field + time_1=[]; + if isfield(Field,'Time') + time_1=Field.Time(1); + end + DataOut=Field;%outcome reproduces the first (projected) field by default + DataOut.Conventions='uvmat'; %suppress Conventions='uvmat/civdata' for civ input files + if isfield(Param,'ProjObject')&& ismember(Param.ProjObject.ProjMode,{'inside','outside'})%case of histograms + for ivar=1:numel(Field.ListVarName)% list of variable names before projection (histogram) + VarName=Field.ListVarName{ivar}; + if isfield(Data{1},VarName) + DataOut.(VarName)=Field.(VarName); + DataOut.([VarName 'Histo'])=zeros(size(DataOut.(VarName))); + VarMesh=DataOut.(VarName)(2)-DataOut.(VarName)(1); + end + end + disp(['mesh for histogram = ' num2str(VarMesh)]) + else + errorvar=zeros(numel(Field.ListVarName));%index of errorflag associated to each variable + if isfield(Field,'VarAttribute') + for ivar=1:numel(Field.ListVarName) + VarName=Field.ListVarName{ivar}; + DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero + NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero + + for iivar=1:length(Field.VarAttribute) + if isequal(Field.VarDimName{iivar},Field.VarDimName{ivar})&& isfield(Field.VarAttribute{iivar},'Role')... + && strcmp(Field.VarAttribute{iivar}.Role,'errorflag') + errorvar(ivar)=iivar; % index of the errorflag variable corresponding to ivar + end + end + end + DataOut.ListVarName(errorvar(errorvar~=0))=[]; %remove errorflag from result + DataOut.VarDimName(errorvar(errorvar~=0))=[]; %remove errorflag from result + DataOut.VarAttribute(errorvar(errorvar~=0))=[]; %remove errorflag from result else - Field=transform_fct(Data{1},XmlData{1}); + for ivar=1:numel(Field.ListVarName) + VarName=Field.ListVarName{ivar}; + DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero + NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero + end end - case 3 - if length(Data)==2 - Field=transform_fct(Data{1},XmlData{1},Data{2}); + + end + end %current field + for ivar=1:length(DataOut.ListVarName) + VarName=DataOut.ListVarName{ivar}; + sizmean=size(DataOut.(VarName)); + siz=size(Field.(VarName)); + if isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'}) + if isfield(Data{1},VarName) + MaxValue=max(DataOut.(VarName));% current max of histogram absissa + MinValue=min(DataOut.(VarName));% current min of histogram absissa + % VarMesh=Field.VarAttribute{ivar}.Mesh; + MaxIndex=round(MaxValue/VarMesh); + MinIndex=round(MinValue/VarMesh); + MaxIndex_new=round(max(Field.(VarName)/VarMesh));% max of the current field + MinIndex_new=round(min(Field.(VarName)/VarMesh)); + if MaxIndex_new>MaxIndex% the variable max for the current field exceeds the previous one + DataOut.(VarName)=[DataOut.(VarName) VarMesh*(MaxIndex+1:MaxIndex_new)];% append the new variable values + DataOut.([VarName 'Histo'])=[DataOut.([VarName 'Histo']) zeros(1,MaxIndex_new-MaxIndex)]; % append the new histo values + end + if MinIndex_new <= MinIndex-1 + DataOut.(VarName)=[VarMesh*(MinIndex_new:MinIndex-1) DataOut.(VarName)];% insert the new variable values + DataOut.([VarName 'Histo'])=[zeros(1,MinIndex-MinIndex_new) DataOut.([VarName 'Histo'])];% insert the new histo values + ind_start=1; + else + ind_start=MinIndex_new-MinIndex+1; + end + DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)=... + DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)+Field.([VarName 'Histo']); + end + elseif ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean) + disp_uvmat('ERROR',['unequal size of input field ' VarName ', need to project on a grid'],checkrun) + return + else + if errorvar(ivar)==0 + check_bad=isnan(Field.(VarName));%=0 for NaN data values, 1 else else - Field=transform_fct(Data{1},XmlData{1}); + check_bad=isnan(Field.(VarName)) | Field.(Field.ListVarName{errorvar(ivar)})~=0;%=0 for NaN or error flagged data values, 1 else end - case 2 - Field=transform_fct(Data{1},XmlData{1}); - case 1 - Field=transform_fct(Data{1}); - end - end - - %% field projection on an object - if Param.CheckObject - if strcmp(Param.ProjObject.ProjMode,'interp_tps') - Field=tps_coeff_field(Field,check_proj_tps);% calculate tps coefficients if needed - end - [Field,errormsg]=proj_field(Field,Param.ProjObject,VarMesh); - if ~isempty(errormsg) - disp_uvmat('ERROR',['error in aver_stat/proj_field:' errormsg],checkrun) - return - end - end - - %%%%%%%%%%%% MAIN RUNNING OPERATIONS %%%%%%%%%%%% - if nbfiles==1 %first field - time_1=[]; - if isfield(Field,'Time') - time_1=Field.Time(1); - end - DataOut=Field;%outcome reproduces the first (projected) field by default - DataOut.Conventions='uvmat'; %suppress Conventions='uvmat/civdata' for civ input files - if isfield(Param,'ProjObject')&& ismember(Param.ProjObject.ProjMode,{'inside','outside'})%case of histograms - for ivar=1:numel(Field.ListVarName)% list of variable names before projection (histogram) - VarName=Field.ListVarName{ivar}; - if isfield(Data{1},VarName) - DataOut.(VarName)=Field.(VarName); - DataOut.([VarName 'Histo'])=zeros(size(DataOut.(VarName))); - VarMesh=DataOut.(VarName)(2)-DataOut.(VarName)(1); - end - end - disp(['mesh for histogram = ' num2str(VarMesh)]) - else - errorvar=zeros(numel(Field.ListVarName));%index of errorflag associated to each variable - if isfield(Field,'VarAttribute') - for ivar=1:numel(Field.ListVarName) - VarName=Field.ListVarName{ivar}; - DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero - NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero - - for iivar=1:length(Field.VarAttribute) - if isequal(Field.VarDimName{iivar},Field.VarDimName{ivar})&& isfield(Field.VarAttribute{iivar},'Role')... - && strcmp(Field.VarAttribute{iivar}.Role,'errorflag') - errorvar(ivar)=iivar; % index of the errorflag variable corresponding to ivar - end - end - end - DataOut.ListVarName(errorvar(errorvar~=0))=[]; %remove errorflag from result - DataOut.VarDimName(errorvar(errorvar~=0))=[]; %remove errorflag from result - DataOut.VarAttribute(errorvar(errorvar~=0))=[]; %remove errorflag from result - else - for ivar=1:numel(Field.ListVarName) - VarName=Field.ListVarName{ivar}; - DataOut.(VarName)=zeros(size(DataOut.(VarName)));% initiate each field to zero - NbData.(VarName)=zeros(size(DataOut.(VarName)));% initiate the nbre of good data to zero - end - end - - end - end %current field - for ivar=1:length(DataOut.ListVarName) - VarName=DataOut.ListVarName{ivar}; - sizmean=size(DataOut.(VarName)); - siz=size(Field.(VarName)); - if isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'}) - if isfield(Data{1},VarName) - MaxValue=max(DataOut.(VarName));% current max of histogram absissa - MinValue=min(DataOut.(VarName));% current min of histogram absissa -% VarMesh=Field.VarAttribute{ivar}.Mesh; - MaxIndex=round(MaxValue/VarMesh); - MinIndex=round(MinValue/VarMesh); - MaxIndex_new=round(max(Field.(VarName)/VarMesh));% max of the current field - MinIndex_new=round(min(Field.(VarName)/VarMesh)); - if MaxIndex_new>MaxIndex% the variable max for the current field exceeds the previous one - DataOut.(VarName)=[DataOut.(VarName) VarMesh*(MaxIndex+1:MaxIndex_new)];% append the new variable values - DataOut.([VarName 'Histo'])=[DataOut.([VarName 'Histo']) zeros(1,MaxIndex_new-MaxIndex)]; % append the new histo values - end - if MinIndex_new <= MinIndex-1 - DataOut.(VarName)=[VarMesh*(MinIndex_new:MinIndex-1) DataOut.(VarName)];% insert the new variable values - DataOut.([VarName 'Histo'])=[zeros(1,MinIndex-MinIndex_new) DataOut.([VarName 'Histo'])];% insert the new histo values - ind_start=1; - else - ind_start=MinIndex_new-MinIndex+1; - end - DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)=... - DataOut.([VarName 'Histo'])(ind_start:ind_start+MaxIndex_new-MinIndex_new)+Field.([VarName 'Histo']); - end - elseif ~isequal(DataOut.(VarName),0)&& ~isequal(siz,sizmean) - disp_uvmat('ERROR',['unequal size of input field ' VarName ', need to project on a grid'],checkrun) - return - else - if errorvar(ivar)==0 - check_bad=isnan(Field.(VarName));%=0 for NaN data values, 1 else - else - check_bad=isnan(Field.(VarName)) | Field.(Field.ListVarName{errorvar(ivar)})~=0;%=0 for NaN or error flagged data values, 1 else - end - Field.(VarName)(check_bad)=0; %set to zero NaN or data marked by error flag - DataOut.(VarName)=DataOut.(VarName)+ double(Field.(VarName)); % update the sum - NbData.(VarName)=NbData.(VarName)+ ~check_bad;% records the number of data for each point - end - end - %%%%%%%%%%%% END MAIN RUNNING OPERATIONS %%%%%%%%%%%% - else - disp(errormsg) - end -end -%%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%% + Field.(VarName)(check_bad)=0; %set to zero NaN or data marked by error flag + DataOut.(VarName)=DataOut.(VarName)+ double(Field.(VarName)); % update the sum + NbData.(VarName)=NbData.(VarName)+ ~check_bad;% records the number of data for each point + end + end + %%%%%%%%%%%% END MAIN RUNNING OPERATIONS %%%%%%%%%%%% + else + disp(errormsg) + end + end + %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%% + if ~(isfield(Param,'ProjObject') && ismember(Param.ProjObject.ProjMode,{'inside','outside'})) for ivar=1:length(Field.ListVarName) @@ -426,5 +462,5 @@ %% writing the result file -OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j,last_j); +OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j_out,last_j_out); if strcmp(FileExtOut,'.png') %case of images if isequal(FileInfo{1}.BitDepth,16)||(numel(FileInfo)==2 &&isequal(FileInfo{2}.BitDepth,16)) @@ -444,4 +480,5 @@ end end % end averaging loop +end %% open the result file with uvmat (in RUN mode) Index: /trunk/src/series/merge_proj.m =================================================================== --- /trunk/src/series/merge_proj.m (revision 958) +++ /trunk/src/series/merge_proj.m (revision 959) @@ -103,5 +103,5 @@ WaitbarHandle=findobj(hseries,'Tag','Waitbar');%handle of waitbar in GUI series end -tild=phys([],[]);% test to provoke the inclusion of the function phys at compilation + %% define the directory for result file (with path=RootPath{1}) OutputDir=[Param.OutputSubDir Param.OutputDirExt];% subdirectory for output files