Index: trunk/src/proj_field.m
===================================================================
--- trunk/src/proj_field.m	(revision 630)
+++ trunk/src/proj_field.m	(revision 644)
@@ -564,5 +564,5 @@
    VarIndex=find(check_proj);
 
-    %identify vector components   
+    %% identify vector components   
     testU=isfield(CellInfo{icell},'VarIndex_vector_x') &&isfield(CellInfo{icell},'VarIndex_vector_y') ;% test for vectors
     if testU
@@ -578,16 +578,17 @@
         errorflag=FieldData.(FFName);
     end   
-    % check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y
+    
+    %% check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y
     if strcmp(CellInfo{icell}.CoordType,'scattered')
-        if  ~isequal(ProjMode,'interp_lin')
+        if  strcmp(ProjMode,'projection') 
             if width==0
                 errormsg='range of the projection object is missing';
-                return      
-            else
-                lambda=2/(width*width); %smoothing factor used for interp_tps: weight exp(-2) at distance width from the line
-            end
-        end
-        if ~isequal(ProjMode,'projection')
-            if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
+                return  
+            end
+%             else
+%                 lambda=2/(width*width); %smoothing factor used for interp_tps: weight exp(-2) at distance width from the line
+%             end
+        else
+            if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX)
                 DX=abs(ObjectData.DX);%mesh of interpolation points along the line
             else
@@ -601,5 +602,6 @@
         coord_y=FieldData.(YName);
     end   
-    %initiate projection
+    
+    %% initiate projection
     for ivar=1:length(VarIndex)
         ProjLine{ivar}=[];
@@ -612,234 +614,282 @@
   %%%%%%%  % A FAIRE CALCULER MEAN DES QUANTITES    %%%%%%
    %case of unstructured coordinates
