Changeset 515

Timestamp:

Aug 15, 2012, 11:36:12 PM (12 years ago)

Author:

sommeria

Message:

improvement of calc-field and combination of two fields

Location:

trunk/src

Files:

• calc_field.m (modified) (6 diffs)
• calc_field_interp.m (added)
• calc_field_tps.m (added)
• calc_tps.m (modified) (2 diffs)
• check_files.m (modified) (2 diffs)
• fill_GUI.m (modified) (2 diffs)
• find_field_cells.m (modified) (9 diffs)
• griddata_uvmat.m (deleted)
• plot_field.m (modified) (2 diffs)
• proj_field.m (modified) (31 diffs)
• read_civdata.m (modified) (8 diffs)
• sub_field.m (modified) (5 diffs)
• uvmat.m (modified) (40 diffs)

trunk/src/calc_field.m

@@ -38 +38 @@
 
 %% list of field options implemented
-FieldOptions={'velocity';...%image correlation corresponding to a vel vector
-    'ima_cor';...%image correlation corresponding to a vel vector
-    'norm_vel';...%norm of the velocity
-    'vort';...%vorticity
-    'div';...%divergence
-    'strain';...%rate of strain
-    'u';... %u velocity component
-    'v';... %v velocity component
+FieldOptions={'vec(U,V)';...%image correlation corresponding to a vel vector
+    'C';...%image correlation corresponding to a vel vector
+    'norm(U,V)';...%norm of the velocity
+    'curl(U,V)';...%vorticity
+    'div(U,V)';...%divergence
+    'strain(U,V)';...%rate of strain
+    'U';... %u velocity component
+    'V';... %v velocity component
     'w';... %w velocity component
     'w_normal';... %w velocity component normal to the plane
@@ -77 +77 @@
 end
 FieldList=FieldList(check_calc==1);
-% if isempty(FieldList)
-%     DataOut=DataIn;
-%     return
-% end
 if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
     nbcoord=3;
@@ -91 +87 @@
 units_cell={};
 
+%% reproduce global attributes
+DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; 
+for ilist=1:numel(DataOut.ListGlobalAttribute)
+    DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
+end
+
 %% interpolation with new civ data
-if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var') || check_grid ||check_der)
-    %reproduce global attributes
-    DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; 
-    for ilist=1:numel(DataOut.ListGlobalAttribute)
-        DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
-    end
-
+[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(DataIn);
+icell_tps=[];
+for icell=1:numel(CellVarIndex)
+    VarType=VarTypeCell{icell};
+    if NbDimVec(icell)>=2 && ~isempty(VarType.subrange_tps)      
+        icell_tps=[icell_tps icell];
+    end
+end
+
+%if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var') || check_grid ||check_der)
+if ~isempty(icell_tps)
     %create a default grid if needed
     if  ~exist('Coord_interp','var')
-            XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates
-            YMax=max(max(DataIn.SubRange(2,:,:)));
-            XMin=min(min(DataIn.SubRange(1,:,:)));
-            YMin=min(min(DataIn.SubRange(2,:,:))); 
+        for ind=1:numel(icell_tps)
+            SubRange=DataIn.(DataIn.ListVarName{VarType{icell_tps(ind)}.subrange_tps});
+            XMax(ind)=max(max(SubRange(1,:,:)));% extrema of the coordinates
+            YMax(ind)=max(max(SubRange(2,:,:)));
+            XMin(ind)=min(min(SubRange(1,:,:)));
+            YMin(ind)=min(min(SubRange(2,:,:))); 
+        end
+        XMax=max(XMax);
+        YMax=max(YMax);
+        XMin=min(XMin);
+        YMin=min(YMin);
         if ~isfield(DataIn,'Mesh')
             DataIn.Mesh=sqrt(2*(XMax-XMin)*(YMax-YMin)/numel(DataIn.Coord_tps));
@@ -119 +132 @@
         coord_x=XMin:DataIn.Mesh:XMax;% increase the recommanded mesh to  visualisation
         coord_y=YMin:DataIn.Mesh:YMax;
-%         npx=length(coord_x);
-%         npy=length(coord_y);
         DataOut.coord_x=[XMin XMax];
         DataOut.coord_y=[YMin YMax];
         [XI,YI]=meshgrid(coord_x,coord_y);
-%         XI=reshape(XI,[],1);
-%         YI=reshape(YI,[],1);
         Coord_interp=cat(3,XI,YI);%[XI YI];
     end
@@ -131 +140 @@
     npy=size(Coord_interp,1);
     Coord_interp=reshape(Coord_interp,npx*npy,size(Coord_interp,3));
-%         npy=length(coord_y);
+
     %initialise output
     nb_sites=size(Coord_interp,1);
@@ -152 +161 @@
     
     % interpolate data in each subdomain
-    for isub=1:NbSubDomain
-        nbvec_sub=DataIn.NbSites(isub);
-        check_range=(Coord_interp >=ones(nb_sites,1)*DataIn.SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*DataIn.SubRange(:,2,isub)');
-        ind_sel=find(sum(check_range,2)==nb_coord);
-        %rho smoothing parameter
-        %                 epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites
-        %                 ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
-        nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
-        if check_grid
-            EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
-        end
-        if check_der
-            [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
-        end
-        ListVar={};
-        for ilist=1:length(FieldList)
-            var_count=numel(ListVar);
-            switch FieldList{ilist}
-                case 'velocity'
-                    ListVar=[ListVar {'U', 'V'}];
-                    VarAttribute{var_count+1}.Role='vector_x';
-                    VarAttribute{var_count+2}.Role='vector_y';
-                    DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
-                    DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
-                case 'u'
-                    ListVar=[ListVar {'u'}];
-                    VarAttribute{var_count+1}.Role='scalar';
-                    DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
-                case 'v'
-                    ListVar=[ListVar {'v'}];
-                    VarAttribute{var_count+1}.Role='scalar';
-                    DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
-                case 'norm_vel'
-                    ListVar=[ListVar {'norm_vel'}];
-                    VarAttribute{var_count+1}.Role='scalar';
-                    U=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
-                    V=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
-                    DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V);
-                case 'vort'
-                    ListVar=[ListVar {'vort'}];
-                    VarAttribute{var_count+1}.Role='scalar';
-                    DataOut.vort(ind_sel)=DataOut.vort(ind_sel)-EMDY *DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
-                case 'div'
-                    ListVar=[ListVar {'div'}];
-                    VarAttribute{var_count+1}.Role='scalar';
-                    DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDY *DataIn.V_tps(1:nbvec_sub+3,isub);
-                case 'strain'
-                    ListVar=[ListVar {'strain'}];
-                    VarAttribute{var_count+1}.Role='scalar';
-                    DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
+    for icell=icell_tps
+        for isub=1:NbSubDomain
+            nbvec_sub=DataIn.(DataIn.ListVarName{VarType{icell}.nbsites_tps});
+            SubRange=DataIn.(DataIn.ListVarName{VarType{icell}.subrange_tps});
+            check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
+            ind_sel=find(sum(check_range,2)==nb_coord);
+            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
+            Coord_tps=DataIn.(DataIn.ListVarName{VarType{icell}.coord_tps});
+            if check_grid
+                EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
             end
-        end
-    end
-    DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
-    nbval(nbval==0)=1;% to avoid division by zero for averaging
-    if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed
-        DataOut.ListVarName=[DataOut.ListVarName {'FF'}];
-        DataOut.VarDimName=[DataOut.VarDimName {{'coord_y','coord_x'}}];
-        DataOut.VarAttribute{length(DataOut.ListVarName)}.Role='errorflag';
-    end
-    DataOut.ListVarName=[DataOut.ListVarName ListVar];
-    for ifield=1:numel(ListVar)
-        VarDimName{ifield}={'coord_y','coord_x'};
-        DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval;
-        DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx);
-    end
-    DataOut.FF=reshape(DataOut.FF,npy,npx);
-    DataOut.VarDimName=[DataOut.VarDimName VarDimName];     
-    DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];
+            if check_der
+                [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
+            end
+            ListVar={};
+            U_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_x_tps});
+            V_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_y_tps});
+            for ilist=1:length(FieldList)
+                var_count=numel(ListVar);
+                switch FieldList{ilist}
+                    case 'velocity'
+                        ListVar=[ListVar {'U', 'V'}];
+                        VarAttribute{var_count+1}.Role='vector_x';
+                        VarAttribute{var_count+2}.Role='vector_y';
+                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
+                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
+                    case 'u'
+                        ListVar=[ListVar {'u'}];
+                        VarAttribute{var_count+1}.Role='scalar';
+                        DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
+                    case 'v'
+                        ListVar=[ListVar {'v'}];
+                        VarAttribute{var_count+1}.Role='scalar';
+                        DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
+                    case 'norm_vel'
+                        ListVar=[ListVar {'norm_vel'}];
+                        VarAttribute{var_count+1}.Role='scalar';
+                        U=DataOut.U(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
+                        V=DataOut.V(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
+                        DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V);
+                    case 'vort'
+                        ListVar=[ListVar {'vort'}];
+                        VarAttribute{var_count+1}.Role='scalar';
+                        DataOut.vort(ind_sel)=DataOut.vort(ind_sel)-EMDY *DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
+                    case 'div'
+                        ListVar=[ListVar {'div'}];
+                        VarAttribute{var_count+1}.Role='scalar';
+                        DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDY *DataIn.V_tps(1:nbvec_sub+3,isub);
+                    case 'strain'
+                        ListVar=[ListVar {'strain'}];
+                        VarAttribute{var_count+1}.Role='scalar';
+                        DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
+                end
+            end
+        end
+        DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
+        nbval(nbval==0)=1;% to avoid division by zero for averaging
+        if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed
+            DataOut.ListVarName=[DataOut.ListVarName {'FF'}];
+            DataOut.VarDimName=[DataOut.VarDimName {{'coord_y','coord_x'}}];
+            DataOut.VarAttribute{length(DataOut.ListVarName)}.Role='errorflag';
+        end
+        DataOut.ListVarName=[DataOut.ListVarName ListVar];
+        for ifield=1:numel(ListVar)
+            VarDimName{ifield}={'coord_y','coord_x'};
+            DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval;
+            DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx);
+        end
+        DataOut.FF=reshape(DataOut.FF,npy,npx);
+        DataOut.VarDimName=[DataOut.VarDimName VarDimName];
+        DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];
+    end
 else
 

trunk/src/calc_tps.m

@@ -1 @@
-function DataOut=calc_tps(DataIn)     
+function DataOut=calc_tps(DataIn,checkall)     
 DataOut=DataIn;%default
 SubDomain=1000; %default, estimated nbre of vectors in a subdomain used for tps
@@ -5 +5 @@
     SubDomain=DataIn.SubDomain;%
 end
