Changeset 515 for trunk/src/proj_field.m


Ignore:
Timestamp:
Aug 15, 2012, 11:36:12 PM (12 years ago)
Author:
sommeria
Message:

improvement of calc-field and combination of two fields

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/proj_field.m

    r512 r515  
    8484ProjData=[];
    8585
    86 %% case of no projection (object is used only as graph display)
    87 if isfield(ObjectData,'ProjMode') && (isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside'))
     86%% in the absence of object Type or projection mode, or object coordinaes, the output is empty
     87if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode')
    8888    return
    8989end
    90 
    91 %% check coincidence of coordinate units
    92 if isfield(FieldData,'CoordUnit') && isfield(ObjectData,'CoordUnit')&&~strcmp(FieldData.CoordUnit,ObjectData.CoordUnit)
    93     errormsg='inconsistent coord units for field and projection object';
     90% case of no projection (object is used only as graph display)
     91if isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside')
    9492    return
    9593end
    96    
    97 %% in the absence of object Type or projection mode, or object coordinaes, the input field is just tranfered without change
    98 if ~isfield(ObjectData,'Type')||~isfield(ObjectData,'ProjMode')
    99     ProjData=FieldData;
    100     return
    101 end
    102 if ~isfield(ObjectData,'Coord')
     94if ~isfield(ObjectData,'Coord')||isempty(ObjectData.Coord)
    10395    if strcmp(ObjectData.Type,'plane')
    10496        ObjectData.Coord=[0 0 0];%default
    10597    else
    106         ProjData=FieldData;
    10798        return
    10899    end
    109100end
    110        
    111 %% OBSOLETE
    112 if isfield(ObjectData,'XMax') && ~isempty(ObjectData.XMax)
    113     ObjectData.RangeX(1)=ObjectData.XMax;
    114 end
    115 if isfield(ObjectData,'XMin') && ~isempty(ObjectData.XMin)
    116     ObjectData.RangeX(2)=ObjectData.XMin;
    117 end
    118 if isfield(ObjectData,'YMax') && ~isempty(ObjectData.YMax)
    119     ObjectData.RangeY(1)=ObjectData.YMax;
    120 end
    121 if isfield(ObjectData,'YMin') && ~isempty(ObjectData.YMin)
    122     ObjectData.RangeY(2)=ObjectData.YMin;
    123 end
    124 if isfield(ObjectData,'ZMax') && ~isempty(ObjectData.ZMax)
    125     ObjectData.RangeZ(1)=ObjectData.ZMax;
    126 end
    127 if isfield(ObjectData,'ZMin') && ~isempty(ObjectData.ZMin)
    128     ObjectData.RangeZ(2)=ObjectData.ZMin;
    129 end
    130 %%%%%%%%%%
    131101
    132102%% apply projection depending on the object type
    … …  
    401371    testproj(VarType.color)=1;
    402372    VarIndex=VarIndex(find(testproj(VarIndex)));%select only the projected variables
    403     if testX %case of unstructured coordinates
     373    if check_unstructured%case of unstructured coordinates
    404374         eval(['nbpoint=numel(FieldData.' FieldData.ListVarName{VarIndex(1)} ');'])
    405375         for ivar=[VarIndex VarType.coord_x VarType.coord_y VarType.errorflag]
    … …  
    547517ProjData.NbDim=1;
    548518%initialisation of the input parameters and defaultoutput
    549 ProjMode='projection';%direct projection on the line by default
    550 if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
     519ProjMode=ObjectData.ProjMode;
    551520% ProjAngle=90; %90 degrees projection by default
    552521
    … …  
    609578        continue
    610579    end
    611     testX=~isempty(VarType.coord_x) && ~isempty(VarType.coord_y);% test for unstructured coordinates
     580    check_unstructured=~isempty(VarType.coord_x) && ~isempty(VarType.coord_y);% test for unstructured coordinates
    612581    test_tps=~isempty(VarType.coord_tps);
    613582    testU=~isempty(VarType.vector_x) && ~isempty(VarType.vector_y);% test for vectors
    … …  
    631600    end   
    632601    % check needed object properties for unstructured positions (position given by the variables with role coord_x, coord_y
    633     if testX
     602    if check_unstructured
    634603        if  ~isequal(ProjMode,'interp')
    635604            if width==0
    … …  
    664633  %%%%%%%  % A FAIRE CALCULER MEAN DES QUANTITES    %%%%%%
    665634   %case of unstructured coordinates
    666     if testX   
     635    if check_unstructured   
    667636        for ip=1:siz_line(1)-1     %Loop on the segments of the polyline
    668637            linelength=sqrt(dlinx(ip)*dlinx(ip)+dliny(ip)*dliny(ip)); 
    … …  
    917886%-----------------------------------------------------------------
    918887
    919 %% initialisation of the input parameters of the projection plane
    920 ProjMode='projection';%direct projection by default
    921 if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
    922 
    923 %% axis origin
    924 if isempty(ObjectData.Coord)
    925     ObjectData.Coord(1,1)=0;%origin of the plane set to [0 0] by default
    926     ObjectData.Coord(1,2)=0;
    927     ObjectData.Coord(1,3)=0;
    928 end
    929 
    930888%% rotation angles
    931889PlaneAngle=[0 0 0];
    … …  
    951909
    952910%% mesh sizes DX and DY
    953 DX=0;
    954 DY=0; %default
     911DX=FieldData.Mesh;
     912DY=FieldData.Mesh; %default
    955913if isfield(ObjectData,'DX') && ~isempty(ObjectData.DX)
    956914     DX=abs(ObjectData.DX);%mesh of interpolation points
    … …  
    959917     DY=abs(ObjectData.DY);%mesh of interpolation points
    960918end
    961 if  ~strcmp(ProjMode,'projection') && (DX==0||DY==0)
     919if  ~strcmp(ObjectData.ProjMode,'projection') && (DX==0||DY==0)
    962920        errormsg='DX or DY missing';
    963921        display(errormsg)
    … …  
    970928testYMin=0;
    971929testYMax=0;
     930XMin=FieldData.XMin;%default
     931XMax=FieldData.XMax;%default
     932YMin=FieldData.YMin;%default
     933YMax=FieldData.YMax;%default
    972934if isfield(ObjectData,'RangeX')
    973935        XMin=min(ObjectData.RangeX);
    … …  
    1003965ListIndex={};
    1004966
    1005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    1006967%% group the variables (fields of 'FieldData') in cells of variables with the same dimensions
    1007 %-----------------------------------------------------------------
    1008 idimvar=0;
    1009 
     968% CellVarIndex=cells of variable index arrays
    1010969[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(FieldData);
    1011970if ~isempty(errormsg)
    … …  
    1014973end
    1015974
     975%% projection modes
     976check_grid=0;
     977ProjMode=cell(size(VarTypeCell));
     978for icell=1:numel(VarTypeCell)% TODO: recalculate coordinates here to get the bounds in the rotated coordinates
     979    ProjMode{icell}=ObjectData.ProjMode;
     980    if isfield(VarTypeCell{icell},'FieldRequest')
     981        switch VarTypeCell{icell}.FieldRequest
     982            case 'interp'
     983                ProjMode{icell}='interp';
     984            case 'derivatives'
     985                ProjMode{icell}='filter';
     986        end
     987    end
     988    if strcmp(ProjMode{icell},'interp')||strcmp(ProjMode{icell},'filter')
     989        check_grid=1;
     990    end
     991end
     992
     993%% define the new coordinates in case of interpolation on a grid
     994if check_grid% TODO: recalculate coordinates to get the bounds in the rotated coordinates
     995    ProjData.ListVarName={'coord_y','coord_x'};
     996    ProjData.VarDimName={'coord_y','coord_x'}; 
     997    ProjData.VarAttribute={[],[]};
     998    ProjData.coord_y=[YMin YMax];
     999    ProjData.coord_x=[XMin XMax];
     1000end
     1001
     1002%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    10161003% LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
    10171004% CellVarIndex=cells of variable index arrays
    … …  
    10201007nbcoord=0;%number of added coordinate variables brought by projection
    10211008nbvar=0;
     1009vector_x_proj=[];
     1010vector_y_proj=[];
    10221011for icell=1:length(CellVarIndex)
    10231012    NbDim=NbDimVec(icell);
    … …  
    10271016    VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list FieldData.ListVarName
    10281017    VarType=VarTypeCell{icell};
    1029 %     ivar_X=VarType.coord_x;
    1030 %     ivar_Y=VarType.coord_y;
    1031 %     ivar_Z=VarType.coord_z;
    10321018    ivar_U=VarType.vector_x;
    10331019    ivar_V=VarType.vector_y;
    … …  
    10371023    end
    10381024    ivar_W=VarType.vector_z;
    1039 %     ivar_Anc=VarType.ancillary;
    1040 %     test_anc=zeros(size(VarIndex));
    1041 %     test_anc(ivar_Anc)=ones(size(ivar_Anc));
    1042 %     ivar_F=VarType.warnflag;
    1043 %     ivar_FF=VarType.errorflag;
    10441025
    10451026    %type of coordinates
    … …  
    10581039    end
    10591040    coord_z=0;%default
    1060    
     1041
    10611042    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    10621043    switch CoordType
    … …  
    10641045        %% case of input fields with unstructured coordinates
    10651046        case 'unstructured'
    1066             if strcmp(ObjectData.ProjMode,'filter')
    1067                 continue %skip for filter (needs tps)
    1068             end
    1069             coord_x=FieldData.(FieldData.ListVarName{VarType.coord_x});
    1070             coord_y=FieldData.(FieldData.ListVarName{VarType.coord_y});
     1047            if strcmp(ProjMode{icell},'filter')
     1048                continue %skip for filter (needs tps field cell)
     1049            end
     1050            coord_x=FieldData.(FieldData.ListVarName{VarType.coord_x});% initial x coordinates
     1051            coord_y=FieldData.(FieldData.ListVarName{VarType.coord_y});% initial y coordinates
    10711052            if ~isempty(VarType.coord_z)
    1072                 coord_z=FieldData.(FieldData.ListVarName{VarType.coord_z});
     1053                coord_z=FieldData.(FieldData.ListVarName{VarType.coord_z});% initial z coordinates
    10731054            end
    10741055           
    1075             % translate  initial coordinates
     1056            % translate  initial coordinates to account for the new origin
    10761057            coord_x=coord_x-ObjectData.Coord(1,1);
    10771058            coord_y=coord_y-ObjectData.Coord(1,2);
    … …  
    10851066                fieldZ=norm_plane(1)*coord_x + norm_plane(2)*coord_y+ norm_plane(3)*coord_z;% distance to the plane
    10861067                indcut=find(abs(fieldZ) <= width);
    1087                 size(indcut)
    10881068                for ivar=VarIndex
    10891069                    VarName=FieldData.ListVarName{ivar};
    … …  
    10961076            end
    10971077           
    1098             %rotate coordinates if needed:
     1078            %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane
    10991079            Psi=PlaneAngle(1);
    11001080            Theta=PlaneAngle(2);
    … …  
    11041084                coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
    11051085                coord_Y=coord_Y+coord_z *sin(Theta);
    1106                 coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
    1107                
     1086                coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER               
    11081087                coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
    11091088            else
    … …  
    11121091            end
    11131092           
    1114             %restriction to the range of x and y if imposed
     1093            %restriction to the range of X and Y if imposed by the projection object
    11151094            testin=ones(size(coord_X)); %default
    11161095            testbound=0;
    … …  
    11391118                for ivar=VarIndex
    11401119                    VarName=FieldData.ListVarName{ivar};
    1141                     eval(['FieldData.' VarName '=FieldData.' VarName '(indcut);'])
     1120                    FieldData.(VarName)=FieldData.(VarName)(indcut);
    11421121                end
    11431122                coord_X=coord_X(indcut);
    … …  
    11491128           
    11501129            % different cases of projection
    1151             switch ObjectData.ProjMode
    1152                 case 'projection'
     1130            switch ProjMode{icell}
     1131                case 'projection' 
     1132                    nbvar=numel(ProjData.ListVarName);
    11531133                    for ivar=VarIndex %transfer variables to the projection plane
    11541134                        VarName=FieldData.ListVarName{ivar};
    … …  
    11721152                    coord_x_proj=XMin:DX:XMax;
    11731153                    coord_y_proj=YMin:DY:YMax;
    1174                     DimCell={'coord_y','coord_x'};
    1175                     ProjData.ListVarName={'coord_y','coord_x'};
    1176                     ProjData.VarDimName={'coord_y','coord_x'};
    1177                     nbcoord=2;
    1178                     ProjData.coord_y=[YMin YMax];
    1179                     ProjData.coord_x=[XMin XMax];
    1180 %                     if isempty(ivar_X), ivar_X=0; end;
    1181 %                     if isempty(ivar_Y), ivar_Y=0; end;
    1182 %                     if isempty(ivar_Z), ivar_Z=0; end;
    1183 %                     if isempty(ivar_U), ivar_U=0; end;
    1184 %                     if isempty(ivar_V), ivar_V=0; end;
    1185 %                     if isempty(ivar_W), ivar_W=0; end;
    1186 %                     if isempty(ivar_F), ivar_F=0; end;
    1187 %                     if isempty(ivar_FF), ivar_FF=0; end;
    1188                   %  if ~isequal(ivar_FF,0)
    1189                    if ~isempty(VarType.errorflag)
     1154                    [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);
     1155                    if ~isempty(VarType.errorflag)
    11901156                        VarName_FF=FieldData.ListVarName{VarType.errorflag};
    11911157                        indsel=find(FieldData.(VarName_FF)==0);
    … …  
    11931159                        coord_Y=coord_Y(indsel);
    11941160                    end
    1195                    
    1196                     FF=zeros(1,length(coord_y_proj)*length(coord_x_proj));
    11971161                    testFF=0;
    1198                     %for ivar=VarIndex % loop on field variables to project
    1199                     for ivar=[VarType.scalar VarType.vector_x VarType.vector_y]
    1200                         VarName=FieldData.ListVarName{ivar};
    1201 %                         if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF ||...
    1202 %                                 ( numel(test_anc)>=ivar && test_anc(ivar)==1))
    1203                             ivar_new=ivar_new+1;
    1204                             ProjData.ListVarName=[ProjData.ListVarName {VarName}];
    1205                             ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1206                             if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
    1207                                 ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar};
    1208                             end
    1209                             %if  ~isequal(ivar_FF,0)
    1210                             if ~isempty(VarType.errorflag)
    1211                                 FieldData.(VarName)=FieldData.(VarName)(indsel);
    1212                             end
    1213                             ProjData.(VarName)=griddata_uvmat(double(coord_X),double(coord_Y),double(FieldData.(VarName)),coord_x_proj,coord_y_proj');
    1214                             varline=reshape(ProjData.(VarName),1,length(coord_y_proj)*length(coord_x_proj));
    1215                             FFlag= isnan(varline); %detect undefined values NaN
    1216                             indnan=find(FFlag);
    1217                             if~isempty(indnan)
    1218                                 varline(indnan)=zeros(size(indnan));
    1219                                 ProjData.(VarName)=reshape(varline,length(coord_y_proj),length(coord_x_proj));
    1220                                 FF(indnan)=ones(size(indnan));
    1221                                 testFF=1;
    1222                             end
    1223                             if ivar==ivar_U
    1224                                 ivar_U=ivar_new;
    1225                             end
    1226                             if ivar==ivar_V
    1227                                 ivar_V=ivar_new;
    1228                             end
    1229                             if ivar==ivar_W
    1230                                 ivar_W=ivar_new;
    1231                             end
    1232 %                         end
    1233                     end
    1234                     if testFF
    1235                         ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
    1236                         ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
    1237                         ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1238                         ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
    1239                     end
    1240 %                     case 'filter'%interpolate data on a regular grid
    1241 %                         errormsg='tps required for filter option'
    1242                        
    1243             end
    1244            
    1245             %% case of tps interpolation (applies only in filter mode)
     1162                    nbvar=numel(ProjData.ListVarName);
     1163                    if isfield(VarType,'vector_x')&&isfield(VarType,'vector_y')&&~isempty(VarType.vector_x)
     1164                        VarName_x=FieldData.ListVarName{VarType.vector_x};
     1165                        VarName_y=FieldData.ListVarName{VarType.vector_y};
     1166                        if ~isempty(VarType.errorflag)
     1167                            FieldData.(VarName_x)=FieldData.(VarName_x)(indsel);
     1168                            FieldData.(VarName_y)=FieldData.(VarName_y)(indsel);
     1169                        end
     1170                        FieldVar=cat(2,FieldData.(VarName_x),FieldData.(VarName_y));
     1171                        if ~isfield(VarType,'CheckSub') || ~VarType.CheckSub
     1172                            vector_x_proj=numel(ProjData.ListVarName)+1;
     1173                            vector_y_proj=numel(ProjData.ListVarName)+2;
     1174                        end
     1175                    end
     1176                    if ~isempty(VarType.scalar)
     1177                        VarName_scalar=FieldData.ListVarName{VarType.scalar};
     1178                        if ~isempty(VarType.errorflag)
     1179                            FieldData.(VarName_scalar)=FieldData.(VarName_scalar)(indsel);
     1180                        end
     1181                        FieldVar=FieldData.(VarName_scalar);
     1182                    end
     1183                    [VarVal,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldVar,VarType.Operation,XI,YI);
     1184                    if isfield(VarType,'CheckSub') && VarType.CheckSub && ~isempty(vector_x_proj)
     1185                        ProjData.(ProjData.ListVarName{vector_x_proj})=ProjData.(ProjData.ListVarName{vector_x_proj})-VarVal{1};
     1186                        ProjData.(ProjData.ListVarName{vector_y_proj})=ProjData.(ProjData.ListVarName{vector_y_proj})-VarVal{2};
     1187                    else
     1188                        VarDimName=cell(size(ListFieldProj));
     1189                        for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
     1190                            ListFieldProj{ilist}=regexprep(ListFieldProj{ilist},'(.+','');
     1191                            if ~isempty(find(strcmp(ListFieldProj{ilist},ProjData.ListVarName)))
     1192                                ListFieldProj{ilist}=[ListFieldProj{ilist} '_1'];
     1193                            end                       
     1194                            ProjData.(ListFieldProj{ilist})=VarVal{ilist};
     1195                            VarDimName{ilist}={'coord_y','coord_x'};
     1196                        end
     1197                        ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
     1198                        ProjData.VarDimName=[ProjData.VarDimName VarDimName];
     1199                        ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
     1200                    end
     1201            end
     1202
     1203            %% case of tps interpolation (applies only in filter mode and for spatial derivatives)
    12461204        case 'tps'
    1247             if strcmp(ObjectData.ProjMode,'filter')
     1205            if strcmp(ProjMode{icell},'filter')
     1206                Coord=FieldData.(FieldData.ListVarName{VarType.coord_tps});
     1207                NbSites=FieldData.(FieldData.ListVarName{VarType.nbsites_tps});
     1208                SubRange=FieldData.(FieldData.ListVarName{VarType.subrange_tps});
     1209                if isfield(VarType,'vector_x_tps')&&isfield(VarType,'vector_y_tps')
     1210                    FieldVar=cat(3,FieldData.(FieldData.ListVarName{VarType.vector_x_tps}),FieldData.(FieldData.ListVarName{VarType.vector_y_tps}));
     1211                end
    12481212                coord_x_proj=XMin:DX:XMax;
    12491213                coord_y_proj=YMin:DY:YMax;
    … …  
    12531217                XI=XI+ObjectData.Coord(1,1);
    12541218                YI=YI+ObjectData.Coord(1,2);
    1255                 DataOut=calc_field(FieldData.FieldList,FieldData,cat(3,XI,YI));
    1256                 ProjData.ListVarName=[ProjData.ListVarName DataOut.ListVarName];
    1257                 ProjData.VarDimName=[ProjData.VarDimName DataOut.VarDimName];
    1258                 ProjData.VarAttribute=[ProjData.VarAttribute DataOut.VarAttribute];   
    1259                 for ilist=1:length(DataOut.ListVarName)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
    1260                     VarName=DataOut.ListVarName{ilist};
    1261 %                     ProjData.(VarName)=DataOut.(VarName);
    1262                     if ilist>=3
    1263                     ProjData.(VarName)=reshape(DataOut.(VarName),np_y,np_x);
    1264                     end
    1265                 end
    1266                 ProjData.coord_x=[XMin XMax];
    1267                 ProjData.coord_y=[YMin YMax];
     1219                [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord,NbSites,SubRange,FieldVar,VarType.Operation,cat(3,XI,YI));   
     1220                ListFieldProj=(fieldnames(DataOut))';
     1221                VarDimName=cell(size(ListFieldProj));
     1222                for ilist=1:numel(ListFieldProj)% reshape data, excluding coordinates (ilist=1-2), TODO: rationalise
     1223                    VarName=ListFieldProj{ilist};
     1224                    ProjData.(VarName)=DataOut.(VarName);
     1225                    VarDimName{ilist}={'coord_y','coord_x'};
     1226                end
     1227%                 if ~isfield(ProjData,'coord_x')
     1228%                 ProjData.coord_x=[XMin XMax];
     1229%                 ProjData.coord_y=[YMin YMax];
     1230%                 ListFieldProj=[{'coord_x','coord_y'} ListFieldProj'];
     1231%                 VarDimName=[{'coord_x','coord_y'} VarDimName];
     1232%                 VarAttribute=[{[],[]} VarAttribute];
     1233%                 end
     1234                ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
     1235                ProjData.VarDimName=[ProjData.VarDimName VarDimName];
     1236                ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
    12681237            end
    12691238           
    … …  
    12721241           
    12731242            VarName=FieldData.ListVarName{VarIndex(1)};%get the first variable of the cell to get the input matrix dimensions
    1274             eval(['DimValue=size(FieldData.' VarName ');'])%input matrix dimensions
     1243            DimValue=size(FieldData.(VarName));%input matrix dimensions
    12751244            DimValue(DimValue==1)=[];%remove singleton dimensions
    12761245            NbDim=numel(DimValue);%update number of space dimensions
    … …  
    14051374            end
    14061375            % case with no  interpolation
    1407             if isequal(ProjMode,'projection') && (~testangle || test90y || test90x)
     1376            if isequal(ProjMode{icell},'projection') && (~testangle || test90y || test90x)
    14081377                if  NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax
    1409                     ProjData=FieldData;% no change by projection
     1378                    ProjData.ListVarName=[ProjData.ListVarName FieldData.ListVarName(VarIndex)];
     1379                    ProjData.VarDimName=[ProjData.VarDimName FieldData.VarDimName(VarIndex)]; 
     1380                    ProjData.VarAttribute=[ProjData.VarAttribute FieldData.VarAttribute(VarIndex)];
     1381                    ProjData.(AYProjName)=FieldData.(AYName);
     1382                    ProjData.(AXProjName)=FieldData.(AXName);
     1383                    for ivar=VarIndex
     1384                        VarName=FieldData.ListVarName{ivar};
     1385                        ProjData.(VarName)=FieldData.(VarName);% no change by projection
     1386                    end
    14101387                else
    14111388                    indY=NbDim-1;
    … …  
    14891466                    YIMA=reshape(round(YIMA),1,npX*npY);
    14901467                    flagin=XIMA>=1 & XIMA<=DimValue(2) & YIMA >=1 & YIMA<=DimValue(1);%flagin=1 inside the original image
    1491                     if isequal(ObjectData.ProjMode,'filter')
     1468                    if isequal(ProjMode{icell},'filter')
    14921469                        npx_filter=ceil(abs(DX/DAX));
    14931470                        npy_filter=ceil(abs(DY/DAY));
    … …  
    15811558end
    15821559
     1560% %prepare substraction in case of two input fields
     1561% SubData.ListVarName={};
     1562% SubData.VarDimName={};
     1563% SubData.VarAttribute={};
     1564% check_remove=zeros(size(ProjData.ListVarName));
     1565% for iproj=1:numel(ProjData.VarAttribute)
     1566%     if isfield(ProjData.VarAttribute{iproj},'CheckSub')&&isequal(ProjData.VarAttribute{iproj}.CheckSub,1)
     1567%         VarName=ProjData.ListVarName{iproj};
     1568%         SubData.ListVarName=[SubData.ListVarName {VarName}];
     1569%         SubData.VarDimName=[SubData.VarDimName ProjData.VarDimName{iproj}];
     1570%         SubData.VarAttribute=[SubData.VarAttribute ProjData.VarAttribute{iproj}];
     1571%         SubData.(VarName)=ProjData.(VarName);
     1572%         check_remove(iproj)=1;       
     1573%     end
     1574% end
     1575% if ~isempty(find(check_remove))
     1576%     ind_remove=find(check_remove);
     1577%     ProjData.ListVarName(ind_remove)=[];
     1578%     ProjData.VarDimName(ind_remove)=[];
     1579%     ProjData.VarAttribute(ind_remove)=[];
     1580%     ProjData=sub_field(ProjData,[],SubData);
     1581% end   
     1582
    15831583%-----------------------------------------------------------------
    15841584%projection in a volume
    … …  
    15861586%-----------------------------------------------------------------
    15871587ProjData=FieldData;%default output
    1588 
    1589 %% initialisation of the input parameters of the projection plane
    1590 ProjMode='projection';%direct projection by default
    1591 if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
    15921588
    15931589%% axis origin
Note: See TracChangeset for help on using the changeset viewer.