Index: /trunk/src/plot_field.m =================================================================== --- /trunk/src/plot_field.m (revision 1093) +++ /trunk/src/plot_field.m (revision 1094) @@ -768,12 +768,12 @@ test_interp_X=0; %default, regularly meshed X coordinate test_interp_Y=0; %default, regularly meshed Y coordinate - if isfield(Data,'VarAttribute') - if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end)},'units') - x_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end)}.units; - end - if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end-1) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)},'units') - y_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)}.units; - end - end +% if isfield(Data,'VarAttribute') +% if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end)},'units') +% x_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end)}.units; +% end +% if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end-1) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)},'units') +% y_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)}.units; +% end +% end if numel(Coord_y)>2 DCoord_y=diff(Coord_y); @@ -790,5 +790,5 @@ DCoord_x_min=min(DCoord_x); DCoord_x_max=max(DCoord_x); - if sign(DCoord_x_min)~=sign(DCoord_x_max);% =1 for increasing values, 0 otherwise + if sign(DCoord_x_min)~=sign(DCoord_x_max)% =1 for increasing values, 0 otherwise errormsg=['errror in plot_field.m: non monotonic dimension variable ' Data.ListVarName{VarRole.coord(2)} ]; return @@ -821,4 +821,7 @@ end %define coordinates as CoordUnits, if not defined as attribute for each variable +% if isfield(Data,'VarAttribute')&& numel(Data.VarAttribute)>=1 && isfield(Data.VarAttribute{1},'unit') +% y_units=Data.VarAttribute{1}.unit; +% end if isfield(Data,'CoordUnit') if isempty(x_units) @@ -827,4 +830,11 @@ if isempty(y_units) y_units=Data.CoordUnit; + end + elseif isfield(Data,'VarAttribute') + if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end)},'units') + x_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end)}.units; + end + if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end-1) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)},'units') + y_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)}.units; end end Index: /trunk/src/proj_field.m =================================================================== --- /trunk/src/proj_field.m (revision 1093) +++ /trunk/src/proj_field.m (revision 1094) @@ -1125,8 +1125,4 @@ ProjMode=num2cell(blanks(numel(CellInfo))); ProjMode=regexprep(ProjMode,' ',ObjectData.ProjMode); -%ProjMode=cell(size(CellInfo)); -% for icell=1:numel(CellInfo) -% ProjMode{icell}=ObjectData.ProjMode;% projection mode of the plane object -% end icell_grid=[];% field cell index which defines the grid if strcmp(ObjectData.ProjMode,'projection') @@ -1192,6 +1188,4 @@ [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates ProjData.VarDimName={AYName,AXName}; - % XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system - % YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3)); else% we use the existing grid from field cell #icell_grid NbDim=NbDimArray(icell_grid); @@ -1204,8 +1198,17 @@ ProjData.(AXName)=FieldData.(AXName); % new (projected ) y coordinates end - ProjData.ListVarName={AYName,AXName}; - + ProjData.ListVarName={AYName,AXName}; ProjData.VarAttribute{1}.Role='coord_y'; ProjData.VarAttribute{2}.Role='coord_x'; + YAttribute=FieldData.VarAttribute{CellInfo{icell_grid}.CoordIndex(NbDim-1)}; + XAttribute=FieldData.VarAttribute{CellInfo{icell_grid}.CoordIndex(NbDim)}; + if ~testangle + if isfield(YAttribute,'units') + ProjData.VarAttribute{1}.units=YAttribute.units; + end + if isfield(XAttribute,'units') + ProjData.VarAttribute{2}.units=XAttribute.units; + end + end end @@ -1244,5 +1247,5 @@ coord_x=FieldData.(CellInfo{icell}.XName);% initial x coordinates coord_y=FieldData.(CellInfo{icell}.YName);% initial y coordinates - + if check3D coord_z=FieldData.(CellInfo{icell}.ZName); @@ -1348,4 +1351,8 @@ VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag}; indsel=find(FieldData.(VarName_FF)==0); + if isempty(indsel) + errormsg='bad projection plane: no data found in the projection domain'; + return + end coord_X=coord_X(indsel); coord_Y=coord_Y(indsel); @@ -1625,13 +1632,5 @@ end [X,Y]=meshgrid(Coord{2},Coord{1});%initial coordinates - %name of error flag variable -% FFName='FF';%default name (if not already used) -% if isfield(ProjData,'FF') -% ind=1; -% while isfield(ProjData,['FF_' num2str(ind)]) -% ind=ind+1; -% end -% FFName=['FF_' num2str(ind)];% append an index to the name of error flag, FF_1,FF_2... -% end + % project all variables in the cell for ivar=VarIndex @@ -1656,11 +1655,5 @@ VarAttribute{length(ListVarName)+nbcoord}=FieldData.VarAttribute{ivar}; end -% ProjData.(FFName)=isnan(ProjData.(VarName));%detact NaN (points outside the interpolation range) -% ProjData.(VarName)(ProjData.(FFName))=0; %set to 0 the NaN data - end - %update list of variables with error flag -% ListVarName=[ListVarName FFName]; -% VarDimName=[VarDimName {DimCell}]; -% VarAttribute{numel(ListVarName)}.Role='errorflag'; + end elseif ~testangle % 3Dcase without change of angle % unstructured z coordinate @@ -1691,5 +1684,4 @@ Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});%initial y coordinates Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});%initial x coordinates - coord_x_proj=ObjectData.RangeX(1):InterpMesh:ObjectData.RangeX(2);% set of coordinates in the projection plane coord_y_proj=ObjectData.RangeY(1):InterpMesh:ObjectData.RangeY(2); @@ -1700,6 +1692,4 @@ YI_proj=M(2,1)*XI+M(2,2)*YI+ObjectData.Coord(1,2); ZI_proj=M(3,1)*XI+M(3,2)*YI+ObjectData.Coord(1,3); - - for ivar=VarIndex VarName=FieldData.ListVarName{ivar}; @@ -1763,26 +1753,4 @@ end end -% %prepare substraction in case of two input fields -% SubData.ListVarName={}; -% SubData.VarDimName={}; -% SubData.VarAttribute={}; -% check_remove=zeros(size(ProjData.ListVarName)); -% for iproj=1:numel(ProjData.VarAttribute) -% if isfield(ProjData.VarAttribute{iproj},'CheckSub')&&isequal(ProjData.VarAttribute{iproj}.CheckSub,1) -% VarName=ProjData.ListVarName{iproj}; -% SubData.ListVarName=[SubData.ListVarName {VarName}]; -% SubData.VarDimName=[SubData.VarDimName ProjData.VarDimName{iproj}]; -% SubData.VarAttribute=[SubData.VarAttribute ProjData.VarAttribute{iproj}]; -% SubData.(VarName)=ProjData.(VarName); -% check_remove(iproj)=1; -% end -% end -% if ~isempty(find(check_remove)) -% ind_remove=find(check_remove); -% ProjData.ListVarName(ind_remove)=[]; -% ProjData.VarDimName(ind_remove)=[]; -% ProjData.VarAttribute(ind_remove)=[]; -% ProjData=sub_field(ProjData,[],SubData); -% end %----------------------------------------------------------------- Index: /trunk/src/series/merge_proj.m =================================================================== --- /trunk/src/series/merge_proj.m (revision 1093) +++ /trunk/src/series/merge_proj.m (revision 1094) @@ -253,4 +253,5 @@ end OutputPath=fullfile(Param.OutputPath,Param.Experiment,Param.Device); + for index=1:NbField update_waitbar(WaitbarHandle,index/NbField) Index: /trunk/src/series/sliding_average.m =================================================================== --- /trunk/src/series/sliding_average.m (revision 1093) +++ /trunk/src/series/sliding_average.m (revision 1094) @@ -25,7 +25,7 @@ % .RUN =0 for GUI input, =1 for function activation % .RunMode='local','background', 'cluster': type of function use -% -% .IndexRange: set the file or frame indices on which the action must be performed -% .FieldTransform: .TransformName: name of the selected transform function +% 900 +% .IndexRange: set the file or frame indices on which the action must be performseriesed +% .FieldTransform: .TransformName: name of the select39ed transform function % .TransformPath: path of the selected transform function % .InputFields: sub structure describing the input fields withfields @@ -34,9 +34,8 @@ % .FieldName_1: name of the second field in case of two input series % .VelType_1: velocity type of the second field in case of two input series -% .Coord_y: name of y coordinate variable +% uvmat .Coord_y: name of y coordinate variable % .Coord_x: name of x coordinate variable % .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% - %======================================================================= % Copyright 2008-2021, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France @@ -62,10 +61,10 @@ if isstruct(Param) && isequal(Param.Action.RUN,0) ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default) - ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default) + ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to mseriesax (options 'off'/'on', 'off' by default) ParamOut.NbSlice=1; %nbre of slices ('off' by default) 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.ProjObject='off';%can use projection object(option 'off'/'on', + 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 @@ -122,5 +121,5 @@ return end - [FileInfo{iview},MovieObject{iview}]=get_file_info(filecell{iview,1}); + [FileInfo{iview},MovieObject{iview}]=get_file_info(filecell{iview,1});900 FileType{iview}=FileInfo{iview}.FileType; CheckImage{iview}=~isempty(find(strcmp(FileType{iview},ImageTypeOptions)));% =1 for images @@ -138,5 +137,5 @@ if size(time,1)>1 diff_time=max(max(diff(time))); - if diff_time>0 + if diff_time>0series msgbox_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time)]) end @@ -177,5 +176,5 @@ %% Set field names and velocity types -InputFields{1}=[];%default (case of images) +InputFields{1}=[];%default (case of images)series if isfield(Param,'InputFields') InputFields{1}=Param.InputFields; @@ -187,12 +186,18 @@ %% initialisation T=24.2; %main wave period +t0=3; % time for motion start (torus at its maximum x) NbPeriod=2; %number of periods for the sliding average omega=2*pi/T; +amplitude=2.5; %oscillation amplitude +Lscale=15;%diameter of the torus, length scale for normalisation +Uscale=amplitude*omega; DataOut.ListGlobalAttribute= {'Conventions','Time'}; DataOut.Conventions='uvmat'; -DataOut.ListVarName={'coord_y','coord_x','Umean','Vmean','Ucos','Vcos','DUDXsin','DUDYsin','DVDXsin','DVDYsin','Ustokes','Vstokes'}; +DataOut.ListVarName={'coord_y','coord_x','Umean','Vmean','Ucos','Vcos','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'},{'coord_y','coord_x'},... + {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}}; %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%% @@ -207,6 +212,9 @@ Time_end=Data.Time; dt=(Time_end-Time_1)/(NbField-1); %time interval -NpTime=round(NbPeriod*T/dt+1) - +NpTime=round(NbPeriod*T/dt+1); + +OutputPath=fullfile(Param.OutputPath,Param.Experiment,Param.Device); +RootFileOut=RootFile{1}; +NomTypeOut='_1'; %%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%% disp('loop for filtering started') @@ -219,8 +227,8 @@ [Field,tild,errormsg] = read_field(filecell{1,index},FileType{iview},InputFields{iview},frame_index{iview}(index)); - %%%%%%%%%%%% MAIN RUNNING OPERATIONS %%%%%%%%%%%% + %%%%%%%%%%% MAIN RUNNING OPERATIONS %%%%%%%%%%%% if index==1 %first field - DataOut.coord_x=Field.coord_x; - DataOut.coord_y=Field.coord_y; + DataOut.coord_x=Field.coord_x/Lscale; + DataOut.coord_y=Field.coord_y/Lscale; npy=numel(DataOut.coord_y); npx=numel(DataOut.coord_x); @@ -229,49 +237,56 @@ Ucos=zeros(NpTime,npy,npx); Vcos=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; + 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 - DUDXsin=circshift(DUDXsin,[-1 0 0]); - DUDYsin=circshift(DUDYsin,[-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)); - DUDXsin(end,:,:)=Field.DUDX*sin(omega*Time(index)); - DUDYsin(end,:,:)=Field.DUDY*sin(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; - DataOut.Umean=squeeze(nanmean(Umean,1)); - DataOut.Vmean=squeeze(nanmean(Vmean,1)); - DataOut.Ucos=2*squeeze(nanmean(Ucos,1)); - DataOut.Vcos=2*squeeze(nanmean(Vcos,1)); - DataOut.DUDXsin=2*squeeze(nanmean(DUDXsin,1)); - DataOut.DUDYsin=2*squeeze(nanmean(DUDYsin,1)); - DataOut.DVDXsin=2*squeeze(nanmean(DVDXsin,1)); - DataOut.DVDYsin=2*squeeze(nanmean(DVDYsin,1)); - DataOut.Ustokes=(1/omega)*(DataOut.Ucos.*DataOut.DUDXsin+DataOut.Vcos.*DataOut.DUDYsin); - DataOut.Vstokes=(1/omega)*(DataOut.Ucos.*DataOut.DVDXsin+DataOut.Vcos.*DataOut.DVDYsin); + 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 + 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)); + 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*squeeze(nanmean(Ucos,1)); + DataOut.Vcos=2*squeeze(nanmean(Vcos,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); + % writing the result file as netcdf file - if index-round(NpTime/2)>=1 - OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomType{1},index-round(NpTime/2)); - %case of netcdf input file , determine global attributes - errormsg=struct2nc(OutputFile,DataOut); %save result file - if isempty(errormsg) - disp([OutputFile ' written']); - else - disp(['error in writting result file: ' errormsg]) - end - end -end -%%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%% - + i1=i1_series{1}(index); + OutputFile=fullfile_uvmat(OutputPath,OutputDir,RootFileOut,'.nc',NomTypeOut,i1); + errormsg=struct2nc(OutputFile, DataOut); + if isempty(errormsg) + disp([OutputFile ' written']) + else + disp(errormsg) + end +end + Index: /trunk/src/transform_field/phys_polar.m =================================================================== --- /trunk/src/transform_field/phys_polar.m (revision 1093) +++ /trunk/src/transform_field/phys_polar.m (revision 1094) @@ -52,5 +52,5 @@ dlg_title = 'set the parameters for the polar coordinates'; num_lines= 2; - def = { '[0 0]';'0';'0';'+'}; + def = { '[0 0]';'';'0';'+'}; if isfield(XmlData,'TransformInput') if isfield(XmlData.TransformInput,'PolarCentre') @@ -87,5 +87,5 @@ DataCell{2}=[];%default checkpixel(1)=0; -if isfield(DataCell{1},'CoorUnit')&& strcmp(DataCell{1}.CoorUnit,'px') +if isfield(DataCell{1},'CoordUnit')&& strcmp(DataCell{1}.CoordUnit,'pixel') checkpixel(1)=1; end @@ -101,5 +101,5 @@ if nargin==4% case of two input fields checkpixel(2)=0; -if isfield(DataCell{2},'CoorUnit')&& strcmp(DataCell{2}.CoorUnit,'px') +if isfield(DataCell{2},'CoordUnit')&& strcmp(DataCell{2}.CoordUnit,'pixel') checkpixel(2)=1; end @@ -124,10 +124,10 @@ end end - if isfield(XmlData.TransformInput,'PolarReferenceRadius') && isnumeric(XmlData.TransformInput.PolarReferenceRadius) + if isfield(XmlData.TransformInput,'PolarReferenceRadius') && ~isempty(XmlData.TransformInput.PolarReferenceRadius) radius_offset=XmlData.TransformInput.PolarReferenceRadius; end if radius_offset > 0 angle_scale=radius_offset; %the azimuth is rescale in terms of the length along the reference radius - check_degree=0; %the output has the same unit asthe input + check_degree=0; %the output has the same unit as the input else angle_scale=180/pi; %polar angle in degrees @@ -144,5 +144,5 @@ nbvar=0;%counter for the number of output variables -nbcoord=0;%counter for the number of variables for radial coordiantes (case of multiple field inputs) +nbcoord=0;%counter for the number of variablecheck_degrees for radial coordiantes (case of multiple field inputs) nbgrid=0;%counter for the number of gridded fields (all linearly interpolated on the same output polar grid) nbscattered=0;%counter of scattered fields @@ -157,5 +157,4 @@ return end - % end %transform of X,Y coordinates for vector fields if isfield(DataCell{ifield},'ZIndex')&& ~isempty(DataCell{ifield}.ZIndex) @@ -180,18 +179,33 @@ Data.VarAttribute{nbvar-1}.Role='coord_x'; check_unit=1; - if isfield(DataCell{ifield},'CoordUnit') - Data=rmfield(Data,'CoordUnit'); - Data.VarAttribute{nbvar-1}.unit=DataCell{ifield}.CoordUnit; - elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit') - Data.VarAttribute{nbvar-1}.unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units + %unit of output field + if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit') + radius_unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units + elseif isfield(DataCell{ifield},'CoordUnit') + radius_unit=DataCell{ifield}.CoordUnit; else - check_unit=0; - end + radius_unit=''; + end + Data.VarAttribute{nbvar-1}.units=radius_unit; + if check_degree + Data.VarAttribute{nbvar}.units='degree'; + else %case of a reference radius + Data.VarAttribute{nbvar}.units=radius_unit; + Data.CoordUnit=radius_unit; + end +% if isfield(DataCell{ifield},'CoordUnit') +% Data=rmfield(Data,'CoordUnit'); +% Data.VarAttribute{nbvar-1}.unit=DataCell{ifield}.CoordUnit; +% elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit') +% Data.VarAttribute{nbvar-1}.unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units +% else +% check_unit=0; +% end Data.VarAttribute{nbvar}.Role='coord_y'; - if check_degree - Data.VarAttribute{nbvar}.unit='degree'; - elseif check_unit - Data.VarAttribute{nbvar}.unit=Data.VarAttribute{nbvar-1}.unit; - end +% if check_degree +% Data.VarAttribute{nbvar}.units='degree'; +% elseif check_unit +% Data.VarAttribute{nbvar}.units=Data.VarAttribute{nbvar-1}.units; +% end %transform u,v into polar coordinates @@ -250,28 +264,29 @@ if nbgrid==0% no gridded data yet, introduce the coordinate variables common to all gridded data nbcoord=nbcoord+1;%add new radial coordinates for the first gridded field - radius_name = rename_indexing(radius_name,Data.ListVarName); - theta_name = rename_indexing(theta_name,Data.ListVarName); - Data.ListVarName = [Data.ListVarName {radius_name} {theta_name}]; + radius_name = rename_indexing(radius_name,Data.ListVarName);% add an index to the name, or increment an existing index, + theta_name = rename_indexing(theta_name,Data.ListVarName);% if the proposed Name already exists in the list + Data.ListVarName = [Data.ListVarName {radius_name} {theta_name}];%add polar coordinates to the list of variables Data.VarDimName=[Data.VarDimName {radius_name} {theta_name}]; nbvar=nbvar+2; if check_reverse - Data.VarAttribute{nbvar-1}.Role='coord_y'; - Data.VarAttribute{nbvar}.Role='coord_x'; + Data.VarAttribute{nbvar-1}.Role='coord_y'; + Data.VarAttribute{nbvar}.Role='coord_x'; else - Data.VarAttribute{nbvar-1}.Role='coord_x'; - Data.VarAttribute{nbvar}.Role='coord_y'; + Data.VarAttribute{nbvar-1}.Role='coord_x'; + Data.VarAttribute{nbvar}.Role='coord_y'; end check_unit=1; - if isfield(DataCell{ifield},'CoordUnit') - Data.VarAttribute{nbvar-1}.unit=DataCell{ifield}.CoordUnit; - elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit') - Data.VarAttribute{nbvar-1}.unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units + + if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit') + Data.VarAttribute{nbvar-1}.units=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units + elseif isfield(DataCell{ifield},'CoordUnit') + Data.VarAttribute{nbvar-1}.units=DataCell{ifield}.CoordUnit;%radius in coord units else check_unit=0; end if check_degree - Data.VarAttribute{nbvar}.unit='degree'; + Data.VarAttribute{nbvar}.units='degree';%angle in degree elseif check_unit - Data.VarAttribute{nbvar}.unit=Data.VarAttribute{nbvar-1}.unit; + Data.VarAttribute{nbvar}.units=Data.VarAttribute{nbvar-1}.units;% angle in coord unit (normalised by reference radiuss) end end @@ -282,6 +297,4 @@ FieldName{nbgrid}=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_scalar}; A{nbgrid}=DataCell{ifield}.(FieldName{nbgrid}); -% Data.ListVarName=[Data.ListVarName {FieldName{nbgrid}}]; -% Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}}]; nbpoint(nbgrid)=numel(A{nbgrid}); check_scalar(nbgrid)=1; @@ -296,5 +309,5 @@ A{nbgrid+1}=DataCell{ifield}.(FieldName{nbgrid+1}); A{nbgrid+2}=DataCell{ifield}.(FieldName{nbgrid+2}); - % Data.ListVarName=[Data.ListVarName {'U_r','U_theta'}]; + % Data.ListVarName=[Data.ListVarName {'U_r','U_theta'}]; %Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}} {{theta_name,radius_name}}]; Data.VarAttribute{nbvar+1}.Role='vector_x'; Index: /trunk/src/uvmat.m =================================================================== --- /trunk/src/uvmat.m (revision 1093) +++ /trunk/src/uvmat.m (revision 1094) @@ -1538,4 +1538,8 @@ XmlData.LIFCalib.Ray1Coord=LineData{1}.Coord; XmlData.LIFCalib.Ray2Coord=LineData{2}.Coord; +if numel(LineData)<3 + msgbox_uvmat('ERROR','draw a reference line of direct laser illumination (without dye absorbsion)'); + return +end XmlData.LIFCalib.RefLineCoord=LineData{3}.Coord; @@ -1545,6 +1549,6 @@ y=linspace(UvData.Field.Coord_y(1),UvData.Field.Coord_y(2),nby)-nby/2; [X,Y]=meshgrid(x,y); -coeff_quad=0.15*4/(nbx*nbx);% image luminosity reduced by 10% at the edge -UvData.Field.A=double(UvData.Field.A).*(1+coeff_quad*(X.*X+Y.*Y)); +%coeff_quad=0.15*4/(nbx*nbx);% image luminosity reduced by 10% at the edge +%UvData.Field.A=double(UvData.Field.A).*(1+coeff_quad*(X.*X+Y.*Y)); %% display the current image in polar axes with origin at the illumination source @@ -1564,5 +1568,5 @@ y_ref=y_ref-y0; [theta_ref,r_ref] = cart2pol(x_ref,y_ref);%theta_ref and r_ref are the polar coordinates of the points on the line - theta_ref=theta_ref*180/pi;% theta_ref in radians + theta_ref=theta_ref*180/pi;% theta_ref in degrees figure(10) plot(theta_ref,r_ref) @@ -1570,5 +1574,5 @@ ylabel('radius from light source') title('ref line in polar coordinates') - azimuth_ima=linspace(DataPol.Coord_y(1),DataPol.Coord_y(2),size(DataPol.A,1));%array of angular index on the transformed image + azimuth_ima=linspace(DataPol.Coord_y(1),DataPol.Coord_y(2),size(DataPol.A,1))-360;%array of angular indices on the transformed image dist_source = interp1(theta_ref,r_ref,azimuth_ima);% get the polar position of the reference line dist_source_pixel=round(size(DataPol.A,2)*(dist_source-DataPol.Coord_x(1))/(DataPol.Coord_x(2)-DataPol.Coord_x(1))); @@ -3719,5 +3723,8 @@ %% choose and read a second field FileName_1 if defined ParamOut_1=[]; -if numel(UvData.FileInfo)>1 +if ~isempty(FileName_1) + if numel(UvData.FileInfo)==1 + UvData.FileInfo{2}=UvData.FileInfo{1}; + end VelType_1=[];%default FieldName_1=[]; @@ -3736,5 +3743,5 @@ NomType_1=get(handles.NomType,'String'); end - if strcmp(UvData.FileInfo{2}.FieldType,'image') + if strcmp(UvData.FileInfo{2}.FileType,'image') FieldName_1='image'; frame_index_1=1;%default @@ -3761,5 +3768,5 @@ end switch UvData.FileInfo{2}.FileType - case {'civx','civdata','netcdf','pivdata_fluidimage'}; + case {'civx','civdata','netcdf','pivdata_fluidimage'} list_fields=get(handles.FieldName_1,'String');% list menu fields FieldName_1= list_fields{get(handles.FieldName_1,'Value')}; % selected field @@ -3849,7 +3856,7 @@ %% update the display menu for the second velocity type (second menuline) test_veltype_1=0; -if isempty(FileName_1) -elseif ~test_keepdata_1 - if strcmp(UvData.FileInfo{2}.FieldType,'civdata')&& ~strcmp(FieldName_1,'get_field...') +if ~isempty(FileName_1) + + if strcmp(UvData.FileInfo{2}.FileType,'civdata')&& ~strcmp(FieldName_1,'get_field...') test_veltype_1=1; set(handles.VelType_1,'Visible','on')