-[DataOut.SubRange,DataOut.NbSites,DataOut.Coord_tps,DataOut.U_tps,DataOut.V_tps] =...
-    filter_tps([DataIn.X(DataIn.FF==0) DataIn.Y(DataIn.FF==0)],DataIn.U(DataIn.FF==0),DataIn.V(DataIn.FF==0),[],SubDomain,0);
-nbvar=numel(DataIn.ListVarName);
-DataOut.ListVarName=[DataIn.ListVarName {'SubRange','NbSites','Coord_tps','U_tps','V_tps'}];
-DataOut.VarDimName=[DataIn.VarDimName {{'nb_coord','nb_bounds','nb_subdomain'},{'nb_subdomain'},...
-    {'nb_tps','nb_coord','nb_subdomain'},{'nb_tps','nb_subdomain'},{'nb_tps','nb_subdomain'}}];
-DataOut.VarAttribute{nbvar+3}.Role='coord_tps';
-DataOut.VarAttribute{nbvar+4}.Role='vector_x';
-DataOut.VarAttribute{nbvar+5}.Role='vector_y';
-if isfield(DataOut,'ListDimName')%cleaning
-    DataOut=rmfield(DataOut,'ListDimName');
+[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(DataIn);
+nbtps=0;
+for icell=1:numel(CellVarIndex);
+    VarType=VarTypeCell{icell};
+    if NbDimVec(icell)>=2 && ~isempty(VarType.coord_x)
+        nbtps=nbtps+1;
+        X=DataIn.(DataIn.ListVarName{VarType.coord_x});
+        Y=DataIn.(DataIn.ListVarName{VarType.coord_y});
+        if ~isempty(VarType.vector_x)&&~isempty(VarType.vector_y)
+            Attr=DataIn.VarAttribute{VarType.vector_x};
+            if ~isfield(Attr,'VarIndex_tps')&& (checkall || (isfield(Attr,'FieldRequest')&&strcmp(Attr.FieldRequest,'derivatives')))               
+                U=DataIn.(DataIn.ListVarName{VarType.vector_x});
+                V=DataIn.(DataIn.ListVarName{VarType.vector_y});
+            else
+                continue
+            end
+        end
+        if ~isempty(VarType.errorflag)
+            FF=DataIn.(DataIn.ListVarName{VarType.errorflag});
+            X=X(FF==0);
+            Y=Y(FF==0);
+            U=U(FF==0);
+            V=V(FF==0);
+        end
+        if nbtps==1
+            ListNewVar={'SubRange','NbSites','Coord_tps','U_tps','V_tps'};
+            ListNewDim={'nb_tps','nb_subdomain'};
+            DataOut.VarDimName=[DataIn.VarDimName {{'nb_coord','nb_bounds','nb_subdomain'},{'nb_subdomain'},...
+                {'nb_tps','nb_coord','nb_subdomain'},{'nb_tps','nb_subdomain'},{'nb_tps','nb_subdomain'}}];
+        else
+            ListNewVar={['SubRange_' num2str(nbtps-1)],['NbSites_' num2str(nbtps-1)],['Coord_tps_' num2str(nbtps-1)],['U_tps_' num2str(nbtps-1)] ,['V_tps_' num2str(nbtps-1)]};
+            ListNewDim={['nb_tps_' num2str(nbtps-1)],['nb_subdomain_' num2str(nbtps-1)]};
+            DataOut.VarDimName=[DataIn.VarDimName {{'nb_coord','nb_bounds',ListNewDim{2}},ListNewDim(2),...
+                {ListNewDim{1},'nb_coord',ListNewDim{2}},ListNewDim,ListNewDim}];
+        end
+        DataOut.ListVarName=[DataIn.ListVarName ListNewVar];
+        
+        [DataOut.(ListNewVar{1}),DataOut.(ListNewVar{2}),DataOut.(ListNewVar{3}),DataOut.(ListNewVar{4}),DataOut.(ListNewVar{5})] =...
+            filter_tps([X Y],U,V,[],SubDomain,0);
+        nbvar=numel(DataIn.ListVarName);
+        
+        DataOut.VarAttribute{nbvar+3}.Role='coord_tps';
+        DataOut.VarAttribute{nbvar+4}=DataIn.VarAttribute{VarType.vector_x};%reproduce attributes of velocity
+         DataOut.VarAttribute{nbvar+4}.Role='vector_x_tps';
+         DataIn.VarAttribute{VarType.vector_x}.VarIndex_tps=nbvar+4;% indicte the correspondance with initial data
+        DataOut.VarAttribute{nbvar+5}=DataIn.VarAttribute{VarType.vector_y};%reproduce attributes of velocity 
+         DataOut.VarAttribute{nbvar+5}.Role='vector_y_tps';
+%          if isfield(DataOut.VarAttribute{VarType.vector_x},'FieldRequest')
+%              DataOut.VarAttribute{VarType.vector_x}=rmfield(DataOut.VarAttribute{VarType.vector_x},'FieldRequest');
+%          end
+%          if isfield(DataOut.VarAttribute{VarType.vector_x},'Operation')
+%              DataOut.VarAttribute{VarType.vector_x}=rmfield(DataOut.VarAttribute{VarType.vector_x},'Operation');
+%          end
+        if isfield(DataOut,'ListDimName')%cleaning'FieldRequest'
+            DataOut=rmfield(DataOut,'ListDimName');
+        end
+        if isfield(DataOut,'DimValue')%cleaning
+            DataOut=rmfield(DataOut,'DimValue');
+        end
+    end
 end
-if isfield(DataOut,'DimValue')%cleaning
-    DataOut=rmfield(DataOut,'DimValue');
-end

trunk/src/check_files.m

@@ -32 +32 @@
 svn_info.status=[];
 list_fct={...
-    'calc_field';...% defines fields (velocity, vort, div...) from civx data and calculate them
+    'calc_field_interp';...% defines fields (velocity, vort, div...) from civx data and calculate them
+    'calc_field_tps';...% defines fields (velocity, vort, div...) and calculate them
     'cell2tab';... %transform a Matlab cell in a character array suitable for display in a table
     'check_field_structure';...% check the validity of the field struture representation consistant with the netcdf format
@@ -64 +65 @@
     'get_file_series';...% determine the list of file names and file indices for functions called by 'series'.
     'get_file_type';...% determine info about a file (image, multimage, civdata,...) .
-    'griddata_uvmat';...%make 2D linear interpolation using griddata, with input appropriate for both Matlab 6.5 and 7
     'hist_update';...%  update of a current global histogram by inclusion of a new field
     'imadoc2struct';...%convert the image documentation file ImaDoc into a Matlab structure

trunk/src/fill_GUI.m

@@ -31 +31 @@
         hh=[];
         input_data=Param.(fields{ifield});
-                    display(fields{ifield})
-                    display(input_data)
+%                     display(fields{ifield})
+%                     display(input_data)
         check_done=0;
         if isfield(handles,fields{ifield})
@@ -68 +68 @@
                     end
                 case 'edit'
+                    input_string='';
                     if isnumeric(input_data)
+                        if numel(input_data)>0
                         input_string=num2str(input_data(ibox));
+                        end
                     else
                         input_string=input_data;

trunk/src/find_field_cells.m

@@ -19 +19 @@
 %      .scalar: scalar field (default)
 %      .coord: vector of indices of coordinate variables corresponding to matrix dimensions
+%
+%      .FieldRequest= 'interp', 'derivatives' indicate whether interpolation  or derivatives (tps) is needed to calculate the requested field
+%      .FieldNames = cell of fields to calculate from the fied cell
+%
 % errormsg: error message
 %   
@@ -49 +53 @@
 %     GNU General Public License (file UVMAT/COPYING.txt) for more details.for input in proj_field and plot_field
 %    group the variables  into 'fields' with common dimensions
-%------------------------------------------------------------------------
-% function  [CellVarIndex,NbDim,CellVarType,errormsg]=find_field_cells(Data)
-%
-% OUTPUT:
-% CellVaxIndex: cell whose elements are arrays of indices in the list data.ListVarName  
-%              CellvarIndex{i} represents a set of variables with the same dimensions
-% NbDim: array with the length of CellVarIndex, giving its  space dimension
-% CellVarType: cell array of structures with fields
-%      .coord_x, y, z: indices (in .ListVarname) of variables representing  unstructured coordinates x, y, z 
-%      .vector_x,_y,_z: indices of variables giving the vector components x, y, z
-%      .warnflag: index of warnflag
-%      .errorflag: index of error flag
-%      .ancillary: indices of ancillary variables
-%      .image   : B/W image, (behaves like scalar)
-%      .color : color image, the last index, which is not a coordinate variable, represent the 3 color components rgb
-%      .discrete: like scalar, but set of data points without continuity, represented as dots in a usual plot, instead of continuous lines otherwise
-%      .scalar: scalar field (default)
-%      .coord: vector of indices of coordinate variables corresponding to matrix dimensions
-% errormsg: error message
-%   
-% INPUT:
-% Data: structure representing fields, output of check_field_structure
-%            .ListVarName: cell array listing the names (cahr strings) of the variables
-%            .VarDimName: cell array of cells containing the set of dimension names for each variable of .ListVarName
-%            .VarAttribute: cell array of structures containing the variable attributes: 
-%                     .VarAttribute{ilist}.key=value, where ilist is the variable
-%                     index, key is the name of the attribute, value its value (char string or number)
-%
-% HELP: 
-% to get the dimensions of arrays common to the field #icell
-%         VarIndex=CellVarIndex{icell}; % list of variable indices
-%         DimIndex=Data.VarDimIndex{VarIndex(1)} % list of dimensions for each variable in the cell #icell
-%
-%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
-%  Copyright Joel Sommeria, 2008, LEGI / CNRS-UJF-INPG, sommeria@coriolis-legi.org.
-%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
-%     This file is part of the toolbox UVMAT.
-% 
-%     UVMAT is free software; you can redistribute it and/or modify
-%     it under the terms of the GNU General Public License as published by
-%     the Free Software Foundation; either version 2 of the License, or
-%     (at your option) any later version.
-% 
-%     UVMAT is distributed in the hope that it will be useful,
-%     but WITHOUT ANY WARRANTY; without even the implied warranty of
-%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-%     GNU General Public License (file UVMAT/COPYING.txt) for more details.
-%AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
 
 function [CellVarIndex,NbDim,CellVarType,errormsg]=find_field_cells(Data)
 CellVarIndex={};
 
-NbDim=[];
 CellVarType=[];
 errormsg=[];
@@ -108 +63 @@
 
 NbDim=[];
-ivardim=0;
 VarDimIndex=[];
 VarDimName={};
@@ -116 +70 @@
 end
 
-%% role of variables
+%% role of variables and list of requested operations
 Role=num2cell(blanks(nbvar));%initialize a cell array of nbvar blanks
+FieldRequest=regexprep(Role,' ',''); % fieldRequest set to '' by default
+Operation=cell(size(Role)); % fieldRequest set to {} by default
+CheckSub=zeros(size(Role));% =1 for fields to substract
 Role=regexprep(Role,' ','scalar'); % Role set to 'scalar' by default
 if isfield(Data,'VarAttribute')
@@ -124 +81 @@
             Role{ivar}=Data.VarAttribute{ivar}.Role;
         end
+        if isfield(Data.VarAttribute{ivar},'FieldRequest')
+            FieldRequest{ivar}=Data.VarAttribute{ivar}.FieldRequest;
+        end
+        if isfield(Data.VarAttribute{ivar},'Operation')
+            Operation{ivar}=Data.VarAttribute{ivar}.Operation;
+        end
+        if isfield(Data.VarAttribute{ivar},'CheckSub')
+            CheckSub(ivar)=Data.VarAttribute{ivar}.CheckSub;
+        end
     end
 end
 
 %% loop on the list of variables, group them by common dimensions
+CellVarType=cell(1,length(CellVarIndex));
 for ivar=1:nbvar
     if ischar(Data.VarDimName{ivar})
@@ -146 +113 @@
         icell=icell+1;
         CellVarIndex{icell}=ivar;%put the current variable index in the new cell 
-        NbDim(icell)=numel(DimCell);%default
-    end
-   
-    %look for dimension variables
-%     if numel(DimCell)==1% if the variable has a single dimension 
-%         if strcmp(DimCell{1},Data.ListVarName{ivar}) %|| strcmp(Role{ivar},'dimvar')
-%             ivardim=ivardim+1;
-%             VarDimIndex(ivardim)=ivar;%index of the variable
-%             VarDimName{ivardim}=DimCell{1};%name of the dimension
-%         end
-%     end
+        NbDim(icell)=numel(DimCell);%default   
+        CellVarType{icell}=[];
+    end
+    if ~isempty(FieldRequest{ivar})
+       CellVarType{icell}.FieldRequest=FieldRequest{ivar};
+    end
+    if ~isempty(Operation{ivar})
+       CellVarType{icell}.Operation=Operation{ivar};
+    end
+    if CheckSub(ivar)
+    CellVarType{icell}.CheckSub=1;
+    end
 end
 
@@ -163 +131 @@
 ind_dim_var_cell=find(checksinglecell);
 %CoordType(ind_dim_var_cell)='dim_var';% to be used in output
+%VarDimIndex=cell(size(ind_dim_var_cell));
+VarDimName=cell(size(ind_dim_var_cell));
 for icoord=1:numel(ind_dim_var_cell)
-VarDimIndex(icoord)=CellVarIndex{ind_dim_var_cell(icoord)};
-VarDimName{icoord}=Data.VarDimName{VarDimIndex(icoord)}{1};
+    VarDimIndex(icoord)=CellVarIndex{ind_dim_var_cell(icoord)};
+    VarDimName{icoord}=Data.VarDimName{VarDimIndex(icoord)}{1};
 end
 
@@ -201 +171 @@
 
 index_remove=[];
-CellVarType=cell(1,length(CellVarIndex));
 for icell=1:length(CellVarIndex)
     if checksinglecell(icell)
@@ -262 +231 @@
                     coord(idim)=VarDimIndex(ind_coord);
                 end
-%                 for ivardim=1:numel(VarDimName)
-%                     if strcmp(VarDimName{ivardim},DimCell{idim})
-%                         coord(idim)=VarDimIndex(ivardim);
-%                         break
-%                     end
-%                 end
             end
             NbDim(icell)=numel(find(coord));

trunk/src/plot_field.m

@@ -122 +122 @@
     return
 end
-index_2D=find(NbDim==2,2);%find 2D fields (at most 2)
+index_2D=find(NbDim==2);%find 2D fields
 index_3D=find(NbDim>2,1);
 if ~isempty(index_3D)
@@ -560 +560 @@
     if ~isempty(ivar_U) && ~isempty(ivar_V)% vector components detected
         if test_vec
-            errormsg='error in plot_field: attempt to plot two vector fields';
+            errormsg='error in plot_field: attempt to plot two vector fields: to get the difference project on a plane with mode interp';
             return
         else

trunk/src/proj_field.m

@@ -84 +84 @@
 ProjData=[];
 
-%% case of no projection (object is used only as graph display)
-if isfield(ObjectData,'ProjMode') && (isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside'))
+%% in the absence of object Type or projection mode, or object coordinaes, the output is empty
+if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode')
     return
 end
-
-%% check coincidence of coordinate units
-if isfield(FieldData,'CoordUnit') && isfield(ObjectData,'CoordUnit')&&~strcmp(FieldData.CoordUnit,ObjectData.CoordUnit)
-    errormsg='inconsistent coord units for field and projection object';
+% case of no projection (object is used only as graph display)
+if isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside')
     return
 end
-    
-%% in the absence of object Type or projection mode, or object coordinaes, the input field is just tranfered without change
-if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode')
-    ProjData=FieldData;
-    return
-end
-if ~isfield(ObjectData,'Coord')
+if ~isfield(ObjectData,'Coord')||isempty(ObjectData.Coord)
     if strcmp(ObjectData.Type,'plane')
         ObjectData.Coord=[0 0 0];%default
     else
-        ProjData=FieldData;
         return
     end
 end
-        
-%% OBSOLETE
-if isfield(ObjectData,'XMax') && ~isempty(ObjectData.XMax)
-    ObjectData.RangeX(1)=ObjectData.XMax;
-end
-if isfield(ObjectData,'XMin') && ~isempty(ObjectData.XMin)
-    ObjectData.RangeX(2)=ObjectData.XMin;
-end
-if isfield(ObjectData,'YMax') && ~isempty(ObjectData.YMax)
-    ObjectData.RangeY(1)=ObjectData.YMax;
-end
-if isfield(ObjectData,'YMin') && ~isempty(ObjectData.YMin)
-    ObjectData.RangeY(2)=ObjectData.YMin;
-end
-if isfield(ObjectData,'ZMax') && ~isempty(ObjectData.ZMax)
-    ObjectData.RangeZ(1)=ObjectData.ZMax;
-end
-if isfield(ObjectData,'ZMin') && ~isempty(ObjectData.ZMin)
-    ObjectData.RangeZ(2)=ObjectData.ZMin;
-end
-%%%%%%%%%%
 
 %% apply projection depending on the object type
@@ -401 +371 @@
     testproj(VarType.color)=1;
     VarIndex=VarIndex(find(testproj(VarIndex)));%select only the projected variables
-    if testX %case of unstructured coordinates
+    if check_unstructured%case of unstructured coordinates
          eval(['nbpoint=numel(FieldData.' FieldData.ListVarName{VarIndex(1)} ');'])
          for ivar=[VarIndex VarType.coord_x VarType.coord_y VarType.errorflag]
@@ -547 +517 @@
 ProjData.NbDim=1;
 %initialisation of the input parameters and defaultoutput
-ProjMode='projection';%direct projection on the line by default
-if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end; 
+ProjMode=ObjectData.ProjMode; 
 % ProjAngle=90; %90 degrees projection by default
 
@@ -609 +578 @@
         continue
     end
-    testX=~isempty(VarType.coord_x) && ~isempty(VarType.coord_y);% test for unstructured coordinates
+    check_unstructured=~isempty(VarType.coord_x) && ~isempty(VarType.coord_y);% test for unstructured coordinates
     test_tps=~isempty(VarType.coord_tps);
     testU=~isempty(VarType.vector_x) && ~isempty(VarType.vector_y);% test for vectors
@@ -631 +600 @@
     end   
     % check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y
-    if testX
+    if check_unstructured
         if  ~isequal(ProjMode,'interp')
             if width==0
@@ -664 +633 @@
   %%%%%%%  % A FAIRE CALCULER MEAN DES QUANTITES    %%%%%%
    %case of unstructured coordinates
-    if testX   
+    if check_unstructured   
         for ip=1:siz_line(1)-1     %Loop on the segments of the polyline
             linelength=sqrt(dlinx(ip)*dlinx(ip)+dliny(ip)*dliny(ip));  
@@ -917 +886 @@
 %-----------------------------------------------------------------
 
-%% initialisation of the input parameters of the projection plane
-ProjMode='projection';%direct projection by default
-if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
-
-%% axis origin
-if isempty(ObjectData.Coord)
-    ObjectData.Coord(1,1)=0;%origin of the plane set to [0 0] by default
-    ObjectData.Coord(1,2)=0;
-    ObjectData.Coord(1,3)=0;
-end
-
 %% rotation angles 
 PlaneAngle=[0 0 0]; 
@@ -951 +909 @@
 
 %% mesh sizes DX and DY
-DX=0;
-DY=0; %default 
+DX=FieldData.Mesh;
+DY=FieldData.Mesh; %default 
 if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX)
      DX=abs(ObjectData.DX);%mesh of interpolation points 
@@ -959 +917 @@
      DY=abs(ObjectData.DY);%mesh of interpolation points 
 end
-if  ~strcmp(ProjMode,'projection') && (DX==0||DY==0)
+if  ~strcmp(ObjectData.ProjMode,'projection') && (DX==0||DY==0)
         errormsg='DX or DY missing';
         display(errormsg)
@@ -970 +928 @@
 testYMin=0;
 testYMax=0;
+XMin=FieldData.XMin;%default
+XMax=FieldData.XMax;%default
+YMin=FieldData.YMin;%default
+YMax=FieldData.YMax;%default
 if isfield(ObjectData,'RangeX')
         XMin=min(ObjectData.RangeX);
@@ -1003 +965 @@
 ListIndex={};
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
-%-----------------------------------------------------------------
-idimvar=0;
-
+% CellVarIndex=cells of variable index arrays
 [CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(FieldData);
 if ~isempty(errormsg)
@@ -1014 +973 @@
 end
 
+%% projection modes
+check_grid=0;
+ProjMode=cell(size(VarTypeCell));
+for icell=1:numel(VarTypeCell)% TODO: recalculate coordinates here to get the bounds in the rotated coordinates
+    ProjMode{icell}=ObjectData.ProjMode;
+    if isfield(VarTypeCell{icell},'FieldRequest')
+        switch VarTypeCell{icell}.FieldRequest
+            case 'interp'
+                ProjMode{icell}='interp';
+            case 'derivatives'
+                ProjMode{icell}='filter';
+        end
+    end
+    if strcmp(ProjMode{icell},'interp')||strcmp(ProjMode{icell},'filter')
+        check_grid=1;
+    end
+end
+
+%% define the new coordinates in case of interpolation on a grid
+if check_grid% TODO: recalculate coordinates to get the bounds in the rotated coordinates
+    ProjData.ListVarName={'coord_y','coord_x'};
+    ProjData.VarDimName={'coord_y','coord_x'};  
+    ProjData.VarAttribute={[],[]};
+    ProjData.coord_y=[YMin YMax];
+    ProjData.coord_x=[XMin XMax];
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
 % CellVarIndex=cells of variable index arrays
@@ -1020 +1007 @@
 nbcoord=0;%number of added coordinate variables brought by projection
 nbvar=0;
+vector_x_proj=[];
+vector_y_proj=[];
 for icell=1:length(CellVarIndex)
     NbDim=NbDimVec(icell);
@@ -1027 +1016 @@
     VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list FieldData.ListVarName
     VarType=VarTypeCell{icell};
-%     ivar_X=VarType.coord_x;
-%     ivar_Y=VarType.coord_y;
-%     ivar_Z=VarType.coord_z;
     ivar_U=VarType.vector_x;
     ivar_V=VarType.vector_y;
@@ -1037 +1023 @@
     end
     ivar_W=VarType.vector_z;
-%     ivar_Anc=VarType.ancillary;
-%     test_anc=zeros(size(VarIndex));
-%     test_anc(ivar_Anc)=ones(size(ivar_Anc));
-%     ivar_F=VarType.warnflag;
-%     ivar_FF=VarType.errorflag;
 
     %type of coordinates
@@ -1058 +1039 @@
     end
     coord_z=0;%default
-    
+
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     switch CoordType
@@ -1064 +1045 @@
         %% case of input fields with unstructured coordinates
         case 'unstructured'
-            if strcmp(ObjectData.ProjMode,'filter')
-                continue %skip for filter (needs tps)
-            end
-            coord_x=FieldData.(FieldData.ListVarName{VarType.coord_x});
-            coord_y=FieldData.(FieldData.ListVarName{VarType.coord_y});
+            if strcmp(ProjMode{icell},'filter')
+                continue %skip for filter (needs tps field cell)
+            end
+            coord_x=FieldData.(FieldData.ListVarName{VarType.coord_x});% initial x coordinates 
+            coord_y=FieldData.(FieldData.ListVarName{VarType.coord_y});% initial y coordinates
             if ~isempty(VarType.coord_z)
-                coord_z=FieldData.(FieldData.ListVarName{VarType.coord_z});
+                coord_z=FieldData.(FieldData.ListVarName{VarType.coord_z});% initial z coordinates
             end
             
-            % translate  initial coordinates
+            % translate  initial coordinates to account for the new origin
             coord_x=coord_x-ObjectData.Coord(1,1);
             coord_y=coord_y-ObjectData.Coord(1,2);
@@ -1085 +1066 @@
                 fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane
                 indcut=find(abs(fieldZ) <= width);
-                size(indcut)
                 for ivar=VarIndex
                     VarName=FieldData.ListVarName{ivar};
@@ -1096 +1076 @@
             end
             
-            %rotate coordinates if needed:
+            %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane
             Psi=PlaneAngle(1);
             Theta=PlaneAngle(2);
@@ -1104 +1084 @@
                 coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
                 coord_Y=coord_Y+coord_z *sin(Theta);
-                coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
-                
+                coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER               
                 coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
             else
@@ -1112 +1091 @@
             end
             
-            %restriction to the range of x and y if imposed
+            %restriction to the range of X and Y if imposed by the projection object
             testin=ones(size(coord_X)); %default
             testbound=0;
@@ -1139 +1118 @@
                 for ivar=VarIndex
                     VarName=FieldData.ListVarName{ivar};
-                    eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
+                    FieldData.(VarName)=FieldData.(VarName)(indcut);
                 end
                 coord_X=coord_X(indcut);
@@ -1149 +1128 @@
             
             % different cases of projection
-            switch ObjectData.ProjMode
-                case 'projection'
+            switch ProjMode{icell}
+                case 'projection'  
+                    nbvar=numel(ProjData.ListVarName);
                     for ivar=VarIndex %transfer variables to the projection plane
                         VarName=FieldData.ListVarName{ivar};
@@ -1172 +1152 @@
                     coord_x_proj=XMin:DX:XMax;
                     coord_y_proj=YMin:DY:YMax;
-                    DimCell={'coord_y','coord_x'};
-                    ProjData.ListVarName={'coord_y','coord_x'};
-                    ProjData.VarDimName={'coord_y','coord_x'};
-                    nbcoord=2;
-                    ProjData.coord_y=[YMin YMax];
-                    ProjData.coord_x=[XMin XMax];
-%                     if isempty(ivar_X), ivar_X=0; end;
-%                     if isempty(ivar_Y), ivar_Y=0; end;
-%                     if isempty(ivar_Z), ivar_Z=0; end;
-%                     if isempty(ivar_U), ivar_U=0; end;
-%                     if isempty(ivar_V), ivar_V=0; end;
-%                     if isempty(ivar_W), ivar_W=0; end;
-%                     if isempty(ivar_F), ivar_F=0; end;
-%                     if isempty(ivar_FF), ivar_FF=0; end;
-                  %  if ~isequal(ivar_FF,0)
-                   if ~isempty(VarType.errorflag)
+                    [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);
+                    if ~isempty(VarType.errorflag)
                         VarName_FF=FieldData.ListVarName{VarType.errorflag};
                         indsel=find(FieldData.(VarName_FF)==0);
@@ -1193 +1159 @@
                         coord_Y=coord_Y(indsel);
                     end
-                    
-                    FF=zeros(1,length(coord_y_proj)*length(coord_x_proj));
                     testFF=0;
-                    %for ivar=VarIndex % loop on field variables to project
-                    for ivar=[VarType.scalar VarType.vector_x VarType.vector_y]
-                        VarName=FieldData.ListVarName{ivar};
-%                         if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF ||...
-%                                 ( numel(test_anc)>=ivar && test_anc(ivar)==1))
-                            ivar_new=ivar_new+1;
-                            ProjData.ListVarName=[ProjData.ListVarName {VarName}];
-                            ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
-                            if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
-                                ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar};
-                            end
-                            %if  ~isequal(ivar_FF,0)
-                            if ~isempty(VarType.errorflag)
-                                FieldData.(VarName)=FieldData.(VarName)(indsel);
-                            end
-                            ProjData.(VarName)=griddata_uvmat(double(coord_X),double(coord_Y),double(FieldData.(VarName)),coord_x_proj,coord_y_proj');
-                            varline=reshape(ProjData.(VarName),1,length(coord_y_proj)*length(coord_x_proj));
-                            FFlag= isnan(varline); %detect undefined values NaN
-                            indnan=find(FFlag);
-                            if~isempty(indnan)
-                                varline(indnan)=zeros(size(indnan));
-                                ProjData.(VarName)=reshape(varline,length(coord_y_proj),length(coord_x_proj));
-                                FF(indnan)=ones(size(indnan));
-                                testFF=1;
-                            end
-                            if ivar==ivar_U
-                                ivar_U=ivar_new;
-                            end
-                            if ivar==ivar_V
-                                ivar_V=ivar_new;
-                            end
-                            if ivar==ivar_W
-                                ivar_W=ivar_new;
-                            end
-%                         end
-                    end
-                    if testFF
-                        ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
-                        ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
-                        ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
-                        ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
-                    end
-%                     case 'filter'%interpolate data on a regular grid
-%                         errormsg='tps required for filter option'
-                        
-            end
-            
-            %% case of tps interpolation (applies only in filter mode)
+                    nbvar=numel(ProjData.ListVarName);
+                    if isfield(VarType,'vector_x')&&isfield(VarType,'vector_y')&&~isempty(VarType.vector_x)
+                        VarName_x=FieldData.ListVarName{VarType.vector_x};
+                        VarName_y=FieldData.ListVarName{VarType.vector_y};
+                        if ~isempty(VarType.errorflag)
+                            FieldData.(VarName_x)=FieldData.(VarName_x)(indsel);
+                            FieldData.(VarName_y)=FieldData.(VarName_y)(indsel);
+                        end
+                        FieldVar=cat(2,FieldData.(VarName_x),FieldData.(VarName_y));
+                        if ~isfield(VarType,'CheckSub') || ~VarType.CheckSub
+                            vector_x_proj=numel(ProjData.ListVarName)+1;
+                            vector_y_proj=numel(ProjData.ListVarName)+2;
+                        end
+                    end
+                    if ~isempty(VarType.scalar)
+                        VarName_scalar=FieldData.ListVarName{VarType.scalar};
+                        if ~isempty(VarType.errorflag)
+                            FieldData.(VarName_scalar)=FieldData.(VarName_scalar)(indsel);
+                        end
+                        FieldVar=FieldData.(VarName_scalar);
+                    end
+                    [VarVal,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldVar,VarType.Operation,XI,YI);
+                    if isfield(VarType,'CheckSub') && VarType.CheckSub && ~isempty(vector_x_proj)
+                        ProjData.(ProjData.ListVarName{vector_x_proj})=ProjData.(ProjData.ListVarName{vector_x_proj})-VarVal{1};
+                        ProjData.(ProjData.ListVarName{vector_y_proj})=ProjData.(ProjData.ListVarName{vector_y_proj})-VarVal{2};
+                    else
+                        VarDimName=cell(size(ListFieldProj));
+                        for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
+                            ListFieldProj{ilist}=regexprep(ListFieldProj{ilist},'(.+','');
+                            if ~isempty(find(strcmp(ListFieldProj{ilist},ProjData.ListVarName)))
+                                ListFieldProj{ilist}=[ListFieldProj{ilist} '_1'];
+                            end                        
+                            ProjData.(ListFieldProj{ilist})=VarVal{ilist};
+                            VarDimName{ilist}={'coord_y','coord_x'};
+                        end
+                        ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
+                        ProjData.VarDimName=[ProjData.VarDimName VarDimName];
+                        ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
+                    end
+            end
+
+            %% case of tps interpolation (applies only in filter mode and for spatial derivatives)
         case 'tps'
-            if strcmp(ObjectData.ProjMode,'filter')
+            if strcmp(ProjMode{icell},'filter')
+                Coord=FieldData.(FieldData.ListVarName{VarType.coord_tps});
+                NbSites=FieldData.(FieldData.ListVarName{VarType.nbsites_tps});
+                SubRange=FieldData.(FieldData.ListVarName{VarType.subrange_tps});
+                if isfield(VarType,'vector_x_tps')&&isfield(VarType,'vector_y_tps')
+                    FieldVar=cat(3,FieldData.(FieldData.ListVarName{VarType.vector_x_tps}),FieldData.(FieldData.ListVarName{VarType.vector_y_tps}));
+                end
                 coord_x_proj=XMin:DX:XMax;
                 coord_y_proj=YMin:DY:YMax;
@@ -1253 +1217 @@
                 XI=XI+ObjectData.Coord(1,1);
                 YI=YI+ObjectData.Coord(1,2);
-                DataOut=calc_field(FieldData.FieldList,FieldData,cat(3,XI,YI));
-                ProjData.ListVarName=[ProjData.ListVarName DataOut.ListVarName];
-                ProjData.VarDimName=[ProjData.VarDimName DataOut.VarDimName];
-                ProjData.VarAttribute=[ProjData.VarAttribute DataOut.VarAttribute];   
-                for ilist=1:length(DataOut.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
-                    VarName=DataOut.ListVarName{ilist};
-%                     ProjData.(VarName)=DataOut.(VarName);
-                    if ilist>=3
-                    ProjData.(VarName)=reshape(DataOut.(VarName),np_y,np_x);
-                    end
-                end
-                ProjData.coord_x=[XMin XMax];
-                ProjData.coord_y=[YMin YMax];
+                [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbSites,SubRange,FieldVar,VarType.Operation,cat(3,XI,YI));   
+                ListFieldProj=(fieldnames(DataOut))';
+                VarDimName=cell(size(ListFieldProj));
+                for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
+                    VarName=ListFieldProj{ilist};
+                    ProjData.(VarName)=DataOut.(VarName);
+                    VarDimName{ilist}={'coord_y','coord_x'};
+                end
+%                 if ~isfield(ProjData,'coord_x')
+%                 ProjData.coord_x=[XMin XMax];
+%                 ProjData.coord_y=[YMin YMax];
+%                 ListFieldProj=[{'coord_x','coord_y'} ListFieldProj'];
+%                 VarDimName=[{'coord_x','coord_y'} VarDimName];
+%                 VarAttribute=[{[],[]} VarAttribute];
+%                 end
+                ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
+                ProjData.VarDimName=[ProjData.VarDimName VarDimName];
+                ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
             end
             
@@ -1272 +1241 @@
             
             VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
-            eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions
+            DimValue=size(FieldData.(VarName));%input matrix dimensions
             DimValue(DimValue==1)=[];%remove singleton dimensions
             NbDim=numel(DimValue);%update number of space dimensions
@@ -1405 +1374 @@
             end
             % case with no  interpolation
-            if isequal(ProjMode,'projection') && (~testangle || test90y || test90x)
+            if isequal(ProjMode{icell},'projection') && (~testangle || test90y || test90x)
                 if  NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax
-                    ProjData=FieldData;% no change by projection
+                    ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName(VarIndex)];
+                    ProjData.VarDimName=[ProjData.VarDimName FieldData.VarDimName(VarIndex)];  
+                    ProjData.VarAttribute=[ProjData.VarAttribute FieldData.VarAttribute(VarIndex)];
+                    ProjData.(AYProjName)=FieldData.(AYName);
+                    ProjData.(AXProjName)=FieldData.(AXName);
+                    for ivar=VarIndex
+                        VarName=FieldData.ListVarName{ivar};
+                        ProjData.(VarName)=FieldData.(VarName);% no change by projection
+                    end
                 else
                     indY=NbDim-1;
@@ -1489 +1466 @@
                     YIMA=reshape(round(YIMA),1,npX*npY);
                     flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image
-                    if isequal(ObjectData.ProjMode,'filter')
+                    if isequal(ProjMode{icell},'filter')
                         npx_filter=ceil(abs(DX/DAX));
                         npy_filter=ceil(abs(DY/DAY));
@@ -1581 +1558 @@
 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   
+
 %-----------------------------------------------------------------
 %projection in a volume 
@@ -1586 +1586 @@
 %-----------------------------------------------------------------
 ProjData=FieldData;%default output
-
-%% initialisation of the input parameters of the projection plane
-ProjMode='projection';%direct projection by default
-if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
 
 %% axis origin

trunk/src/read_civdata.m

@@ -10 +10 @@
 %                    .NbCoord: number of vector components
 %                    .NbDim: number of dimensions (=2 or 3)
-%                    .dt: time interval for the corresponding image pair
+%                    .Dt: time interval for the corresponding image pair
 %                    .Time: absolute time (average of the initial image pair)
 %                    .CivStage: =0, 
@@ -63 +63 @@
 end
 errormsg='';
+if ischar(FieldNames), FieldNames={FieldNames}; end;
+FieldRequest='vec';
+for ilist=1:length(FieldNames)
+    if ~isempty(FieldNames{ilist})
+        switch FieldNames{ilist}
+            case{'U','V','norm(U,V)'}
+                FieldRequest='interp';
+            case {'curl(U,V)','div(U,V)','strain(U,V)'}
+                FieldRequest='derivatives';
+        end
+    end
+end
 
 %% reading data
-[varlist,role,units,VelTypeOut]=varcivx_generator(FieldNames,VelType,CivStage);
+[varlist,role,VelTypeOut]=varcivx_generator(FieldRequest,VelType,CivStage);
 if isempty(varlist)
     erromsg=['error in read_civdata: unknow velocity type ' VelType];
@@ -97 +109 @@
 end
 Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'Dt','Time'}];
+ivar_U_tps=[];
 var_ind=find(vardetect);
 for ivar=1:numel(var_ind)
     Field.VarAttribute{ivar}.Role=role{var_ind(ivar)};
-    Field.VarAttribute{ivar}.Unit=units{var_ind(ivar)};
+    Field.VarAttribute{ivar}.FieldRequest=FieldRequest;
+    if strcmp(role{var_ind(ivar)},'vector_x')
+        Field.VarAttribute{ivar}.Operation=FieldNames;
+        ivar_U=ivar;
+    end
+    if strcmp(role{var_ind(ivar)},'vector_x_tps')
+        Field.VarAttribute{ivar}.Operation=FieldNames;
+        ivar_U_tps=ivar;
+    end
+%     Field.VarAttribute{ivar}.Unit=units{var_ind(ivar)};
     Field.VarAttribute{ivar}.Mesh=0.1;%typical mesh for histograms O.1 pixel
 end
+if ~isempty(ivar_U_tps)
+    Field.VarAttribute{ivar_U}.VarIndex_tps=ivar_U_tps;
+end
+% if strcmp(FieldRequest,'derivatives')% fields will be calculated from the tps
+%     Field.VarAttribute{ivar_U_tps}.FieldRequest=FieldRequest;
+%     Field.VarAttribute{ivar_U_tps}.Operation=FieldNames;
+% else% fields will be calculated from the initial fields
+%     Field.VarAttribute{ivar_U}.FieldRequest=FieldRequest;%
+%     Field.VarAttribute{ivar_U}.Operation=FieldNames;
+% end
+
 Field.ListGlobalAttribute=[Field.ListGlobalAttribute {'NbCoord','NbDim','TimeUnit','CoordUnit'}];
 % %% update list of global attributes
@@ -121 +154 @@
 %            if vel_type=[] or'*', a  priority choice is done, civ2 considered better than civ1 )
 
