Ignore:
Timestamp:
Dec 19, 2018, 10:16:08 AM (2 years ago)
Author:
sommeria
Message:

projection repaired

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/proj_field.m

    r1058 r1060  
    969969
    970970%-----------------------------------------------------------------
    971 %project on a plane 
     971%project on a plane
    972972% AJOUTER flux,circul,error
    973973function  [ProjData,errormsg] = proj_plane(FieldData, ObjectData)
     
    988988%     ObjectData.Angle(1)=90*Delta_x/Delta_mod;
    989989%     ObjectData.0(2)=90*Delta_y/Delta_mod;
    990 % end   
    991 if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 2])&& ~isequal(ObjectData.Angle,[0 0])
    992     test90y=0;%isequal(ObjectData.Angle,[0 90 0]);
     990% end
     991if isfield(ObjectData,'Angle')
     992    checkM2=0;
     993    if numel(ObjectData.Angle)==2 && ~isequal(ObjectData.Angle,[0 0])
     994        checkM2=1;
     995        test90y=0;%isequal(ObjectData.Angle,[0 90 0]);
     996    end
    993997    PlaneAngle=(pi/180)*ObjectData.Angle;
    994998    %     om=norm(PlaneAngle);%norm of rotation angle in radians
     
    10031007   
    10041008    M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1];
    1005     M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))];
    1006     M=M2*M1;% first rotate in the x,y plane with angle PlaneAngle(1), then slant around the new x axis0 with angle PlaneAngle(2)
     1009    M=M1;
     1010    if checkM2
     1011        M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))];
     1012        M=M2*M1;% first rotate in the x,y plane with angle PlaneAngle(1), then slant around the new x axis0 with angle PlaneAngle(2)
     1013    end
    10071014    norm_plane=M*[0 0 1]';
    10081015   
     
    10291036InterpMesh=min(DX,DY);%mesh used for interpolation in a slanted plane
    10301037% if strcmp(ObjectData.Type,'plane_z')
    1031 %     InterpMesh=10*InterpMesh;%TODO: temporary, to shorten computation 
     1038%     InterpMesh=10*InterpMesh;%TODO: temporary, to shorten computation
    10321039% end
    10331040
     
    11181125if ~strcmp(ObjectData.ProjMode,'projection')&& ~strcmp(ObjectData.Type,'plane_z')% TODO:rationalize
    11191126    %% define the new coordinates in case of interpolation on a imposed grid
    1120     if ~testYMin 
     1127    if ~testYMin
    11211128        errormsg='min Y value not defined for the projection grid';return
    11221129    end
     
    11771184        [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
    11781185        ProjData.VarDimName={AYName,AXName};
    1179 %         XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system
    1180 %         YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3));
     1186        %         XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system
     1187        %         YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3));
    11811188    else% we use the existing grid from field cell #icell_grid
    11821189        NbDim=NbDimArray(icell_grid);
     
    11851192        AYDimName=FieldData.VarDimName{CellInfo{icell_grid}.CoordIndex(NbDim-1)};%
    11861193        AXDimName=FieldData.VarDimName{CellInfo{icell_grid}.CoordIndex(NbDim)};%
    1187          ProjData.VarDimName={AYDimName,AXDimName};
     1194        ProjData.VarDimName={AYDimName,AXDimName};
    11881195        ProjData.(AYName)=FieldData.(AYName); % new (projected ) y coordinates
    11891196        ProjData.(AXName)=FieldData.(AXName); % new (projected ) y coordinates
     
    11941201    ProjData.VarAttribute{2}.Role='coord_x';
    11951202end
    1196    
     1203
    11971204%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    11981205% LOOP ON FIELD CELLS, PROJECT VARIABLES
     
    12561263           
    12571264            %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane
    1258             Psi=PlaneAngle(1);
    1259             Theta=PlaneAngle(2);
    1260            % Phi=PlaneAngle(3);
     1265            Phi=PlaneAngle(1);
     1266            %             Theta=PlaneAngle(2);
     1267            % Phi=PlaneAngle(3);
    12611268            if testangle && ~test90y && ~test90x;%=1 for slanted plane
    1262                 coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
    1263                 coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
    1264                 coord_Y=coord_Y+coord_z *sin(Theta);
    1265                 coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
    1266                 coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
     1269                coord_X=coord_x *cos(Phi) + coord_y* sin(Phi);
     1270                coord_Y=-coord_x *sin(Phi) + coord_y *cos(Phi);
     1271                if checkM2
     1272                    coord_Y=coord_Y*cos(PlaneAngle(2));
     1273                    coord_Y=coord_Y+coord_z *sin(PlaneAngle(2));
     1274                    coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
     1275                    coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
     1276                end
    12671277            else
    12681278                coord_X=coord_x;
     
    13771387                    end
    13781388                    if isfield (CellInfo{icell},'VarIndex_vector_x')&& isfield (CellInfo{icell},'VarIndex_vector_y')
    1379                     vector_x_proj=CellInfo{icell}.VarIndex_vector_x; %preserve for next cell
    1380                     vector_y_proj=CellInfo{icell}.VarIndex_vector_y; %preserve for next cell
     1389                        vector_x_proj=CellInfo{icell}.VarIndex_vector_x; %preserve for next cell
     1390                        vector_y_proj=CellInfo{icell}.VarIndex_vector_y; %preserve for next cell
    13811391                    end
    13821392            end
     
    13971407               
    13981408                %rotate coordinates if needed: coord_X,coord_Y= = coordinates in the new plane
    1399                 Psi=PlaneAngle(1);
    1400                 Theta=PlaneAngle(2);
    1401                % Phi=PlaneAngle(3);
     1409                Phi=PlaneAngle(1);
    14021410                if testangle && ~test90y && ~test90x;%=1 for slanted plane
    1403                     new_XI=XI*cos(Phi) - YI*sin(Phi)+ObjectData.Coord(1);
     1411                    new_XI=XI *cos(Phi) - YI* sin(Phi)+ObjectData.Coord(1);
    14041412                    YI=XI *sin(Phi) + YI *cos(Phi)+ObjectData.Coord(2);
    14051413                    XI=new_XI;
    1406                     %                 if checkUV
    1407                     %                     UValue=cos(Phi)*FieldVar(:,:,1)+ sin(Phi)*FieldVar(:,:,2);
    1408                     %                     FieldVar(:,:,2)=-sin(Phi)*FieldVar(:,:,1)+ cos(Phi)*FieldVar(:,:,2);
    1409                     %                     FieldVar(:,:,1)=UValue;
    1410                     %                 end
    14111414                end
    14121415               
     
    14451448            DimValue(DimValue==1)=[];%remove singleton dimensions
    14461449            NbDim=numel(DimValue);%update number of space dimensions
    1447            % nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
     1450            % nbcolor=1; %default number of 'color' components: third matrix index without corresponding coordinate
    14481451            if NbDim>=3
    14491452                if NbDim>3
     
    14921495                    end
    14931496                    if testYMax
    1494                          YIndexMax=(YMax-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the max y value for the new field
     1497                        YIndexMax=(YMax-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the max y value for the new field
    14951498                        if testYMin%test_direct(indY)
    14961499                            YIndexMin=(YMin-Coord{NbDim-1}(1))/DY+1;% matrix index corresponding to the min x value for the new field
    14971500                        else
    14981501                            YIndexMin=1;
    1499                         end         
     1502                        end
    15001503                    else
    15011504                        YIndexMax=numel(Coord{NbDim-1});
     
    15031506                    end
    15041507                    if testXMax
    1505                          XIndexMax=(XMax-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the max y value for the new field
     1508                        XIndexMax=(XMax-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the max y value for the new field
    15061509                        if testYMin%test_direct(indY)
    15071510                            XIndexMin=(XMin-Coord{NbDim}(1))/DX+1;% matrix index corresponding to the min x value for the new field
    15081511                        else
    15091512                            XIndexMin=1;
    1510                         end         
     1513                        end
    15111514                    else
    15121515                        XIndexMax=numel(Coord{NbDim});
     
    15611564                        end
    15621565                        if testXMax
    1563                          ProjData.(AXName)=Coord{NbDim}(1)+DX*(XIndexRange-1); %record the new (projected ) x coordinates
     1566                            ProjData.(AXName)=Coord{NbDim}(1)+DX*(XIndexRange-1); %record the new (projected ) x coordinates
    15641567                        else
    1565                           ProjData.(AXName)=FieldData.(AXName);
     1568                            ProjData.(AXName)=FieldData.(AXName);
    15661569                        end
    15671570                        if testYMax
    15681571                            ProjData.(AYName)=Coord{NbDim-1}(1)+DY*(YIndexRange-1); %record the new (projected ) x coordinates
    15691572                        else
    1570                           ProjData.(AYName)=FieldData.(AYName);
    1571                         end                           
     1573                            ProjData.(AYName)=FieldData.(AYName);
     1574                        end
    15721575                    end
    15731576                end
     
    16091612                    XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(2))-YI*sin(PlaneAngle(1));%corresponding coordinates in the original system
    16101613                    YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(2))+YI*cos(PlaneAngle(1));
    1611 
     1614                   
    16121615                    if numel(Coord{1})==2% x coordinate defined by its bounds, get the whole set
    16131616                        Coord{1}=linspace(Coord{1}(1),Coord{1}(2),CellInfo{icell}.CoordSize(1));
     
    16341637                            ProjData.(VarName)=interp2(X,Y,double(FieldData.(VarName)(:,:,1)),XI,YI,'*linear');
    16351638                            for icolor=2:size(FieldData.(VarName),3)% project 'color' components
    1636                                 ProjData.(VarName)=cat(3,ProjData.(VarName),interp2(X,Y,double(FieldData.(VarName)(:,:,icolor)),XI,YI,'*linear')); 
     1639                                ProjData.(VarName)=cat(3,ProjData.(VarName),interp2(X,Y,double(FieldData.(VarName)(:,:,icolor)),XI,YI,'*linear'));
    16371640                            end
    16381641                        end
     
    16741677                        end
    16751678                    end
    1676                 else   %projection of structured coordinates on oblique plane 
     1679                else   %projection of structured coordinates on oblique plane
    16771680                    % determine the boundaries of the projected field,
    16781681                    % first find the 8 summits of the initial volume in the
     
    16831686                    Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});%initial x coordinates
    16841687                    summit=zeros(3,8);% initialize summit coordinates
    1685                     summit(1,1:4)=[Coord{3}(1) Coord{3}(end) Coord{3}(1) Coord{3}(end)];%square 
     1688                    summit(1,1:4)=[Coord{3}(1) Coord{3}(end) Coord{3}(1) Coord{3}(end)];%square
    16861689                    summit(2,1:4)=[Coord{2}(1) Coord{2}(1) Coord{2}(end) Coord{2}(end)];% square at z= Coord{1}(1)
    16871690                    summit(1:2,5:8)=summit(1:2,1:4);
     
    16911694                    ObjectData.Coord=ObjectData.Coord';% set ObjectData.Coord as a vertical vector
    16921695                    if size(ObjectData.Coord,1)<3
    1693                     ObjectData.Coord=[ObjectData.Coord; 0];%add z origin at z=0 by default
    1694                     end
    1695                
     1696                        ObjectData.Coord=[ObjectData.Coord; 0];%add z origin at z=0 by default
     1697                    end
     1698                   
    16961699                    M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1];
    16971700                    M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))];
     
    17121715                   
    17131716                    %modangle=sqrt(PlaneAngle(1)*PlaneAngle(1)+PlaneAngle(2)*PlaneAngle(2));
    1714 %                     cosphi=PlaneAngle(1)/modangle;
    1715 %                     sinphi=PlaneAngle(2)/modangle;
     1717                    %                     cosphi=PlaneAngle(1)/modangle;
     1718                    %                     sinphi=PlaneAngle(2)/modangle;
    17161719                    iX=[coord_x_proj(end)-coord_x_proj(1);0;0]/(npx-1);
    17171720                    iY=[0;coord_y_proj(end)-coord_y_proj(1);0]/(npy-1);
    17181721                    iZ=[0;0;coord_z_proj(end)-coord_z_proj(1)]/(npz-1);
    1719 %                     iX(1:2)=[cosphi -sinphi;sinphi cosphi]*iX(1:2);
    1720 %                     iY(1:2)=[-cosphi -sinphi;sinphi cosphi]*iY(1:2);
     1722                    %                     iX(1:2)=[cosphi -sinphi;sinphi cosphi]*iX(1:2);
     1723                    %                     iY(1:2)=[-cosphi -sinphi;sinphi cosphi]*iY(1:2);
    17211724                   
    17221725                    ix=M_inv*iX;%  vector along the new x coordinates transformed into old coordinates
    17231726                    iy=M_inv*iY;% vector along y coordinates
    17241727                    iz=M_inv*iZ;% vector along z coordinates
    1725 
     1728                   
    17261729                    [Grid_x,Grid_y,Grid_z]=meshgrid(0:npx-1,0:npy-1,0:npz-1);
    17271730                    if ismatrix(Grid_x)% add a singleton in case of a single z value
     
    17331736                    YI=Origin(2)+ix(2)*Grid_x+iy(2)*Grid_y+iz(2)*Grid_z;
    17341737                    ZI=Origin(3)+ix(3)*Grid_x+iy(3)*Grid_y+iz(3)*Grid_z;
    1735                    [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1});% mesh in the initial coordinates
     1738                    [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1});% mesh in the initial coordinates
    17361739                    for ivar=VarIndex
    1737                             VarName=FieldData.ListVarName{ivar};
    1738                             ListVarName=[ListVarName VarName];
    1739                             VarDimName=[VarDimName {{'coord_y','coord_x'}}];
    1740                             VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
    1741                             FieldData.(VarName)=permute(FieldData.(VarName),[2 3 1]);
    1742                             ProjData.coord_x=coord_x_proj;
    1743                             ProjData.coord_y=coord_y_proj;
    1744                             ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI,'*linear');
    1745                             ProjData.(VarName)=nanmean(ProjData.(VarName),3);
    1746                             ProjData.(VarName)=squeeze(ProjData.(VarName));
     1740                        VarName=FieldData.ListVarName{ivar};
     1741                        ListVarName=[ListVarName VarName];
     1742                        VarDimName=[VarDimName {{'coord_y','coord_x'}}];
     1743                        VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
     1744                        FieldData.(VarName)=permute(FieldData.(VarName),[2 3 1]);
     1745                        ProjData.coord_x=coord_x_proj;
     1746                        ProjData.coord_y=coord_y_proj;
     1747                        ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI,'*linear');
     1748                        ProjData.(VarName)=nanmean(ProjData.(VarName),3);
     1749                        ProjData.(VarName)=squeeze(ProjData.(VarName));
    17471750                    end
    17481751                end
     
    17751778                UName=ListVarName{ivar_U};
    17761779                VName=ListVarName{ivar_V};
    1777                 UValue=cos(PlaneAngle(3))*ProjData.(UName)+ sin(PlaneAngle(3))*ProjData.(VName);
    1778                 ProjData.(VName)=(-sin(PlaneAngle(3))*ProjData.(UName)+ cos(PlaneAngle(3))*ProjData.(VName));
     1780                if checkM2
     1781                    UValue=cos(PlaneAngle(1))*ProjData.(UName)+ sin(PlaneAngle(1))*ProjData.(VName);
     1782                    ProjData.(VName)=(-sin(PlaneAngle(1))*ProjData.(UName)+ cos(PlaneAngle(1))*ProjData.(VName));
     1783                else
     1784                    UValue=cos(PlaneAngle(1))*ProjData.(UName)+ sin(PlaneAngle(1))*ProjData.(VName);
     1785                    ProjData.(VName)=(-sin(PlaneAngle(1))*ProjData.(UName)+ cos(PlaneAngle(1))*ProjData.(VName));
     1786                end
    17791787                ProjData.(UName)=UValue;
    17801788                if ~isempty(ivar_W)
Note: See TracChangeset for help on using the changeset viewer.