-    if strcmp(CellInfo{icell}.CoordType,'scattered')  
-        for ip=1:siz_line(1)-1     %Loop on the segments of the polyline
-            linelength=sqrt(dlinx(ip)*dlinx(ip)+dliny(ip)*dliny(ip));  
-            %select the vector indices in the range of action
-            if testfalse
-                flagsel=(errorflag==0); % keep only non false vectors
-            else
-                flagsel=ones(size(coord_x));
-            end
-            if isequal(ProjMode,'projection') %|| isequal(ProjMode,'interp_tps')
-                flagsel=flagsel & ((coord_y -yinf(ip))*(xinf(ip+1)-xinf(ip))>(coord_x-xinf(ip))*(yinf(ip+1)-yinf(ip))) ...
-                & ((coord_y -ysup(ip))*(xsup(ip+1)-xsup(ip))<(coord_x-xsup(ip))*(ysup(ip+1)-ysup(ip))) ...
-                & ((coord_y -yinf(ip+1))*(xsup(ip+1)-xinf(ip+1))>(coord_x-xinf(ip+1))*(ysup(ip+1)-yinf(ip+1))) ...
-                & ((coord_y -yinf(ip))*(xsup(ip)-xinf(ip))<(coord_x-xinf(ip))*(ysup(ip)-yinf(ip)));
-            end
-            indsel=find(flagsel);%indsel =indices of good vectors 
-            X_sel=coord_x(indsel);
-            Y_sel=coord_y(indsel);
-            nbvar=0;
-            for iselect=1:numel(VarIndex)-2*testU
-                VarName=FieldData.ListVarName{VarIndex(iselect)};
-                ProjVar{iselect}=FieldData.(VarName)(indsel);%scalar value
-            end   
-            if testU
-                ProjVar{numel(VarIndex)-1}=cos(theta(ip))*vector_x(indsel)+sin(theta(ip))*vector_y(indsel);% longitudinal component
-                ProjVar{numel(VarIndex)}=-sin(theta(ip))*vector_x(indsel)+cos(theta(ip))*vector_y(indsel);%transverse component         
-            end
-            if isequal(ProjMode,'projection')
-                sintheta=sin(theta(ip));
-                costheta=cos(theta(ip));
-                Xproj=(X_sel-ObjectData.Coord(ip,1))*costheta + (Y_sel-ObjectData.Coord(ip,2))*sintheta; %projection on the line
-                [Xproj,indsort]=sort(Xproj);
-                for ivar=1:numel(ProjVar)
-                    if ~isempty(ProjVar{ivar})
-                        ProjVar{ivar}=ProjVar{ivar}(indsort);
-                     end
-                end
-            elseif isequal(ProjMode,'interp_lin') %linear interpolation:
-                npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX)
-                Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
-                xreg=cos(theta(ip))*Xproj+ObjectData.Coord(ip,1);
-                yreg=sin(theta(ip))*Xproj+ObjectData.Coord(ip,2);
-                for ivar=1:numel(ProjVar)
-                     if ~isempty(ProjVar{ivar})
-                        ProjVar{ivar}=griddata_uvmat(X_sel,Y_sel,ProjVar{ivar},xreg,yreg);
-                     end
-                end
-            elseif isequal(ProjMode,'interp_tps') %filtering
-                npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX)
-                Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
-                siz=size(X_sel);
-                xregij=cos(theta(ip))*ones(siz(1),1)*Xproj+ObjectData.Coord(ip,1);
-                yregij=sin(theta(ip))*ones(siz(1),1)*Xproj+ObjectData.Coord(ip,2);
-                xij=X_sel*ones(1,npoint);
-                yij=Y_sel*ones(1,npoint);
-                Aij=exp(-lambda*((xij-xregij).*(xij-xregij)+(yij-yregij).*(yij-yregij)));
-                norm=Aij'*ones(siz(1),1);
-                for ivar=1:numel(ProjVar)
-                     if ~isempty(ProjVar{ivar})
-                        ProjVar{ivar}=Aij'*ProjVar{ivar}./norm;
-                     end
-                end              
-            end
-            %prolongate the total record
-            for ivar=1:numel(ProjVar)
-                  if ~isempty(ProjVar{ivar})
-                     ProjLine{ivar}=[ProjLine{ivar}; ProjVar{ivar}];
-                  end
-            end
-            XLine=[XLine ;(Xproj+linelengthtot)];%along line abscissa
-            linelengthtot=linelengthtot+linelength;
-            %     circul=circul+(sum(U_sel))*linelength/npoint;
-            %     flux=flux+(sum(V_sel))*linelength/npoint;
-        end
-        ProjData.X=XLine';
-        ProjData.ListVarName=[ProjData.ListVarName {XName}];
-        ProjData.VarDimName=[ProjData.VarDimName {XName}];
-        ProjData.VarAttribute{1}.long_name='abscissa along line';
-        for iselect=1:numel(VarIndex)
-            VarName=FieldData.ListVarName{VarIndex(iselect)};
-            eval(['ProjData.' VarName '=ProjLine{iselect};'])
-            ProjData.ListVarName=[ProjData.ListVarName {VarName}];
-            ProjData.VarDimName=[ProjData.VarDimName {XName}];
-            ProjData.VarAttribute{iselect}=FieldData.VarAttribute{VarIndex(iselect)};
-            if strcmp(ProjMode,'projection')
-                ProjData.VarAttribute{iselect}.Role='discrete';
-            else
-                 ProjData.VarAttribute{iselect}.Role='continuous';
-            end
-        end
-    
-    %case of structured coordinates
-    elseif strcmp(CellInfo{icell}.CoordType,'grid')
-        if ~isequal(ObjectData.Type,'line')% exclude polyline
-            errormsg=['no  projection available on ' ObjectData.Type 'for structured coordinates']; % 
-        else
-            test_Amat=1;%image or 2D matrix
-            test_interp2=0;%default 
-            AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
-            AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
-            eval(['AX=FieldData.' AXName ';']);% set of x positions
-            eval(['AY=FieldData.' AYName ';']);% set of y positions  
-            AName=FieldData.ListVarName{VarIndex(1)};
-            eval(['A=FieldData.' AName ';']);% scalar
-            npxy=size(A);
-            npx=npxy(2);
-            npy=npxy(1); 
-            if numel(AX)==2
-                DX=(AX(2)-AX(1))/(npx-1);
-            else
-                DX_vec=diff(AX);
-                DX=max(DX_vec);
-                DX_min=min(DX_vec);
-                if (DX-DX_min)>0.0001*abs(DX) 
-                    test_interp2=1;
-                    DX=DX_min;
-                end    
-            end
-            if numel(AY)==2
-                DY=(AY(2)-AY(1))/(npy-1);
-            else
-                DY_vec=diff(AY);
-                DY=max(DY_vec);
-                DY_min=min(DY_vec);
-                if (DY-DY_min)>0.0001*abs(DY)
+   if strcmp(CellInfo{icell}.CoordType,'scattered')
+       for ip=1:siz_line(1)-1     %Loop on the segments of the polyline
+           linelength=sqrt(dlinx(ip)*dlinx(ip)+dliny(ip)*dliny(ip));
+           %select the vector indices in the range of action
+           if testfalse
+               flagsel=(errorflag==0); % keep only non false vectors
+           else
+               flagsel=ones(size(coord_x));
+           end
+           if isequal(ProjMode,'projection') %|| isequal(ProjMode,'interp_tps')
+               flagsel=flagsel & ((coord_y -yinf(ip))*(xinf(ip+1)-xinf(ip))>(coord_x-xinf(ip))*(yinf(ip+1)-yinf(ip))) ...
+                   & ((coord_y -ysup(ip))*(xsup(ip+1)-xsup(ip))<(coord_x-xsup(ip))*(ysup(ip+1)-ysup(ip))) ...
+                   & ((coord_y -yinf(ip+1))*(xsup(ip+1)-xinf(ip+1))>(coord_x-xinf(ip+1))*(ysup(ip+1)-yinf(ip+1))) ...
+                   & ((coord_y -yinf(ip))*(xsup(ip)-xinf(ip))<(coord_x-xinf(ip))*(ysup(ip)-yinf(ip)));
+           end
+           indsel=find(flagsel);%indsel =indices of good vectors
+           X_sel=coord_x(indsel);
+           Y_sel=coord_y(indsel);
+           nbvar=0;
+           for iselect=1:numel(VarIndex)-2*testU
+               VarName=FieldData.ListVarName{VarIndex(iselect)};
+               ProjVar{iselect}=FieldData.(VarName)(indsel);%scalar value
+           end
+           if testU
+               ProjVar{numel(VarIndex)-1}=cos(theta(ip))*vector_x(indsel)+sin(theta(ip))*vector_y(indsel);% longitudinal component
+               ProjVar{numel(VarIndex)}=-sin(theta(ip))*vector_x(indsel)+cos(theta(ip))*vector_y(indsel);%transverse component
+           end
+           if isequal(ProjMode,'projection')
+               sintheta=sin(theta(ip));
+               costheta=cos(theta(ip));
+               Xproj=(X_sel-ObjectData.Coord(ip,1))*costheta + (Y_sel-ObjectData.Coord(ip,2))*sintheta; %projection on the line
+               [Xproj,indsort]=sort(Xproj);
+               for ivar=1:numel(ProjVar)
+                   if ~isempty(ProjVar{ivar})
+                       ProjVar{ivar}=ProjVar{ivar}(indsort);
+                   end
+               end
+           elseif isequal(ProjMode,'interp_lin')||isequal(ProjMode,'interp_tps') %filtering %linear interpolation:
+               npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX)
+               Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
+               xreg=cos(theta(ip))*Xproj+ObjectData.Coord(ip,1);
+               yreg=sin(theta(ip))*Xproj+ObjectData.Coord(ip,2);
+               if isfield(CellInfo{icell},'VarIndex_vector_x')&&isfield(CellInfo{icell},'VarIndex_vector_y')
+                   VarName_x=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x};
+                   VarName_y=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y};
+                   if isfield(CellInfo{icell},'VarIndex_errorflag')
+                       FieldData.(VarName_x)=FieldData.(VarName_x)(indsel);
+                       FieldData.(VarName_y)=FieldData.(VarName_y)(indsel);
+                   end
+                   if ~isfield(CellInfo{icell},'CheckSub') || ~CellInfo{icell}.CheckSub
+                       vector_x_proj=numel(ProjData.ListVarName)+1;
+                       vector_y_proj=numel(ProjData.ListVarName)+2;
+                   end
+               end
+               if isfield(CellInfo{icell},'VarIndex_scalar')
+                   VarName_scalar=FieldData.ListVarName{CellInfo{icell}.VarIndex_scalar};
+                   if isfield(CellInfo{icell},'errorflag') && ~isempty(CellInfo{icell}.errorflag)
+                       FieldData.(VarName_scalar)=FieldData.(VarName_scalar)(indsel);
+                   end
+               end
+               if isfield(CellInfo{icell},'VarIndex_ancillary')% do not project ancillary data with interp
+                   FieldData=rmfield(FieldData,FieldData.ListVarName{CellInfo{icell}.VarIndex_ancillary});
+               end
+               if isfield(CellInfo{icell},'VarIndex_warnflag')% do not project ancillary data with interp
+                   FieldData=rmfield(FieldData,FieldData.ListVarName{CellInfo{icell}.VarIndex_warnflag});
+               end
+               if isfield(CellInfo{icell},'VarIndex_errorflag')% do not project ancillary data with interp
+                   FieldData=rmfield(FieldData,FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag});
+               end
+               if isequal(ProjMode,'interp_lin')
+               [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([X_sel Y_sel],FieldData,CellInfo{icell}.FieldName,xreg',yreg');
+               else
+                  [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_tps([X_sel Y_sel],FieldData,CellInfo{icell}.FieldName,xreg',yreg');
+               end
+               ivar_vector_x=[];
+               ivar_vector_y=[];
+               for ivar=1:numel(VarAttribute)
+                   if isfield(VarAttribute{ivar},'Role')
+                       if strcmp(VarAttribute{ivar}.Role,'vector_x')
+                       ivar_vector_x=ivar;
+                       elseif strcmp(VarAttribute{ivar}.Role,'vector_y')
+                   ivar_vector_y=ivar;
+                       end
+                   end
+               end
+               if ~isempty(ivar_vector_x)&&~isempty(ivar_vector_y)
+                                  ProjVar{ivar_vector_x}=cos(theta(ip))*ProjVar{ivar_vector_x}+sin(theta(ip))*ProjVar{ivar_vector_y};% longitudinal component
+               ProjVar{ivar_vector_y}=-sin(theta(ip))*ProjVar{ivar_vector_x}+cos(theta(ip))*ProjVar{ivar_vector_y};%transverse component
+               end
+           elseif isequal(ProjMode,'interp_tps') %filtering
+               %   TODO 
+               %                 npoint=floor(linelength/DX)+1;% nbre of points in the profile (interval DX)
+               %                 Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
+               %                 siz=size(X_sel);
+               %                 xregij=cos(theta(ip))*ones(siz(1),1)*Xproj+ObjectData.Coord(ip,1);
+               %                 yregij=sin(theta(ip))*ones(siz(1),1)*Xproj+ObjectData.Coord(ip,2);
+               %                 xij=X_sel*ones(1,npoint);
+               %                 yij=Y_sel*ones(1,npoint);
+               %                 Aij=exp(-lambda*((xij-xregij).*(xij-xregij)+(yij-yregij).*(yij-yregij)));
+               %                 norm=Aij'*ones(siz(1),1);
+               %                 for ivar=1:numel(ProjVar)
+               %                      if ~isempty(ProjVar{ivar})
+               %                         ProjVar{ivar}=Aij'*ProjVar{ivar}./norm;
+               %
+               %                      end
+               %                 end
+           end
+           %prolongate the total record
+           for ivar=1:numel(ProjVar)
+               if ~isempty(ProjVar{ivar})
+                   if numel(ProjLine)>=ivar
+                       ProjLine{ivar}=[ProjLine{ivar}; ProjVar{ivar}];
+                   else
+                       ProjLine{ivar}=ProjVar{ivar};
+                   end
+               end
+           end
+           XLine=[XLine ;(Xproj+linelengthtot)];%along line abscissa
+           linelengthtot=linelengthtot+linelength;
+           %     circul=circul+(sum(U_sel))*linelength/npoint;
+           %     flux=flux+(sum(V_sel))*linelength/npoint;
+       end
+       ProjData.X=XLine';
+       ProjData.ListVarName=[ProjData.ListVarName {XName}];
+       ProjData.VarDimName=[ProjData.VarDimName {XName}];
+       ProjData.VarAttribute{1}.long_name='abscissa along line';
+       for iselect=1:numel(VarIndex)
+           VarName=FieldData.ListVarName{VarIndex(iselect)};
+           eval(['ProjData.' VarName '=ProjLine{iselect};'])
+           ProjData.ListVarName=[ProjData.ListVarName {VarName}];
+           ProjData.VarDimName=[ProjData.VarDimName {XName}];
+           ProjData.VarAttribute{iselect}=FieldData.VarAttribute{VarIndex(iselect)};
+           if strcmp(ProjMode,'projection')
+               ProjData.VarAttribute{iselect}.Role='discrete';
+           else
+               ProjData.VarAttribute{iselect}.Role='continuous';
+           end
+       end
+       
+       %case of structured coordinates
+   elseif strcmp(CellInfo{icell}.CoordType,'grid')
+       if ~isequal(ObjectData.Type,'line')% exclude polyline
+           errormsg=['no  projection available on ' ObjectData.Type 'for structured coordinates']; %
+       else
+           test_Amat=1;%image or 2D matrix
+           test_interp2=0;%default
+           AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
+           AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
+           eval(['AX=FieldData.' AXName ';']);% set of x positions
+           eval(['AY=FieldData.' AYName ';']);% set of y positions
+           AName=FieldData.ListVarName{VarIndex(1)};
+           eval(['A=FieldData.' AName ';']);% scalar
+           npxy=size(A);
+           npx=npxy(2);
+           npy=npxy(1);
+           if numel(AX)==2
+               DX=(AX(2)-AX(1))/(npx-1);
+           else
+               DX_vec=diff(AX);
+               DX=max(DX_vec);
+               DX_min=min(DX_vec);
+               if (DX-DX_min)>0.0001*abs(DX)
                    test_interp2=1;
-                    DY=DY_min;
-                end     
-            end              
-            AXI=linspace(AX(1),AX(end), npx);%set of  x  positions for the interpolated input data
-            AYI=linspace(AY(1),AY(end), npy);%set of  x  positions for the interpolated input data
-            if isfield(ObjectData,'DX')
-                DXY_line=ObjectData.DX;%mesh on the projection line
-            else
-                DXY_line=sqrt(abs(DX*DY));% mesh on the projection line
-            end
-            dlinx=ObjectData.Coord(2,1)-ObjectData.Coord(1,1);
-            dliny=ObjectData.Coord(2,2)-ObjectData.Coord(1,2);
-            linelength=sqrt(dlinx*dlinx+dliny*dliny);
-            theta=angle(dlinx+i*dliny);%angle of the line   
-            if isfield(FieldData,'RangeX')
-                XMin=min(FieldData.RangeX);%shift of the origin on the line
-            else
-                XMin=0;
-            end
-            eval(['ProjData.' AXName '=linspace(XMin,XMin+linelength,linelength/DXY_line+1);'])%abscissa of the new pixels along the line
-            y=linspace(-width,width,2*width/DXY_line+1);%ordintes of the new pixels (coordinate across the line)
-            eval(['npX=length(ProjData.' AXName ');'])
-            npY=length(y); %TODO: utiliser proj_grid
-            eval(['[X,Y]=meshgrid(ProjData.' AXName ',y);'])%grid in the line coordinates
-            XIMA=ObjectData.Coord(1,1)+(X-XMin)*cos(theta)-Y*sin(theta);
-            YIMA=ObjectData.Coord(1,2)+(X-XMin)*sin(theta)+Y*cos(theta);
-            XIMA=(XIMA-AX(1))/DX+1;%  index of the original image along x
-            YIMA=(YIMA-AY(1))/DY+1;% index of the original image along y
-            XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
-            YIMA=reshape(round(YIMA),1,npX*npY);
-            flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
-            ind_in=find(flagin);
-            ind_out=find(~flagin);
-            ICOMB=(XIMA-1)*npy+YIMA;
-            ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
-            nbcolor=1; %color images
-            if numel(npxy)==2
-                nbcolor=1;
-            elseif length(npxy)==3
-                nbcolor=npxy(3);
-            else
-                errormsg='multicomponent field not projected';
-                display(errormsg)
-                return
-            end 
-            nbvar=length(ProjData.ListVarName);% number of var from previous cells
-            ProjData.ListVarName=[ProjData.ListVarName {AXName}];
-            ProjData.VarDimName=[ProjData.VarDimName {AXName}];
-            for ivar=VarIndex
-                %VarName{ivar}=FieldData.ListVarName{ivar};
-                if test_interp2% interpolate on new grid
-                    FieldData.(FieldData.ListVarName{ivar})=interp2(FieldData.(AXName),FieldData.(AYName),FieldData.(FieldData.ListVarName{ivar}),AXI,AYI);%TO TEST
-                end
-                vec_A=reshape(squeeze(FieldData.(FieldData.ListVarName{ivar})),npx*npy,nbcolor); %put the original image in colum
-                if nbcolor==1
-                    vec_B(ind_in)=vec_A(ICOMB);
-                    vec_B(ind_out)=zeros(size(ind_out));
-                    A_out=reshape(vec_B,npY,npX);
-                    ProjData.(FieldData.ListVarName{ivar}) =sum(A_out,1)/npY;
-                elseif nbcolor==3
-                    vec_B(ind_in,1:3)=vec_A(ICOMB,:);
-                    vec_B(ind_out,1)=zeros(size(ind_out));
-                    vec_B(ind_out,2)=zeros(size(ind_out));
-                    vec_B(ind_out,3)=zeros(size(ind_out));
-                    A_out=reshape(vec_B,npY,npX,nbcolor);
-                    ProjData.(FieldData.ListVarName{ivar})=squeeze(sum(A_out,1)/npY);
-                end  
-                ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName{ivar}];
-                ProjData.VarDimName=[ProjData.VarDimName {AXName}];%to generalize with the initial name of the x coordinate
-                ProjData.VarAttribute{ivar}.Role='continuous';% for plot with continuous line
-            end
-            if testU
-                 vector_x =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x});
-                 vector_y =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y});
-                 ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x}) =cos(theta)*vector_x+sin(theta)*vector_y;
-                 ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y}) =-sin(theta)*vector_x+cos(theta)*vector_y;
-            end
-            ProjData.VarAttribute{nbvar+1}.long_name='abscissa along line';
-            if nbcolor==3
-                ProjData.VarDimName{end}={AXName,'rgb'};
-            end
-        end
-    elseif strcmp(CellInfo{icell}.CoordType,'tps')
-         if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
-                DX=abs(ObjectData.DX);%mesh of interpolation points along the line
-                Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
-                xreg=cos(theta(ip))*Xproj+ObjectData.Coord(ip,1)
-                yreg=sin(theta(ip))*Xproj+ObjectData.Coord(ip,2)
-%                 coord_x_proj=XMin:DX:XMax;
-%                 coord_y_proj=YMin:DY:YMax;
-                DataOut=calc_field_tps(FieldData.FieldList,FieldData,cat(3,xreg,yreg));
-                ProjData.ListVarName=[ProjData.ListVarName DataOut.ListVarName];
-                ProjData.VarDimName=[ProjData.VarDimName DataOut.VarDimName];
-                ProjData.VarAttribute=[ProjData.VarAttribute DataOut.VarAttribute];   
-                DataOut.ListVarName(1)=[];
-                DataOut.VarDimName(1)=[];
-                DataOut.VarAttribute(1)=[];
-                for ilist=2:length(DataOut.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
-                    VarName=DataOut.ListVarName{ilist};
-                     ProjData.(VarName)=DataOut.(VarName);
-                end
-                ProjData.coord_x=Xproj;
-         end
-    end
+                   DX=DX_min;
+               end
+           end
+           if numel(AY)==2
+               DY=(AY(2)-AY(1))/(npy-1);
+           else
+               DY_vec=diff(AY);
+               DY=max(DY_vec);
+               DY_min=min(DY_vec);
+               if (DY-DY_min)>0.0001*abs(DY)
+                   test_interp2=1;
+                   DY=DY_min;
+               end
+           end
+           AXI=linspace(AX(1),AX(end), npx);%set of  x  positions for the interpolated input data
+           AYI=linspace(AY(1),AY(end), npy);%set of  x  positions for the interpolated input data
+           if isfield(ObjectData,'DX')
+               DXY_line=ObjectData.DX;%mesh on the projection line
+           else
+               DXY_line=sqrt(abs(DX*DY));% mesh on the projection line
+           end
+           dlinx=ObjectData.Coord(2,1)-ObjectData.Coord(1,1);
+           dliny=ObjectData.Coord(2,2)-ObjectData.Coord(1,2);
+           linelength=sqrt(dlinx*dlinx+dliny*dliny);
+           theta=angle(dlinx+i*dliny);%angle of the line
+           if isfield(FieldData,'RangeX')
+               XMin=min(FieldData.RangeX);%shift of the origin on the line
+           else
+               XMin=0;
+           end
+           eval(['ProjData.' AXName '=linspace(XMin,XMin+linelength,linelength/DXY_line+1);'])%abscissa of the new pixels along the line
+           y=linspace(-width,width,2*width/DXY_line+1);%ordintes of the new pixels (coordinate across the line)
+           eval(['npX=length(ProjData.' AXName ');'])
+           npY=length(y); %TODO: utiliser proj_grid
+           eval(['[X,Y]=meshgrid(ProjData.' AXName ',y);'])%grid in the line coordinates
+           XIMA=ObjectData.Coord(1,1)+(X-XMin)*cos(theta)-Y*sin(theta);
+           YIMA=ObjectData.Coord(1,2)+(X-XMin)*sin(theta)+Y*cos(theta);
+           XIMA=(XIMA-AX(1))/DX+1;%  index of the original image along x
+           YIMA=(YIMA-AY(1))/DY+1;% index of the original image along y
+           XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
+           YIMA=reshape(round(YIMA),1,npX*npY);
+           flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
+           ind_in=find(flagin);
+           ind_out=find(~flagin);
+           ICOMB=(XIMA-1)*npy+YIMA;
+           ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
+           nbcolor=1; %color images
+           if numel(npxy)==2
+               nbcolor=1;
+           elseif length(npxy)==3
+               nbcolor=npxy(3);
+           else
+               errormsg='multicomponent field not projected';
+               display(errormsg)
+               return
+           end
+           nbvar=length(ProjData.ListVarName);% number of var from previous cells
+           ProjData.ListVarName=[ProjData.ListVarName {AXName}];
+           ProjData.VarDimName=[ProjData.VarDimName {AXName}];
+           for ivar=VarIndex
+               %VarName{ivar}=FieldData.ListVarName{ivar};
+               if test_interp2% interpolate on new grid
+                   FieldData.(FieldData.ListVarName{ivar})=interp2(FieldData.(AXName),FieldData.(AYName),FieldData.(FieldData.ListVarName{ivar}),AXI,AYI);%TO TEST
+               end
+               vec_A=reshape(squeeze(FieldData.(FieldData.ListVarName{ivar})),npx*npy,nbcolor); %put the original image in colum
+               if nbcolor==1
+                   vec_B(ind_in)=vec_A(ICOMB);
+                   vec_B(ind_out)=zeros(size(ind_out));
+                   A_out=reshape(vec_B,npY,npX);
+                   ProjData.(FieldData.ListVarName{ivar}) =sum(A_out,1)/npY;
+               elseif nbcolor==3
+                   vec_B(ind_in,1:3)=vec_A(ICOMB,:);
+                   vec_B(ind_out,1)=zeros(size(ind_out));
+                   vec_B(ind_out,2)=zeros(size(ind_out));
+                   vec_B(ind_out,3)=zeros(size(ind_out));
+                   A_out=reshape(vec_B,npY,npX,nbcolor);
+                   ProjData.(FieldData.ListVarName{ivar})=squeeze(sum(A_out,1)/npY);
+               end
+               ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName{ivar}];
+               ProjData.VarDimName=[ProjData.VarDimName {AXName}];%to generalize with the initial name of the x coordinate
+               ProjData.VarAttribute{ivar}.Role='continuous';% for plot with continuous line
+           end
+           if testU
+               vector_x =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x});
+               vector_y =ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y});
+               ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x}) =cos(theta)*vector_x+sin(theta)*vector_y;
+               ProjData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y}) =-sin(theta)*vector_x+cos(theta)*vector_y;
+           end
+           ProjData.VarAttribute{nbvar+1}.long_name='abscissa along line';
+           if nbcolor==3
+               ProjData.VarDimName{end}={AXName,'rgb'};
+           end
+       end
+   elseif strcmp(CellInfo{icell}.CoordType,'tps')
+       if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
+           DX=abs(ObjectData.DX);%mesh of interpolation points along the line
+           Xproj=linelength/(2*npoint):linelength/npoint:linelength-linelength/(2*npoint);
+           xreg=cos(theta(ip))*Xproj+ObjectData.Coord(ip,1)
+           yreg=sin(theta(ip))*Xproj+ObjectData.Coord(ip,2)
+           %                 coord_x_proj=XMin:DX:XMax;
+           %                 coord_y_proj=YMin:DY:YMax;
+           DataOut=calc_field_tps(FieldData.FieldList,FieldData,cat(3,xreg,yreg));
+           ProjData.ListVarName=[ProjData.ListVarName DataOut.ListVarName];
+           ProjData.VarDimName=[ProjData.VarDimName DataOut.VarDimName];
+           ProjData.VarAttribute=[ProjData.VarAttribute DataOut.VarAttribute];
+           DataOut.ListVarName(1)=[];
+           DataOut.VarDimName(1)=[];
+           DataOut.VarAttribute(1)=[];
+           for ilist=2:length(DataOut.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
+               VarName=DataOut.ListVarName{ilist};
+               ProjData.(VarName)=DataOut.(VarName);
+           end
+           ProjData.coord_x=Xproj;
+       end
+   end
 end
 