-function [var,role,units,vel_type_out,errormsg]=varcivx_generator(FieldNames,vel_type,CivStage) 
+function [var,role,vel_type_out,errormsg]=varcivx_generator(FieldRequest,vel_type,CivStage) 
 
 %% default input values
 if ~exist('vel_type','var'),vel_type='';end;
 if iscell(vel_type),vel_type=vel_type{1}; end;%transform cell to string if needed
-if ~exist('FieldNames','var'),FieldNames={'ima_cor'};end;%default scalar 
-if ischar(FieldNames), FieldNames={FieldNames}; end;
 errormsg='';
 
 %% select the priority order for automatic vel_type selection
-testder=0;
-testpatch=0;
-for ilist=1:length(FieldNames)
-    if ~isempty(FieldNames{ilist})
-        switch FieldNames{ilist}
-            case{'u','v'}
-                testpatch=1;
-            case {'vort','div','strain'}
-                testder=1;
-        end
-    end
-end
-if strcmp(vel_type,'civ2') && testder
+if strcmp(vel_type,'civ2') && strcmp(FieldRequest,'derivatives')
     vel_type='filter2';
-elseif strcmp(vel_type,'civ1') && testder
+elseif strcmp(vel_type,'civ1') && strcmp(FieldRequest,'derivatives')
     vel_type='filter1';
 end
