Changeset 684 for trunk/src/proj_field.m


Ignore:
Timestamp:
Sep 8, 2013, 10:24:51 PM (11 years ago)
Author:
sommeria
Message:

corrections in proj_field: reducing the field of an image, and read_field: reading norm of vectors on gridded mesh

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/proj_field.m

    r675 r684  
    249249                if numel(Coord{idim})==2
    250250                    DCoord_min(idim)= (Coord{idim}(2)-Coord{idim}(1))/(npxy(idim)-1);
     251                    test_direct(idim)=DCoord_min(idim)>0;% =1 for increasing values, 0 otherwise
    251252                else
    252253                    DCoord=diff(Coord{idim});
    … …  
    924925    XMin=min(ObjectData.RangeX);
    925926    XMax=max(ObjectData.RangeX);
    926     testXMin=XMax>XMin;
     927    testXMin=XMax>XMin;%=1 if XMin defined (i.e. RangeY has tow distinct elements)
    927928    testXMax=1;% range restriction along X
    928929end
    … …  
    930931    YMin=min(ObjectData.RangeY);
    931932    YMax=max(ObjectData.RangeY);
    932     testYMin=YMax>YMin;
     933    testYMin=YMax>YMin;%=1 if YMin defined (i.e. RangeY has tow distinct elements)
    933934    testYMax=1;
    934935end
    … …  
    10411042        AYName='coord_y';
    10421043        AXName='coord_x';
    1043         if ~strcmp(ObjectData.ProjMode,'projection')
     1044        if strcmp(ObjectData.ProjMode,'projection')
     1045            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
     1046            ProjData.coord_x=[FieldData.XMin FieldData.XMax];
     1047            coord_x_proj=FieldData.XMin:FieldData.CoordMesh:FieldData.XMax;
     1048            coord_y_proj=FieldData.YMin:FieldData.CoordMesh:FieldData.YMax;
     1049        else
    10441050            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
    1045         ProjData.coord_x=[ObjectData.RangeX(1) ObjectData.RangeX(2)];
    1046         coord_x_proj=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2);
    1047         coord_y_proj=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2);
    1048         else
    1049         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
    1050         ProjData.coord_x=[FieldData.XMin FieldData.XMax];
    1051         coord_x_proj=FieldData.XMin:FieldData.CoordMesh:FieldData.XMax;
    1052         coord_y_proj=FieldData.YMin:FieldData.CoordMesh:FieldData.YMax;
    1053         end
    1054         [X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
    1055         XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system
    1056         YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3));
     1051            ProjData.coord_x=[ObjectData.RangeX(1) ObjectData.RangeX(2)];
     1052            coord_x_proj=ObjectData.RangeX(1):ObjectData.DX:ObjectData.RangeX(2);
     1053            coord_y_proj=ObjectData.RangeY(1):ObjectData.DY:ObjectData.RangeY(2);
     1054        end
     1055        [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
     1056%         XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system
     1057%         YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3));
    10571058    else% we use the existing grid from field cell #icell_grid
    10581059        NbDim=NbDimArray(icell_grid);
    … …  
    10621063        ProjData.(AXName)=FieldData.(AXName); % new (projected ) y coordinates
    10631064    end
    1064             ProjData.ListVarName={AYName,AXName};
    1065         ProjData.VarDimName={AYName,AXName};
    1066         ProjData.VarAttribute={[],[]};
     1065    ProjData.ListVarName={AYName,AXName};
     1066    ProjData.VarDimName={AYName,AXName};
     1067    ProjData.VarAttribute={[],[]};
    10671068end
    10681069   
    … …  
    12761277                ProjMode{icell}='interp_lin'; %request linear interpolation for projection on a tilted plane
    12771278            end
    1278 
     1279           
    12791280            if isequal(ProjMode{icell},'projection')% && (~testangle || test90y || test90x)
    12801281                if  NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax% no range restriction
    … …  
    12841285                        VarAttribute=[VarAttribute FieldData.VarAttribute(VarIndex)];
    12851286                    end
    1286                    
    12871287                    ProjData.(AYName)=FieldData.(AYName);
    12881288                    ProjData.(AXName)=FieldData.(AXName);
    … …  
    12921292                    end
    12931293                else
    1294                     indY=NbDim-1;
    1295                     if test_direct(indY)
    1296                         min_indy=ceil((YMin-Coord{indY}(1))/DYinit)+1;
    1297                         max_indy=floor((YMax-Coord{indY}(1))/DYinit)+1;
    1298                         Ybound(1)=Coord{indY}(1)+DYinit*(min_indy-1);
    1299                         Ybound(2)=Coord{indY}(1)+DYinit*(max_indy-1);
     1294                    Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});
     1295                    Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});
     1296                    if numel(Coord{NbDim-1})==2
     1297                        DY=(Coord{NbDim-1}(2)-Coord{NbDim-1}(1))/(DimValue(1)-1);
     1298                    end
     1299                    if numel(Coord{NbDim})==2
     1300                        DX=(Coord{NbDim}(2)-Coord{NbDim}(1))/(DimValue(2)-1);
     1301                    end
     1302                    if testYMin%test_direct(indY)
     1303                        YIndexMin=(YMin-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the min y value for the new field
     1304                        YIndexMax=(YMax-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the max y value for the new field
    13001305                    else
    1301                         min_indy=ceil((Coord{indY}(1)-YMax)/DYinit)+1;
    1302                         max_indy=floor((Coord{indY}(1)-YMin)/DYinit)+1;
    1303                         Ybound(2)=Coord{indY}(1)-DYinit*(max_indy-1);
    1304                         Ybound(1)=Coord{indY}(1)-DYinit*(min_indy-1);
    1305                     end
    1306                     if test_direct(NbDim)==1
    1307                         min_indx=ceil((XMin-Coord{NbDim}(1))/DXinit)+1;
    1308                         max_indx=floor((XMax-Coord{NbDim}(1))/DXinit)+1;
    1309                         Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
    1310                         Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
     1306                        YIndexMin=(Coord{NbDim-1}(1)-YMax)/DY+1;
     1307                        YIndexMax=(Coord{NbDim-1}(1)-YMin)/DY+1;
     1308                        Ybound(2)=Coord{NbDim-1}(1)-DY*(YIndexMax-1);
     1309                        Ybound(1)=Coord{NbDim-1}(1)-DY*(YIndexMin-1);
     1310                    end
     1311                    if testXMin%test_direct(NbDim)==1
     1312                        XIndexMin=(XMin-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the min x value for the new field
     1313                        XIndexMax=(XMax-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the max x value for the new field
     1314                        Xbound(1)=Coord{NbDim}(1)+DX*(XIndexMin-1);%  x value corresponding to XIndexMin
     1315                        Xbound(2)=Coord{NbDim}(1)+DX*(XIndexMax-1);%  x value corresponding to XIndexMax
    13111316                    else
    1312                         min_indx=ceil((Coord{NbDim}(1)-XMax)/DXinit)+1;
    1313                         max_indx=floor((Coord{NbDim}(1)-XMin)/DXinit)+1;
    1314                         Xbound(2)=Coord{NbDim}(1)+DXinit*(max_indx-1);
    1315                         Xbound(1)=Coord{NbDim}(1)+DXinit*(min_indx-1);
    1316                     end
    1317                     min_indy=max(min_indy,1);% deals with margin (bound lower than the first index)
    1318                     min_indx=max(min_indx,1);
     1317                        XIndexMin=(Coord{NbDim}(1)-XMax)/DX+1;
     1318                        XIndexMax=(Coord{NbDim}(1)-XMin)/DX+1;
     1319                        Xbound(2)=Coord{NbDim}(1)+DX*(XIndexMax-1);
     1320                        Xbound(1)=Coord{NbDim}(1)+DX*(XIndexMin-1);
     1321                    end
     1322                    YIndexRange(1)=ceil(min(YIndexMin,YIndexMax));%first y index to select from the previous field
     1323                    YIndexRange(1)=max(YIndexRange(1),1);% avoid bound lower than the first index
     1324                    YIndexRange(2)=floor(max(YIndexMin,YIndexMax));%last y index to select from the previous field
     1325                    YIndexRange(2)=min(YIndexRange(2),DimValue(NbDim-1));% limit to the last available index
     1326                    XIndexRange(1)=ceil(min(XIndexMin,XIndexMax));%first x index to select from the previous field
     1327                    XIndexRange(1)=max(XIndexRange(1),1);% avoid bound lower than the first index
     1328                    XIndexRange(2)=floor(max(XIndexMin,XIndexMax));%last x index to select from the previous field
     1329                    XIndexRange(2)=min(XIndexRange(2),DimValue(NbDim));% limit to the last available index
    13191330                    if test90y
    13201331                        ind_new=[3 2 1];
    … …  
    13291340                            eval(['ProjData.' VarName '=squeeze(ProjData.' VarName '(iz,:,:));'])% select the z index iz
    13301341                        end
    1331                         eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
    1332                         eval(['ProjData.' AXName '=[Coord{1}(end),Coord{1}(1)];']) %record the new (projected ) x coordinates
     1342                        ProjData.(AYName)=[Ybound(1) Ybound(2)]; %record the new (projected ) y coordinates
     1343                        ProjData.(AXName)=[Coord{1}(end),Coord{1}(1)]; %record the new (projected ) x coordinates
    13331344                    else
    13341345                        if NbDim==3
    … …  
    13411352                            end
    13421353                        end
    1343                         max_indy=min(max_indy,DimValue(1));%introduce bounds in y and x indices
    1344                         max_indx=min(max_indx,DimValue(2));
    13451354                        for ivar=VarIndex% loop on non coordinate variables
    13461355                            VarName=FieldData.ListVarName{ivar};
    … …  
    13511360                            end
    13521361                            if NbDim==3
    1353                                 eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,min_indy:max_indy,min_indx:max_indx));']);
     1362                                ProjData.(VarName)=squeeze(FieldData.(VarName)(iz,YIndexRange(1):YIndexRange(end),XIndexRange(1):XIndexRange(end)));
    13541363                            else
    1355                                 eval(['ProjData.' VarName '=FieldData.' VarName '(min_indy:max_indy,min_indx:max_indx,:);']);
     1364                                ProjData.(VarName)=FieldData.(VarName)(YIndexRange(1):YIndexRange(end),XIndexRange(1):XIndexRange(end),:);
    13561365                            end
    13571366                        end
    1358                         eval(['ProjData.' AYName '=[Ybound(1) Ybound(2)];']) %record the new (projected ) y coordinates
    1359                         eval(['ProjData.' AXName '=[Xbound(1) Xbound(2)];']) %record the new (projected ) x coordinates
     1367                        ProjData.(AYName)=Coord{NbDim-1}(1)+DY*(YIndexRange-1); %record the new (projected ) y coordinates
     1368                        ProjData.(AXName)=Coord{NbDim}(1)+DX*(XIndexRange-1); %record the new (projected ) x coordinates
    13601369                    end
    13611370                end
    … …  
    13701379                        test_interp_tps=0;
    13711380                    end
    1372                     coord_x_proj=XMin:DX:XMax;
    1373                     coord_y_proj=YMin:DY:YMax;
     1381                    Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});
     1382                    Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});
     1383                    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
     1384                    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
     1385                    xcor_new=xcorner*cos(PlaneAngle(3))+ycorner*sin(PlaneAngle(3));%coordinates of the corners in new frame
     1386                    ycor_new=-xcorner*sin(PlaneAngle(3))+ycorner*cos(PlaneAngle(3));
     1387                    if testXMin
     1388                        xcor_new=max(xcor_new,XMin);
     1389                    end
     1390                    if testXMax
     1391                        xcor_new=min(xcor_new,XMax);
     1392                    end
     1393                    if testYMin
     1394                        ycor_new=max(ycor_new,YMin);
     1395                    end
     1396                    if testYMax
     1397                        ycor_new=min(ycor_new,YMax);
     1398                    end
     1399                    coord_x_proj=min(xcor_new):DX:max(xcor_new);
     1400                    coord_y_proj=min(ycor_new):DY:max(ycor_new);
     1401                    ProjData.(AYName)=[coord_y_proj(1) coord_y_proj(end)]; %record the new (projected ) y coordinates
     1402                    ProjData.(AXName)=[coord_x_proj(1) coord_x_proj(end)]; %record the new (projected ) x coordinates
    13741403                    [X,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
    13751404                    XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system
    13761405                    YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3));
    1377                     Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});
    1378                     Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});
    13791406                    if numel(Coord{1})==2% x coordiante defiend by its bounds, get the whole set
    13801407                        Coord{1}=linspace(Coord{1}(1),Coord{1}(2),CellInfo{icell}.CoordSize(1));
    … …  
    13921419                        end
    13931420                        FFName=['FF_' num2str(ind)];% append an index to the name of error flag, FF_1,FF_2...
    1394                     end 
     1421                    end
    13951422                    % project all variables in the cell
    13961423                    for ivar=VarIndex
    … …  
    19471974                if test_direct(indY)
    19481975                    min_indy=ceil((YMin-Coord{indY}(1))/DYinit)+1;
    1949                     max_indy=floor((YMax-Coord{indY}(1))/DYinit)+1;
     1976                    YIndexFirst=floor((YMax-Coord{indY}(1))/DYinit)+1;
    19501977                    Ybound(1)=Coord{indY}(1)+DYinit*(min_indy-1);
    1951                     Ybound(2)=Coord{indY}(1)+DYinit*(max_indy-1);
     1978                    Ybound(2)=Coord{indY}(1)+DYinit*(YIndexFirst-1);
    19521979                else
    19531980                    min_indy=ceil((Coord{indY}(1)-YMax)/DYinit)+1;
Note: See TracChangeset for help on using the changeset viewer.