source: trunk/src/proj_field.m @ 869

Last change on this file since 869 was 869, checked in by sommeria, 6 years ago

bug_repaired

File size: 109.2 KB
Line 
1%'proj_field': projects the field on a projection object
2%--------------------------------------------------------------------------
3%  function [ProjData,errormsg]=proj_field(FieldData,ObjectData)
4%
5% OUTPUT:
6% ProjData structure containing the fields of the input field FieldData,
7% transmitted or projected on the object, plus the additional fields
8%    .UMax, .UMin, .VMax, .VMin: min and max of velocity components in a domain
9%    .UMean,VMean: mean of the velocity components in a domain
10%    .AMin, AMax: min and max of a scalar
11%    .AMean: mean of a scalar in a domain 
12%  .NbPix;
13%  .DimName=  names of the matrix dimensions (matlab cell)
14%  .VarName= names of the variables [ProjData.VarName {'A','AMean','AMin','AMax'}];
15%  .VarDimNameIndex= dimensions of the variables, indicated by indices in the list .DimName;
16%
17%INPUT
18% ObjectData: structure characterizing the projection object
19%    .Type : type of projection object
20%    .ProjMode=mode of projection ;
21%    .CoordUnit: 'px', 'cm' units for the coordinates defining the object
22%    .Angle (  angles of rotation (=[0 0 0] by default)
23%    .ProjAngle=angle of projection;
24%    .DX,.DY,.DZ=increments along each coordinate
25%    .Coord(nbpoints,3): set of coordinates defining the object position;
26
27%FieldData: data of the field to be projected on the projection object, with optional fields
28%    .Txt: error message, transmitted to the projection
29%    .FieldList: cell array of strings representing the fields to calculate
30%    .CoordMesh: typical distance between data points (used for mouse action or display), transmitted
31%    .CoordUnit, .TimeUnit, .dt: transmitted
32% standardised description of fields, nc-formated Matlab structure with fields:
33%         .ListGlobalAttribute: cell listing the names of the global attributes
34%        .Att_1,Att_2... : values of the global attributes
35%            .ListVarName: cell listing the names of the variables
36%           .VarAttribute: cell of structures s containing names and values of variable attributes (s.name=value) for each variable of .ListVarName
37%        .Var1, .Var2....: variables (Matlab arrays) with names listed in .ListVarName
38% The variables are grouped in 'fields', made of a set of variables with common dimensions (using the function find_field_cells)
39% The variable attribute 'Role' is used to define the role for plotting:
40%       Role = 'scalar':  (default) represents a scalar field
41%            = 'coord':  represents a set of unstructured coordinates, whose
42%                     space dimension is given by the last array dimension (called 'NbDim').
43%            = 'coord_x', 'coord_y',  'coord_z': represents a separate set of
44%                        unstructured coordinate x, y  or z
45%            = 'vector': represents a vector field whose number of components
46%                is given by the last dimension (called 'NbDim')
47%            = 'vector_x', 'vector_y', 'vector_z'  :represents the x, y or z  component of a vector 
48%            = 'warnflag' : provides a warning flag about the quality of data in a 'Field', default=0, no warning
49%            = 'errorflag': provides an error flag marking false data,
50%                   default=0, no error. Different non zero values can represent different criteria of elimination.
51%
52% Default role of variables (by name)
53%  vector field:
54%    .X,.Y: position of the velocity vectors, projected on the object
55%    .U, .V, .W: velocity components, projected on the object
56%    .C, .CName: scalar associated to the vector
57%    .F : equivalent to 'warnflag'
58%    .FF: equivalent to 'errorflag'
59%  scalar field or image:
60%    .AName: name of a scalar (to be calculated from velocity fields after projection), transmitted
61%    .A: scalar, projected on the object
62%    .AX, .AY: positions for the scalar
63%     case of a structured grid: A is a dim 2 matrix and .AX=[first last] (length 2 vector) represents the first and last abscissa of the grid
64%     case of an unstructured scalar: A is a vector, AX and AY the corresponding coordinates
65
66%=======================================================================
67% Copyright 2008-2014, LEGI UMR 5519 / CNRS UJF G-INP, Grenoble, France
68%   http://www.legi.grenoble-inp.fr
69%   Joel.Sommeria - Joel.Sommeria (A) legi.cnrs.fr
70%
71%     This file is part of the toolbox UVMAT.
72%
73%     UVMAT is free software; you can redistribute it and/or modify
74%     it under the terms of the GNU General Public License as published
75%     by the Free Software Foundation; either version 2 of the license,
76%     or (at your option) any later version.
77%
78%     UVMAT is distributed in the hope that it will be useful,
79%     but WITHOUT ANY WARRANTY; without even the implied warranty of
80%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
81%     GNU General Public License (see LICENSE.txt) for more details.
82%=======================================================================
83
84function [ProjData,errormsg]=proj_field(FieldData,ObjectData)
85errormsg='';%default
86ProjData=[];
87
88%% check input projection object: type, projection mode and Coord:
89if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode')
90    return
91end
92ListProjMode={'projection','interp_lin','interp_tps','inside','outside'};%list of effective projection modes
93if isempty(find(strcmp(ObjectData.ProjMode,ListProjMode), 1))% no projection in case
94    return
95end
96if ~isfield(ObjectData,'Coord')||isempty(ObjectData.Coord)
97    if strcmp(ObjectData.Type,'plane')
98        ObjectData.Coord=[0 0];%default
99    else
100        return
101    end
102end
103
104%% apply projection depending on the object type
105switch ObjectData.Type
106    case 'points'
107        [ProjData,errormsg]=proj_points(FieldData,ObjectData);
108    case {'line','polyline'}
109        [ProjData,errormsg] = proj_line(FieldData,ObjectData);
110    case {'polygon','rectangle','ellipse'}
111        if isequal(ObjectData.ProjMode,'inside')||isequal(ObjectData.ProjMode,'outside')
112            [ProjData,errormsg] = proj_patch(FieldData,ObjectData);
113        else
114            [ProjData,errormsg] = proj_line(FieldData,ObjectData);
115        end
116    case 'plane'
117        [ProjData,errormsg] = proj_plane(FieldData,ObjectData);
118    case 'volume'
119        [ProjData,errormsg] = proj_volume(FieldData,ObjectData);
120end
121
122% %% take the difference of projected input fields if relevant
123% [CellInfo,NbDim,errormsg]=find_field_cells(ProjData);
124% ind_remove=zeros(size(ProjData.ListVarName));
125% ivar=[];
126% ivar_1=[];
127% for icell=1:numel(CellInfo)
128%     if ~isempty(CellInfo{icell})
129%         % if two scalar are in the same cell
130%         if isfield(CellInfo{icell},'VarIndex_scalar') && numel(CellInfo{icell}.VarIndex_scalar)==2 && ProjData.VarAttribute{CellInfo{icell}.VarIndex_scalar(2)}.CheckSub;
131%             ivar=[ivar CellInfo{icell}.VarIndex_scalar(1)];
132%             ivar_1=[ivar_1 CellInfo{icell}.VarIndex_scalar(2)];
133%         end
134%         % if two vector u components are in the same cell
135%         if isfield(CellInfo{icell},'VarIndex_vector_x') && numel(CellInfo{icell}.VarIndex_vector_x)==2 && ProjData.VarAttribute{CellInfo{icell}.VarIndex_vector_x(2)}.CheckSub;
136%             ivar=[ivar CellInfo{icell}.VarIndex_vector_x(1)];
137%             ivar_1=[ivar_1 CellInfo{icell}.VarIndex_vector_x(2)];
138%         end
139%          % if two vector v components are in the same cell
140%         if isfield(CellInfo{icell},'VarIndex_vector_y') && numel(CellInfo{icell}.VarIndex_vector_y)==2 && ProjData.VarAttribute{CellInfo{icell}.VarIndex_vector_y(2)}.CheckSub;
141%             ivar=[ivar CellInfo{icell}.VarIndex_vector_y(1)];
142%             ivar_1=[ivar_1 CellInfo{icell}.VarIndex_vector_y(2)];
143%         end
144%     end
145% end
146% % subtract fields if relevant
147% for imod=1:numel(ivar)
148%         VarName=ProjData.ListVarName{ivar(imod)};
149%         VarName_1=ProjData.ListVarName{ivar_1(imod)};
150%         ProjData.(VarName)=double(ProjData.(VarName))-double(ProjData.(VarName_1));
151%         ind_remove(ivar_1(imod))=1;
152% end
153% ProjData.ListVarName(find(ind_remove))=[];
154% ProjData.VarDimName(find(ind_remove))=[];
155% ProjData.VarAttribute(find(ind_remove))=[];
156
157%-----------------------------------------------------------------
158%project on a set of points
159function  [ProjData,errormsg]=proj_points(FieldData,ObjectData)%%
160%-------------------------------------------------------------------
161
162siz=size(ObjectData.Coord);
163width=0;
164if isfield(ObjectData,'Range')
165    width=ObjectData.Range(1,2);
166end
167if isfield(ObjectData,'RangeX')&&~isempty(ObjectData.RangeX)
168    width=max(ObjectData.RangeX);
169end
170if isfield(ObjectData,'RangeY')&&~isempty(ObjectData.RangeY)
171    width=max(width,max(ObjectData.RangeY));
172end
173if isfield(ObjectData,'RangeZ')&&~isempty(ObjectData.RangeZ)
174    width=max(width,max(ObjectData.RangeZ));
175end
176if isequal(ObjectData.ProjMode,'projection')
177    if width==0
178        errormsg='projection range around points needed';
179        return
180    end
181elseif  ~isequal(ObjectData.ProjMode,'interp_lin')
182    errormsg=(['ProjMode option ' ObjectData.ProjMode ' not available in proj_field']);
183        return
184end
185[ProjData,errormsg]=proj_heading(FieldData,ObjectData);
186if ~isempty(errormsg)
187    return
188end
189ProjData.NbDim=0;
190[CellInfo,NbDimArray,errormsg]=find_field_cells(FieldData);
191if ~isempty(errormsg)
192    errormsg=['error in proj_field/proj_points:' errormsg];
193    return
194end
195%LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
196for icell=1:length(CellInfo)
197    if NbDimArray(icell)<=1
198        continue %projection only for multidimensional fields
199    end
200    VarIndex=CellInfo{icell}.VarIndex;%  indices of the selected variables in the list FieldData.ListVarName
201    ivar_X=CellInfo{icell}.CoordIndex(end);
202    ivar_Y=CellInfo{icell}.CoordIndex(end-1);
203    ivar_Z=[];
204    if NbDimArray(icell)==3
205        ivar_Z=CellInfo{icell}.CoordIndex(1);
206    end
207    ivar_FF=[];
208    if isfield(CellInfo{icell},'VarIndex_errorflag')
209        ivar_FF=CellInfo{icell}.VarIndex_errorflag;
210        if numel(ivar_FF)>1
211            errormsg='multiple error flag input';
212            return
213        end
214    end   
215    % select types of  variables to be projected
216   ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'};
217      check_proj=false(size(FieldData.ListVarName));
218   for ilist=1:numel(ListProj)
219       if isfield(CellInfo{icell},ListProj{ilist})
220           check_proj(CellInfo{icell}.(ListProj{ilist}))=1;
221       end
222   end
223   VarIndex=find(check_proj);
224    ProjData.ListVarName={'Y','X','NbVal'};
225    ProjData.VarDimName={'nb_points','nb_points','nb_points'};
226    ProjData.VarAttribute{1}.Role='ancillary';
227    ProjData.VarAttribute{2}.Role='ancillary';
228    ProjData.VarAttribute{3}.Role='ancillary';
229    for ivar=VarIndex       
230        VarName=FieldData.ListVarName{ivar};
231        ProjData.ListVarName=[ProjData.ListVarName {VarName}];% add the current variable to the list of projected variables
232        ProjData.VarDimName=[ProjData.VarDimName {'nb_points'}]; % projected VarName has a single dimension called 'nb_points' (set of projection points)
233
234    end
235    if strcmp( CellInfo{icell}.CoordType,'scattered')
236        coord_x=FieldData.(FieldData.ListVarName{ivar_X});
237        coord_y=FieldData.(FieldData.ListVarName{ivar_Y});
238        test3D=0;% TEST 3D CASE : NOT COMPLETED ,  3D CASE : NOT COMPLETED
239        if length(ivar_Z)==1
240            coord_z=FieldData.(FieldData.ListVarName{ivar_Z});
241            test3D=1;
242        end
243   
244        for ipoint=1:siz(1)
245           Xpoint=ObjectData.Coord(ipoint,:);
246           distX=coord_x-Xpoint(1);
247           distY=coord_y-Xpoint(2);         
248           dist=distX.*distX+distY.*distY;
249           indsel=find(dist<width*width);
250           ProjData.X(ipoint,1)=Xpoint(1);
251           ProjData.Y(ipoint,1)=Xpoint(2);
252           if isequal(length(ivar_FF),1)
253               FFName=FieldData.ListVarName{ivar_FF};
254               FF=FieldData.(FFName)(indsel);
255               indsel=indsel(~FF);
256           end
257           ProjData.NbVal(ipoint,1)=length(indsel);
258            for ivar=VarIndex
259               VarName=FieldData.ListVarName{ivar};
260               if isempty(indsel)
261                    ProjData.(VarName)(ipoint,1)=NaN;
262               else
263                    Var=FieldData.(VarName)(indsel);
264                    ProjData.(VarName)(ipoint,1)=mean(Var);
265                    if isequal(ObjectData.ProjMode,'interp_lin')
266                         ProjData.(VarName)(ipoint,1)=griddata_uvmat(coord_x(indsel),coord_y(indsel),Var,Xpoint(1),Xpoint(2));
267                    end
268               end
269            end
270        end
271    else    %case of structured coordinates
272        if  strcmp( CellInfo{icell}.CoordType,'grid')
273            AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
274            AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
275            eval(['AX=FieldData.' AXName ';']);% set of x positions
276            eval(['AY=FieldData.' AYName ';']);% set of y positions 
277            AName=FieldData.ListVarName{VarIndex(1)};% a single variable assumed in the current cell
278            eval(['A=FieldData.' AName ';']);% scalar
279            npxy=size(A);         
280            %update VarDimName in case of components (non coordinate dimensions e;g. color components)
281            if numel(npxy)>NbDimArray(icell)
282                ProjData.VarDimName{end}={'nb_points','component'};
283            end
284            for idim=1:NbDimArray(icell) %loop on space dimensions
285                test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
286                test_coord(idim)=0;%test for defined coordinates, =0 by default
287                ivar=CellInfo{icell}.CoordIndex(idim);
288                Coord{idim}=FieldData.(FieldData.ListVarName{ivar}); % position for the first index
289                if numel(Coord{idim})==2
290                    DCoord_min(idim)= (Coord{idim}(2)-Coord{idim}(1))/(npxy(idim)-1);
291                    test_direct(idim)=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise
292                else
293                    DCoord=diff(Coord{idim});
294                    DCoord_min(idim)=min(DCoord);
295                    DCoord_max=max(DCoord);
296                    test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
297                    test_direct_min=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise
298                    if ~isequal(test_direct(idim),test_direct_min)
299                        errormsg=['non monotonic dimension variable # ' num2str(idim)  ' in proj_field.m'];
300                        return
301                    end
302                    test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
303                    test_coord(idim)=1;
304                end
305            end
306            DX=DCoord_min(2);
307            DY=DCoord_min(1);
308            for ipoint=1:siz(1)
309                xwidth=width/(abs(DX));
310                ywidth=width/(abs(DY));
311                i_min=round((ObjectData.Coord(ipoint,1)-Coord{2}(1))/DX+0.5-xwidth); %minimum index of the selected region
312                i_min=max(1,i_min);%restrict to field limit
313                i_plus=round((ObjectData.Coord(ipoint,1)-Coord{2}(1))/DX+0.5+xwidth);
314                i_plus=min(npxy(2),i_plus); %restrict to field limit
315                j_min=round((ObjectData.Coord(ipoint,2)-Coord{1}(1))/DY-ywidth+0.5);
316                j_min=max(1,j_min);
317                j_plus=round((ObjectData.Coord(ipoint,2)-Coord{1}(1))/DY+ywidth+0.5);
318                j_plus=min(npxy(1),j_plus);
319                ProjData.X(ipoint,1)=ObjectData.Coord(ipoint,1);
320                ProjData.Y(ipoint,1)=ObjectData.Coord(ipoint,2);
321                i_int=(i_min:i_plus);
322                j_int=(j_min:j_plus);
323                ProjData.NbVal(ipoint,1)=length(j_int)*length(i_int);
324                if isempty(i_int) || isempty(j_int)
325                   for ivar=VarIndex   
326                        eval(['ProjData.' FieldData.ListVarName{ivar} '(ipoint,:)=NaN;']);
327                   end
328                   errormsg=['no data points in the selected projection range ' num2str(width) ];
329                else
330                    %TODO: introduce circle in the selected subregion
331                    %[I,J]=meshgrid([1:j_int],[1:i_int]);
332                    for ivar=VarIndex   
333                        Avalue=FieldData.(FieldData.ListVarName{ivar})(j_int,i_int,:);
334                        ProjData.(FieldData.ListVarName{ivar})(ipoint,:)=mean(mean(Avalue));
335                    end
336                end
337            end
338        end
339   end
340end
341
342%-----------------------------------------------------------------
343%project in a patch
344function  [ProjData,errormsg]=proj_patch(FieldData,ObjectData)%%
345%-------------------------------------------------------------------
346[ProjData,errormsg]=proj_heading(FieldData,ObjectData);
347if ~isempty(errormsg)
348    return
349end
350%objectfield=fieldnames(ObjectData);
351widthx=0;
352widthy=0;
353if isfield(ObjectData,'RangeX') && ~isempty(ObjectData.RangeX)
354    widthx=max(ObjectData.RangeX);
355end
356if isfield(ObjectData,'RangeY') && ~isempty(ObjectData.RangeY)
357    widthy=max(ObjectData.RangeY);
358end
359
360%A REVOIR, GENERALISER: UTILISER proj_line
361ProjData.NbDim=1;
362ProjData.ListVarName={};
363ProjData.VarDimName={};
364ProjData.VarAttribute={};
365
366CoordMesh=zeros(1,numel(FieldData.ListVarName));
367if isfield (FieldData,'VarAttribute')
368    for iattr=1:length(FieldData.VarAttribute)%initialization of variable attribute values
369        if isfield(FieldData.VarAttribute{iattr},'Unit')
370            unit{iattr}=FieldData.VarAttribute{iattr}.Unit;
371        end
372        if isfield(FieldData.VarAttribute{iattr},'CoordMesh')
373            CoordMesh(iattr)=FieldData.VarAttribute{iattr}.CoordMesh;
374        end
375    end
376end
377
378%group the variables (fields of 'FieldData') in cells of variables with the same dimensions
379[CellInfo,NbDim,errormsg]=find_field_cells(FieldData);
380if ~isempty(errormsg)
381    errormsg=['error in proj_field/proj_patch:' errormsg];
382    return
383end
384
385%LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
386for icell=1:length(CellInfo)
387    CoordType=CellInfo{icell}.CoordType;
388    test_Amat=0;
389    if NbDim(icell)~=2% proj_patch acts only on fields of space dimension 2
390        continue
391    end
392    ivar_FF=[];
393    testfalse=isfield(CellInfo{icell},'VarIndex_errorflag');
394    if testfalse
395        ivar_FF=CellInfo{icell}.VarIndex_errorflag;
396        FFName=FieldData.ListVarName{ivar_FF};
397        errorflag=FieldData.(FFName);
398    end
399    % select types of  variables to be projected
400    ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'};
401    check_proj=false(size(FieldData.ListVarName));
402    for ilist=1:numel(ListProj)
403        if isfield(CellInfo{icell},ListProj{ilist})
404            check_proj(CellInfo{icell}.(ListProj{ilist}))=1;
405        end
406    end
407    VarIndex=find(check_proj);
408   
409    ivar_X=CellInfo{icell}.CoordIndex(end);
410    ivar_Y=CellInfo{icell}.CoordIndex(end-1);
411    ivar_Z=[];
412    if NbDim(icell)==3
413        ivar_Z=CellInfo{icell}.CoordIndex(1);
414    end
415    switch CellInfo{icell}.CoordType
416        case 'scattered' %case of unstructured coordinates
417            for ivar=[VarIndex ivar_X ivar_Y ivar_FF]
418                VarName=FieldData.ListVarName{ivar};
419                FieldData.(VarName)=reshape(FieldData.(VarName),[],1);
420            end
421            XName=FieldData.ListVarName{ivar_X};
422            YName=FieldData.ListVarName{ivar_Y};
423            coord_x=FieldData.(FieldData.ListVarName{ivar_X});
424            coord_y=FieldData.(FieldData.ListVarName{ivar_Y});
425            % image or 2D matrix
426        case 'grid' %case of structured coordinates
427            test_Amat=1;% test for image or 2D matrix
428            AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
429            AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
430            AX=FieldData.(AXName);% x coordinate
431            AY=FieldData.(AYName);% y coordinate
432            VarName=FieldData.ListVarName{VarIndex(1)};
433            DimValue=size(FieldData.(VarName));
434            if length(AX)==2
435                AX=linspace(AX(1),AX(end),DimValue(2));
436            end
437            if length(AY)==2
438                AY=linspace(AY(1),AY(end),DimValue(1));
439            end
440            if length(DimValue)==3
441                testcolor=1;
442                npxy(3)=3;
443            else
444                testcolor=0;
445                npxy(3)=1;
446            end
447            [Xi,Yi]=meshgrid(AX,AY);
448            npxy(1)=length(AY);
449            npxy(2)=length(AX);
450            Xi=reshape(Xi,npxy(1)*npxy(2),1);
451            Yi=reshape(Yi,npxy(1)*npxy(2),1);
452            for ivar=1:length(VarIndex)
453                VarName=FieldData.ListVarName{VarIndex(ivar)};
454                FieldData.(VarName)=reshape(FieldData.(VarName),npxy(1)*npxy(2),npxy(3)); % keep only non false vectors
455            end
456    end
457    %select the indices in the range of action
458    testin=[];%default
459    if isequal(ObjectData.Type,'rectangle')
460        if strcmp(CellInfo{icell}.CoordType,'scattered')
461            distX=abs(coord_x-ObjectData.Coord(1,1));
462            distY=abs(coord_y-ObjectData.Coord(1,2));
463            testin=distX<widthx & distY<widthy;
464        elseif test_Amat
465            distX=abs(Xi-ObjectData.Coord(1,1));
466            distY=abs(Yi-ObjectData.Coord(1,2));
467            testin=distX<widthx & distY<widthy;
468        end
469    elseif isequal(ObjectData.Type,'polygon')
470        if strcmp(CoordType,'scattered')
471            testin=inpolygon(coord_x,coord_y,ObjectData.Coord(:,1),ObjectData.Coord(:,2));
472        elseif strcmp(CoordType,'grid')
473            testin=inpolygon(Xi,Yi,ObjectData.Coord(:,1),ObjectData.Coord(:,2));
474        else%calculate the scalar
475            testin=[]; %A REVOIR
476        end
477    elseif isequal(ObjectData.Type,'ellipse')
478        X2Max=widthx*widthx;
479        Y2Max=(widthy)*(widthy);
480        if strcmp(CoordType,'scattered')
481            distX=(coord_x-ObjectData.Coord(1,1));
482            distY=(coord_y-ObjectData.Coord(1,2));
483            testin=(distX.*distX/X2Max+distY.*distY/Y2Max)<1;
484        elseif strcmp(CoordType,'grid') %case of usual 2x2 matrix
485            distX=(Xi-ObjectData.Coord(1,1));
486            distY=(Yi-ObjectData.Coord(1,2));
487            testin=(distX.*distX/X2Max+distY.*distY/Y2Max)<1;
488        end
489    end
490    %selected indices
491    if isequal(ObjectData.ProjMode,'outside')
492        testin=~testin;
493    end
494    if testfalse
495        testin=testin & (errorflag==0); % keep only non false vectors
496    end
497    indsel=find(testin);
498    for ivar=VarIndex
499        VarName=FieldData.ListVarName{ivar};
500        ProjData.([VarName 'Mean'])=mean(double(FieldData.(VarName)(indsel,:))); % take the mean in the selected region, for each color component
501        ProjData.([VarName 'Min'])=min(double(FieldData.(VarName)(indsel,:))); % take the min in the selected region , for each color component
502        ProjData.([VarName 'Max'])=max(double(FieldData.(VarName)(indsel,:))); % take the max in the selected region , for each color component
503        if isequal(CoordMesh(ivar),0)
504            [ProjData.([VarName 'Histo']),ProjData.(VarName)]=hist(double(FieldData.(VarName)(indsel,:,:)),100); % default histogram with 100 bins
505        else
506            ProjData.(VarName)=ProjData.([VarName 'Min'])+CoordMesh(ivar)/2:CoordMesh(ivar):ProjData.([VarName 'Max']); % list of bin values
507            ProjData.([VarName 'Histo'])=hist(double(FieldData.(VarName)(indsel,:)),ProjData.(VarName)); % histogram at predefined bin positions
508        end
509        ProjData.ListVarName=[ProjData.ListVarName {VarName} {[VarName 'Histo']} {[VarName 'Mean']} {[VarName 'Min']} {[VarName 'Max']}];
510        if test_Amat && testcolor
511            ProjData.VarDimName=[ProjData.VarDimName  {VarName} {{VarName,'rgb'}} {'rgb'} {'rgb'} {'rgb'}];%{{'nb_point','rgb'}};
512        else
513            ProjData.VarDimName=[ProjData.VarDimName {VarName} {VarName} {'one'} {'one'} {'one'}];
514        end
515        if isfield(FieldData,'VarAttribute')&& numel(FieldData.VarAttribute)>=ivar
516            ProjData.VarAttribute=[ProjData.VarAttribute FieldData.VarAttribute{ivar} {[]} {[]} {[]} {[]}];
517        end
518    end
519end
520
521%-----------------------------------------------------------------
522%project on a line
523% AJOUTER flux,circul,error
524% OUTPUT:
525% ProjData: projected field
526%
527function  [ProjData,errormsg] = proj_line(FieldData, ObjectData)
528%-----------------------------------------------------------------
529[ProjData,errormsg]=proj_heading(FieldData,ObjectData);%transfer global attributes
530if ~isempty(errormsg)
531    return
532end
533ProjData.NbDim=1;
534%initialisation of the input parameters and defaultoutput
535ProjMode=ObjectData.ProjMode; %rmq: ProjMode always defined from input={'projection','interp_lin','interp_tps'}
536% ProjAngle=90; %90 degrees projection by default
537width=0;
538if isfield(ObjectData,'RangeY')
539    width=max(ObjectData.RangeY);%Rangey needed bfor mode 'projection'
540end
541% default output
542errormsg='';%default
543Xline=[];
544flux=0;
545circul=0;
546liny=ObjectData.Coord(:,2);
547NbPoints=size(ObjectData.Coord,1);
548testfalse=0;
549ListIndex={};
550
551
552%% projection line: object types selected from  proj_field='line','polyline','polygon','rectangle','ellipse':
553LineCoord=ObjectData.Coord;
554switch ObjectData.Type
555    case 'ellipse'
556        LineLength=2*pi*ObjectData.RangeX*ObjectData.RangeY;
557        NbSegment=0;
558    case 'rectangle'
559        LineCoord([1 4],1)=ObjectData.Coord(1,1)-ObjectData.RangeX;
560        LineCoord([1 2],2)=ObjectData.Coord(1,2)-ObjectData.RangeY;
561        LineCoord([2 3],1)=ObjectData.Coord(1,1)+ObjectData.RangeX;
562        LineCoord([4 1],2)=ObjectData.Coord(1,2)+ObjectData.RangeY;
563    case 'polygon'
564        LineCoord(NbPoints+1)=LineCoord(1);
565end
566if ~strcmp(ObjectData.Type,'ellipse')
567    if ~strcmp(ObjectData.Type,'rectangle') && NbPoints<2
568        return% line needs at least 2 points to be defined
569    end
570    dlinx=diff(LineCoord(:,1));
571    dliny=diff(LineCoord(:,2));
572    [theta,dlength]=cart2pol(dlinx,dliny);%angle and length of each segment
573    LineLength=sum(dlength);
574    NbSegment=numel(LineLength);
575end
576CheckClosedLine=~isempty(find(strcmp(ObjectData.Type,{'rectangle','ellipse','polygon'})));
577
578%     x = a \ \cosh \mu \ \cos \nu
579%
580%     y = a \ \sinh \mu \ \sin \nu
581
582%% angles of the polyline and boundaries of action for mode 'projection'
583
584% determine a rectangles at +-width from the line (only used for the ProjMode='projection or 'interp_tps')
585xsup=zeros(1,NbPoints); xinf=zeros(1,NbPoints); ysup=zeros(1,NbPoints); yinf=zeros(1,NbPoints);
586if isequal(ProjMode,'projection')
587    if strcmp(ObjectData.Type,'line')
588        xsup=ObjectData.Coord(:,1)-width*sin(theta);
589        xinf=ObjectData.Coord(:,1)+width*sin(theta);
590        ysup=ObjectData.Coord(:,2)+width*cos(theta);
591        yinf=ObjectData.Coord(:,2)-width*cos(theta);
592    else
593        errormsg='mode projection only available for simple line, use interpolation otherwise';
594        return
595    end
596else % need to define the set of interpolation points
597    if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX)
598        DX=abs(ObjectData.DX);%mesh of interpolation points along the line
599        if CheckClosedLine
600            NbPoint=ceil(LineLength/DX);
601            DX=LineLength/NbPoint;%adjust DX to get an integer nbre of intervals in a closed line
602            DX_edge=DX/2;
603        else
604            DX_edge=(LineLength-DX*floor(LineLength/DX))/2;%margin from the first point and first interpolation point, the same for the end point
605        end
606        XI=[];
607        YI=[];
608        ThetaI=[];
609        dlengthI=[];
610        if strcmp(ObjectData.Type,'ellipse')
611            phi=(DX_edge:DX:LineLength)*2*pi/LineLength;
612            XI=ObjectData.RangeX*cos(phi);
613            YI=ObjectData.RangeY*sin(phi);
614            dphi=2*pi*DX/LineLength;
615            [ThetaI,dlengthI]=cart2pol(-ObjectData.RangeX*sin(phi)*dphi,ObjectData.RangeY*cos(phi)*dphi);
616        else
617            for isegment=1:NbSegment
618                costheta=cos(theta(isegment));
619                sintheta=sin(theta(isegment));
620                %                 XIsegment=LineCoord(isegment,1)+DX_edge*costheta:DX*costheta:LineCoord(isegment+1,1));
621                %                 YIsegment=(LineCoord(isegment,2)+DX_edge*sintheta:DX*sintheta:LineCoord(isegment+1,2));
622                NbInterval=floor((dlength(isegment)-DX_edge)/DX);
623                LastX=DX_edge+DX*NbInterval;
624                NbPoint=NbInterval+1;
625                XIsegment=linspace(LineCoord(isegment,1)+DX_edge*costheta,LineCoord(isegment,1)+LastX*costheta,NbPoint);
626                YIsegment=linspace(LineCoord(isegment,2)+DX_edge*sintheta,LineCoord(isegment,2)+LastX*sintheta,NbPoint);
627                XI=[XI XIsegment];
628                YI=[YI YIsegment];
629                ThetaI=[ThetaI theta(isegment)*ones(1,numel(XIsegment))];
630                dlengthI=[dlengthI DX*ones(1,numel(XIsegment))];
631                DX_edge=DX-(dlength(isegment)-LastX);%edge for the next segment set to keep DX=DX_end+DX_edge between two segments
632            end
633        end
634        Xproj=cumsum(dlengthI);
635    else
636        errormsg='abscissa mesh along line DX needed for interpolation';
637        return
638    end
639end
640
641
642%% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
643[CellInfo,NbDim,errormsg]=find_field_cells(FieldData);
644if ~isempty(errormsg)
645    errormsg=['error in proj_field/proj_line:' errormsg];
646    return
647end
648CellInfo=CellInfo(NbDim==2); %keep only the 2D cells
649%%%%%% TODO: treat 1D fields: project as identity so that P o P=P for projection operation
650
651%% loop on variable cells with the same space dimension 2
652ProjData.ListVarName={};
653ProjData.VarDimName={};
654check_abscissa=0;
655for icell=1:length(CellInfo)
656    % list of variable types to be projected
657    ListProj={'VarIndex_scalar','VarIndex_image','VarIndex_color','VarIndex_vector_x','VarIndex_vector_y'};
658    check_proj=false(size(FieldData.ListVarName));
659    for ilist=1:numel(ListProj)
660        if isfield(CellInfo{icell},ListProj{ilist})
661            check_proj(CellInfo{icell}.(ListProj{ilist}))=1;
662        end
663    end
664    VarIndex=find(check_proj);% indices of the variables to be projected
665   
666    %% identify vector components
667    %testU=isfield(CellInfo{icell},'VarIndex_vector_x') &&isfield(CellInfo{icell},'VarIndex_vector_y') ;% test for vectors
668    %     if testU
669    %         UName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x};
670    %         VName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y};
671    %         vector_x=FieldData.(UName);
672    %         vector_y=FieldData.(VName);
673    %     end
674    %identify error flag
675    errorflag=0; %default, no error flag
676    if isfield(CellInfo{icell},'VarIndex_errorflag');% test for error flag
677        FFName=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};
678        errorflag=FieldData.(FFName);
679    end
680    VarName=FieldData.ListVarName(VarIndex);% cell array of the names of variables to pje
681    ivar_U=[];
682    ivar_V=[];
683    %% check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y
684   
685    %         circul=0;
686    %         flux=0;
687    %%%%%%%  % A FAIRE CALCULER MEAN DES QUANTITES    %%%%%%
688    switch CellInfo{icell}.CoordType
689        %case of unstructured coordinates
690        case 'scattered'
691%             XName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
692%             YName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
693            coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});
694            coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)});
695           
696            if isequal(ProjMode,'projection')
697                if width==0
698                    errormsg='range of the projection object is missing';
699                    return
700                end
701                ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName(CellInfo{icell}.CoordIndex(end))];
702                ProjData.VarDimName=[ProjData.VarDimName FieldData.ListVarName(CellInfo{icell}.CoordIndex(end))];
703                nbvar=numel(ProjData.ListVarName);
704                ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
705                % select the (non false) input data located in the band of projection
706                flagsel=(errorflag==0) & ((coord_y -yinf(1))*(xinf(2)-xinf(1))>(coord_x-xinf(1))*(yinf(2)-yinf(1))) ...
707                    & ((coord_y -ysup(1))*(xsup(2)-xsup(1))<(coord_x-xsup(1))*(ysup(2)-ysup(1))) ...
708                    & ((coord_y -yinf(2))*(xsup(2)-xinf(2))>(coord_x-xinf(2))*(ysup(2)-yinf(2))) ...
709                    & ((coord_y -yinf(1))*(xsup(1)-xinf(1))<(coord_x-xinf(1))*(ysup(1)-yinf(1)));
710                coord_x=coord_x(flagsel);
711                coord_y=coord_y(flagsel);
712                costheta=cos(theta);
713                sintheta=sin(theta);
714                Xproj=(coord_x-ObjectData.Coord(1,1))*costheta + (coord_y-ObjectData.Coord(1,2))*sintheta; %projection on the line
715                [Xproj,indsort]=sort(Xproj);% sort points by increasing absissa along the projection line
716                ProjData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)})=Xproj;
717                for ivar=1:numel(VarIndex)
718                    ProjData.(VarName{ivar})=FieldData.(VarName{ivar})(flagsel);% restrict variables to the projection band
719                    ProjData.(VarName{ivar})=ProjData.(VarName{ivar})(indsort);% sort by absissa
720                    ProjData.ListVarName=[ProjData.ListVarName VarName{ivar}];
721                    ProjData.VarDimName=[ProjData.VarDimName FieldData.ListVarName(CellInfo{icell}.CoordIndex(end))];
722                    ProjData.VarAttribute{nbvar+ivar}=FieldData.VarAttribute{VarIndex(ivar)};%reproduce var attribute
723                    if isfield(ProjData.VarAttribute{nbvar+ivar},'Role')
724                        if  strcmp(ProjData.VarAttribute{nbvar+ivar}.Role,'vector_x');
725                            ivar_U=nbvar+ivar;
726                        elseif strcmp(ProjData.VarAttribute{nbvar+ivar}.Role,'vector_y');
727                            ivar_V=nbvar+ivar;
728                        end
729                    end
730                    ProjData.VarAttribute{ivar+nbvar}.Role='discrete';% will promote plots of the profiles with continuous lines
731                end
732            elseif isequal(ProjMode,'interp_lin')  %filtering %linear interpolation:
733                if ~check_abscissa
734                    XName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
735                    ProjData.ListVarName=[ProjData.ListVarName {XName}];
736                    ProjData.VarDimName=[ProjData.VarDimName {XName}];
737                    nbvar=numel(ProjData.ListVarName);
738                    ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
739                    check_abscissa=1; % define abcissa only once
740                end
741                if ~isequal(errorflag,0)
742                    VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};
743                    indsel=find(FieldData.(VarName_FF)==0);
744                    coord_x=coord_x(indsel);
745                    coord_y=coord_y(indsel);
746                    for ivar=1:numel(CellInfo{icell}.VarIndex)
747                        VarName=FieldData.ListVarName{CellInfo{icell}.VarIndex(ivar)};
748                        FieldData.(VarName)=FieldData.(VarName)(indsel);
749                    end
750                end
751                [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_x coord_y],FieldData,CellInfo{icell}.FieldName,XI,YI);
752                ProjData.X=Xproj;
753                nbvar=numel(ProjData.ListVarName);
754                ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
755                ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
756                for ivar=1:numel(VarAttribute)
757                    ProjData.VarDimName=[ProjData.VarDimName {XName}];
758                    if isfield(VarAttribute{ivar},'Role')
759                        if  strcmp(VarAttribute{ivar}.Role,'vector_x');
760                            ivar_U=ivar+nbvar;
761                        elseif strcmp(VarAttribute{ivar}.Role,'vector_y');
762                            ivar_V=ivar+nbvar;
763                        end
764                    end
765                    ProjData.VarAttribute{ivar+nbvar}.Role='continuous';% will promote plots of the profiles with continuous lines
766                    ProjData.(ListFieldProj{ivar})=ProjVar{ivar};
767                end
768            end
769        case 'tps'%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
770            if strcmp(ProjMode,'interp_tps')
771                Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex});
772                NbCentres=FieldData.(FieldData.ListVarName{CellInfo{icell}.NbCentres_tps});
773                SubRange=FieldData.(FieldData.ListVarName{CellInfo{icell}.SubRange_tps});
774                if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps')
775                    FieldVar=cat(3,FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x_tps}),FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y_tps}));
776                end
777                [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI));
778                ProjData.ListVarName=[ProjData.ListVarName {'X'}];
779                ProjData.VarDimName=[ProjData.VarDimName {'X'}];
780                ProjData.X=Xproj;
781                nbvar=numel(ProjData.ListVarName);
782                ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
783                ProjVarName=(fieldnames(DataOut))';
784                ProjData.ListVarName=[ProjData.ListVarName ProjVarName];
785                ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
786                for ivar=1:numel(VarAttribute)
787                    ProjData.VarDimName=[ProjData.VarDimName {'X'}];
788                    if isfield(VarAttribute{ivar},'Role')
789                        if  strcmp(VarAttribute{ivar}.Role,'vector_x');
790                            ivar_U=ivar+nbvar;
791                        elseif strcmp(VarAttribute{ivar}.Role,'vector_y');
792                            ivar_V=ivar+nbvar;
793                        end
794                    end
795                    ProjData.VarAttribute{ivar+nbvar}.Role='continuous';% will promote plots of the profiles with continuous lines
796                    ProjData.(ProjVarName{ivar})=DataOut.(ProjVarName{ivar});
797                end
798            end
799            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
800           
801        case 'grid'   %case of structured coordinates
802            if ~isequal(ObjectData.Type,'line')% exclude polyline
803                errormsg=['no  projection available on ' ObjectData.Type 'for structured coordinates']; %
804            else
805                test_Amat=1;%image or 2D matrix
806                test_interp2=0;%default
807                AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
808                AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
809                eval(['AX=FieldData.' AXName ';']);% set of x positions
810                eval(['AY=FieldData.' AYName ';']);% set of y positions
811                AName=FieldData.ListVarName{VarIndex(1)};
812                eval(['A=FieldData.' AName ';']);% scalar
813                npxy=size(A);
814                npx=npxy(2);
815                npy=npxy(1);
816                if numel(AX)==2
817                    DX=(AX(2)-AX(1))/(npx-1);
818                else
819                    DX_vec=diff(AX);
820                    DX=max(DX_vec);
821                    DX_min=min(DX_vec);
822                    if (DX-DX_min)>0.0001*abs(DX)
823                        test_interp2=1;
824                        DX=DX_min;
825                    end
826                end
827                if numel(AY)==2
828                    DY=(AY(2)-AY(1))/(npy-1);
829                else
830                    DY_vec=diff(AY);
831                    DY=max(DY_vec);
832                    DY_min=min(DY_vec);
833                    if (DY-DY_min)>0.0001*abs(DY)
834                        test_interp2=1;
835                        DY=DY_min;
836                    end
837                end
838                AXI=linspace(AX(1),AX(end), npx);%set of  x  positions for the interpolated input data
839                AYI=linspace(AY(1),AY(end), npy);%set of  x  positions for the interpolated input data
840                if isfield(ObjectData,'DX')
841                    DXY_line=ObjectData.DX;%mesh on the projection line
842                else
843                    DXY_line=sqrt(abs(DX*DY));% mesh on the projection line
844                end
845                dlinx=ObjectData.Coord(2,1)-ObjectData.Coord(1,1);
846                dliny=ObjectData.Coord(2,2)-ObjectData.Coord(1,2);
847                linelength=sqrt(dlinx*dlinx+dliny*dliny);
848                theta=angle(dlinx+i*dliny);%angle of the line
849                if isfield(FieldData,'RangeX')
850                    XMin=min(FieldData.RangeX);%shift of the origin on the line
851                else
852                    XMin=0;
853                end
854                eval(['ProjData.' AXName '=linspace(XMin,XMin+linelength,linelength/DXY_line+1);'])%abscissa of the new pixels along the line
855                y=linspace(-width,width,2*width/DXY_line+1);%ordintes of the new pixels (coordinate across the line)
856                eval(['npX=length(ProjData.' AXName ');'])
857                npY=length(y); %TODO: utiliser proj_grid
858                eval(['[X,Y]=meshgrid(ProjData.' AXName ',y);'])%grid in the line coordinates
859                XIMA=ObjectData.Coord(1,1)+(X-XMin)*cos(theta)-Y*sin(theta);
860                YIMA=ObjectData.Coord(1,2)+(X-XMin)*sin(theta)+Y*cos(theta);
861                XIMA=(XIMA-AX(1))/DX+1;%  index of the original image along x
862                YIMA=(YIMA-AY(1))/DY+1;% index of the original image along y
863                XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
864                YIMA=reshape(round(YIMA),1,npX*npY);
865                flagin=XIMA>=1 & XIMA<=npx & YIMA >=1 & YIMA<=npy;%flagin=1 inside the original image
866                ind_in=find(flagin);
867                ind_out=find(~flagin);
868                ICOMB=(XIMA-1)*npy+YIMA;
869                ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
870                nbcolor=1; %color images
871                if numel(npxy)==2
872                    nbcolor=1;
873                elseif length(npxy)==3
874                    nbcolor=npxy(3);
875                else
876                    errormsg='multicomponent field not projected';
877                    display(errormsg)
878                    return
879                end
880                nbvar=length(ProjData.ListVarName);% number of var from previous cells
881                ProjData.ListVarName=[ProjData.ListVarName {AXName}];
882                ProjData.VarDimName=[ProjData.VarDimName {AXName}];
883                for ivar=VarIndex
884                    %VarName{ivar}=FieldData.ListVarName{ivar};
885                    if test_interp2% interpolate on new grid
886                        FieldData.(FieldData.ListVarName{ivar})=interp2(FieldData.(AXName),FieldData.(AYName),FieldData.(FieldData.ListVarName{ivar}),AXI,AYI);%TO TEST
887                    end
888                    vec_A=reshape(squeeze(FieldData.(FieldData.ListVarName{ivar})),npx*npy,nbcolor); %put the original image in colum
889                    if nbcolor==1
890                        vec_B(ind_in)=vec_A(ICOMB);
891                        vec_B(ind_out)=zeros(size(ind_out));
892                        A_out=reshape(vec_B,npY,npX);
893                        ProjData.(FieldData.ListVarName{ivar}) =sum(A_out,1)/npY;
894                    elseif nbcolor==3
895                        vec_B(ind_in,1:3)=vec_A(ICOMB,:);
896                        vec_B(ind_out,1)=zeros(size(ind_out));
897                        vec_B(ind_out,2)=zeros(size(ind_out));
898                        vec_B(ind_out,3)=zeros(size(ind_out));
899                        A_out=reshape(vec_B,npY,npX,nbcolor);
900                        ProjData.(FieldData.ListVarName{ivar})=squeeze(sum(A_out,1)/npY);
901                    end
902                    ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName{ivar}];
903                    ProjData.VarDimName=[ProjData.VarDimName {AXName}];%to generalize with the initial name of the x coordinate
904                    ProjData.VarAttribute{ivar}.Role='continuous';% for plot with continuous line
905                end
906                if nbcolor==3
907                    ProjData.VarDimName{end}={AXName,'rgb'};
908                end
909            end
910    end
911    if ~isempty(ivar_U) && ~isempty(ivar_V)
912        vector_x =ProjData.(ProjData.ListVarName{ivar_U});
913        ProjData.(ProjData.ListVarName{ivar_U}) =cos(theta)*vector_x+sin(theta)*ProjData.(ProjData.ListVarName{ivar_V});
914        ProjData.(ProjData.ListVarName{ivar_V}) =-sin(theta)*vector_x+cos(theta)*ProjData.(ProjData.ListVarName{ivar_V});
915    end
916end
917
918% %shotarter case for horizontal or vertical line (A FAIRE
919% %     Rangx=[0.5 npx-0.5];%image coordiantes of corners
920% %     Rangy=[npy-0.5 0.5];
921% %     if isfield(Calib,'Pxcmx')&isfield(Calib,'Pxcmy')%old calib
922% %         Rangx=Rangx/Calib.Pxcmx;
923% %         Rangy=Rangy/Calib.Pxcmy;
924% %     else
925% %         [Rangx]=phys_XYZ(Calib,Rangx,[0.5 0.5],[0 0]);%case of translations without rotation and quadratic deformation
926% %         [xx,Rangy]=phys_XYZ(Calib,[0.5 0.5],Rangy,[0 0]);
927% %     end
928%
929% %     test_scal=0;%default% 3- 'UserData':(get(handles.Tag,'UserData')
930
931
932%-----------------------------------------------------------------
933%project on a plane
934% AJOUTER flux,circul,error
935function  [ProjData,errormsg] = proj_plane(FieldData, ObjectData)
936%-----------------------------------------------------------------
937
938%% rotation angles
939PlaneAngle=[0 0 0];
940norm_plane=[0 0 1];
941cos_om=1;
942sin_om=0;
943test90x=0;%=1 for 90 degree rotation alround x axis
944test90y=0;%=1 for 90 degree rotation alround y axis
945if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0])
946    test90y=isequal(ObjectData.Angle,[0 90 0]);
947    PlaneAngle=(pi/180)*ObjectData.Angle;
948    om=norm(PlaneAngle);%norm of rotation angle in radians
949    OmAxis=PlaneAngle/om; %unit vector marking the rotation axis
950    cos_om=cos(om);
951    sin_om=sin(om);
952    coeff=OmAxis(3)*(1-cos_om);
953    %components of the unity vector norm_plane normal to the projection plane
954    norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
955    norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
956    norm_plane(3)=OmAxis(3)*coeff+cos_om;
957end
958testangle=~isequal(PlaneAngle,[0 0 0])||~isequal(ObjectData.Coord(1:2),[0 0 ]) ;% && ~test90y && ~test90x;%=1 for slanted plane
959
960%% mesh sizes DX and DY
961DX=[];
962DY=[];%default
963if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX)
964    DX=abs(ObjectData.DX);%mesh of interpolation points
965elseif isfield(FieldData,'CoordMesh')
966    DX=FieldData.CoordMesh;
967end
968if isfield(ObjectData,'DY') && ~isempty(ObjectData.DY)
969    DY=abs(ObjectData.DY);%mesh of interpolation points
970elseif isfield(FieldData,'CoordMesh')
971    DY=FieldData.CoordMesh;
972end
973if  ~strcmp(ObjectData.ProjMode,'projection') && (isempty(DX)||isempty(DY))
974    errormsg='DX or DY not defined';
975    return
976end
977
978%% extrema along each axis
979testXMin=0;% test if min of X coordinates defined on the projection object, =0 by default
980testXMax=0;% test if max of X coordinates defined on the projection object, =0 by default
981testYMin=0;% test if min of Y coordinates defined on the projection object, =0 by default
982testYMax=0;% test if max of Y coordinates defined on the projection object, =0 by default
983if isfield(ObjectData,'RangeX') % rangeX defined by the projection object
984    XMin=min(ObjectData.RangeX);
985    XMax=max(ObjectData.RangeX);
986    testXMin=XMax>XMin;%=1 if XMin defined (i.e. RangeY has two distinct elements)
987    testXMax=1;% max of X coordinates defined on the projection object
988end
989if isfield(ObjectData,'RangeY') % rangeY defined by the projection object
990    YMin=min(ObjectData.RangeY);
991    YMax=max(ObjectData.RangeY);
992    testYMin=YMax>YMin;%=1 if YMin defined (i.e. RangeY has tow distinct elements)
993    testYMax=1;% max of Y coordinates defined on the projection object
994end
995width=0;%default width of the projection band
996if isfield(ObjectData,'RangeZ')
997    width=max(ObjectData.RangeZ);
998end
999
1000%% initiate Matlab  structure for physical field
1001[ProjData,errormsg]=proj_heading(FieldData,ObjectData);
1002if ~isempty(errormsg)
1003    return
1004end
1005
1006%% reproduce initial plane position and angle
1007if isfield(FieldData,'PlaneCoord')&&length(FieldData.PlaneCoord)==3&& isfield(ProjData,'ProjObjectCoord')
1008    if length(ProjData.ProjObjectCoord)==3% if the projection plane has a z coordinate
1009        if isfield(ProjData,'.PlaneCoord') && ~isequal(ProjData.PlaneCoord(3),ProjData.ProjObjectCoord) %check the consistency with the z coordinate of the field plane (set by calibration)
1010            errormsg='inconsistent z position for field and projection plane';
1011            return
1012        end
1013    else % the z coordinate is set only by the field plane (by calibration)
1014        ProjData.ProjObjectCoord(3)=FieldData.PlaneCoord(3);
1015    end
1016    if isfield(FieldData,'PlaneAngle')
1017        if isfield(ProjData,'ProjObjectAngle')
1018            if ~isequal(FieldData.PlaneAngle,ProjData.ProjObjectAngle) %check the consistency with the z coordinate of the field plane (set by calibration)
1019                errormsg='inconsistent plane angle for field and projection plane';
1020                return
1021            end
1022        else
1023            ProjData.ProjObjectAngle=FieldData.PlaneAngle;
1024        end
1025    end
1026end
1027ProjData.NbDim=2;
1028ProjData.ListVarName={};
1029ProjData.VarDimName={};
1030ProjData.VarAttribute={};
1031if ~isempty(DX) && ~isempty(DY)
1032    ProjData.CoordMesh=sqrt(DX*DY);%define typical data mesh, useful for mouse selection in plots
1033elseif isfield(FieldData,'CoordMesh')
1034    ProjData.CoordMesh=FieldData.CoordMesh;
1035end
1036error=0;%default
1037flux=0;
1038testfalse=0;
1039ListIndex={};
1040
1041%% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
1042[CellInfo,NbDimArray,errormsg]=find_field_cells(FieldData);
1043
1044if ~isempty(errormsg)
1045    errormsg=['error in proj_field/proj_plane:' errormsg];
1046    return
1047end
1048check_grid=zeros(size(CellInfo));% =1 if a grid is needed , =0 otherwise, for each field cell
1049
1050ProjMode=cell(size(CellInfo));
1051for icell=1:numel(CellInfo)
1052    ProjMode{icell}=ObjectData.ProjMode;% projection mode of the plane object
1053end
1054    icell_grid=[];% field cell index which defines the grid
1055if ~strcmp(ObjectData.ProjMode,'projection')
1056    %% define the new coordinates in case of interpolation on a imposed grid
1057    if ~testYMin
1058        errormsg='min Y value not defined for the projection grid';return
1059    end
1060    if ~testYMax
1061        errormsg='max Y value not defined for the projection grid';return
1062    end
1063    if ~testXMin
1064        errormsg='min X value not defined for the projection grid';return
1065    end
1066    if ~testXMax
1067        errormsg='max X value not defined for the projection grid';return
1068    end
1069else
1070    %% case of a grid requested by the input field
1071    for icell=1:numel(CellInfo)% TODO: recalculate coordinates here to get the bounds in the rotated coordinates
1072        if isfield(CellInfo{icell},'ProjModeRequest')
1073            switch CellInfo{icell}.ProjModeRequest
1074                case 'interp_lin'
1075                    ProjMode{icell}='interp_lin';
1076                case 'interp_tps'
1077                    ProjMode{icell}='interp_tps';
1078            end
1079        end
1080        if strcmp(ProjMode{icell},'interp_lin')||strcmp(ProjMode{icell},'interp_tps')
1081            check_grid(icell)=1;
1082        end
1083        if strcmp(CellInfo{icell}.CoordType,'grid')&&NbDimArray(icell)>=2
1084            if ~testangle && isempty(icell_grid)% if the input gridded data is not modified, choose the first one in case of multiple gridded field cells
1085                icell_grid=icell;
1086                ProjMode{icell}='projection';
1087            end
1088            check_grid(icell)=1;
1089        end
1090    end
1091    if ~isempty(find(check_grid))% if a grid is requested by the input field
1092        if isempty(icell_grid)%  if the grid is not given by cell #icell_grid
1093            if ~isfield(FieldData,'XMax')
1094                FieldData=find_field_bounds(FieldData);
1095            end
1096        end
1097    end
1098end
1099if ~isempty(find(check_grid))||~strcmp(ObjectData.ProjMode,'projection')%no existing gridded data used
1100    if isempty(icell_grid)||~strcmp(ObjectData.ProjMode,'projection')%no existing gridded data used
1101        AYName='coord_y';
1102        AXName='coord_x';
1103        if strcmp(ObjectData.ProjMode,'projection')
1104            ProjData.coord_y=[FieldData.YMin FieldData.YMax];%note that if projection is done on a grid, the Min and Max along each direction must have been defined
1105            ProjData.coord_x=[FieldData.XMin FieldData.XMax];
1106            coord_x_proj=FieldData.XMin:FieldData.CoordMesh:FieldData.XMax;
1107            coord_y_proj=FieldData.YMin:FieldData.CoordMesh:FieldData.YMax;
1108        else
1109            ProjData.coord_y=[ObjectData.RangeY(1) ObjectData.RangeY(2)];%note that if projection is done on a grid, the Min and Max along each direction must have been defined
1110            ProjData.coord_x=[ObjectData.RangeX(1) ObjectData.RangeX(2)];
1111            coord_x_proj=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2);
1112            coord_y_proj=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2);
1113        end
1114        [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
1115%         XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system
1116%         YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3));
1117    else% we use the existing grid from field cell #icell_grid
1118        NbDim=NbDimArray(icell_grid);
1119        AYName=FieldData.ListVarName{CellInfo{icell_grid}.CoordIndex(NbDim-1)};%name of input x coordinate (name preserved on projection)
1120        AXName=FieldData.ListVarName{CellInfo{icell_grid}.CoordIndex(NbDim)};%name of input y coordinate (name preserved on projection)
1121        ProjData.(AYName)=FieldData.(AYName); % new (projected ) y coordinates
1122        ProjData.(AXName)=FieldData.(AXName); % new (projected ) y coordinates
1123    end
1124    ProjData.ListVarName={AYName,AXName};
1125    ProjData.VarDimName={AYName,AXName};
1126    ProjData.VarAttribute={[],[]};
1127end
1128   
1129%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1130% LOOP ON FIELD CELLS, PROJECT VARIABLES
1131% CellVarIndex=cells of variable index arrays
1132%ivar_new=0; % index of the current variable in the projected field
1133% icoord=0;
1134nbcoord=0;%number of added coordinate variables brought by projection
1135%nbvar=0;
1136vector_x_proj=[];
1137vector_y_proj=[];
1138for icell=1:length(CellInfo)
1139    NbDim=NbDimArray(icell);
1140    if NbDim<2
1141        continue % only cells represnting 2D or 3D fields are involved
1142    end
1143    VarIndex=CellInfo{icell}.VarIndex;%  indices of the selected variables in the list FieldData.ListVarName
1144    %dimensions
1145    DimCell=FieldData.VarDimName{VarIndex(1)};
1146    if ischar(DimCell)
1147        DimCell={DimCell};%name of dimensions
1148    end
1149    coord_z=0;%default
1150    ListVarName={};% initiate list of projected variables for cell # icell
1151    VarDimName={};% initiate coresponding list of dimensions for cell # icell
1152    VarAttribute={};% initiate coresponding list of var attributes  for cell # icell
1153    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1154    switch CellInfo{icell}.CoordType
1155       
1156        case 'scattered'
1157            %% case of input fields with unstructured coordinates (applies for projMode ='projection' or 'interp_lin')
1158            if strcmp(ProjMode{icell},'interp_tps')
1159                continue %skip for next cell (needs tps field cell)
1160            end
1161            coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});% initial x coordinates
1162            coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)});% initial y coordinates
1163            check3D=(numel(CellInfo{icell}.CoordIndex)==3);
1164            if check3D
1165                coord_z=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});
1166            end
1167           
1168            % translate  initial coordinates to account for the new origin
1169            coord_x=coord_x-ObjectData.Coord(1,1);
1170            coord_y=coord_y-ObjectData.Coord(1,2);
1171            if check3D
1172                coord_z=coord_z-ObjectData.Coord(1,3);
1173            end
1174           
1175            % selection of the vectors in the projection range (3D case)
1176            if check3D &&  width > 0
1177                %components of the unitiy vector normal to the projection plane
1178                fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane
1179                indcut=find(abs(fieldZ) <= width);
1180                for ivar=VarIndex
1181                    VarName=FieldData.ListVarName{ivar};
1182                    FieldData.(VarName)=FieldData.(VarName)(indcut);
1183                end
1184                coord_x=coord_x(indcut);
1185                coord_y=coord_y(indcut);
1186                coord_z=coord_z(indcut);
1187            end
1188           
1189            %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane
1190            Psi=PlaneAngle(1);
1191            Theta=PlaneAngle(2);
1192            Phi=PlaneAngle(3);
1193            if testangle && ~test90y && ~test90x;%=1 for slanted plane
1194                coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
1195                coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
1196                coord_Y=coord_Y+coord_z *sin(Theta);
1197                coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
1198                coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
1199            else
1200                coord_X=coord_x;
1201                coord_Y=coord_y;
1202            end
1203           
1204            %restriction to the range of X and Y if imposed by the projection object
1205            testin=ones(size(coord_X)); %default
1206            testbound=0;
1207            if testXMin
1208                testin=testin & (coord_X >= XMin);
1209                testbound=1;
1210            end
1211            if testXMax
1212                testin=testin & (coord_X <= XMax);
1213                testbound=1;
1214            end
1215            if testYMin
1216                testin=testin & (coord_Y >= YMin);
1217                testbound=1;
1218            end
1219            if testYMin
1220                testin=testin & (coord_Y <= YMax);
1221                testbound=1;
1222            end
1223            if testbound
1224                indcut=find(testin);
1225                if isempty(indcut)
1226                    errormsg='data outside the bounds of the projection object';
1227                    return
1228                end
1229                for ivar=VarIndex
1230                    VarName=FieldData.ListVarName{ivar};
1231                    FieldData.(VarName)=FieldData.(VarName)(indcut);
1232                end
1233                coord_X=coord_X(indcut);
1234                coord_Y=coord_Y(indcut);
1235                if check3D
1236                    coord_Z=coord_Z(indcut);
1237                end
1238            end
1239           
1240            % two cases of projection for scattered coordinates
1241            switch ProjMode{icell}
1242                case 'projection'
1243                    nbvar=0;
1244                    %nbvar=numel(ProjData.ListVarName);
1245                    for ivar=VarIndex %transfer variables to the projection plane
1246                        VarName=FieldData.ListVarName{ivar};
1247                        if ivar==CellInfo{icell}.CoordIndex(end)
1248                            ProjData.(VarName)=coord_X;
1249                        elseif ivar==CellInfo{icell}.CoordIndex(end-1)  % y coordinate
1250                            ProjData.(VarName)=coord_Y;
1251                        elseif ~(check3D && ivar==CellInfo{icell}.CoordIndex(1)) % other variables (except Z coordinate wyhich is not reproduced)
1252                            ProjData.(VarName)=FieldData.(VarName);
1253                        end
1254                        if ~(check3D && ivar==CellInfo{icell}.CoordIndex(1))
1255                            ListVarName=[ListVarName VarName];
1256                            VarDimName=[VarDimName DimCell];
1257                            nbvar=nbvar+1;
1258                            if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
1259                                VarAttribute{nbvar}=FieldData.VarAttribute{ivar};
1260                            end
1261                        end
1262                    end
1263                case 'interp_lin'%interpolate data on a regular grid
1264                    if isfield(CellInfo{icell},'VarIndex_errorflag')
1265                        VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};
1266                        indsel=find(FieldData.(VarName_FF)==0);
1267                        coord_X=coord_X(indsel);
1268                        coord_Y=coord_Y(indsel);
1269                        for ivar=1:numel(CellInfo{icell}.VarIndex)
1270                            VarName=FieldData.ListVarName{CellInfo{icell}.VarIndex(ivar)};
1271                            FieldData.(VarName)=FieldData.(VarName)(indsel);
1272                        end
1273                    end
1274                    % interpolate and calculate field on the grid
1275                    [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldData,CellInfo{icell}.FieldName,XI,YI);
1276                   
1277                    if isfield(CellInfo{icell},'CheckSub') && CellInfo{icell}.CheckSub && ~isempty(vector_x_proj)
1278                        ProjData.(FieldData.ListVarName{vector_x_proj})=ProjData.(FieldData.ListVarName{vector_x_proj})-VarVal{1};
1279                        ProjData.(FieldData.ListVarName{vector_y_proj})=ProjData.(FieldData.ListVarName{vector_y_proj})-VarVal{2};
1280                        ListVarName={};% no new variable
1281                        VarAttribute={};
1282                    else
1283                        VarDimName=cell(size(ListVarName));
1284                        for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
1285                            ListVarName{ilist}=regexprep(ListVarName{ilist},'(.+','');
1286                            if ~isempty(find(strcmp(ListVarName{ilist},ProjData.ListVarName)))
1287                                ListVarName{ilist}=[ListVarName{ilist} '_1'];
1288                            end
1289                            ProjData.(ListVarName{ilist})=VarVal{ilist};
1290                            VarDimName{ilist}={'coord_y','coord_x'};
1291                        end
1292                    end
1293                    if isfield (CellInfo{icell},'VarIndex_vector_x')&& isfield (CellInfo{icell},'VarIndex_vector_y')
1294                    vector_x_proj=CellInfo{icell}.VarIndex_vector_x; %preserve for next cell
1295                    vector_y_proj=CellInfo{icell}.VarIndex_vector_y; %preserve for next cell
1296                    end
1297            end
1298           
1299        case 'tps'
1300            %% case of tps data (applies only in interp_tps mode)
1301            if strcmp(ProjMode{icell},'interp_tps')
1302                Coord=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex});
1303                NbCentres=FieldData.(FieldData.ListVarName{CellInfo{icell}.NbCentres_tps});
1304                SubRange=FieldData.(FieldData.ListVarName{CellInfo{icell}.SubRange_tps});
1305                if isfield(CellInfo{icell},'VarIndex_vector_x_tps')&&isfield(CellInfo{icell},'VarIndex_vector_y_tps')
1306                    FieldVar=cat(3,FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x_tps}),FieldData.(FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y_tps}));
1307                end
1308                [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbCentres,SubRange,FieldVar,CellInfo{icell}.FieldName,cat(3,XI,YI));
1309                ListVarName=(fieldnames(DataOut))';
1310                VarDimName=cell(size(ListVarName));
1311                for ilist=1:numel(ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
1312                    VarName=ListVarName{ilist};
1313                    ProjData.(VarName)=DataOut.(VarName);
1314                    VarDimName{ilist}={'coord_y','coord_x'};
1315                end
1316            end
1317           
1318        case 'grid'
1319            %% case of input fields defined on a structured  grid
1320            VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
1321            DimValue=size(FieldData.(VarName));%input matrix dimensions
1322            DimValue(DimValue==1)=[];%remove singleton dimensions
1323            NbDim=numel(DimValue);%update number of space dimensions
1324            nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
1325            if NbDim>=3
1326                if NbDim>3
1327                    errormsg='matrices with more than 3 dimensions not handled';
1328                    return
1329                else
1330                    if numel(CellInfo{icell}.CoordIndex)==2% the third matrix dimension does not correspond to a space coordinate
1331                        nbcolor=DimValue(3);
1332                        DimValue(3)=[]; %number of 'color' components updated
1333                        NbDim=2;% space dimension set to 2
1334                    end
1335                end
1336            end
1337            Coord_z=[];
1338            Coord_y=[];
1339            Coord_x=[];
1340           
1341            if testangle
1342                ProjMode{icell}='interp_lin'; %request linear interpolation for projection on a tilted plane
1343            end
1344           
1345            if isequal(ProjMode{icell},'projection')% && (~testangle || test90y || test90x)
1346                if  NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax% no range restriction
1347                    ListVarName=[ListVarName FieldData.ListVarName(VarIndex)];
1348                    VarDimName=[VarDimName FieldData.VarDimName(VarIndex)];
1349                    if isfield(FieldData,'VarAttribute')
1350                        VarAttribute=[VarAttribute FieldData.VarAttribute(VarIndex)];
1351                    end
1352                    ProjData.(AYName)=FieldData.(AYName);
1353                    ProjData.(AXName)=FieldData.(AXName);
1354                    for ivar=VarIndex
1355                        VarName=FieldData.ListVarName{ivar};
1356                        ProjData.(VarName)=FieldData.(VarName);% no change by projection
1357                    end
1358                else
1359                    Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});
1360                    Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});
1361                    if NbDim==3
1362                        Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});
1363                    end
1364                    if numel(Coord{NbDim-1})==2
1365                        DY=(Coord{NbDim-1}(2)-Coord{NbDim-1}(1))/(DimValue(1)-1);
1366                    end
1367                    if numel(Coord{NbDim})==2
1368                        DX=(Coord{NbDim}(2)-Coord{NbDim}(1))/(DimValue(2)-1);
1369                    end
1370                    if testYMax
1371                         YIndexMax=(YMax-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the max y value for the new field
1372                        if testYMin%test_direct(indY)
1373                            YIndexMin=(YMin-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the min x value for the new field
1374                        else
1375                            YIndexMin=1;
1376                        end         
1377                    else
1378                        YIndexMax=Coord{NbDim-1}(end)/DY;
1379                        YIndexMin=1;
1380                    end
1381                    if testXMax
1382                         XIndexMax=(XMax-Coord{NbDim}(1))/DY+1;% matrix index corresponding to the max y value for the new field
1383                        if testYMin%test_direct(indY)
1384                            XIndexMin=(XMin-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the min x value for the new field
1385                        else
1386                            XIndexMin=1;
1387                        end         
1388                    else
1389                        XIndexMax=Coord{NbDim}(end)/DX;
1390                        XIndexMin=1;
1391                    end
1392                    YIndexRange(1)=ceil(min(YIndexMin,YIndexMax));%first y index to select from the previous field
1393                    YIndexRange(1)=max(YIndexRange(1),1);% avoid bound lower than the first index
1394                    YIndexRange(2)=floor(max(YIndexMin,YIndexMax));%last y index to select from the previous field
1395                    YIndexRange(2)=min(YIndexRange(2),DimValue(NbDim-1));% limit to the last available index
1396                    XIndexRange(1)=ceil(min(XIndexMin,XIndexMax));%first x index to select from the previous field
1397                    XIndexRange(1)=max(XIndexRange(1),1);% avoid bound lower than the first index
1398                    XIndexRange(2)=floor(max(XIndexMin,XIndexMax));%last x index to select from the previous field
1399                    XIndexRange(2)=min(XIndexRange(2),DimValue(NbDim));% limit to the last available index
1400                    if test90y
1401                        ind_new=[3 2 1];
1402                        DimCell={AYProjName,AXProjName};
1403                        iz=ceil((ObjectData.Coord(1,1)-Coord{3}(1))/DX)+1;
1404                        for ivar=VarIndex
1405                            VarName=FieldData.ListVarName{ivar};
1406                            ListVarName=[ListVarName VarName];
1407                            VarDimName=[VarDimName {DimCell}];
1408                            VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
1409                            ProjData.(VarName)=permute(FieldData.(VarName),ind_new);% permute x and z indices for 90 degree rotation
1410                            ProjData.(VarName)=squeeze(ProjData.(VarName)(iz,:,:));% select the z index iz
1411                        end
1412                        ProjData.(AYName)=[Ybound(1) Ybound(2)]; %record the new (projected ) y coordinates
1413                        ProjData.(AXName)=[Coord{1}(end),Coord{1}(1)]; %record the new (projected ) x coordinates
1414                    else
1415                        if NbDim==3
1416                            DZ=(Coord{1}(end)-Coord{1}(1))/(numel(Coord{1})-1);
1417                            DimCell(1)=[]; %suppress z variable
1418                            DimValue(1)=[];
1419                            test_direct=1;%TOdo; GENERALIZE, SEE CASE OF points
1420                            if test_direct(1)
1421                                iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1;
1422                            else
1423                                iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1;
1424                            end
1425                        end
1426                        for ivar=VarIndex% loop on non coordinate variables
1427                            VarName=FieldData.ListVarName{ivar};
1428                            ListVarName=[ListVarName VarName];
1429                            VarDimName=[VarDimName {DimCell}];
1430                            if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar
1431                                VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar};
1432                            end
1433                            if NbDim==3
1434                                ProjData.(VarName)=squeeze(FieldData.(VarName)(iz,YIndexRange(1):YIndexRange(end),XIndexRange(1):XIndexRange(end)));
1435                            else
1436                                ProjData.(VarName)=FieldData.(VarName)(YIndexRange(1):YIndexRange(end),XIndexRange(1):XIndexRange(end),:);
1437                            end
1438                        end
1439                        if testXMax
1440                         ProjData.(AXName)=Coord{NbDim}(1)+DX*(XIndexRange-1); %record the new (projected ) x coordinates
1441                        else
1442                          ProjData.(AXName)=FieldData.(AXName);
1443                        end
1444                        if testYMax
1445                            ProjData.(AYName)=Coord{NbDim-1}(1)+DY*(YIndexRange-1); %record the new (projected ) x coordinates
1446                        else
1447                          ProjData.(AYName)=FieldData.(AYName);
1448                        end                           
1449                    end
1450                end
1451            else       % case with interpolation on a grid
1452                if NbDim==2 %2D case
1453                    if isequal(ProjMode{icell},'interp_tps')
1454                        npx_interp_tps=ceil(abs(DX/DAX));
1455                        npy_interp_tps=ceil(abs(DY/DAY));
1456                        Minterp_tps=ones(npy_interp_tps,npx_interp_tps)/(npx_interp_tps*npy_interp_tps);
1457                        test_interp_tps=1;
1458                    else
1459                        test_interp_tps=0;
1460                    end
1461                    Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});
1462                    Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});
1463                    if ~(testXMin && testYMin)% % if the range of the projected coordinates is not fully defined by the projection object, find the extrema of the projected field
1464                        xcorner=[min(Coord{NbDim}) max(Coord{NbDim}) max(Coord{NbDim}) min(Coord{NbDim})]-ObjectData.Coord(1,1);% corner absissa of the original grid with respect to the new origin
1465                        ycorner=[min(Coord{NbDim-1}) min(Coord{NbDim-1}) max(Coord{NbDim-1}) max(Coord{NbDim-1})]-ObjectData.Coord(1,2);% corner ordinates of the original grid
1466                        xcor_new=xcorner*cos(PlaneAngle(3))+ycorner*sin(PlaneAngle(3));%coordinates of the corners in new frame
1467                        ycor_new=-xcorner*sin(PlaneAngle(3))+ycorner*cos(PlaneAngle(3));
1468                        if ~testXMin
1469                            XMin=min(xcor_new);
1470                        end
1471                        if ~testXMax
1472                            XMax=max(xcor_new);
1473                        end
1474                        if ~testYMin
1475                            YMin=min(ycor_new);
1476                        end
1477                        if ~testYMax
1478                            YMax=max(ycor_new);
1479                        end
1480                    end
1481                    coord_x_proj=XMin:DX:XMax;
1482                    coord_y_proj=YMin:DY:YMax;
1483                    ProjData.(AYName)=[coord_y_proj(1) coord_y_proj(end)]; %record the new (projected ) y coordinates
1484                    ProjData.(AXName)=[coord_x_proj(1) coord_x_proj(end)]; %record the new (projected ) x coordinates
1485                    [X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
1486                    XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system
1487                    YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3));
1488                    if numel(Coord{1})==2% x coordiante defiend by its bounds, get the whole set
1489                        Coord{1}=linspace(Coord{1}(1),Coord{1}(2),CellInfo{icell}.CoordSize(1));
1490                    end
1491                    if numel(Coord{2})==2% y coordiante defiend by its bounds, get the whole set
1492                        Coord{2}=linspace(Coord{2}(1),Coord{2}(2),CellInfo{icell}.CoordSize(2));
1493                    end
1494                    [X,Y]=meshgrid(Coord{2},Coord{1});%initial coordinates
1495                    %name of error flag variable
1496                    FFName='FF';%default name (if not already used)
1497                    if isfield(ProjData,'FF')
1498                        ind=1;
1499                        while isfield(ProjData,['FF_' num2str(ind)])
1500                            ind=ind+1;
1501                        end
1502                        FFName=['FF_' num2str(ind)];% append an index to the name of error flag, FF_1,FF_2...
1503                    end
1504                    % project all variables in the cell
1505                    for ivar=VarIndex
1506                        VarName=FieldData.ListVarName{ivar};
1507                        if size(FieldData.(VarName),3)==1
1508                            ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)),XI,YI,'*linear');%interpolation fct
1509                        else
1510                            ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)(:,:,1)),XI,YI,'*linear');
1511                            for icolor=2:size(FieldData.(VarName),3)% project 'color' components
1512                                ProjData.(VarName)=cat(3,ProjData.(VarName),interp2(X,Y,double(FieldData.(VarName)(:,:,icolor)),XI,YI,'*linear'));
1513                            end
1514                        end
1515                        if isa(FieldData.(VarName),'uint8')
1516                            ProjData.(VarName)=uint8(ProjData.(VarName));%put result to integer 8 bits if the initial field is integer (image)
1517                        elseif isa(FieldData.(VarName),'uint16')
1518                            ProjData.(VarName)=uint16(ProjData.(VarName));%put result to integer 16 bits if the initial field is integer (image)
1519                        end
1520                        ListVarName=[ListVarName VarName];
1521                        DimCell(1:2)={AYName,AXName};
1522                        VarDimName=[VarDimName {DimCell}];
1523                        if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar
1524                            VarAttribute{length(ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
1525                        end;
1526                        ProjData.(FFName)=isnan(ProjData.(VarName));%detact NaN (points outside the interpolation range)
1527                        ProjData.(VarName)(ProjData.(FFName))=0; %set to 0 the NaN data
1528                    end
1529                    %update list of variables with error flag
1530                    ListVarName=[ListVarName FFName];
1531                    VarDimName=[VarDimName {DimCell}];
1532                    VarAttribute{numel(ListVarName)}.Role='errorflag';
1533                elseif ~testangle
1534                    % unstructured z coordinate
1535                    test_sup=(Coord{1}>=ObjectData.Coord(1,3));
1536                    iz_sup=find(test_sup);
1537                    iz=iz_sup(1);
1538                    if iz>=1 & iz<=npz
1539                        %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
1540                        %ProjData.DimValue=[ProjData.DimValue npY npX];
1541                        for ivar=VarIndex
1542                            VarName=FieldData.ListVarName{ivar};
1543                            ListVarName=[ListVarName VarName];
1544                            VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
1545                            ProjData.(VarName)=squeeze(FieldData.(VarName)(iz,:,:));% select the z index iz
1546                            %TODO : do a vertical average for a thick plane
1547                            if test_interp(2) || test_interp(3)
1548                                ProjData.(VarName)=interp2(Coord{3},Coord{2},ProjData.(VarName),Coord_x,Coord_y);
1549                            end
1550                        end
1551                    end
1552                else
1553                    errormsg='projection of structured coordinates on oblique plane not yet implemented';
1554                    %TODO: use interp_lin3
1555                    return
1556                end
1557            end
1558    end
1559    % update the global list of projected variables:
1560    ProjData.ListVarName=[ProjData.ListVarName ListVarName];
1561    ProjData.VarDimName=[ProjData.VarDimName VarDimName];
1562    ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
1563   
1564    %% projection of  velocity components in the rotated coordinates
1565    if testangle
1566        ivar_U=[];ivar_V=[];ivar_W=[];
1567        for ivar=1:numel(VarAttribute)
1568            if isfield(VarAttribute{ivar},'Role')
1569                if strcmp(VarAttribute{ivar}.Role,'vector_x')
1570                    ivar_U=ivar;
1571                elseif strcmp(VarAttribute{ivar}.Role,'vector_y')
1572                    ivar_V=ivar;
1573                elseif strcmp(VarAttribute{ivar}.Role,'vector_z')
1574                    ivar_W=ivar;
1575                end
1576            end
1577        end
1578        if ~isempty(ivar_U)
1579            if isempty(ivar_V)
1580                msgbox_uvmat('ERROR','v velocity component missing in proj_field.m')
1581                return
1582            else
1583                UName=ListVarName{ivar_U};
1584                VName=ListVarName{ivar_V};
1585                ProjData.(UName)=cos(PlaneAngle(3))*ProjData.(UName)+ sin(PlaneAngle(3))*ProjData.(VName);
1586                ProjData.(VName)=(-sin(PlaneAngle(3))*ProjData.(UName)+ cos(PlaneAngle(3))*ProjData.(VName));
1587                if ~isempty(ivar_W)
1588                    WName=FieldData.ListVarName{ivar_W};
1589                    ProjData.(VName)=ProjData.(VName)+ ProjData.(WName)*sin(Theta);%
1590                    ProjData.(WName)=NormVec_X*ProjData.(UName)+ NormVec_Y*ProjData.(VName)+ NormVec_Z* ProjData.(WName);
1591                end
1592            end
1593        end
1594    end
1595end
1596% %prepare substraction in case of two input fields
1597% SubData.ListVarName={};
1598% SubData.VarDimName={};
1599% SubData.VarAttribute={};
1600% check_remove=zeros(size(ProjData.ListVarName));
1601% for iproj=1:numel(ProjData.VarAttribute)
1602%     if isfield(ProjData.VarAttribute{iproj},'CheckSub')&&isequal(ProjData.VarAttribute{iproj}.CheckSub,1)
1603%         VarName=ProjData.ListVarName{iproj};
1604%         SubData.ListVarName=[SubData.ListVarName {VarName}];
1605%         SubData.VarDimName=[SubData.VarDimName ProjData.VarDimName{iproj}];
1606%         SubData.VarAttribute=[SubData.VarAttribute ProjData.VarAttribute{iproj}];
1607%         SubData.(VarName)=ProjData.(VarName);
1608%         check_remove(iproj)=1;
1609%     end
1610% end
1611% if ~isempty(find(check_remove))
1612%     ind_remove=find(check_remove);
1613%     ProjData.ListVarName(ind_remove)=[];
1614%     ProjData.VarDimName(ind_remove)=[];
1615%     ProjData.VarAttribute(ind_remove)=[];
1616%     ProjData=sub_field(ProjData,[],SubData);
1617% end
1618
1619%-----------------------------------------------------------------
1620%projection in a volume
1621function  [ProjData,errormsg] = proj_volume(FieldData, ObjectData)
1622
1623%-----------------------------------------------------------------
1624ProjData=FieldData;%default output
1625
1626%% axis origin
1627if isempty(ObjectData.Coord)
1628    ObjectData.Coord(1,1)=0;%origin of the plane set to [0 0] by default
1629    ObjectData.Coord(1,2)=0;
1630    ObjectData.Coord(1,3)=0;
1631end
1632
1633%% rotation angles
1634VolumeAngle=[0 0 0];
1635norm_plane=[0 0 1];
1636if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 3])&& ~isequal(ObjectData.Angle,[0 0 0])
1637    PlaneAngle=ObjectData.Angle;
1638    VolumeAngle=ObjectData.Angle;
1639    om=norm(VolumeAngle);%norm of rotation angle in radians
1640    OmAxis=VolumeAngle/om; %unit vector marking the rotation axis
1641    cos_om=cos(pi*om/180);
1642    sin_om=sin(pi*om/180);
1643    coeff=OmAxis(3)*(1-cos_om);
1644    %components of the unity vector norm_plane normal to the projection plane
1645    norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
1646    norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
1647    norm_plane(3)=OmAxis(3)*coeff+cos_om;
1648end
1649testangle=~isequal(VolumeAngle,[0 0 0]);
1650
1651%% mesh sizes DX, DY, DZ
1652DX=0;
1653DY=0; %default
1654DZ=0;
1655if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
1656     DX=abs(ObjectData.DX);%mesh of interpolation points
1657end
1658if isfield(ObjectData,'DY')&~isempty(ObjectData.DY)
1659     DY=abs(ObjectData.DY);%mesh of interpolation points
1660end
1661if isfield(ObjectData,'DZ')&~isempty(ObjectData.DZ)
1662     DZ=abs(ObjectData.DZ);%mesh of interpolation points
1663end
1664if  ~strcmp(ProjMode,'projection') && (DX==0||DY==0||DZ==0)
1665        errormsg='grid mesh DX , DY or DZ is missing';
1666        return
1667end
1668
1669%% extrema along each axis
1670testXMin=0;
1671testXMax=0;
1672testYMin=0;
1673testYMax=0;
1674if isfield(ObjectData,'RangeX')
1675        XMin=min(ObjectData.RangeX);
1676        XMax=max(ObjectData.RangeX);
1677        testXMin=XMax>XMin;
1678        testXMax=1;
1679end
1680if isfield(ObjectData,'RangeY')
1681        YMin=min(ObjectData.RangeY);
1682        YMax=max(ObjectData.RangeY);
1683        testYMin=YMax>YMin;
1684        testYMax=1;
1685end
1686width=0;%default width of the projection band
1687if isfield(ObjectData,'RangeZ')
1688        ZMin=min(ObjectData.RangeZ);
1689        ZMax=max(ObjectData.RangeZ);
1690        testZMin=ZMax>ZMin;
1691        testZMax=1;
1692end
1693
1694%% initiate Matlab  structure for physical field
1695[ProjData,errormsg]=proj_heading(FieldData,ObjectData);
1696if ~isempty(errormsg)
1697    return
1698end
1699
1700ProjData.NbDim=3;
1701ProjData.ListVarName={};
1702ProjData.VarDimName={};
1703if ~isequal(DX,0)&& ~isequal(DY,0)
1704    ProjData.CoordMesh=sqrt(DX*DY);%define typical data mesh, useful for mouse selection in plots
1705elseif isfield(FieldData,'CoordMesh')
1706    ProjData.CoordMesh=FieldData.CoordMesh;
1707end
1708
1709error=0;%default
1710flux=0;
1711testfalse=0;
1712ListIndex={};
1713
1714%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1715%% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
1716%-----------------------------------------------------------------
1717idimvar=0;
1718% LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
1719% CellVarIndex=cells of variable index arrays
1720ivar_new=0; % index of the current variable in the projected field
1721icoord=0;
1722nbcoord=0;%number of added coordinate variables brought by projection
1723nbvar=0;
1724for icell=1:length(CellVarIndex)
1725    NbDim=NbDimVec(icell);
1726    if NbDim<3
1727        continue
1728    end
1729    VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list FieldData.ListVarName
1730    VarType=VarTypeCell{icell};
1731    ivar_X=VarType.coord_x;
1732    ivar_Y=VarType.coord_y;
1733    ivar_Z=VarType.coord_z;
1734    ivar_U=VarType.vector_x;
1735    ivar_V=VarType.vector_y;
1736    ivar_W=VarType.vector_z;
1737    ivar_C=VarType.scalar ;
1738    ivar_Anc=VarType.ancillary;
1739    test_anc=zeros(size(VarIndex));
1740    test_anc(ivar_Anc)=ones(size(ivar_Anc));
1741    ivar_F=VarType.warnflag;
1742    ivar_FF=VarType.errorflag;
1743    check_unstructured_coord=~isempty(ivar_X) && ~isempty(ivar_Y);
1744    DimCell=FieldData.VarDimName{VarIndex(1)};
1745    if ischar(DimCell)
1746        DimCell={DimCell};%name of dimensions
1747    end
1748
1749%% case of input fields with unstructured coordinates
1750    if check_unstructured_coord
1751        XName=FieldData.ListVarName{ivar_X};
1752        YName=FieldData.ListVarName{ivar_Y};
1753        eval(['coord_x=FieldData.' XName ';'])
1754        eval(['coord_y=FieldData.' YName ';'])
1755        if length(ivar_Z)==1
1756            ZName=FieldData.ListVarName{ivar_Z};
1757            eval(['coord_z=FieldData.' ZName ';'])
1758        end
1759
1760        % translate  initial coordinates
1761        coord_x=coord_x-ObjectData.Coord(1,1);
1762        coord_y=coord_y-ObjectData.Coord(1,2);
1763        if ~isempty(ivar_Z)
1764            coord_z=coord_z-ObjectData.Coord(1,3);
1765        end
1766       
1767        % selection of the vectors in the projection range
1768%         if length(ivar_Z)==1 &&  width > 0
1769%             %components of the unitiy vector normal to the projection plane
1770%             fieldZ=NormVec_X*coord_x + NormVec_Y*coord_y+ NormVec_Z*coord_z;% distance to the plane           
1771%             indcut=find(abs(fieldZ) <= width);
1772%             for ivar=VarIndex
1773%                 VarName=FieldData.ListVarName{ivar};
1774%                 eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);']) 
1775%                     % A VOIR : CAS DE VAR STRUCTUREE MAIS PAS GRILLE REGULIERE : INTERPOLER SUR GRILLE REGULIERE             
1776%             end
1777%             coord_x=coord_x(indcut);
1778%             coord_y=coord_y(indcut);
1779%             coord_z=coord_z(indcut);
1780%         end
1781
1782       %rotate coordinates if needed: TODO modify
1783       if testangle
1784           coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
1785           coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
1786           if ~isempty(ivar_Z)
1787               coord_Y=coord_Y+coord_z *sin(Theta);
1788           end
1789           
1790           coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
1791           coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
1792           
1793       else
1794           coord_X=coord_x;
1795           coord_Y=coord_y;
1796           coord_Z=coord_z;
1797       end
1798        %restriction to the range of x and y if imposed
1799        testin=ones(size(coord_X)); %default
1800        testbound=0;
1801        if testXMin
1802            testin=testin & (coord_X >= XMin);
1803            testbound=1;
1804        end
1805        if testXMax
1806            testin=testin & (coord_X <= XMax);
1807            testbound=1;
1808        end
1809        if testYMin
1810            testin=testin & (coord_Y >= YMin);
1811            testbound=1;
1812        end
1813        if testYMax
1814            testin=testin & (coord_Y <= YMax);
1815            testbound=1;
1816        end
1817        if testbound
1818            indcut=find(testin);
1819            for ivar=VarIndex
1820                VarName=FieldData.ListVarName{ivar};
1821                eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])           
1822            end
1823            coord_X=coord_X(indcut);
1824            coord_Y=coord_Y(indcut);
1825            if length(ivar_Z)==1
1826                coord_Z=coord_Z(indcut);
1827            end
1828        end
1829        % different cases of projection
1830        if isequal(ObjectData.ProjMode,'projection')%%%%%%%   NOT USED %%%%%%%%%%
1831            for ivar=VarIndex %transfer variables to the projection plane
1832                VarName=FieldData.ListVarName{ivar};
1833                if ivar==ivar_X %x coordinate
1834                    eval(['ProjData.' VarName '=coord_X;'])
1835                elseif ivar==ivar_Y % y coordinate
1836                    eval(['ProjData.' VarName '=coord_Y;'])
1837                elseif isempty(ivar_Z) || ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced)
1838                    eval(['ProjData.' VarName '=FieldData.' VarName ';'])
1839                end
1840                if isempty(ivar_Z) || ivar~=ivar_Z
1841                    ProjData.ListVarName=[ProjData.ListVarName VarName];
1842                    ProjData.VarDimName=[ProjData.VarDimName DimCell];
1843                    nbvar=nbvar+1;
1844                    if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
1845                        ProjData.VarAttribute{nbvar}=FieldData.VarAttribute{ivar};
1846                    end
1847                end
1848            end 
1849        elseif isequal(ObjectData.ProjMode,'interp_lin')||isequal(ObjectData.ProjMode,'interp_tps')%interpolate data on a regular grid
1850            coord_x_proj=XMin:DX:XMax;
1851            coord_y_proj=YMin:DY:YMax;
1852            coord_z_proj=ZMin:DZ:ZMax;
1853            DimCell={'coord_z','coord_y','coord_x'};
1854            ProjData.ListVarName={'coord_z','coord_y','coord_x'};
1855            ProjData.VarDimName={'coord_z','coord_y','coord_x'};   
1856            nbcoord=2; 
1857            ProjData.coord_z=[ZMin ZMax];
1858            ProjData.coord_y=[YMin YMax];
1859            ProjData.coord_x=[XMin XMax];
1860            if isempty(ivar_X), ivar_X=0; end;
1861            if isempty(ivar_Y), ivar_Y=0; end;
1862            if isempty(ivar_Z), ivar_Z=0; end;
1863            if isempty(ivar_U), ivar_U=0; end;
1864            if isempty(ivar_V), ivar_V=0; end;
1865            if isempty(ivar_W), ivar_W=0; end;
1866            if isempty(ivar_F), ivar_F=0; end;
1867            if isempty(ivar_FF), ivar_FF=0; end;
1868            if ~isequal(ivar_FF,0)
1869                VarName_FF=FieldData.ListVarName{ivar_FF};
1870                eval(['indsel=find(FieldData.' VarName_FF '==0);'])
1871                coord_X=coord_X(indsel);
1872                coord_Y=coord_Y(indsel);
1873            end
1874            FF=zeros(1,length(coord_y_proj)*length(coord_x_proj));
1875            testFF=0;
1876            [X,Y,Z]=meshgrid(coord_y_proj,coord_z_proj,coord_x_proj);%grid in the new coordinates
1877            for ivar=VarIndex
1878                VarName=FieldData.ListVarName{ivar};
1879                if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1)                 
1880                    ivar_new=ivar_new+1;
1881                    ProjData.ListVarName=[ProjData.ListVarName {VarName}];
1882                    ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
1883                    if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
1884                        ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar};
1885                    end
1886                    if  ~isequal(ivar_FF,0)
1887                        eval(['FieldData.' VarName '=FieldData.' VarName '(indsel);'])
1888                    end
1889                    % linear interpolation
1890                    InterpFct=TriScatteredInterp(double(coord_X),double(coord_Y),double(coord_Z),double(FieldData.(VarName)));
1891                    ProjData.(VarName)=InterpFct(X,Y,Z);
1892%                     eval(['varline=reshape(ProjData.' VarName ',1,length(coord_y_proj)*length(coord_x_proj));'])
1893%                     FFlag= isnan(varline); %detect undefined values NaN
1894%                     indnan=find(FFlag);
1895%                     if ~isempty(indnan)
1896%                         varline(indnan)=zeros(size(indnan));
1897%                         eval(['ProjData.' VarName '=reshape(varline,length(coord_y_proj),length(coord_x_proj));'])
1898%                         FF(indnan)=ones(size(indnan));
1899%                         testFF=1;
1900%                     end
1901                    if ivar==ivar_U
1902                        ivar_U=ivar_new;
1903                    end
1904                    if ivar==ivar_V
1905                        ivar_V=ivar_new;
1906                    end
1907                    if ivar==ivar_W
1908                        ivar_W=ivar_new;
1909                    end
1910                end
1911            end
1912            if testFF
1913                ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
1914                ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
1915               ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
1916                ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
1917            end
1918        end
1919       
1920%% case of input fields defined on a structured  grid
1921    else
1922        VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
1923        eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions
1924        DimValue(DimValue==1)=[];%remove singleton dimensions       
1925        NbDim=numel(DimValue);%update number of space dimensions
1926        nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
1927        if NbDim>=3
1928            if NbDim>3
1929                errormsg='matrices with more than 3 dimensions not handled';
1930                return
1931            else
1932                if numel(find(VarType.coord))==2% the third matrix dimension does not correspond to a space coordinate
1933                    nbcolor=DimValue(3);
1934                    DimValue(3)=[]; %number of 'color' components updated
1935                    NbDim=2;% space dimension set to 2
1936                end
1937            end
1938        end
1939        AYName=FieldData.ListVarName{VarType.coord(NbDim-1)};%name of input x coordinate (name preserved on projection)
1940        AXName=FieldData.ListVarName{VarType.coord(NbDim)};%name of input y coordinate (name preserved on projection)   
1941        eval(['AX=FieldData.' AXName ';'])
1942        eval(['AY=FieldData.' AYName ';'])
1943        ListDimName=FieldData.VarDimName{VarIndex(1)};
1944        ProjData.ListVarName=[ProjData.ListVarName {AYName} {AXName}]; %TODO: check if it already exists in Projdata (several cells)
1945        ProjData.VarDimName=[ProjData.VarDimName {AYName} {AXName}];
1946
1947%         for idim=1:length(ListDimName)
1948%             DimName=ListDimName{idim};
1949%             if strcmp(DimName,'rgb')||strcmp(DimName,'nb_coord')||strcmp(DimName,'nb_coord_i')
1950%                nbcolor=DimValue(idim);
1951%                DimValue(idim)=[];
1952%             end
1953%             if isequal(DimName,'nb_coord_j')% NOTE: CASE OF TENSOR NOT TREATED
1954%                 DimValue(idim)=[];
1955%             end
1956%         end 
1957        Coord_z=[];
1958        Coord_y=[];
1959        Coord_x=[];   
1960
1961        for idim=1:NbDim %loop on space dimensions
1962            test_interp(idim)=0;%test for coordiate interpolation (non regular grid), =0 by default
1963            ivar=VarType.coord(idim);% index of the variable corresponding to the current dimension
1964            if ~isequal(ivar,0)%  a variable corresponds to the dimension #idim
1965                eval(['Coord{idim}=FieldData.' FieldData.ListVarName{ivar} ';']) ;% coord values for the input field
1966                if numel(Coord{idim})==2 %input array defined on a regular grid
1967                   DCoord_min(idim)=(Coord{idim}(2)-Coord{idim}(1))/DimValue(idim);
1968                else
1969                    DCoord=diff(Coord{idim});%array of coordinate derivatives for the input field
1970                    DCoord_min(idim)=min(DCoord);
1971                    DCoord_max=max(DCoord);
1972                %    test_direct(idim)=DCoord_max>0;% =1 for increasing values, 0 otherwise
1973                    if abs(DCoord_max-DCoord_min(idim))>abs(DCoord_max/1000)
1974                        msgbox_uvmat('ERROR',['non monotonic dimension variable # ' num2str(idim)  ' in proj_field.m'])
1975                                return
1976                    end               
1977                    test_interp(idim)=(DCoord_max-DCoord_min(idim))> 0.0001*abs(DCoord_max);% test grid regularity
1978                end
1979                test_direct(idim)=(DCoord_min(idim)>0);
1980            else  % no variable associated with the  dimension #idim, the coordinate value is set equal to the matrix index by default
1981                Coord_i_str=['Coord_' num2str(idim)];
1982                DCoord_min(idim)=1;%default
1983                Coord{idim}=[0.5 DimValue(idim)-0.5];
1984                test_direct(idim)=1;
1985            end
1986        end
1987        if DY==0
1988            DY=abs(DCoord_min(NbDim-1));
1989        end
1990        npY=1+round(abs(Coord{NbDim-1}(end)-Coord{NbDim-1}(1))/DY);%nbre of points after interpol
1991        if DX==0
1992            DX=abs(DCoord_min(NbDim));
1993        end
1994        npX=1+round(abs(Coord{NbDim}(end)-Coord{NbDim}(1))/DX);%nbre of points after interpol
1995        for idim=1:NbDim
1996            if test_interp(idim)
1997                DimValue(idim)=1+round(abs(Coord{idim}(end)-Coord{idim}(1))/abs(DCoord_min(idim)));%nbre of points after possible interpolation on a regular gri
1998            end
1999        end       
2000        Coord_y=linspace(Coord{NbDim-1}(1),Coord{NbDim-1}(end),npY);
2001        test_direct_y=test_direct(NbDim-1);
2002        Coord_x=linspace(Coord{NbDim}(1),Coord{NbDim}(end),npX);
2003        test_direct_x=test_direct(NbDim);
2004        DAX=DCoord_min(NbDim);
2005        DAY=DCoord_min(NbDim-1); 
2006        minAX=min(Coord_x);
2007        maxAX=max(Coord_x);
2008        minAY=min(Coord_y);
2009        maxAY=max(Coord_y);
2010        xcorner=[minAX maxAX minAX maxAX]-ObjectData.Coord(1,1);
2011        ycorner=[maxAY maxAY minAY minAY]-ObjectData.Coord(1,2);
2012        xcor_new=xcorner*cos(Phi)+ycorner*sin(Phi);%coord new frame
2013        ycor_new=-xcorner*sin(Phi)+ycorner*cos(Phi);
2014        if ~testXMax
2015            XMax=max(xcor_new);
2016        end
2017        if ~testXMin
2018            XMin=min(xcor_new);
2019        end
2020        if ~testYMax
2021            YMax=max(ycor_new);
2022        end
2023        if ~testYMin
2024            YMin=min(ycor_new);
2025        end
2026        DXinit=(maxAX-minAX)/(DimValue(NbDim)-1);
2027        DYinit=(maxAY-minAY)/(DimValue(NbDim-1)-1);
2028        if DX==0
2029            DX=DXinit;
2030        end
2031        if DY==0
2032            DY=DYinit;
2033        end
2034        if NbDim==3
2035            DZ=(Coord{1}(end)-Coord{1}(1))/(DimValue(1)-1);
2036            if ~test_direct(1)
2037                DZ=-DZ;
2038            end
2039            Coord_z=linspace(Coord{1}(1),Coord{1}(end),DimValue(1));
2040            test_direct_z=test_direct(1);
2041        end
2042        npX=floor((XMax-XMin)/DX+1);
2043        npY=floor((YMax-YMin)/DY+1);   
2044        if test_direct_y
2045            coord_y_proj=linspace(YMin,YMax,npY);%abscissa of the new pixels along the line
2046        else
2047            coord_y_proj=linspace(YMax,YMin,npY);%abscissa of the new pixels along the line
2048        end
2049        if test_direct_x
2050            coord_x_proj=linspace(XMin,XMax,npX);%abscissa of the new pixels along the line
2051        else
2052            coord_x_proj=linspace(XMax,XMin,npX);%abscissa of the new pixels along the line
2053        end
2054       
2055        % case with no rotation and interpolation
2056        if isequal(ProjMode,'projection') && isequal(Phi,0) && isequal(Theta,0) && isequal(Psi,0)
2057            if ~testXMin && ~testXMax && ~testYMin && ~testYMax && NbDim==2
2058                ProjData=FieldData;
2059            else
2060                indY=NbDim-1;
2061                if test_direct(indY)
2062                    min_indy=ceil((YMin-Coord{indY}(1))/DYinit)+1;
2063                    YIndexFirst=floor((YMax-Coord{indY}(1))/DYinit)+1;
2064                    Ybound(1)=Coord{indY}(1)+DYinit*(min_indy-1);
2065                    Ybound(2)=Coord{indY}(1)+DYinit*(YIndexFirst-1);
2066                else
2067                    min_indy=ceil((Coord{indY}(1)-YMax)/DYinit)+1;
2068                    max_indy=floor((Coord{indY}(1)-YMin)/DYinit)+1;
2069                    Ybound(2)=Coord{indY}(1)-DYinit*(max_indy-1);
2070                    Ybound(1)=Coord{indY}(1)-DYinit*(min_indy-1);
2071                end   
2072                if test_direct(NbDim)==1
2073                    min_indx=ceil((XMin-Coord{NbDim}(1))/DXinit)+1;
2074                    max_indx=floor((XMax-Coord{NbDim}(1))/DXinit)+1;
2075                    Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
2076                    Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
2077                else
2078                    min_indx=ceil((Coord{NbDim}(1)-XMax)/DXinit)+1;
2079                    max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1;
2080                    Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
2081                    Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
2082                end
2083                if NbDim==3
2084                    DimCell(1)=[]; %suppress z variable
2085                    DimValue(1)=[];
2086                                        %structured coordinates
2087                    if test_direct(1)
2088                        iz=ceil((ObjectData.Coord(1,3)-Coord{1}(1))/DZ)+1;
2089                    else
2090                        iz=ceil((Coord{1}(1)-ObjectData.Coord(1,3))/DZ)+1;
2091                    end
2092                end
2093                min_indy=max(min_indy,1);% deals with margin (bound lower than the first index)
2094                min_indx=max(min_indx,1);
2095                max_indy=min(max_indy,DimValue(1));
2096                max_indx=min(max_indx,DimValue(2));
2097                for ivar=VarIndex% loop on non coordinate variables
2098                    VarName=FieldData.ListVarName{ivar};
2099                    ProjData.ListVarName=[ProjData.ListVarName VarName];
2100                    ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
2101                    if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute)>=ivar
2102                        ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar};
2103                    end
2104                    if NbDim==3
2105                        eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']);
2106                    else
2107                        eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']);
2108                    end
2109                end 
2110                eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
2111                eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates
2112            end
2113        elseif isfield(FieldData,'A') %TO GENERALISE       % case with rotation and/or interpolation
2114            if NbDim==2 %2D case
2115                [X,Y]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
2116                XIMA=ObjectData.Coord(1,1)+(X)*cos(Phi)-Y*sin(Phi);%corresponding coordinates in the original image
2117                YIMA=ObjectData.Coord(1,2)+(X)*sin(Phi)+Y*cos(Phi);
2118                XIMA=(XIMA-minAX)/DXinit+1;% image index along x
2119                YIMA=(-YIMA+maxAY)/DYinit+1;% image index along y
2120                XIMA=reshape(round(XIMA),1,npX*npY);%indices reorganized in 'line'
2121                YIMA=reshape(round(YIMA),1,npX*npY);
2122                flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image
2123                if isequal(ObjectData.ProjMode,'interp_tps')
2124                    npx_interp_tps=ceil(abs(DX/DAX));
2125                    npy_interp_tps=ceil(abs(DY/DAY));
2126                    Minterp_tps=ones(npy_interp_tps,npx_interp_tps)/(npx_interp_tps*npy_interp_tps);
2127                    test_interp_tps=1;
2128                else
2129                    test_interp_tps=0;
2130                end
2131                eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates
2132                eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates
2133                for ivar=VarIndex
2134                    VarName=FieldData.ListVarName{ivar};
2135                    if test_interp(1) || test_interp(2)%interpolate on a regular grid       
2136                          eval(['ProjData.' VarName '=interp2(Coord{2},Coord{1},FieldData.' VarName ',Coord_x,Coord_y'');']) %TO TEST
2137                    end
2138                    %filter the field (image) if option 'interp_tps' is used
2139                    if test_interp_tps 
2140                         Aclass=class(FieldData.A);
2141                         ProjData.(VarName)=interp_tps2(Minterp_tps,FieldData.(VarName),'valid');
2142                         if ~isequal(Aclass,'double')
2143                             ProjData.(VarName)=Aclass(FieldData.(VarName));%revert to integer values
2144                         end
2145                    end
2146                    eval(['vec_A=reshape(FieldData.' VarName ',[],nbcolor);'])%put the original image in line             
2147                    %ind_in=find(flagin);
2148                    ind_out=find(~flagin);
2149                    ICOMB=(XIMA-1)*DimValue(1)+YIMA;
2150                    ICOMB=ICOMB(flagin);%index corresponding to XIMA and YIMA in the aligned original image vec_A
2151                    vec_B(flagin,1:nbcolor)=vec_A(ICOMB,:);
2152                    for icolor=1:nbcolor
2153                        vec_B(ind_out,icolor)=zeros(size(ind_out));
2154                    end
2155                    ProjData.ListVarName=[ProjData.ListVarName VarName];
2156                    ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
2157                    if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar
2158                        ProjData.VarAttribute{length(ProjData.ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
2159                    end     
2160                    eval(['ProjData.' VarName '=reshape(vec_B,npY,npX,nbcolor);']);
2161                end
2162                ProjData.FF=reshape(~flagin,npY,npX);%false flag A FAIRE: tenir compte d'un flga antï¿œrieur 
2163                ProjData.ListVarName=[ProjData.ListVarName 'FF'];
2164                ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
2165                ProjData.VarAttribute{length(ProjData.ListVarName)}.Role='errorflag';
2166            else %3D case
2167                if ~testangle     
2168                    % unstructured z coordinate
2169                    test_sup=(Coord{1}>=ObjectData.Coord(1,3));
2170                    iz_sup=find(test_sup);
2171                    iz=iz_sup(1);
2172                    if iz>=1 & iz<=npz
2173                        %ProjData.ListDimName=[ProjData.ListDimName ListDimName(2:end)];
2174                        %ProjData.DimValue=[ProjData.DimValue npY npX];
2175                        for ivar=VarIndex
2176                            VarName=FieldData.ListVarName{ivar};
2177                            ProjData.ListVarName=[ProjData.ListVarName VarName];
2178                            ProjData.VarAttribute{length(ProjData.ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes 
2179                            eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
2180                            %TODO : do a vertical average for a thick plane
2181                            if test_interp(2) || test_interp(3)
2182                                eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');'])
2183                            end
2184                        end
2185                    end
2186                else
2187                    errormsg='projection of structured coordinates on oblique plane not yet implemented';
2188                    %TODO: use interp3
2189                    return
2190                end
2191            end
2192        end
2193    end
2194
2195    %% projection of  velocity components in the rotated coordinates
2196    if testangle
2197        if isempty(ivar_V)
2198            msgbox_uvmat('ERROR','v velocity component missing in proj_field.m')
2199            return
2200        end
2201        UName=FieldData.ListVarName{ivar_U};
2202        VName=FieldData.ListVarName{ivar_V};   
2203        eval(['ProjData.' UName  '=cos(Phi)*ProjData.' UName '+ sin(Phi)*ProjData.' VName ';'])
2204        eval(['ProjData.' VName  '=cos(Theta)*(-sin(Phi)*ProjData.' UName '+ cos(Phi)*ProjData.' VName ');'])
2205        if ~isempty(ivar_W)
2206            WName=FieldData.ListVarName{ivar_W};
2207            eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])%
2208            eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']);
2209        end
2210        if ~isequal(Psi,0)
2211            eval(['ProjData.' UName '=cos(Psi)* ProjData.' UName '- sin(Psi)*ProjData.' VName ';']);
2212            eval(['ProjData.' VName '=sin(Psi)* ProjData.' UName '+ cos(Psi)*ProjData.' VName ';']);
2213        end
2214    end
2215end
2216
2217%------------------------------------------------------------------------
2218%--- transfer the global attributes
2219function [ProjData,errormsg]=proj_heading(FieldData,ObjectData)
2220%------------------------------------------------------------------------
2221ProjData=[];%default
2222errormsg='';%default
2223
2224%% transfer error
2225if isfield(FieldData,'Txt')
2226    errormsg=FieldData.Txt; %transmit erreur message
2227    return;
2228end
2229
2230%% transfer global attributes
2231if ~isfield(FieldData,'ListGlobalAttribute')
2232    ProjData.ListGlobalAttribute={};
2233else
2234    ProjData.ListGlobalAttribute=FieldData.ListGlobalAttribute;
2235end
2236for iattr=1:length(ProjData.ListGlobalAttribute)
2237    AttrName=ProjData.ListGlobalAttribute{iattr};
2238    if isfield(FieldData,AttrName)
2239        ProjData.(AttrName)=FieldData.(AttrName);
2240    end
2241end
2242
2243%% transfer coordinate unit
2244if isfield(ProjData,'CoordUnit')
2245    ProjData=rmfield(ProjData,'CoordUnit');% do not transfer by default (to avoid x/y=1 for profiles)
2246end
2247if isfield(FieldData,'CoordUnit')
2248    if isfield(ObjectData,'CoordUnit') && ~strcmp(FieldData.CoordUnit,ObjectData.CoordUnit)
2249        errormsg=[ObjectData.Type ' in ' ObjectData.CoordUnit ' coordinates, while field in ' FieldData.CoordUnit ];
2250        return
2251    elseif strcmp(ObjectData.Type,'plane')|| strcmp(ObjectData.Type,'volume')
2252         ProjData.CoordUnit=FieldData.CoordUnit;
2253    end
2254end
2255
2256%% store the properties of the projection object
2257ListObject={'Name','Type','ProjMode','angle','RangeX','RangeY','RangeZ','DX','DY','DZ','Coord'};
2258for ilist=1:length(ListObject)
2259    if isfield(ObjectData,ListObject{ilist})
2260        val=ObjectData.(ListObject{ilist});
2261        if ~isempty(val)
2262            ProjData.(['ProjObject' ListObject{ilist}])=val;
2263            ProjData.ListGlobalAttribute=[ProjData.ListGlobalAttribute {['ProjObject' ListObject{ilist}]}];
2264        end
2265    end   
2266end
2267
Note: See TracBrowser for help on using the repository browser.