Index: trunk/src/proj_field.m
===================================================================
--- trunk/src/proj_field.m	(revision 864)
+++ trunk/src/proj_field.m	(revision 866)
@@ -120,4 +120,39 @@
 end
 
+% %% take the difference of projected input fields if relevant
+% [CellInfo,NbDim,errormsg]=find_field_cells(ProjData);
+% ind_remove=zeros(size(ProjData.ListVarName));
+% ivar=[];
+% ivar_1=[];
+% for icell=1:numel(CellInfo)
+%     if ~isempty(CellInfo{icell})
+%         % if two scalar are in the same cell
+%         if isfield(CellInfo{icell},'VarIndex_scalar') && numel(CellInfo{icell}.VarIndex_scalar)==2 && ProjData.VarAttribute{CellInfo{icell}.VarIndex_scalar(2)}.CheckSub;
+%             ivar=[ivar CellInfo{icell}.VarIndex_scalar(1)];
+%             ivar_1=[ivar_1 CellInfo{icell}.VarIndex_scalar(2)];
+%         end
+%         % if two vector u components are in the same cell
+%         if isfield(CellInfo{icell},'VarIndex_vector_x') && numel(CellInfo{icell}.VarIndex_vector_x)==2 && ProjData.VarAttribute{CellInfo{icell}.VarIndex_vector_x(2)}.CheckSub;
+%             ivar=[ivar CellInfo{icell}.VarIndex_vector_x(1)];
+%             ivar_1=[ivar_1 CellInfo{icell}.VarIndex_vector_x(2)];
+%         end
+%          % if two vector v components are in the same cell
+%         if isfield(CellInfo{icell},'VarIndex_vector_y') && numel(CellInfo{icell}.VarIndex_vector_y)==2 && ProjData.VarAttribute{CellInfo{icell}.VarIndex_vector_y(2)}.CheckSub;
+%             ivar=[ivar CellInfo{icell}.VarIndex_vector_y(1)];
+%             ivar_1=[ivar_1 CellInfo{icell}.VarIndex_vector_y(2)];
+%         end
+%     end
+% end
+% % subtract fields if relevant
+% for imod=1:numel(ivar)
+%         VarName=ProjData.ListVarName{ivar(imod)};
+%         VarName_1=ProjData.ListVarName{ivar_1(imod)};
+%         ProjData.(VarName)=double(ProjData.(VarName))-double(ProjData.(VarName_1));
+%         ind_remove(ivar_1(imod))=1;
+% end
+% ProjData.ListVarName(find(ind_remove))=[];
+% ProjData.VarDimName(find(ind_remove))=[];
+% ProjData.VarAttribute(find(ind_remove))=[];
+
 %-----------------------------------------------------------------
 %project on a set of points
@@ -387,6 +422,6 @@
             XName=FieldData.ListVarName{ivar_X};
             YName=FieldData.ListVarName{ivar_Y};
-            coord_x=FieldData.(XName);
-            coord_y=FieldData.(YName);
+            coord_x=FieldData.(FieldData.ListVarName{ivar_X});
+            coord_y=FieldData.(FieldData.ListVarName{ivar_Y});
             % image or 2D matrix
         case 'grid' %case of structured coordinates
@@ -584,6 +619,6 @@
                 costheta=cos(theta(isegment));
                 sintheta=sin(theta(isegment));
-%                 XIsegment=LineCoord(isegment,1)+DX_edge*costheta:DX*costheta:LineCoord(isegment+1,1));
-%                 YIsegment=(LineCoord(isegment,2)+DX_edge*sintheta:DX*sintheta:LineCoord(isegment+1,2));
+                %                 XIsegment=LineCoord(isegment,1)+DX_edge*costheta:DX*costheta:LineCoord(isegment+1,1));
+                %                 YIsegment=(LineCoord(isegment,2)+DX_edge*sintheta:DX*sintheta:LineCoord(isegment+1,2));
                 NbInterval=floor((dlength(isegment)-DX_edge)/DX);
                 LastX=DX_edge+DX*NbInterval;
