Changeset 966 for trunk/src/proj_field.m


Ignore:
Timestamp:
Jul 9, 2016, 10:02:57 PM (5 years ago)
Author:
sommeria
Message:

3D projection improved

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/proj_field.m

    r965 r966  
    552552    return
    553553end
    554 CellInfo=CellInfo(NbDim==2); %keep only the 2D cells
     554CellInfo=CellInfo(NbDim>=2); %keep only the 2D cells
    555555%%%%%% TODO: treat 1D fields: project as identity so that P o P=P for projection operation
    556556cell_select=true(size(CellInfo));
     
    676676    VarIndex=find(check_proj);% indices of the variables to be projected
    677677   
    678     %% identify vector components
    679     %testU=isfield(CellInfo{icell},'VarIndex_vector_x') &&isfield(CellInfo{icell},'VarIndex_vector_y') ;% test for vectors
    680     %     if testU
    681     %         UName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x};
    682     %         VName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y};
    683     %         vector_x=FieldData.(UName);
    684     %         vector_y=FieldData.(VName);
    685     %     end
    686678    %identify error flag
    687679    errorflag=0; %default, no error flag
     
    701693        %case of unstructured coordinates
    702694        case 'scattered'
    703 %             XName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
    704 %             YName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
     695            %             XName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
     696            %             YName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
    705697            coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});
    706698            coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)});
     
    813805        case 'grid'   %case of structured coordinates
    814806            if ~isequal(ObjectData.Type,'line')% exclude polyline
    815                 errormsg=['no  projection available on ' ObjectData.Type 'for structured coordinates']; %
     807                errormsg=['no  projection available on ' ObjectData.Type 'for structured coordinates'];
     808                return
     809            end%
     810            test_interp2=0;%default
     811            AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
     812            AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
     813            AX=FieldData.(AXName);% set of x positions
     814            AY=FieldData.(AYName);% set of y positions
     815            AName=FieldData.ListVarName{VarIndex(1)};
     816            npxy=size(FieldData.(AName));
     817            if max(NbDim)==3 % 3D case
     818               
    816819            else
    817                 test_Amat=1;%image or 2D matrix
    818                 test_interp2=0;%default
    819                 AYName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
    820                 AXName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
    821                 eval(['AX=FieldData.' AXName ';']);% set of x positions
    822                 eval(['AY=FieldData.' AYName ';']);% set of y positions
    823                 AName=FieldData.ListVarName{VarIndex(1)};
    824                 eval(['A=FieldData.' AName ';']);% scalar
    825                 npxy=size(A);
    826820                npx=npxy(2);
    827821                npy=npxy(1);
     
    921915            end
    922916    end
    923     if ~isempty(ivar_U) && ~isempty(ivar_V)
    924         vector_x =ProjData.(ProjData.ListVarName{ivar_U});
    925         ProjData.(ProjData.ListVarName{ivar_U}) =cos(theta)*vector_x+sin(theta)*ProjData.(ProjData.ListVarName{ivar_V});
    926         ProjData.(ProjData.ListVarName{ivar_V}) =-sin(theta)*vector_x+cos(theta)*ProjData.(ProjData.ListVarName{ivar_V});
    927     end
     917end
     918if ~isempty(ivar_U) && ~isempty(ivar_V)
     919    vector_x =ProjData.(ProjData.ListVarName{ivar_U});
     920    ProjData.(ProjData.ListVarName{ivar_U}) =cos(theta)*vector_x+sin(theta)*ProjData.(ProjData.ListVarName{ivar_V});
     921    ProjData.(ProjData.ListVarName{ivar_V}) =-sin(theta)*vector_x+cos(theta)*ProjData.(ProjData.ListVarName{ivar_V});
     922end
    928923end
    929924
     
    961956%     ObjectData.Angle=[0 0 0];
    962957%     ObjectData.Angle(1)=90*Delta_x/Delta_mod;
    963 %     ObjectData.Angle(2)=90*Delta_y/Delta_mod;
     958%     ObjectData.0(2)=90*Delta_y/Delta_mod;
    964959% end   
    965960if isfield(ObjectData,'Angle')&& isequal(size(ObjectData.Angle),[1 2])&& ~isequal(ObjectData.Angle,[0 0])
    966961    test90y=0;%isequal(ObjectData.Angle,[0 90 0]);
    967962    PlaneAngle=(pi/180)*ObjectData.Angle;
    968 %     om=norm(PlaneAngle);%norm of rotation angle in radians
    969 %     OmAxis=PlaneAngle/om; %unit vector marking the rotation axis
    970 %     cos_om=cos(om);
    971 %     sin_om=sin(om);
    972 %     coeff=OmAxis(3)*(1-cos_om);
    973 %     %components of the unity vector norm_plane normal to the projection plane
    974 %     norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
    975 %     norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
    976 %     norm_plane(3)=OmAxis(3)*coeff+cos_om;
     963    %     om=norm(PlaneAngle);%norm of rotation angle in radians
     964    %     OmAxis=PlaneAngle/om; %unit vector marking the rotation axis
     965    %     cos_om=cos(om);
     966    %     sin_om=sin(om);
     967    %     coeff=OmAxis(3)*(1-cos_om);
     968    %     %components of the unity vector norm_plane normal to the projection plane
     969    %     norm_plane(1)=OmAxis(1)*coeff+OmAxis(2)*sin_om;
     970    %     norm_plane(2)=OmAxis(2)*coeff-OmAxis(1)*sin_om;
     971    %     norm_plane(3)=OmAxis(3)*coeff+cos_om;
    977972   
    978 M2=[cos(PlaneAngle(2)) sin(PlaneAngle(2)) 0;-sin(PlaneAngle(2)) cos(PlaneAngle(2)) 0;0 0 1];
    979 M1=[1 0 0;0 cos(PlaneAngle(1)) sin(PlaneAngle(1));0 -sin(PlaneAngle(1)) cos(PlaneAngle(1))];
    980 M=M1*M2;
    981 norm_plane=M*[0 0 1]';
     973    M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1];
     974    M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))];
     975    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)
     976    norm_plane=M*[0 0 1]';
    982977   
    983978end
     
    16481643                    % determine the boundaries of the projected field,
    16491644                    % first find the 8 summits of the initial volume in the
    1650                     Angle=ObjectData.Angle*pi/180;
     1645                    PlaneAngle=ObjectData.Angle*pi/180;
    16511646                    % new coordinates
    16521647                    Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});%initial z coordinates
     
    16601655                    %Mrot_inv=rodrigues(-PlaneAngle);
    16611656                    newsummit=zeros(3,8);% initialize the rotated summit coordinates
    1662                     ObjectData.Coord=ObjectData.Coord';
    1663                     if size(ObjectData.Coord,2)<3
     1657                    ObjectData.Coord=ObjectData.Coord';% set ObjectData.Coord as a vertical vector
     1658                    if size(ObjectData.Coord,1)<3
    16641659                    ObjectData.Coord=[ObjectData.Coord; 0];%add z origin at z=0 by default
    16651660                    end
    1666                     M2=[cos(Angle(2)) sin(Angle(2)) 0;-sin(Angle(2)) cos(Angle(2)) 0;0 0 1];
    1667                     M1=[1 0 0;0 cos(Angle(1)) sin(Angle(1));0 -sin(Angle(1)) cos(Angle(1))];
    1668                     M=M1*M2;
     1661               
     1662                    M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1];
     1663                    M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))];
     1664                    M=M2*M1;
    16691665                    M_inv=inv(M);
    16701666                   
Note: See TracChangeset for help on using the changeset viewer.