Index: /trunk/src/get_field.m =================================================================== --- /trunk/src/get_field.m (revision 1081) +++ /trunk/src/get_field.m (revision 1082) @@ -200,6 +200,6 @@ set(handles.scalar,'Value',1)% fill the menu of y vector components set(handles.scalar,'String',Field.Display.ListVarName)% fill the menu for scalar -set(handles.ordinate,'Value',1)% fill the menu of y vector components -set(handles.ordinate,'String',Field.Display.ListVarName)% fill the menu of y coordinate for 1D plots +%set(handles.ordinate,'Value',1)% fill the menu of y vector components +%set(handles.ordinate,'String',Field.Display.ListVarName)% fill the menu of y coordinate for 1D plots checkseries=0; if isfield(ParamIn,'SeriesInput') && ParamIn.SeriesInput% case of call by series @@ -405,10 +405,10 @@ case '1D plot' set(handles.Coordinates,'Visible','on') - set(handles.PanelOrdinate,'Visible','on') - pos=get(handles.PanelOrdinate,'Position'); - pos(1)=2; - pos_coord=get(handles.Coordinates,'Position'); - pos(2)=pos_coord(2)-pos(4)-2; - set(handles.PanelOrdinate,'Position',pos) + %set(handles.PanelOrdinate,'Visible','on') + %pos=get(handles.PanelOrdinate,'Position'); +% pos(1)=2; +% pos_coord=get(handles.Coordinates,'Position'); +% pos(2)=pos_coord(2)-pos(4)-2; + %set(handles.PanelOrdinate,'Position',pos) set(handles.PanelScalar,'Visible','off') set(handles.PanelVectors,'Visible','off') @@ -420,5 +420,5 @@ case {'scalar'} set(handles.Coordinates,'Visible','on') - set(handles.PanelOrdinate,'Visible','off') + %set(handles.PanelOrdinate,'Visible','off') set(handles.PanelScalar,'Visible','on') set(handles.PanelVectors,'Visible','off') @@ -456,5 +456,5 @@ set(handles.PanelVectors,'Visible','on') set(handles.Coordinates,'Visible','on') - set(handles.PanelOrdinate,'Visible','off') + %set(handles.PanelOrdinate,'Visible','off') set(handles.PanelScalar,'Visible','off') pos=get(handles.PanelVectors,'Position'); @@ -496,5 +496,5 @@ vector_Callback(handles) case 'civdata...' - set(handles.PanelOrdinate,'Visible','off') + %set(handles.PanelOrdinate,'Visible','off') set(handles.PanelScalar,'Visible','off') set(handles.PanelVectors,'Visible','off') @@ -502,5 +502,5 @@ end -%------------------------------------------------------------------------ +%NOT USED : TO DELETE------------------------------------------------------------------------ function ordinate_Callback(hObject, DimCell, handles) %------------------------------------------------------------------------ @@ -629,10 +629,10 @@ %% set default selection for grid coordinates -if ~isempty(var_coord) -coord_val(1)=var_coord(end); -coord_val(2)=var_coord(end-1); -if numel(var_coord)>=3 - coord_val(3)=var_coord(end-2); -end +if numel(var_coord)>=2 + coord_val(1)=var_coord(end); + coord_val(2)=var_coord(end-1); + if numel(var_coord)>=3 + coord_val(3)=var_coord(end-2); + end end % if numel(find(test_coord))>3 Index: /trunk/src/nc2struct.m =================================================================== --- /trunk/src/nc2struct.m (revision 1081) +++ /trunk/src/nc2struct.m (revision 1082) @@ -62,5 +62,5 @@ hhh=which('netcdf.open');% look for built-in matlab NetCDF library -if ~isequal(hhh,'') +if ~isempty(hhh) %% default output Data=[];%default @@ -113,5 +113,5 @@ CheckTimeVar=1; TimeIndex=varargin{3}; - input_index=4;% list of varibles to read is at fourth argument + input_index=4;% list of variables to read is at fourth argument elseif isequal(varargin{1},'TimeDimName') TimeDimName=varargin{2}; Index: /trunk/src/read_GUI.m =================================================================== --- /trunk/src/read_GUI.m (revision 1081) +++ /trunk/src/read_GUI.m (revision 1082) @@ -71,6 +71,6 @@ end case {'listbox','popupmenu'} - listinput=get(hchild(ichild),'String') - value=get(hchild(ichild),'Value') + listinput=get(hchild(ichild),'String'); + value=get(hchild(ichild),'Value'); if ~isempty(listinput) if numel(value)==1% single selection Index: /trunk/src/rodrigues.m =================================================================== --- /trunk/src/rodrigues.m (revision 1082) +++ /trunk/src/rodrigues.m (revision 1082) @@ -0,0 +1,273 @@ +function [out,dout]=rodrigues(in) + +% RODRIGUES Transform rotation matrix into rotation vector and viceversa. +% +% Sintax: [OUT]=RODRIGUES(IN) +% If IN is a 3x3 rotation matrix then OUT is the +% corresponding 3x1 rotation vector +% if IN is a rotation 3-vector then OUT is the +% corresponding 3x3 rotation matrix +% + +%% +%% Copyright (c) March 1993 -- Pietro Perona +%% California Institute of Technology +%% + +%% ALL CHECKED BY JEAN-YVES BOUGUET, October 1995. +%% FOR ALL JACOBIAN MATRICES !!! LOOK AT THE TEST AT THE END !! + +%% BUG when norm(om)=pi fixed -- April 6th, 1997; +%% Jean-Yves Bouguet + +%% Add projection of the 3x3 matrix onto the set of special ortogonal matrices SO(3) by SVD -- February 7th, 2003; +%% Jean-Yves Bouguet + +% BUG FOR THE CASE norm(om)=pi fixed by Mike Burl on Feb 27, 2007 + + +[m,n] = size(in); +%bigeps = 10e+4*eps; +bigeps = 10e+20*eps; + +if ((m==1) & (n==3)) | ((m==3) & (n==1)) %% it is a rotation vector + theta = norm(in); + if theta < eps + R = eye(3); + + %if nargout > 1, + + dRdin = [0 0 0; + 0 0 1; + 0 -1 0; + 0 0 -1; + 0 0 0; + 1 0 0; + 0 1 0; + -1 0 0; + 0 0 0]; + + %end; + + else + if n==length(in) in=in'; end; %% make it a column vec. if necess. + + %m3 = [in,theta] + + dm3din = [eye(3);in'/theta]; + + omega = in/theta; + + %m2 = [omega;theta] + + dm2dm3 = [eye(3)/theta -in/theta^2; zeros(1,3) 1]; + + alpha = cos(theta); + beta = sin(theta); + gamma = 1-cos(theta); + omegav=[[0 -omega(3) omega(2)];[omega(3) 0 -omega(1)];[-omega(2) omega(1) 0 ]]; + A = omega*omega'; + + %m1 = [alpha;beta;gamma;omegav;A]; + + dm1dm2 = zeros(21,4); + dm1dm2(1,4) = -sin(theta); + dm1dm2(2,4) = cos(theta); + dm1dm2(3,4) = sin(theta); + dm1dm2(4:12,1:3) = [0 0 0 0 0 1 0 -1 0; + 0 0 -1 0 0 0 1 0 0; + 0 1 0 -1 0 0 0 0 0]'; + + w1 = omega(1); + w2 = omega(2); + w3 = omega(3); + + dm1dm2(13:21,1) = [2*w1;w2;w3;w2;0;0;w3;0;0]; + dm1dm2(13: 21,2) = [0;w1;0;w1;2*w2;w3;0;w3;0]; + dm1dm2(13:21,3) = [0;0;w1;0;0;w2;w1;w2;2*w3]; + + R = eye(3)*alpha + omegav*beta + A*gamma; + + dRdm1 = zeros(9,21); + + dRdm1([1 5 9],1) = ones(3,1); + dRdm1(:,2) = omegav(:); + dRdm1(:,4:12) = beta*eye(9); + dRdm1(:,3) = A(:); + dRdm1(:,13:21) = gamma*eye(9); + + dRdin = dRdm1 * dm1dm2 * dm2dm3 * dm3din; + + + end; + out = R; + dout = dRdin; + + %% it is prob. a rot matr. +elseif ((m==n) & (m==3) & (norm(in' * in - eye(3)) < bigeps)... + & (abs(det(in)-1) < bigeps)) + R = in; + + % project the rotation matrix to SO(3); + [U,S,V] = svd(R); + R = U*V'; + + tr = (trace(R)-1)/2; + dtrdR = [1 0 0 0 1 0 0 0 1]/2; + theta = real(acos(tr)); + + + if sin(theta) >= 1e-4, + + dthetadtr = -1/sqrt(1-tr^2); + + dthetadR = dthetadtr * dtrdR; + % var1 = [vth;theta]; + vth = 1/(2*sin(theta)); + dvthdtheta = -vth*cos(theta)/sin(theta); + dvar1dtheta = [dvthdtheta;1]; + + dvar1dR = dvar1dtheta * dthetadR; + + + om1 = [R(3,2)-R(2,3), R(1,3)-R(3,1), R(2,1)-R(1,2)]'; + + dom1dR = [0 0 0 0 0 1 0 -1 0; + 0 0 -1 0 0 0 1 0 0; + 0 1 0 -1 0 0 0 0 0]; + + % var = [om1;vth;theta]; + dvardR = [dom1dR;dvar1dR]; + + % var2 = [om;theta]; + om = vth*om1; + domdvar = [vth*eye(3) om1 zeros(3,1)]; + dthetadvar = [0 0 0 0 1]; + dvar2dvar = [domdvar;dthetadvar]; + + + out = om*theta; + domegadvar2 = [theta*eye(3) om]; + + dout = domegadvar2 * dvar2dvar * dvardR; + + + else + if tr > 0; % case norm(om)=0; + + out = [0 0 0]'; + + dout = [0 0 0 0 0 1/2 0 -1/2 0; + 0 0 -1/2 0 0 0 1/2 0 0; + 0 1/2 0 -1/2 0 0 0 0 0]; + else + + % case norm(om)=pi; + if(0) + + %% fixed April 6th by Bouguet -- not working in all cases! + out = theta * (sqrt((diag(R)+1)/2).*[1;2*(R(1,2:3)>=0)'-1]); + %keyboard; + + else + + % Solution by Mike Burl on Feb 27, 2007 + % This is a better way to determine the signs of the + % entries of the rotation vector using a hash table on all + % the combinations of signs of a pairs of products (in the + % rotation matrix) + + % Define hashvec and Smat + hashvec = [0; -1; -3; -9; 9; 3; 1; 13; 5; -7; -11]; + Smat = [1,1,1; 1,0,-1; 0,1,-1; 1,-1,0; 1,1,0; 0,1,1; 1,0,1; 1,1,1; 1,1,-1; + 1,-1,-1; 1,-1,1]; + + M = (R+eye(3,3))/2; + uabs = sqrt(M(1,1)); + vabs = sqrt(M(2,2)); + wabs = sqrt(M(3,3)); + + mvec = ([M(1,2), M(2,3), M(1,3)] + [M(2,1), M(3,2), M(3,1)])/2; + syn = ((mvec > eps) - (mvec < -eps)); % robust sign() function + hash = syn * [9; 3; 1]; + idx = find(hash == hashvec); + svec = Smat(idx,:)'; + + out = theta * [uabs; vabs; wabs] .* svec; + + end; + + if nargout > 1, + fprintf(1,'WARNING!!!! Jacobian domdR undefined!!!\n'); + dout = NaN*ones(3,9); + end; + end; + end; + +else + error('Neither a rotation matrix nor a rotation vector were provided'); +end; + +return; + +%% test of the Jacobians: + +%% TEST OF dRdom: +om = randn(3,1); +dom = randn(3,1)/1000000; +[R1,dR1] = rodrigues(om); +R2 = rodrigues(om+dom); +R2a = R1 + reshape(dR1 * dom,3,3); +gain = norm(R2 - R1)/norm(R2 - R2a) + +%% TEST OF dOmdR: +om = randn(3,1); +R = rodrigues(om); +dom = randn(3,1)/10000; +dR = rodrigues(om+dom) - R; + +[omc,domdR] = rodrigues(R); +[om2] = rodrigues(R+dR); +om_app = omc + domdR*dR(:); +gain = norm(om2 - omc)/norm(om2 - om_app) + + +%% OTHER BUG: (FIXED NOW!!!) +omu = randn(3,1); +omu = omu/norm(omu) +om = pi*omu; +[R,dR]= rodrigues(om); +[om2] = rodrigues(R); +[om om2] + +%% NORMAL OPERATION +om = randn(3,1); +[R,dR]= rodrigues(om); +[om2] = rodrigues(R); +[om om2] + +%% Test: norm(om) = pi +u = randn(3,1); +u = u / sqrt(sum(u.^2)); +om = pi*u; +R = rodrigues(om); +R2 = rodrigues(rodrigues(R)); +norm(R - R2) + +%% Another test case where norm(om)=pi from Chen Feng (June 27th, 2014) +R = [-0.950146567583153 -6.41765854280073e-05 0.311803617668748; ... + -6.41765854277654e-05 -0.999999917385145 -0.000401386434914383; ... + 0.311803617668748 -0.000401386434914345 0.950146484968298]; +om = rodrigues(R) +norm(om) - pi + +%% Another test case where norm(om)=pi from 䜙成义 (July 1st, 2014) +R = [-0.999920129411407 -6.68593208347372e-05 -0.0126384464118876; ... + 9.53007036072085e-05 -0.999997464662094 -0.00224979713751896; ... + -0.0126382639492467 -0.00225082189773293 0.999917600647740]; +om = rodrigues(R) +norm(om) - pi + + + + Index: /trunk/src/series/merge_proj_polar_multifile.m =================================================================== --- /trunk/src/series/merge_proj_polar_multifile.m (revision 1081) +++ /trunk/src/series/merge_proj_polar_multifile.m (revision 1082) @@ -270,5 +270,5 @@ for index=1:NbField disp(['index=' num2str(index)]) - disp(['ellapsed time ' num2str(toc(tstart)/60,4) ' minutes']) + %disp(['ellapsed time ' num2str(toc(tstart)/60,4) ' minutes']) update_waitbar(WaitbarHandle,index/NbField) if ~isempty(RUNHandle) && ~strcmp(get(RUNHandle,'BusyAction'),'queue') @@ -374,12 +374,12 @@ ProjData{iview}.VarDimName={'radius','azimuth'}; ProjData{iview}.VarAttribute=[{'coord_x'} {'coord_y'} VarAttribute]; - for ivar=1:numel(ListVarName) - ProjData{iview}.VarDimName{ivar+2}={'radius','azimuth'}; - VarName=ListVarName{ivar}; - if ~isempty(thresh2) - DataOut.(VarName)(Dist>thresh2)=NaN;% put to NaN interpolated positions further than RangeInterp from initial data - end - ProjData{iview}.(VarName)=(DataOut.(VarName))'; - end +% for ivar=1:numel(ListVarName) +% ProjData{iview}.VarDimName{ivar+2}={'radius','azimuth'}; +% VarName=ListVarName{ivar}; +% if ~isempty(thresh2) +% DataOut.(VarName)(Dist>thresh2)=NaN;% put to NaN interpolated positions further than RangeInterp from initial data +% end +% ProjData{iview}.(VarName)=(DataOut.(VarName))'; +% end end Index: /trunk/src/series/turb_stat.m =================================================================== --- /trunk/src/series/turb_stat.m (revision 1081) +++ /trunk/src/series/turb_stat.m (revision 1082) @@ -63,5 +63,5 @@ 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.NbSlice=1; %nbre of slices ('off' by default) + ParamOut.NbSlice='on'; %nbre of slices ('on' if needed as input, fixed value e.g. 1, 'off' by default) ParamOut.VelType='one';% menu for selecting the velocity type (options 'off'/'one'/'two', 'off' by default) ParamOut.FieldName='off';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default) @@ -179,7 +179,6 @@ DataOut.ListGlobalAttribute= {'Conventions'}; DataOut.Conventions= 'uvmat'; -DataOut.ListVarName={'coord_y', 'coord_x' ,'UMean' , 'VMean','u2Mean','v2Mean','u2Mean_1','v2Mean_1','uvMean','CurlMean','DivMean','Curl2Mean','Div2Mean','Counter'}; -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'}}; +DataOut.ListVarName={}; +DataOut.VarDimName={}; DataOut.UMean=0; DataOut.VMean=0; @@ -213,7 +212,19 @@ %%%%%%%%%%%% MAIN RUNNING OPERATIONS %%%%%%%%%%%% if index==1 %initiate the output data structure in the first field - - DataOut.coord_y=Field.coord_y; - DataOut.coord_x=Field.coord_x; + [CellInfo,NbDim,errormsg]=find_field_cells(Field); + YName='coord_y';%default + XName='coord_x';%default + for icell=1:numel(NbDim) + if NbDim(icell)==2 && strcmp(CellInfo{icell}.CoordType,'grid') + YName=CellInfo{icell}.YName; + XName=CellInfo{icell}.XName; + break + end + end + DataOut.ListVarName={YName, XName ,'UMean' , 'VMean','u2Mean','v2Mean','u2Mean_1','v2Mean_1','uvMean','CurlMean','DivMean','Curl2Mean','Div2Mean','Counter'}; + DataOut.VarDimName={YName,XName,{YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName},... + {YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName},{YName,XName}}; + DataOut.(YName)=Field.(YName); + DataOut.(XName)=Field.(XName); Uprev=Field.U;% store the current field for next iteration Vprev=Field.V; @@ -276,5 +287,5 @@ % VarName=DataOut.ListVarName{ivar};% name of the variable % DataOut.ListVarName=[DataOut.ListVarName {[VarName 'Profile']}];%append the name of the profile variable -% DataOut.VarDimName=[DataOut.VarDimName {'coord_y'}]; +% DataOut.VarDimName=[DataOut.VarDimName {'(YName)'}]; % DataOut.([VarName 'Profile'])=mean(DataOut.(VarName)(:,band),2); %take the mean profile of U, excluding the edges % end @@ -283,10 +294,10 @@ OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomTypeOut,first_i,last_i,first_j,last_j); %case of netcdf input file , determine global attributes - errormsg=struct2nc(OutputFile,DataOut); %save result file - if isempty(errormsg) +errormsg=struct2nc(OutputFile,DataOut); %save result file +if isempty(errormsg) disp([OutputFile ' written']); - else +else disp(['error in writting result file: ' errormsg]) - end +end Index: /trunk/src/uvmat.m =================================================================== --- /trunk/src/uvmat.m (revision 1081) +++ /trunk/src/uvmat.m (revision 1082) @@ -4687,5 +4687,5 @@ FieldList={AName}; case '1D plot' - YName=GetFieldData.PanelOrdinate.ordinate; + YName=GetFieldData.Coordinates.Coord_y; FieldList={''}; set(handles.uvmat,'ToolBar','figure')