@@ -153 +172 @@
             vel_type='filter2';
         case {4,5}% civ2 available but not filter2
-            if testder% derivatives needed
+            if strcmp(FieldRequest,'derivatives')% derivatives needed
                 vel_type='filter1';
             else
@@ -171 +190 @@
             'Civ1_X','Civ1_Y','Civ1_Z','Civ1_U','Civ1_V','Civ1_W','Civ1_C','Civ1_F','Civ1_FF'};
         role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
-        units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''};
+    %    units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''};
     case 'filter1'
         var={'X','Y','Z','U','V','W','C','F','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbSites';...
@@ -178 +197 @@
         role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag','coord_tps','vector_x_tps',...
             'vector_y_tps','vector_z_tps','ancillary','ancillary'};
-        units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''};
+     %   units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''};
     case 'civ2'
         var={'X','Y','Z','U','V','W','C','F','FF';...
             'Civ2_X','Civ2_Y','Civ2_Z','Civ2_U','Civ2_V','Civ2_W','Civ2_C','Civ2_F','Civ2_FF'};
         role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag'};
-        units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''};
+      %  units={'pixel','pixel','pixel','pixel','pixel','pixel','','',''};
     case 'filter2'
         var={'X','Y','Z','U','V','W','C','F','FF','Coord_tps','U_tps','V_tps','W_tps','SubRange','NbSites';...
@@ -190 +209 @@
         role={'coord_x','coord_y','coord_z','vector_x','vector_y','vector_z','ancillary','warnflag','errorflag','coord_tps','vector_x_tps',...
             'vector_y_tps','vector_z_tps','ancillary','ancillary'};
-        units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''};
+       % units={'pixel','pixel','pixel','pixel','pixel','pixel','','','','pixel','pixel','pixel','pixel','pixel',''};
+end
+if ~strcmp(FieldRequest,'derivatives')
+    var=var(:,1:9);%suppress tps if not needed
 end
 vel_type_out=vel_type;

trunk/src/sub_field.m

@@ -7 +7 @@
 % when scalar and vectors are combined, the fields are just merged in a single matlab structure for common visualisation
 %-----------------------------------------------------------------------
-% function SubData=sub_field(Field,Field_1)
+% function SubData=sub_field(Field,XmlData,Field_1)
 %
 % OUPUT: 
@@ -16 +16 @@
 % Field_1:matlab structure representing the second field
 
-function [SubData,errormsg]=sub_field(Field,Field_1)
+function SubData=sub_field(Field,XmlData,Field_1)
 
+SubData=[];
+if strcmp(Field,'*')
+    return
+end
+if nargin<3
+    SubData=Field;
+    return
+end
 %% global attributes
-test_attr=0;
+SubData.ListGlobalAttribute={};%default
+%transfer global attributes of Field
 if isfield(Field,'ListGlobalAttribute')
     SubData.ListGlobalAttribute=Field.ListGlobalAttribute;
@@ -26 +35 @@
         SubData.(AttrName)=Field.(AttrName);
     end
-    test_attr=1;
 end
+%transfer global attributes of Field_1
 if isfield(Field_1,'ListGlobalAttribute')
     for ilist=1:numel(Field_1.ListGlobalAttribute)
-        test_1=1;
         AttrName=Field_1.ListGlobalAttribute{ilist};
-        if test_attr
-            for i_prev=1:numel(Field.ListGlobalAttribute)
-                if isequal(Field.ListGlobalAttribute{i_prev},AttrName)
-                    test_1=0; %attribute already written
-                    eval(['Val=Field.' AttrName ';'])                  
-                    eval(['Val_1=Field_1.' AttrName ';'])
-                    if isequal(Val,Val_1)           
-                        break% data already written
-                    else
-                        eval(['SubData.' AttrName '_1=Field_1.' AttrName ';']) 
-                    end
-                end 
-            end
+        AttrNameNew=AttrName;
+        while ~isempty(find(strcmp(AttrNameNew,SubData.ListGlobalAttribute)))&&~isequal(Field_1.(AttrNameNew),Field.(AttrNameNew))
+            AttrNameNew=[AttrNameNew '_1'];
         end
-        if test_1
-            eval(['SubData.' AttrName '=Field_1.' AttrName ';']) 
+        if ~isfield(Field,AttrName) || ~isequal(Field_1.(AttrName),Field.(AttrName))
+        SubData.ListGlobalAttribute=[SubData.ListGlobalAttribute {AttrNameNew}];
+        SubData.(AttrNameNew)=Field_1.(AttrName);
         end
     end
@@ -53 +52 @@
 
 %% variables
+%reproduce variables of the first field and list its dimensions
 SubData.ListVarName=Field.ListVarName;
 SubData.VarDimName=Field.VarDimName;
@@ -58 +58 @@
     SubData.VarAttribute=Field.VarAttribute;
 end
-%reproduce the first field Field by default
-for ivar=1:numel(Field.ListVarName)
-   VarName=Field.ListVarName{ivar};
-   SubData.(VarName)=Field.(VarName);
+ListDimName={};
+for ilist=1:numel(Field.ListVarName)
+    VarName=Field.ListVarName{ilist};
+    SubData.(VarName)=Field.(VarName);
+    SubData.VarAttribute{ilist}.CheckSub=0;
+    DimCell=Field.VarDimName{ilist};
+    if ischar(DimCell)
+        DimCell={DimCell};
+    end
+    for idim=1:numel(DimCell)
+        if isempty(find(strcmp(DimCell{idim},ListDimName)))
+            ListDimName=[ListDimName DimCell(idim)];
+        end
+    end
 end
 
-%% check the two input fields     
-[CellVarIndex,NbDim,VarTypeCell,errormsg]=find_field_cells(Field);
-if ~isempty(errormsg)
-    errormsg=['invalid  first input to sub_field:' errormsg];
-    return
-end
-[CellVarIndex_1,NbDim_1,VarTypeCell_1,errormsg]=find_field_cells(Field_1);
-if ~isempty(errormsg)
-    errormsg=['invalid second input to sub_field:' errormsg];
-    return
-end
-iselect=find(NbDim==2);
-for icell=iselect
-    if ~isempty(VarTypeCell{icell}.coord_tps)
-        NbDim(icell)=0;
+%% field request
+FieldRequest=cell(size(Field.ListVarName));
+for ilist=1:numel(Field.VarAttribute)
+    if isfield(Field.VarAttribute{ilist},'FieldRequest')
+        FieldRequest{ilist}=Field.VarAttribute{ilist}.FieldRequest;
     end
 end
-iselect=find(NbDim==2);
-if ~isequal(numel(iselect),1)
-    errormsg='invalid  first input to sub_field: it must  contain a single 2D field cell';
-    return
-end
-iselect_1=find(NbDim_1==2);
-for icell=iselect_1
-    if ~isempty(VarTypeCell_1{icell}.coord_tps)
-        NbDim_1(icell)=0;
+FieldRequest_1=cell(size(Field_1.ListVarName));
+for ilist=1:numel(Field_1.VarAttribute)
+    if isfield(Field_1.VarAttribute{ilist},'FieldRequest')
+        FieldRequest_1{ilist}=Field_1.VarAttribute{ilist}.FieldRequest;
     end
 end