@@ -618,4 +653,5 @@
 ProjData.ListVarName={};
 ProjData.VarDimName={};
+check_abscissa=0;
 for icell=1:length(CellInfo)
     % list of variable types to be projected
@@ -631,10 +667,10 @@
     %% identify vector components
     testU=isfield(CellInfo{icell},'VarIndex_vector_x') &&isfield(CellInfo{icell},'VarIndex_vector_y') ;% test for vectors
-%     if testU
-%         UName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x};
-%         VName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y};
-%         vector_x=FieldData.(UName);
-%         vector_y=FieldData.(VName);
-%     end
+    %     if testU
+    %         UName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x};
+    %         VName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y};
+    %         vector_x=FieldData.(UName);
+    %         vector_y=FieldData.(VName);
+    %     end
     %identify error flag
     errorflag=0; %default, no error flag
@@ -654,12 +690,9 @@
         %case of unstructured coordinates
         case 'scattered'
-            XName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
-            YName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
-            coord_x=FieldData.(XName);
-            coord_y=FieldData.(YName);
-            ProjData.ListVarName=[ProjData.ListVarName {XName}];
-            ProjData.VarDimName=[ProjData.VarDimName {XName}];
-            nbvar=numel(ProjData.ListVarName);
-            ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
+%             XName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
+%             YName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
+            coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});
+            coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)});
+            
             if isequal(ProjMode,'projection')
                 if width==0
@@ -667,4 +700,8 @@
                     return
                 end
+                ProjData.ListVarName=[ProjData.ListVarName {FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}}];
+                ProjData.VarDimName=[ProjData.VarDimName {FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}}];
+                nbvar=numel(ProjData.ListVarName);
+                ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
                 % select the (non false) input data located in the band of projection
                 flagsel=(errorflag==0) & ((coord_y -yinf(1))*(xinf(2)-xinf(1))>(coord_x-xinf(1))*(yinf(2)-yinf(1))) ...
@@ -678,10 +715,10 @@
                 Xproj=(coord_x-ObjectData.Coord(1,1))*costheta + (coord_y-ObjectData.Coord(1,2))*sintheta; %projection on the line
                 [Xproj,indsort]=sort(Xproj);% sort points by increasing absissa along the projection line
-                ProjData.(XName)=Xproj;
+                ProjData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)})=Xproj;
                 for ivar=1:numel(VarIndex)
                     ProjData.(VarName{ivar})=FieldData.(VarName{ivar})(flagsel);% restrict vrtibles to the projection band
                     ProjData.(VarName{ivar})=ProjData.(VarName{ivar})(indsort);% sort by absissa
                     ProjData.ListVarName=[ProjData.ListVarName VarName{ivar}];
-                    ProjData.VarDimName=[ProjData.VarDimName {XName}];
+                    ProjData.VarDimName=[ProjData.VarDimName {FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}}];
                     ProjData.VarAttribute{nbvar+ivar}=FieldData.VarAttribute{VarIndex(ivar)};%reproduce var attribute
                     if isfield(ProjData.VarAttribute{nbvar+ivar},'Role')
@@ -695,4 +732,12 @@
                 end
             elseif isequal(ProjMode,'interp_lin')  %filtering %linear interpolation:
+                if ~check_abscissa
+                    XName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
+                    ProjData.ListVarName=[ProjData.ListVarName {XName}];
+                    ProjData.VarDimName=[ProjData.VarDimName {XName}];
+                    nbvar=numel(ProjData.ListVarName);
+                    ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
+                    check_abscissa=1; % define abcissa only once
+                end
                 if ~isequal(errorflag,0)
                     VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};
@@ -704,7 +749,8 @@
                         FieldData.(VarName)=FieldData.(VarName)(indsel);
                     end
-                end                    
+                end
                 [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_x coord_y],FieldData,CellInfo{icell}.FieldName,XI,YI);
                 ProjData.X=Xproj;
+                nbvar=numel(ProjData.ListVarName);
                 ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
                 ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
@@ -713,7 +759,7 @@
                     if isfield(VarAttribute{ivar},'Role')
                         if  strcmp(VarAttribute{ivar}.Role,'vector_x');
