Changeset 961 for trunk/src/proj_field.m


Ignore:
Timestamp:
Jun 28, 2016, 9:05:33 AM (5 years ago)
Author:
sommeria
Message:

NomType? level introduced

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/proj_field.m

    r959 r961  
    16421642                    % determine the boundaries of the projected field,
    16431643                    % first find the 8 summits of the initial volume in the
     1644                    Angle=ObjectData.Angle*pi/180;
    16441645                    % new coordinates
    16451646                    Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});%initial z coordinates
     
    16511652                    summit(1:2,5:8)=summit(1:2,1:4);
    16521653                    summit(3,:)=[Coord{1}(1)*ones(1,4) Coord{1}(end)*ones(1,4)];
    1653                     Mrot=rodrigues(PlaneAngle);
     1654                    %Mrot_inv=rodrigues(-PlaneAngle);
    16541655                    newsummit=zeros(3,8);% initialize the rotated summit coordinates
    16551656                    ObjectData.Coord=ObjectData.Coord';
     
    16571658                    ObjectData.Coord=[ObjectData.Coord; 0];%add z origin at z=0 by default
    16581659                    end
     1660                    M2=[cos(Angle(2)) sin(Angle(2)) 0;-sin(Angle(2)) cos(Angle(2)) 0;0 0 1];
     1661                    M1=[1 0 0;0 cos(Angle(1)) sin(Angle(1));0 -sin(Angle(1)) cos(Angle(1))];
     1662                    M=M1*M2;
     1663                    M_inv=inv(M);
     1664                   
    16591665                    for isummit=1:8% TODO: introduce a function for rotation of n points (to use also for scattered data)
    1660                         newsummit(:,isummit)=Mrot*(summit(:,isummit)-(ObjectData.Coord));
     1666                        newsummit(:,isummit)=M*(summit(:,isummit)-(ObjectData.Coord));
    16611667                    end
    16621668                    coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:));% set of coordinqtes in the projection plane
    16631669                    coord_y_proj=min(newsummit(2,:)):InterpMesh: max(newsummit(2,:));
    16641670                    coord_z_proj=-width:width;
    1665                     Mrot_inverse=rodrigues(-PlaneAngle);% inverse rotation matrix
    1666                     Origin=Mrot_inverse*[coord_x_proj(1);coord_y_proj(1);coord_z_proj(1)]+ObjectData.Coord;
     1671                    %Mrot=rodrigues(PlaneAngle);% inverse rotation matrix
     1672                    Origin=M_inv*[coord_x_proj(1);coord_y_proj(1);coord_z_proj(1)]+ObjectData.Coord;
    16671673                    npx=numel(coord_x_proj);
    16681674                    npy=numel(coord_y_proj);
    16691675                    npz=numel(coord_z_proj);
    1670                     ix=Mrot_inverse*[coord_x_proj(end)-coord_x_proj(1);0;0]/(npx-1);% unit vector along the new x coordinates
    1671                     iy=Mrot_inverse*[0;coord_y_proj(end)-coord_y_proj(1);0]/(npy-1);% unit vector y coordinates
    1672                     iz=Mrot_inverse*[0;0;coord_z_proj(end)-coord_z_proj(1)]/(npz-1);% x unit vector z coordinates
     1676                   
     1677                    %modangle=sqrt(PlaneAngle(1)*PlaneAngle(1)+PlaneAngle(2)*PlaneAngle(2));
     1678%                     cosphi=PlaneAngle(1)/modangle;
     1679%                     sinphi=PlaneAngle(2)/modangle;
     1680                    iX=[coord_x_proj(end)-coord_x_proj(1);0;0]/(npx-1);
     1681                    iY=[0;coord_y_proj(end)-coord_y_proj(1);0]/(npy-1);
     1682                    iZ=[0;0;coord_z_proj(end)-coord_z_proj(1)]/(npz-1);
     1683%                     iX(1:2)=[cosphi -sinphi;sinphi cosphi]*iX(1:2);
     1684%                     iY(1:2)=[-cosphi -sinphi;sinphi cosphi]*iY(1:2);
     1685                   
     1686                    ix=M_inv*iX;%  vector along the new x coordinates transformed into old coordinates
     1687                    iy=M_inv*iY;% vector along y coordinates
     1688                    iz=M_inv*iZ;% vector along z coordinates
     1689
    16731690                    [Grid_x,Grid_y,Grid_z]=meshgrid(0:npx-1,0:npy-1,0:npz-1);
    16741691                    if ismatrix(Grid_x)% add a singleton in case of a single z value
     
    16801697                    YI=Origin(2)+ix(2)*Grid_x+iy(2)*Grid_y+iz(2)*Grid_z;
    16811698                    ZI=Origin(3)+ix(3)*Grid_x+iy(3)*Grid_y+iz(3)*Grid_z;
    1682                    [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1});
    1683 %                     X=permute(X,[3 1 2]);
    1684 %                     Y=permute(Y,[3 1 2]);
    1685 %                     Z=permute(Z,[3 1 2]);
     1699                   [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1});% mesh in the initial coordinates
    16861700                    for ivar=VarIndex
    16871701                            VarName=FieldData.ListVarName{ivar};
Note: See TracChangeset for help on using the changeset viewer.