-iselect_1=find(NbDim_1==2);
-if ~isequal(numel(iselect_1),1)
-    errormsg='invalid  second input to sub_field: it must  contain a single 2D field cell';
-    return
-end
-VarType=VarTypeCell{iselect};
-VarType_1=VarTypeCell_1{iselect_1};
-testX=~isempty(VarType.coord_x)&& ~isempty(VarType.coord_y);%unstructured coordiantes
-testX_1=~isempty(VarType_1.coord_x)&& ~isempty(VarType_1.coord_y);%unstructured coordiantes
-testU=~isempty(VarType.vector_x)&& ~isempty(VarType.vector_y);%vector field
-testU_1=~isempty(VarType_1.vector_x)&& ~isempty(VarType_1.vector_y);%vector field
-testfalse_1=~isempty(VarType_1.errorflag);
-ivar_C=[VarType.scalar VarType.image VarType.color VarType.ancillary]; %defines index (indices) for the scalar or ancillary fields
-if numel(ivar_C)>1
-    errormsg='too many scalar fields in the first input of sub_field.m';
-    return
-end
-ivar_C_1=[VarType_1.scalar VarType_1.image VarType_1.color VarType_1.ancillary]; %defines index (indices) for the scalar or ancillary fields
-if numel(ivar_C_1)>1
-    errormsg='too many scalar fields in the second input of sub_field.m';
-    return
+
+%% rename the dimensions of the second field if identical to those of the first
+for ilist=1:numel(Field_1.VarDimName)
+    DimCell=Field_1.VarDimName{ilist};
+    if ischar(DimCell)
+        DimCell={DimCell};
+    end
+    for idim=1:numel(DimCell)
+        ind_dim=find(strcmp(DimCell{idim},ListDimName));
+        if ~isempty(ind_dim)
+            if ischar(Field_1.VarDimName{ilist})
+                Field_1.VarDimName{ilist}=Field_1.VarDimName(ilist);
+            end
+            Field_1.VarDimName{ilist}{idim}=[ListDimName{ind_dim} '_1'];
+        end
+    end
 end
 
-%% substract two vector fields or two scalars
-if (testU && testU_1) || (~testU && ~testU_1)
-   %check coincidence in positions
-   %unstructured coordinates for the first field
-   if testX  
-       XName=Field.ListVarName{VarType.coord_x};
-       YName=Field.ListVarName{VarType.coord_y};
-       vec_X=Field.(XName);
-       vec_Y=Field.(YName);
-       nbpoints=numel(vec_X);
-       vec_X=reshape(vec_X,nbpoints,1);
-       vec_Y=reshape(vec_Y,nbpoints,1);
-       if testX_1 %unstructured coordinates for the second field
-            X_1_Name=Field_1.ListVarName{VarType_1.coord_x};
-            Y_1_Name=Field_1.ListVarName{VarType_1.coord_y};
-            eval(['vec_X_1=Field_1.' X_1_Name ';']) 
-            eval(['vec_Y_1=Field_1.' Y_1_Name ';'])
+%look for coordinates common to Field in Field_1
+ind_remove=zeros(size(Field_1.ListVarName));
+for ilist=1:numel(Field_1.ListVarName)
+    if ischar(Field_1.VarDimName{ilist})||numel(Field_1.VarDimName{ilist})==1
+        OldDim=Field_1.VarDimName{ilist};
+        if ischar(OldDim)
+            OldDim=Field_1.VarDimName(ilist);
+        end
+        VarVal=Field_1.(Field_1.ListVarName{ilist});
+        for i1=1:numel(Field.ListVarName)
+            if (isempty(FieldRequest{i1})&&isempty(FieldRequest_1{ilist})||strcmp(FieldRequest{i1},FieldRequest_1{ilist})) && isequal(Field.(Field.ListVarName{i1}),VarVal)
+               ind_remove(ilist)=1;
+               NewDim=Field.VarDimName{i1};
+               if ischar(NewDim)
+                   NewDim={NewDim};
+               end
+               Field_1.VarDimName=regexprep_r(Field_1.VarDimName,OldDim{1},NewDim{1});
+            end
+        end
+    end
+end
+Field_1.ListVarName(find(ind_remove))=[];%removes these redondent coordinates
+Field_1.VarDimName(find(ind_remove))=[];
+Field_1.VarAttribute(find(ind_remove))=[];
 