-                            ivar_U=ivar;
+                            ivar_U=ivar+nbvar;
                         elseif strcmp(VarAttribute{ivar}.Role,'vector_y');
-                            ivar_V=ivar;
+                            ivar_V=ivar+nbvar;
                         end
                     end
@@ -741,7 +787,7 @@
                     if isfield(VarAttribute{ivar},'Role')
                         if  strcmp(VarAttribute{ivar}.Role,'vector_x');
-                            ivar_U=ivar;
+                            ivar_U=ivar+nbvar;
                         elseif strcmp(VarAttribute{ivar}.Role,'vector_y');
-                            ivar_V=ivar;
+                            ivar_V=ivar+nbvar;
                         end
                     end
@@ -861,11 +907,10 @@
                 end
             end
-    end 
+    end
     if ~isempty(ivar_U) && ~isempty(ivar_V)
-        vector_x =ProjData.(ProjData.ListVarName{nbvar+ivar_U});
-         ProjData.(ProjData.ListVarName{nbvar+ivar_U}) =cos(theta)*vector_x+sin(theta)*ProjData.(ProjData.ListVarName{nbvar+ivar_V});
-          ProjData.(ProjData.ListVarName{nbvar+ivar_V}) =-sin(theta)*vector_x+cos(theta)*ProjData.(ProjData.ListVarName{nbvar+ivar_V});
-    end
-        
+        vector_x =ProjData.(ProjData.ListVarName{ivar_U});
+        ProjData.(ProjData.ListVarName{ivar_U}) =cos(theta)*vector_x+sin(theta)*ProjData.(ProjData.ListVarName{ivar_V});
+        ProjData.(ProjData.ListVarName{ivar_V}) =-sin(theta)*vector_x+cos(theta)*ProjData.(ProjData.ListVarName{ivar_V});
+    end
 end
 
@@ -1230,6 +1275,8 @@
                     
                     if isfield(CellInfo{icell},'CheckSub') && CellInfo{icell}.CheckSub && ~isempty(vector_x_proj)
-                        ProjData.(ListVarName{vector_x_proj})=ProjData.(ListVarName{vector_x_proj})-VarVal{1};
-                        ProjData.(ListVarName{vector_y_proj})=ProjData.(ListVarName{vector_y_proj})-VarVal{2};
+                        ProjData.(FieldData.ListVarName{vector_x_proj})=ProjData.(FieldData.ListVarName{vector_x_proj})-VarVal{1};
+                        ProjData.(FieldData.ListVarName{vector_y_proj})=ProjData.(FieldData.ListVarName{vector_y_proj})-VarVal{2};
+                        ListVarName={};% no new variable
+                        VarAttribute={};
                     else
                         VarDimName=cell(size(ListVarName));
@@ -1243,4 +1290,6 @@
                         end
                     end
+                    vector_x_proj=CellInfo{icell}.VarIndex_vector_x; %preserve for next cell
+                    vector_y_proj=CellInfo{icell}.VarIndex_vector_y; %preserve for next cell
             end
             
@@ -1491,8 +1540,8 @@
                             ListVarName=[ListVarName VarName];
                             VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
-                            eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
+                            ProjData.(VarName)=squeeze(FieldData.(VarName)(iz,:,:));% select the z index iz
                             %TODO : do a vertical average for a thick plane
                             if test_interp(2) || test_interp(3)
-                                eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');'])
+                                ProjData.(VarName)=interp2(Coord{3},Coord{2},ProjData.(VarName),Coord_x,Coord_y);
                             end
                         end
@@ -1535,6 +1584,6 @@
                 if ~isempty(ivar_W)
                     WName=FieldData.ListVarName{ivar_W};
-                    eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])%
-                    eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']);
+                    ProjData.(VName)=ProjData.(VName)+ ProjData.(WName)*sin(Theta);%
+                    ProjData.(WName)=NormVec_X*ProjData.(UName)+ NormVec_Y*ProjData.(VName)+ NormVec_Z* ProjData.(WName);
                 end
             end