-       else   %structured coordinates for the second field
-           y_1_Name=Field_1.ListVarName{VarType_1.coord(1)};
-           x_1_Name=Field_1.ListVarName{VarType_1.coord(2)};
-           eval(['y_1=Field_1.' y_1_Name ';']) 
-           eval(['x_1=Field_1.' x_1_Name ';'])
-           if isequal(numel(x_1),2)  
-               x_1=linspace(x_1(1),x_1(2),nbpoints_x_1);
-           end
-           if isequal(numel(y_1),2)  
-               y_1=linspace(y_1(1),y_1(2),nbpoints_y_1);
-           end
-           [vec_X_1,vec_Y_1]=meshgrid(x_1,y_1);
-       end
-       vec_X_1=reshape(vec_X_1,[],1);
-       vec_Y_1=reshape(vec_Y_1,[],1);
-       if testfalse_1
-           FFName_1=Field_1.ListVarName{VarType_1.errorflag};          
-           eval(['vec_FF_1=Field_1.' FFName_1 ';']) 
-           vec_FF_1=reshape(vec_FF_1,[],1);
-           indsel=find(~vec_FF_1);
-           vec_X_1=vec_X_1(indsel);
-           vec_Y_1=vec_Y_1(indsel);
-       end
-       if testU % vector fields
-            U_1_Name=Field_1.ListVarName{VarType_1.vector_x};
-            V_1_Name=Field_1.ListVarName{VarType_1.vector_y};
-            eval(['vec_U_1=Field_1.' U_1_Name ';']) 
-            eval(['vec_V_1=Field_1.' V_1_Name ';'])
-            nbpoints_x_1=size(vec_U_1,2);
-            nbpoints_y_1=size(vec_U_1,1);
-            vec_U_1=reshape(vec_U_1,nbpoints_x_1*nbpoints_y_1,1);
-            vec_V_1=reshape(vec_V_1,nbpoints_x_1*nbpoints_y_1,1);
-            if testfalse_1
-                vec_U_1=vec_U_1(indsel);
-                vec_V_1=vec_V_1(indsel);
-            end            
-       else %(~testU && ~testU_1)
-           A_1_Name=Field_1.ListVarName{ivar_C_1};
-           eval(['vec_A_1=Field_1.' A_1_Name ';'])   
-           nbpoints_x_1=size(vec_A_1,2);
-           nbpoints_y_1=size(vec_A_1,1);%TODO: use a faster interpolation method for a regular grid (size(x)=2)
-           vec_A_1=reshape(vec_A_1,nbpoints_x_1*nbpoints_y_1,1);
-           if testfalse_1
-                vec_A_1=vec_A_1(indsel);
-           end
-       end
-
-       if ~isequal(vec_X_1,vec_X) && ~isequal(vec_Y_1,vec_Y) % if the unstructured positions are not the same
-           if testU
-               vec_U_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_U_1,vec_X,vec_Y);  %interpolate vectors in the second field
-               vec_V_1=griddata_uvmat(vec_X_1,vec_Y_1,vec_V_1,vec_X,vec_Y);  %interpolate vectors in the second field   
-           else
-               vec_A_1=griddata_uvmat(vec_X_1,vec_Y_1,double(vec_A_1),vec_X,vec_Y);  %interpolate vectors in the second field
-           end
-       end 
-       if testU
-           UName=Field.ListVarName{VarType.vector_x};
-           VName=Field.ListVarName{VarType.vector_y};  
-           eval(['vec_U=Field.' UName ';']) 
-           eval(['vec_V=Field.' VName ';'])       
-           vec_U=reshape(vec_U,numel(vec_U),1);
-           vec_V=reshape(vec_V,numel(vec_V),1);
-           eval(['SubData.' UName '=vec_U-vec_U_1;'])
-           eval(['SubData.' VName '=vec_V-vec_V_1;'])
-       else
-           AName=Field.ListVarName{ivar_C};
-          % size(Field.vort)
-           eval(['SubData.' AName '=Field.' AName '-vec_A_1;'])
-       end
-   else  %structured coordiantes
-       XName=Field.ListVarName{VarType.coord(2)};
-       YName=Field.ListVarName{VarType.coord(1)};
-       eval(['x=Field.' XName ';']) 
-       eval(['y=Field.' YName ';'])
-       if testX_1 %unstructured coordinates for the second field
-           errormsg='the second input scalar is not on a regular grid: comparison option not implemented';
-           return
-       else
-           XName_1=Field.ListVarName{VarType_1.coord(2)};
-           YName_1=Field.ListVarName{VarType_1.coord(1)};
-           eval(['x_1=Field_1.' XName_1 ';']) 
-           eval(['y_1=Field_1.' YName_1 ';'])
-       end
-       if testU % vector fields
-           UName=Field.ListVarName{VarType.vector_x};
-           VName=Field.ListVarName{VarType.vector_y};
-           U_1_Name=Field_1.ListVarName{VarType_1.vector_x};
-           V_1_Name=Field_1.ListVarName{VarType_1.vector_y};
-           eval(['U_1=Field_1.' U_1_Name ';']) 
-           eval(['V_1=Field_1.' V_1_Name ';'])
-           if ~isequal(x_1,x)||~isequal(y_1,y)
-                [X_1,Y_1]=meshgrid(x_1,y_1);
-                U_1 =interp2(X_1,Y_1,U_1,x,y');
-                V_1 =interp2(X_1,Y_1,V_1,x,y');
-           end
-           eval(['SubData.' UName '=Field.' UName '-U_1;'])
-           eval(['SubData.' VName '=Field.' VName '-V_1;'])
-       else
-           AName=Field.ListVarName{ivar_C};
-           A_1_Name=Field_1.ListVarName{ivar_C_1};
-           eval(['A_1=double(Field_1.' A_1_Name ');'])
-           if ~isequal(x_1,x)||~isequal(y_1,y)
-                [X_1,Y_1]=meshgrid(x_1,y_1);
-                A_1 =interp2(X_1,Y_1,A_1,x,y');
-           end
-           eval(['SubData.' AName '=double(Field.' AName ')-A_1;'])
-       end
-   end
+%append the other variables of the second field, modifying their name if needed
+for ilist=1:numel(Field_1.ListVarName)
+    VarName=Field_1.ListVarName{ilist};
+    ind_prev=find(strcmp(VarName,Field.ListVarName));
+    if isempty(ind_prev)% variable name does not exist in Field
+        VarNameNew=VarName;
+    else  % variable name exists in Field     
+            VarNameNew=[VarName '_1'];   
+    end
+        SubData.ListVarName=[SubData.ListVarName {VarNameNew}];
+        SubData.VarDimName=[SubData.VarDimName Field_1.VarDimName(ilist)];
+        SubData.(VarNameNew)=Field_1.(VarName);
+        SubData.VarAttribute=[SubData.VarAttribute Field_1.VarAttribute(ilist)];
+        SubData.VarAttribute{end}.CheckSub=1;% mark that the field needs to be substracted
 end
 
-% merge a vector field and a scalar as second input
-if testU && ~testU_1
-    AName_1=Field_1.ListVarName{ivar_C_1};
-    if isfield(Field_1,'VarAttribute') & numel(Field_1.VarAttribute)>=ivar_C_1
-        AAttr=Field_1.VarAttribute{ivar_C_1} ;
-    else
-        AAttr=[];
+%append the other variables of the second field, modifying their name if needed
+
+%% substrat fields when possible
+[CellVarIndex,NbDim,CellVarType,errormsg]=find_field_cells(SubData);
+ind_remove=zeros(size(SubData.ListVarName));
+ivar=[];
+ivar_1=[];
+for icell=1:numel(CellVarIndex)
+    if ~isempty(CellVarIndex{icell})
+        if isfield(CellVarType{icell},'scalar') && numel(CellVarType{icell}.scalar)==2 && SubData.VarAttribute{CellVarType{icell}.scalar(2)}.CheckSub;
+            ivar=[ivar CellVarType{icell}.scalar(1)];
+            ivar_1=[ivar_1 CellVarType{icell}.scalar(2)];
+        end
+        if isfield(CellVarType{icell},'vector_x') && numel(CellVarType{icell}.vector_x)==2 && SubData.VarAttribute{CellVarType{icell}.vector_x(2)}.CheckSub;
+            ivar=[ivar CellVarType{icell}.vector_x(1)];
+            ivar_1=[ivar_1 CellVarType{icell}.vector_x(2)];
+        end
+        if isfield(CellVarType{icell},'vector_y') && numel(CellVarType{icell}.vector_y)==2 && SubData.VarAttribute{CellVarType{icell}.vector_y(2)}.CheckSub;
+            ivar=[ivar CellVarType{icell}.vector_y(1)];
+            ivar_1=[ivar_1 CellVarType{icell}.vector_y(2)];
+        end
     end
-    if testX_1 %unstructured coordinate
-       XName_1=Field_1.ListVarName{VarType_1.coord_x};
-       YName_1=Field_1.ListVarName{VarType_1.coord_y};
-       DimCell=Field_1.VarDimName([VarType_1.coord_x VarType_1.coord_y ]);
-       if isfield(Field_1,'VarAttribute') 
-           if numel(Field_1.VarAttribute)>=VarType_1.coord_x
-                XAttr=Field_1.VarAttribute{VarType_1.coord_x}; 
-           else
-                XAttr=[];
-           end
-           if numel(Field_1.VarAttribute)>=VarType_1.coord_y
-               YAttr=Field_1.VarAttribute{VarType_1.coord_y}; 
-           else
-               YAttr=[];
-           end
-           SubData.VarAttribute=[SubData.VarAttribute {XAttr} {YAttr}];
-       end
-    else
-       XName_1=Field_1.ListVarName{VarType_1.coord(2)};
-       YName_1=Field_1.ListVarName{VarType_1.coord(1)};
-%        DimCell=[{YName_1} {XName_1}];
-       if isfield(Field_1,'VarAttribute') 
-           if numel(Field_1.VarAttribute)>=VarType_1.coord(2)
-                XAttr=Field_1.VarAttribute{VarType_1.coord(2)} ;
-           else
-                XAttr=[];
-           end
-           if numel(Field_1.VarAttribute)>=VarType_1.coord(1)
-               YAttr=Field_1.VarAttribute{VarType_1.coord(1)} ;
-           else
-               YAttr=[];
-           end
-           SubData.VarAttribute=[SubData.VarAttribute {YAttr} {XAttr}];
-       end
-    end  
-    %look for previously used variable names
-    XName_1_1=XName_1;%default
-    YName_1_1=YName_1;%default
-    AName_1_1=AName_1;%default
-    for iprev=1:numel(SubData.ListVarName)
-        switch SubData.ListVarName{iprev}
-            case XName_1
-                XName_1_1=[XName_1 '_1'];
-            case YName_1
-                YName_1_1=[YName_1 '_1'];
-            case AName_1
-                AName_1_1=[AName_1 '_1']; 
-        end
-    end     
-    if ~testX_1% if the second field has structured coordinates
-        ilist_1=find(strcmp(AName_1,Field_1.ListVarName));
-        DimCell_1=Field_1.VarDimName{ilist_1};
-        if numel(DimCell_1)>2
-             DimCell={XName_1_1,YName_1_1, [{YName_1_1,XName_1_1} DimCell_1(3:end)]};
-        else
-          DimCell={XName_1_1,YName_1_1, {YName_1_1,XName_1_1}};
-        end
-    else
-        DimCell=[DimCell Field_1.VarDimName(ivar_C_1)];
-    end
-    SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {AName_1_1}];
-    % check that a  different dimension name is used for the two fields
-     if testX_1% if the second field has unstructured coordinates
-        for icell=1:numel(DimCell)
-            if isequal(DimCell{icell}{1},SubData.VarDimName{1}{1})
-                DimCell{icell}{1}=[DimCell{icell}{1} '_1'];
-            end
-        end
-     end
-    SubData.VarDimName=[SubData.VarDimName DimCell];
-    if isfield(Field_1,'VarAttribute')
-        SubData.VarAttribute=[SubData.VarAttribute {AAttr}];
-    end
-    eval(['SubData.' XName_1_1 '=Field_1.' XName_1 ';'])
-    eval(['SubData.' YName_1_1 '=Field_1.' YName_1 ';'])
-    eval(['SubData.' AName_1_1 '=Field_1.' AName_1 ';'])
 end
-
-%merge a scalar as the first input and a vector field as second input
-if ~testU && testU_1
-    UName_1=Field_1.ListVarName{VarType_1.vector_x};
-    VName_1=Field_1.ListVarName{VarType_1.vector_y};
-    UAttr=Field_1.VarAttribute{VarType_1.vector_x};
-    VAttr=Field_1.VarAttribute{VarType_1.vector_y};
-    if testX_1 %unstructured coordinate for the second field
-       XName_1=Field_1.ListVarName{VarType_1.coord_x};
-       YName_1=Field_1.ListVarName{VarType_1.coord_y};
-       
-       XAttr=Field_1.VarAttribute{VarType_1.coord_x};
-       YAttr=Field_1.VarAttribute{VarType_1.coord_y};
-%        SubData.ListVarName=[SubData.ListVarName {XName_1} {YName_1}];
-       DimCell=Field_1.VarDimName([VarType_1.coord_x VarType_1.coord_y ]);
-    else
-       XName_1=Field_1.ListVarName{VarType_1.coord(2)};
-       YName_1=Field_1.ListVarName{VarType_1.coord(1)};
-       if numel(Field_1.VarAttribute)>=VarType_1.coord(2)
-           XAttr=Field_1.VarAttribute{VarType_1.coord(2)};
-       else
-           XAttr=[];
-       end
-       if numel(Field_1.VarAttribute)>=VarType_1.coord(1)
-           YAttr=Field_1.VarAttribute{VarType_1.coord(1)};
-       else
-           YAttr=[];
-       end     
-    end  
-    %check for the existence of the same  variable name
-    XName_1_1=XName_1; %default
-    YName_1_1=YName_1; %default
-    UName_1_1=UName_1; %default
-    VName_1_1=VName_1; %default
-    for iprev=1:numel(SubData.ListVarName)
-        switch SubData.ListVarName{iprev}
-            case XName_1
-                XName_1_1=[XName_1 '_1'];
-            case YName_1
-                YName_1_1=[YName_1 '_1'];
-            case UName_1
-                UName_1_1=[UName_1 '_1'];
-            case VName_1
-                VName_1_1=[VName_1 '_1']; 
-        end
-    end    
+for imod=1:numel(ivar)
+        VarName=SubData.ListVarName{ivar(imod)};
+        VarName_1=SubData.ListVarName{ivar_1(imod)};
+        SubData.(VarName)=SubData.(VarName)-SubData.(VarName_1);
+        ind_remove(ivar_1(imod))=1;
+end
+SubData.ListVarName(find(ind_remove))=[];
+SubData.VarDimName(find(ind_remove))=[];
+SubData.VarAttribute(find(ind_remove))=[];
+        
     
-    if ~testX_1
-          DimCell=[{XName_1_1} {YName_1_1}];
-    end
-    SubData.ListVarName=[SubData.ListVarName {XName_1_1} {YName_1_1} {UName_1_1} {VName_1_1}];
-    DimCell=[DimCell Field_1.VarDimName([VarType_1.vector_x VarType_1.vector_y ])];
-    SubData.VarDimName=[SubData.VarDimName DimCell];
-    if isfield(SubData,'VarAttribute')
-        if ~(numel(SubData.VarAttribute)==numel(SubData.ListVarName))
-            for ivar=numel(SubData.VarAttribute)+1:numel(SubData.ListVarName)-4
-                SubData.VarAttribute{ivar}=[];
-            end
-        end
-        SubData.VarAttribute=[SubData.VarAttribute {XAttr} {YAttr} {UAttr} {VAttr}];
-    end
-    SubData.(XName_1_1)=Field_1.(XName_1);
-    SubData.(YName_1_1)=Field_1.(YName_1);
-    SubData.(UName_1_1)=Field_1.(UName_1);
-    SubData.(VName_1_1)=Field_1.(VName_1); 
-end
-  

trunk/src/uvmat.m

@@ -70 +70 @@
 %                        (calc_tps.m)               calculate tps coefficients (for filter projection or spatial derivatives).
 %                            |
-%                       (calc_field.m)               calculate field 
-%                            |
 %                       UvData.Field-------------->histogram
 %               _____________|____________
 %              |                          |                    
-%        proj_field.m               proj_field.m       project the field on the projection objects              
+%        proj_field.m               proj_field.m       project the field on the projection objects (use calc_field.m)           
 %              |                          |
 %         UvData.PlotAxes          ViewData.PlotAxes (on view_field)
 %              |                          |
 %       plot_field.m (uvmat)       plot_field.m (view_field)      plot the projected fields
-%
-% rmq: calc_field can be performed instead at the level of proj_field when needed 
 %
 %
@@ -209 +205 @@
 
 %% TRANSFORM menu: builtin fcts
-transform_menu={'';'phys';'px';'phys_polar'};
+transform_menu={'';'sub_field';'phys';'phys_polar'};
 UvData.OpenParam.NbBuiltin=numel(transform_menu); %number of functions
 path_list=(num2cell(blanks(UvData.OpenParam.NbBuiltin)))';%initialize a cell array of nbvar blanks
@@ -299 +295 @@
 
 %% plot input field if exists
-if testinputfield
-    %delete drawn objects
-    hother=findobj(handles.PlotAxes,'Tag','proj_object');%find all the proj objects
-    for iobj=1:length(hother)
-        delete_object(hother(iobj))
-    end  
-    if isempty(inputfile)
-        errormsg=refresh_field(handles,[],[],[],[],[],[],{Field});
-        set(handles.MenuTools,'Enable','on')
-        set(handles.OBJECT_txt,'Visible','on')
-        set(handles.edit_object,'Visible','on')
-%         set(handles.ListObject_1,'Visible','on')
-        set(handles.frame_object,'Visible','on')
-        if ~isempty(errormsg)
-            msgbox_uvmat('ERROR',errormsg)
-        end
-    end
-end
+% if testinputfield
+%     %delete drawn objects
+%     hother=findobj(handles.PlotAxes,'Tag','proj_object');%find all the proj objects
+%     for iobj=1:length(hother)
+%         delete_object(hother(iobj))
+%     end  
+%     if isempty(inputfile)
+%         errormsg=refresh_field(handles,[],[],[],[],[],[],{Field});
+%         set(handles.MenuTools,'Enable','on')
+%         set(handles.OBJECT_txt,'Visible','on')
+%         set(handles.edit_object,'Visible','on')
+% %         set(handles.ListObject_1,'Visible','on')
+%         set(handles.frame_object,'Visible','on')
+%         if ~isempty(errormsg)
+%             msgbox_uvmat('ERROR',errormsg)
+%         end
+%     end
+% end
 
 set_vec_col_bar(handles) %update the display of color code for vectors
@@ -679 +675 @@
             set(handles.j1,'String',num2stra(j1,NomType));
             set(handles.j2,'String',num2stra(j2,NomType));
-        else %read the current field index if the second file series is opened
-            i1=str2num(get(handles.i1,'String'));
-            i2=str2num(get(handles.i2,'String'));
-            j1=stra2num(get(handles.j1,'String'));
-            j2=stra2num(get(handles.j2,'String'));
+        else %read the current field index to synchronise with the first series
+            i1_s=str2num(get(handles.i1,'String'));
+            i2_0=str2num(get(handles.i2,'String'));
+            if ~isempty(i2_0)
+                i2_s=i2_0;
+            else
+               i2_s=i2; 
+            end
+            j1_0=stra2num(get(handles.j1,'String'));
+            if ~isempty(j1_0)
+                j1_s=j1_0;
+            else
+                j1_s=j1;
+            end
+            j2_0=stra2num(get(handles.j2,'String'));
+            if ~isempty(j2_0)
+                j2_s=j2_0;
+            else
+                j2_s=j2;
+            end
         end
         
@@ -705 +716 @@
             end
             %updtate the indices of the second field series to correspond to the newly opened one
-            FileName_1=fullfile_uvmat(Input.RootPath_1,Input.SubDir_1,Input.RootFile_1,Input.FileExt_1,Input.NomType_1,i1,i2,j1,j2);
+            FileName_1=fullfile_uvmat(Input.RootPath_1,Input.SubDir_1,Input.RootFile_1,Input.FileExt_1,Input.NomType_1,i1_s,i2_s,j1_s,j2_s);
             if exist(FileName_1,'file')
+                FileIndex_1=fullfile_uvmat('','','','',Input.NomType_1,i1_s,i2_s,j1_s,j2_s);
+            else
                 FileIndex_1=fullfile_uvmat('','','','',Input.NomType_1,i1,i2,j1,j2);
-                set(handles.FileIndex_1,'String',FileIndex_1)
-            else
                 msgbox_uvmat('WARNING','unable to synchronise the indices of the two series')
             end
+            set(handles.FileIndex_1,'String',FileIndex_1)
         end
         
@@ -757 +769 @@
 % --- Update information about a new field series (indices to scan, timing,
 %     calibration from an xml file, then refresh current plots
+
 function update_rootinfo(handles,i1_series,i2_series,j1_series,j2_series,FileType,VideoObject,index)
 %------------------------------------------------------------------------
@@ -943 +956 @@
         end
         if ~get(handles.CheckFixLimits,'Value')
-            set(handles.transform_fct,'Value',2); % phys transform by default if fixedLimits is off
+            set(handles.transform_fct,'Value',3); % phys transform by default if fixedLimits is off
         end
         if isfield(GeometryCalib,'SliceCoord')            
@@ -1014 +1027 @@
     end
 end
-[ref_j,ref_i]=find(i1_series);
+[ref_j,ref_i]=find(squeeze(i1_series(1,:,:)));
 if ~isempty(j1_series) 
         state_j='on';
         if index==1
-            if isequal(ref_i,ref_i(1)*ones(size(ref_j)))% if ref_i is always equal to its first vzlue
+            if isequal(ref_i,ref_i(1)*ones(size(ref_j)))% if ref_i is always equal to its first value
                 scan_option='j'; %scan j indext               
             end 
@@ -1057 +1070 @@
 
 %% apply the effect of the transform fct and view the field  
+transform=get(handles.path_transform,'UserData');
+if index==2 && (~isa(transform,'function_handle')||nargin(transform)<3)
+    set(handles.transform_fct,'value',2); % set transform to sub_field if the current fct doe not accept two input fields
+end
 transform_fct_Callback([],[],handles)
-%run0_Callback([],[], handles); %view field
 mask_test=get(handles.CheckMask,'value');
 if mask_test
@@ -1106 +1122 @@
 function i1_Callback(hObject, eventdata, handles)
 %------------------------------------------------------------------------
-set(handles.i1,'BackgroundColor',[0.7 0.7 0.7])
+update_ij(handles,1)
+
+%------------------------------------------------------------------------
+function i2_Callback(hObject, eventdata, handles)
+%------------------------------------------------------------------------
+update_ij(handles,2)
+
+%------------------------------------------------------------------------
+function j1_Callback(hObject, eventdata, handles)
+%------------------------------------------------------------------------
+update_ij(handles,3)
+
+%------------------------------------------------------------------------
+function j2_Callback(hObject, eventdata, handles)
+%------------------------------------------------------------------------
+update_ij(handles,4)
+
+%------------------------------------------------------------------------
+%--- update the index display after action on edit boxes i1, i2, j1 or j2
+function update_ij(handles,index_rank)
 NomType=get(handles.NomType,'String');
-i1=stra2num(get(handles.i1,'String'));
-i2=stra2num(get(handles.i2,'String'));
-j1=stra2num(get(handles.j1,'String'));
-j2=stra2num(get(handles.j2,'String'));
-indices=fullfile_uvmat('','','','',NomType,i1,i2,j1,j2);
-%indices=name_generator('',num1,num_a,'',NomType,1,num2,num_b,'');
+indices=get(handles.FileIndex,'String');
+[tild,tild,tild,i1,i2,j1,j2]=fileparts_uvmat(indices);% the indices for the second series taken from FileIndex
+switch index_rank
+    case 1
+        indices=fullfile_uvmat('','','','',NomType,stra2num(get(handles.i1,'String')),i2,j1,j2);
+        set(handles.i1,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation
+    case 2
+        indices=fullfile_uvmat('','','','',NomType,i1,stra2num(get(handles.i2,'String')),j1,j2);
+        set(handles.i2,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation
+    case 3
+        indices=fullfile_uvmat('','','','',NomType,i1,i2,stra2num(get(handles.j1,'String')),j2);
+        set(handles.j1,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation
+    case 4
+        indices=fullfile_uvmat('','','','',NomType,i1,i2,j1,stra2num(get(handles.j2,'String')));
+        set(handles.j2,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation
+end
 set(handles.FileIndex,'String',indices)
-set(handles.FileIndex,'BackgroundColor',[0.7 0.7 0.7])
-if get(handles.SubField,'Value')==1
+set(handles.FileIndex,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation
+% update the second index if relevant
+if strcmp(get(handles.FileIndex_1,'Visible'),'on')
     NomType_1=get(handles.NomType_1,'String');
-    indices=fullfile_uvmat('','','','',NomType,i1,i2,j1,j2);
-     %indices=name_generator('',num1,num_a,'',NomType_1,1,num2,num_b,'');
-     set(handles.FileIndex_1,'String',indices)
-     set(handles.FileIndex_1,'BackgroundColor',[0.7 0.7 0.7])
-end
-
-%------------------------------------------------------------------------
-function i2_Callback(hObject, eventdata, handles)
-%------------------------------------------------------------------------
-set(handles.i2,'BackgroundColor',[0.7 0.7 0.7])
-i1_Callback(hObject, eventdata, handles)
-
-%------------------------------------------------------------------------
-function j1_Callback(hObject, eventdata, handles)
-%------------------------------------------------------------------------
-set(handles.j1,'BackgroundColor',[0.7 0.7 0.7])
-i1_Callback(hObject, eventdata, handles)
-
-%------------------------------------------------------------------------
-function j2_Callback(hObject, eventdata, handles)
-%------------------------------------------------------------------------
-set(handles.j2,'BackgroundColor',[0.7 0.7 0.7])
-i1_Callback(hObject, eventdata, handles)
+    indices_1=get(handles.FileIndex_1,'String');
+    [tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(indices_1);% the indices for the second series taken from FileIndex_1
+    switch index_rank
+        case 1
+            indices_1=fullfile_uvmat('','','','',NomType_1,stra2num(get(handles.i1,'String')),i2_1,j1_1,j2_1);
+        case 2
+            indices_1=fullfile_uvmat('','','','',NomType_1,i1_1,stra2num(get(handles.i2,'String')),j1_1,j2_1);
+        case 3
+            indices_1=fullfile_uvmat('','','','',NomType_1,i1_1,i2_1,stra2num(get(handles.j1,'String')),j2_1);
+        case 4
+            indices_1=fullfile_uvmat('','','','',NomType_1,i1_1,i2_1,j1_1,stra2num(get(handles.j2,'String')));
+    end
+    set(handles.FileIndex_1,'String',indices_1)
+    set(handles.FileIndex_1,'BackgroundColor',[0.7 0.7 0.7])% mark the edit box in grey, then RUN0 will mark it in white for confirmation
+end
+    
+%------------------------------------------------------------------------
 
 %------------------------------------------------------------------------
@@ -1523 +1563 @@
 if sub_value % a second input file has been entered 
      [InputFile.RootPath_1,InputFile.SubDir_1,InputFile.RootFile_1,InputFile.FileIndex_1,InputFile.FileExt_1,InputFile.NomType_1]=read_file_boxes_1(handles);   
-    [tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(InputFile.FileIndex_1);
+    [tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(InputFile.FileIndex_1);% the indices for the second series taken from FileIndex_1
 else
     filename_1=[];
@@ -1529 +1569 @@
 
 %% increment (or decrement) the field indices and update the input filename(s)
-CheckSearch=0;
 if isempty(increment)
-    CheckSearch=1;% search for the next available file
     set(handles.CheckFixPair,'Value',0)
 end
@@ -1545 +1583 @@
             i2_1=i2_1+increment;
         end
+        
     else % case of scanning along index j (burst numbers)
         j1=j1+increment;
@@ -1604 +1643 @@
     if ~isempty(errormsg),return,end
     siz=size(UvData.i1_series{1});
-    ref_indices=ref_i*siz(1)*siz(2)+ref_j*siz(1)+1:ref_i*siz(1)*siz(2)+(ref_j+1)*siz(1);   
+    ref_indices=ref_i*siz(1)*siz(2)+ref_j*siz(1)+1:ref_i*siz(1)*siz(2)+(ref_j+1)*siz(1);
     i1_subseries=UvData.i1_series{1}(ref_indices);
     ref_indices=ref_indices(i1_subseries>0);
     if isempty(ref_indices)% case of pairs (free index i)
-        ref_indices=ref_i*siz(1)*siz(2)+1:(ref_i+1)*siz(1)*siz(2); 
+        ref_indices=ref_i*siz(1)*siz(2)+1:(ref_i+1)*siz(1)*siz(2);
         i1_subseries=UvData.i1_series{1}(ref_indices);
         ref_indices=ref_indices(i1_subseries>0);
@@ -1623 +1662 @@
     if ~isempty(UvData.j2_series{1})
         j2=UvData.j2_series{1}(ref_indices(end));
-    end  
-    % case of a second file series (TODO:revise)
-    if numel(UvData.i1_series)>=2
-        i1_subseries=UvData.i1_series{2}(ref_i+1,ref_j+1,:);
-        i1_subseries=i1_subseries(i1_subseries>0);
-        i1_1=i1_subseries(end);
+    end
+    
+    % case of a second file series 
+    if sub_value
+        ref_i_1=i1_1;
+        if ~isempty(i2_1)
+            ref_i_1=floor((i1_1+i2_1)/2);% current reference index i
+        end
+        ref_j_1=1;
+        if ~isempty(j1_1)
+            ref_j_1=j1_1;
+            if ~isempty(j2_1)
+                ref_j_1=floor((j1_1+j2_1)/2);% current reference index j
+            end
+        end
+        if ~isempty(increment)
+            if get(handles.scan_i,'Value')==1% case of scanning along index i
+                ref_i_1=ref_i_1+increment;% increment the current reference index i
+            else % case of scanning along index j (burst numbers)
+                ref_j_1=ref_j_1+increment;% increment the current reference index j if scan_j option is used
+            end
+        else % free increment, synchronise the ref indices with the first series
+            ref_i_1=ref_i;
+            ref_j_1=ref_j; 
+        end
+        if numel(UvData.i1_series)==1
+            UvData.i1_series{2}=UvData.i1_series{1};
+            UvData.j1_series{2}=UvData.j1_series{1};
+            UvData.i2_series{2}=UvData.i2_series{1};
+            UvData.j2_series{2}=UvData.j2_series{1};
+        end
+        if ref_i_1<0
+            errormsg='minimum i index reached';
+        elseif ref_j_1<0
+            errormsg='minimum j index reached';
+        elseif ref_i_1+1>size(UvData.i1_series{2},3)
+            errormsg='maximum i index reached for the second series (reload the input file to update the index bound)';
+        elseif ref_j_1+1>size(UvData.i1_series{2},2)
+            errormsg='maximum j index reached for the second series(reload the input file to update the index bound)';
+        end
+        if ~isempty(errormsg),return,end
+        siz=size(UvData.i1_series{2});
+        ref_indices=ref_i_1*siz(1)*siz(2)+ref_j_1*siz(1)+1:ref_i_1*siz(1)*siz(2)+(ref_j_1+1)*siz(1);
+        i1_subseries=UvData.i1_series{2}(ref_indices);
+        ref_indices=ref_indices(i1_subseries>0);
+        if isempty(ref_indices)% case of pairs (free index i)
+            ref_indices=ref_i_1*siz(1)*siz(2)+1:(ref_i_1+1)*siz(1)*siz(2);
+            i1_subseries=UvData.i1_series{2}(ref_indices);
+            ref_indices=ref_indices(i1_subseries>0);
+        end
+        i1_1=UvData.i1_series{2}(ref_indices(end));
         if ~isempty(UvData.i2_series{2})
-            i2_subseries=UvData.i2_series{2}(ref_i+1,ref_j+1,:);
-            i2_subseries=i2_subseries(i2_subseries>0);
-            i2_1=i2_subseries(end);
+            i2_1=UvData.i2_series{2}(ref_indices(end));
         end
         if ~isempty(UvData.j1_series{2})
-            j1_subseries=UvData.j1_series{2}(ref_i+1,ref_j+1,:);
-            j1_subseries=j1_subseries(j1_subseries>0);
-            j1_1=j1_subseries(end);
+            j1_1=UvData.j1_series{2}(ref_indices(end));
         end
         if ~isempty(UvData.j2_series{2})
-            j2_subseries=UvData.j2_series{2}(ref_i+1,ref_j+1,:);
-            j2_subseries=j2_subseries(j2_subseries>0);
-            j2_1=j2_subseries(end);
-        end
+            j2_1=UvData.j2_series{1}(ref_indices(end));
+        end   
+    else% the second series (if needed) is the same file as the first
+        i1_1=i1;
+        i2_1=i2;
+        j1_1=j1;
+        j2_1=j2;
     end
 end
@@ -1652 +1735 @@
 
 %% refresh plots
-errormsg=refresh_field(handles,filename,filename_1,i1,i2,j1,j2);
+errormsg=refresh_field(handles,filename,filename_1,i1,i2,j1,j2,i1_1,i2_1,j1_1,j2_1);
 
 %% update the index counters if the index move is successfull
@@ -1846 +1929 @@
 set(handles.run0,'BackgroundColor',[1 1 0])%paint the command button in yellow
 drawnow
-[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
-filename=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
-filename_1=[];%default
+[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
+filename=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
+filename_1='';%default
+FileIndex_1='';
 if get(handles.SubField,'Value')
-    [RootPath_1,SubDir_1,RootFile_1,FileIndices_1,FileExt_1]=read_file_boxes_1(handles);
-    filename_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndices_1 FileExt_1];
+    [RootPath_1,SubDir_1,RootFile_1,FileIndex_1,FileExt_1]=read_file_boxes_1(handles);
+    filename_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1];
 end
 num_i1=stra2num(get(handles.i1,'String'));
@@ -1857 +1941 @@
 num_j1=stra2num(get(handles.j1,'String'));
 num_j2=stra2num(get(handles.j2,'String'));
-
-errormsg=refresh_field(handles,filename,filename_1,num_i1,num_i2,num_j1,num_j2);
+[tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(FileIndex_1);
+
+errormsg=refresh_field(handles,filename,filename_1,num_i1,num_i2,num_j1,num_j2,i1_1,i2_1,j1_1,j2_1);
 
 if ~isempty(errormsg)
@@ -1883 +1968 @@
 % Field: structure describing an optional input field (then replace the input file)
 
-function errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2,Field)
+function errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2,i1_1,i2_1,j1_1,j2_1)
 %------------------------------------------------------------------------
 
@@ -1923 +2008 @@
 if length(masknumber)>=z_index
     set(handles.masklevel,'Value',z_index)
+end
+
+%% test for need of tps
+check_proj_tps=0;
+if  (strcmp(UvData.FileType{1},'civdata')||strcmp(UvData.FileType{1},'civx'))
+    for iobj=1:numel(UvData.Object)
+        if isfield(UvData.Object{iobj},'ProjMode')&& strcmp(UvData.Object{iobj}.ProjMode,'filter')
+            check_proj_tps=1;
+            break
+        end
+    end
 end
 
@@ -1977 +2073 @@
 ParamIn.VelType=VelType;
 ParamIn.GUIName='get_field';
+end
+check_tps=0;         
+if strcmp(UvData.FileType{1},'civdata')&&~strcmp(ParamIn.FieldName,'velocity')&&~strcmp(ParamIn.FieldName,'get_field...') 
+       check_tps=1;%tps needed to get the requested field
 end
 [Field{1},ParamOut,errormsg] = read_field(FileName,UvData.FileType{1},ParamIn,frame_index);
@@ -2030 +2130 @@
             ParamIn_1=UvData.MovieObject{2};
                         if ~strcmp(NomType_1,'*')
-                frame_index_1=num_j1;%frame index for movies or multimage
+                frame_index_1=j1_1;%frame index for movies or multimage
             else
-                frame_index_1=num_i1;
+                frame_index_1=i1_1;
             end  
          case 'multimage'
             if ~strcmp(NomType_1,'*')
-                frame_index_1=num_j1;%frame index for movies or multimage
+                frame_index_1=j1_1;%frame index for movies or multimage
             else
-                frame_index_1=num_i1;
+                frame_index_1=i1_1;
             end   
         case 'vol' %TODO: update
@@ -2055 +2155 @@
     end
     test_keepdata_1=0;% test for keeping the previous stored data if the input files are unchanged
-    if ~isequal(NomType_1,'*')%in case of a series of files (not avi movie)
-        if isfield(UvData,'FileName_1')%&& isfield(UvData,'VelType_1') && isfield(UvData,'FieldName_1')
-            test_keepdata_1= strcmp(FileName_1,UvData.FileName_1) ;%&& strcmp(FieldName_1,UvData.FieldName_1);
-        end
+    if ~isequal(NomType_1,'*')&& isfield(UvData,'FileName_1')
+           test_keepdata_1= strcmp(FileName_1,UvData.FileName_1) ;%&& strcmp(FieldName_1,UvData.FieldName_1);
     end
     if test_keepdata_1
@@ -2068 +2166 @@
         ParamIn_1.VelType=VelType_1;
         ParamIn_1.GUIName='get_field_1';
-        end
+        end  
         [Field{2},ParamOut_1,errormsg] = read_field(Name,UvData.FileType{2},ParamIn_1,frame_index_1);
         if ~isempty(errormsg)
             errormsg=['error in reading ' FieldName_1 ' in ' FileName_1 ': ' errormsg];
             return
+        end
+        if isstruct(ParamOut_1)&&~strcmp(ParamOut_1.FieldName,'get_field...')&& (strcmp(UvData.FileType{2},'civdata')||strcmp(UvData.FileType{2},'civx'))...
+                &&~strcmp(ParamOut_1.FieldName,'velocity') && ~strcmp(ParamOut_1.FieldName,'get_field...')
+            if ~check_proj_tps
+             %   Field{2}=calc_field([{ParamOut_1.FieldName} {ParamOut_1.ColorVar}],Field{2});
+            end
         end
     end
@@ -2130 +2234 @@
         test_veltype_1=1;
         set(handles.VelType_1,'Visible','on')
-        menu=set_veltype_display(ParamOut_1.CivStage);
+        menu=set_veltype_display(ParamOut_1.CivStage,UvData.FileType{2});
         index_menu=strcmp(ParamOut_1.VelType,menu);
         set(handles.VelType_1,'Value',1+find(index_menu,1))
@@ -2219 +2323 @@
     end
     if numel(UvData.XmlData)==2
-        [tild,tild,tild,num_i1,num_i2,num_j1,num_j2]=fileparts_uvmat(['xx' get(handles.FileIndex_1,'String') get(handles.FileExt_1,'String')]);
-        if isempty(num_i2)
-            num_i2=num_i1;
-        end
-        if isempty(num_j1)
-            num_j1=1;
-        end
-        if isempty(num_j2)
-            num_j2=num_j1;
+        if isempty(i2_1)
+            i2_1=num_i1;
+        end
+        if isempty(j1_1)
+            j1_1=1;
+        end
+        if isempty(j2_1)
+            j2_1=j1_1;
         end
         siz=size(UvData.XmlData{2}.Time);
-        if ~isempty(num_i1) && siz(1)>=max(num_i1,num_i2) && siz(2)>=max(num_j1,num_j2)
-            abstime_1=(UvData.XmlData{2}.Time(num_i1,num_j1)+UvData.XmlData{2}.Time(num_i2,num_j2))/2;%overset the time read from files
-            Field{2}.Dt=(UvData.XmlData{2}.Time(num_i2,num_j2)-UvData.XmlData{2}.Time(num_i1,num_j1));
+        if ~isempty(i1_1) && siz(1)>=max(i1_1,i2_1) && siz(2)>=max(j1_1,j2_1)
+            abstime_1=(UvData.XmlData{2}.Time(i1_1,j1_1)+UvData.XmlData{2}.Time(i2_1,j2_1))/2;%overset the time read from files
+            Field{2}.Dt=(UvData.XmlData{2}.Time(i2_1,j2_1)-UvData.XmlData{2}.Time(i1_1,j1_1));
         end
     end
@@ -2264 +2367 @@
 %% apply coordinate transform or other user fct
 transform=get(handles.path_transform,'UserData');
-if ~isempty(transform)
+if isempty(transform)
+    UvData.Field=Field{1};
+else
     XmlData=[];%default
     XmlData_1=[];%default
@@ -2294 +2399 @@
 end
 
-%% check whether tps is needed, then calculate tps coefficients if needed
-check_proj_tps=0;
-if  (strcmp(UvData.FileType{1},'civdata')||strcmp(UvData.FileType{1},'civx'))
-    for iobj=1:numel(UvData.Object)
-        if isfield(UvData.Object{iobj},'ProjMode')&& strcmp(UvData.Object{iobj}.ProjMode,'filter')
-            check_proj_tps=1;
-            break
-        end
-    end
-end
-check_tps=0;         
-if strcmp(UvData.FileType{1},'civdata')&&~strcmp(ParamOut.FieldName,'velocity')&&~strcmp(ParamOut.FieldName,'get_field...')
-       check_tps=1;%tps needed to get the requested field
-end
-if (check_tps ||check_proj_tps)&&~isfield(UvData.Field,'Coord_tps')
-    UvData.Field=calc_tps(UvData.Field);
-end
-UvData.Field.FieldList=[{ParamOut.FieldName} {ParamOut.ColorVar}];
-
-%% calculate scalar
-if isstruct(ParamOut)&&~strcmp(ParamOut.FieldName,'get_field...')&& (strcmp(UvData.FileType{1},'civdata')||strcmp(UvData.FileType{1},'civx'))...
-         &&~strcmp(ParamOut.FieldName,'velocity') && ~strcmp(ParamOut.FieldName,'get_field...') 
-    if ~check_proj_tps
-        UvData.Field=calc_field([{ParamOut.FieldName} {ParamOut.ColorVar}],UvData.Field);
-    end
-end
-if isstruct(ParamOut_1)&& numel(Field)==2 && ~strcmp(ParamOut_1.FieldName,'get_field...')&& ~test_keepdata_1 && (strcmp(UvData.FileType{2},'civdata')||strcmp(UvData.FileType{2},'civx'))...
-        &&~strcmp(ParamOut_1.FieldName,'velocity') && ~strcmp(ParamOut_1.FieldName,'get_field...') 
-    if check_proj_tps
-        Field{2}.FieldList=[{ParamOut_1.FieldName} {ParamOut_1.ColorVar}];
-    else
-     Field{2}=calc_field([{ParamOut_1.FieldName} {ParamOut_1.ColorVar}],Field{2});
-    end
-end
-
+%% calculate tps coefficients if needed
+UvData.Field=calc_tps(UvData.Field,check_proj_tps);
 
 %% analyse input field
@@ -2942 +3014 @@
     end 
     set(handles.uvmat,'UserData',UvData);
-    run0_Callback(hObject, eventdata, handles); %run
+    transform_fct_list=get(handles.transform_fct,'String');
+    transform_fct=transform_fct_list(get(handles.transform_fct,'Value'));
+    if strcmp(transform_fct,'sub_field')
+        set(handles.transform_fct,'Value',1)%suppress the sub_field transform
+        transform_fct_Callback(hObject, eventdata, handles); 
+    else
+        run0_Callback(hObject, eventdata, handles)
+    end  
 else
     MenuBrowse_1_Callback(hObject, eventdata, handles)
@@ -3114 +3193 @@
 FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1];
 [tild,tild,tild,i1,i2,j1,j2]=fileparts_uvmat(get(handles.FileIndex,'String'));
-% set(handles.FileIndex_1,'Visible','on')
-% set(handles.FileExt_1,'Visible','on')
 switch field_1
     case 'get_field...'
@@ -3128 +3205 @@
         set(hhget_field.list_fig,'Value',1)
         set(hhget_field.list_fig,'String',{'uvmat'})
-%         set(handles.transform_fct,'Value',1)% no transform by default
-%         set(handles.path_transform,'String','')
         if check_new
             UvData.FileType{2}=UvData.FileType{1};
             set(handles.FileIndex_1,'String',get(handles.FileIndex,'String'))
-%             set(handles.FileExt_1,'String',get(handles.FileExt,'String'))
               set(handles.uvmat,'UserData',UvData)
         end
@@ -3239 +3313 @@
 
 %------------------------------------------------------------------------
-% --- Executes on button press in VelType_1.
+% --- Executes on choice selection in VelType_1.
 function VelType_1_Callback(hObject, eventdata, handles)
 %------------------------------------------------------------------------
 set(handles.FixVelType,'Value',1)% the velocity type is now imposed by the GUI (not automatic)
 UvData=get(handles.uvmat,'UserData');
-%refresh field with a second FileName=first file name
-set(handles.run0,'BackgroundColor',[1 1 0])%paint the command button in yellow
+set(handles.run0,'BackgroundColor',[1 1 0])%paint run0 button in yellow to indicate its activation
 drawnow   
-InputFile=read_GUI(handles.InputFile);
-[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
-FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
-
+InputFile=read_GUI(handles.InputFile);% read the input file parameters
+[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
+[RootPath_1,SubDir_1,RootFile_1,FileIndex_1,FileExt_1]=read_file_boxes_1(handles);
+FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];% name of the first input file
+
+check_refresh=0;
 if isempty(InputFile.VelType_1)
-        FileName_1='';% we plot the current field without the second field
+        FileName_1='';% we plot the first input field without the second field
         set(handles.SubField,'Value',0)
-        SubField_Callback(hObject, eventdata, handles)
+        SubField_Callback(hObject, eventdata, handles)% activate SubField_Callback and refresh current display, removing the second field
 elseif get(handles.SubField,'Value')% if subfield is already 'on'
-      [RootPath_1,SubDir_1,RootFile_1,FileIndices_1,FileExt_1]=read_file_boxes_1(handles);
-     FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndices_1 FileExt_1];
-else
-     FileName_1=FileName;% we compare two fields in the same file by default
+     FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1];% name of the second input file
+     check_refresh=1;%will refresh the current display
+else% we introduce the same file (with a different field) for the second series
+     FileName_1=FileName;% we compare two fields in the same file 
      UvData.FileType{2}=UvData.FileType{1};
+     UvData.XmlData{2}= UvData.XmlData{1};
      set(handles.SubField,'Value',1)
-end
-if isfield(UvData,'Field_1')
-    UvData=rmfield(UvData,'Field_1');% removes the stored second field if it exists
-end
-UvData.FileName_1='';% desactivate the use of a constant second file
-set(handles.uvmat,'UserData',UvData)
-num_i1=stra2num(get(handles.i1,'String'));
-num_i2=stra2num(get(handles.i2,'String'));
-num_j1=stra2num(get(handles.j1,'String'));
-num_j2=stra2num(get(handles.j2,'String'));
-
-errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2);
-
-if ~isempty(errormsg)
-    msgbox_uvmat('ERROR',errormsg);
-else
-    set(handles.i1,'BackgroundColor',[1 1 1])
-    set(handles.i2,'BackgroundColor',[1 1 1])
-    set(handles.j1,'BackgroundColor',[1 1 1])
-    set(handles.j2,'BackgroundColor',[1 1 1])
-    set(handles.FileIndex,'BackgroundColor',[1 1 1])
-    set(handles.FileIndex_1,'BackgroundColor',[1 1 1])
-end
-set(handles.run0,'BackgroundColor',[1 0 0])
+     transform=get(handles.path_transform,'UserData');
+     if (~isa(transform,'function_handle')||nargin(transform)<3)
+        set(handles.transform_fct,'value',2); % set transform to sub_field if the current fct doe not accept two input fields
+        transform_fct_Callback(hObject, eventdata, handles)% activate transform_fct_Callback and refresh current display
+     else
+         check_refresh=1;
+     end  
+end
+
+% refresh the current display if it has not been done previously
+if check_refresh
+    UvData.FileName_1='';% desactivate the use of a constant second file
+    set(handles.uvmat,'UserData',UvData)
+    num_i1=stra2num(get(handles.i1,'String'));
+    num_i2=stra2num(get(handles.i2,'String'));
+    num_j1=stra2num(get(handles.j1,'String'));
+    num_j2=stra2num(get(handles.j2,'String'));
+    [tild,tild,tild,i1_1,i2_1,j1_1,j2_1]=fileparts_uvmat(['xx' FileIndex_1]);
+    errormsg=refresh_field(handles,FileName,FileName_1,num_i1,num_i2,num_j1,num_j2,i1_1,i2_1,j1_1,j2_1);
+    if ~isempty(errormsg)
+        msgbox_uvmat('ERROR',errormsg);
+    else
+        set(handles.i1,'BackgroundColor',[1 1 1])
+        set(handles.i2,'BackgroundColor',[1 1 1])
+        set(handles.j1,'BackgroundColor',[1 1 1])
+        set(handles.j2,'BackgroundColor',[1 1 1])
+        set(handles.FileIndex,'BackgroundColor',[1 1 1])
+        set(handles.FileIndex_1,'BackgroundColor',[1 1 1])
+    end
+    set(handles.run0,'BackgroundColor',[1 0 0])
+end
 
 %-----------------------------------------------
@@ -3605 +3688 @@
 
 %% refresh the current plot
-run0_Callback(hObject, eventdata, handles)
-
-
+if isempty(list_path{ichoice}) || nargin(transform_handle)<3
+    set(handles.SubField,'Value',0)
+    SubField_Callback(hObject, eventdata, handles)
+else
+    run0_Callback(hObject, eventdata, handles)
+end
 
 %------------------------------------------------
@@ -3961 +4047 @@
 %replot all the objects within the new projected field
 for IndexObj=1:numel(list_str)
-    IndexObj
-        hobject=UvData.Object{IndexObj}.DisplayHandle.uvmat
+        hobject=UvData.Object{IndexObj}.DisplayHandle.uvmat;
         if isempty(hobject) || ~ishandle(hobject)
-            hobject=handles.PlotAxes
+            hobject=handles.PlotAxes;
         end
         if isequal(IndexObj,get(handles.ListObject,'Value'))
@@ -4151 +4236 @@
         data.Type='plane';
     end
-    if isfield(UvData,'Field')
-        Field=UvData.Field;
-        if isfield(UvData.Field,'Mesh')&&~isempty(UvData.Field.Mesh)
-            data.RangeX=[UvData.Field.XMin UvData.Field.XMax];
-            if strcmp(data.Type,'line')||strcmp(data.Type,'polyline')
-                data.RangeY=UvData.Field.Mesh;
-            else
-                data.RangeY=[UvData.Field.YMin UvData.Field.YMax];
-            end
-            data.DX=UvData.Field.Mesh;
-            data.DY=UvData.Field.Mesh;
-        end
-        if isfield(Field,'NbDim')&& isequal(Field.NbDim,3)
-            data.Coord=[0 0 0]; %default
-        end
-        if isfield(Field,'CoordUnit')
-            data.CoordUnit=Field.CoordUnit;
-        end
-    end
+%     if isfield(UvData,'Field')
+%         Field=UvData.Field;
+%         if isfield(UvData.Field,'Mesh')&&~isempty(UvData.Field.Mesh)
+%             data.RangeX=[UvData.Field.XMin UvData.Field.XMax];
+%             if strcmp(data.Type,'line')||strcmp(data.Type,'polyline')
+%                 data.RangeY=UvData.Field.Mesh;
+%             else
+%                 data.RangeY=[UvData.Field.YMin UvData.Field.YMax];
+%             end
+%             data.DX=UvData.Field.Mesh;
+%             data.DY=UvData.Field.Mesh;
+%         end
+%         if isfield(Field,'NbDim')&& isequal(Field.NbDim,3)
+%             data.Coord=[0 0 0]; %default
+%         end
+%         if isfield(Field,'CoordUnit')
+%             data.CoordUnit=Field.CoordUnit;
+%         end
+%     end
     hset_object=set_object(data,[],ZBounds);
     hhset_object=guidata(hset_object);
@@ -4297 +4382 @@
 UvData=get(huvmat,'UserData');
 %[xx,xx,FileBase]=read_file_boxes(handles);
-[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
+[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
 FileBase=fullfile(RootPath,RootFile);
  %read the current input file name
@@ -4587 +4672 @@
 pos(1)=pos(1)+pos(3)-0.311+0.04; %0.311= width of the geometry_calib interface (units relative to the srcreen)
 pos(2)=pos(2)-0.02;
-[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
-FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
+[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
+FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
 set(handles.view_xml,'Backgroundcolor',[1 1 0])%indicate the reading of the current xml file by geometry_calib
 if isfield(UvData.OpenParam,'CalOrigin')
@@ -4720 +4805 @@
 
 %prepare display of the set_grid GUI
-[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
-FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
+[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
+FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
 CoordList=get(handles.transform_fct,'String');
 val=get(handles.transform_fct,'Value');
@@ -4746 +4831 @@
 %------------------------------------------------------------------------
 series; %first display of the GUI to fill waiting time
-[RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
-param.FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
+[RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
+param.FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
 if isequal(get(handles.SubField,'Value'),1)
-    [RootPath_1,SubDir_1,RootFile_1,FileIndices_1,FileExt_1]=read_file_boxes_1(handles);
-    FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndices_1 FileExt_1];
+    [RootPath_1,SubDir_1,RootFile_1,FileIndex_1,FileExt_1]=read_file_boxes_1(handles);
+    FileName_1=[fullfile(RootPath_1,SubDir_1,RootFile_1) FileIndex_1 FileExt_1];
     if ~isequal(FileName_1,param.FileName)
         param.FileName_1=FileName_1;
@@ -4788 +4873 @@
 function MenuPIV_Callback(hObject, eventdata, handles)
 %------------------------------------------------------------------------
- [RootPath,SubDir,RootFile,FileIndices,FileExt]=read_file_boxes(handles);
- FileName=[fullfile(RootPath,SubDir,RootFile) FileIndices FileExt];
+ [RootPath,SubDir,RootFile,FileIndex,FileExt]=read_file_boxes(handles);
+ FileName=[fullfile(RootPath,SubDir,RootFile) FileIndex FileExt];
 civ(FileName);% interface de civ(not in the uvmat file)