Ignore:
Timestamp:
Jan 17, 2020, 8:13:53 PM (3 years ago)
Author:
sommeria
Message:

LIF updated and bug corrections

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/proj_field.m

    r1071 r1072  
    999999    end
    10001000    PlaneAngle=(pi/180)*ObjectData.Angle;
     1001    if PlaneAngle==0
     1002        PlaneAngle=[0 0];
     1003    end
    10011004    %     om=norm(PlaneAngle);%norm of rotation angle in radians
    10021005    %     OmAxis=PlaneAngle/om; %unit vector marking the rotation axis
     
    10091012    %     norm_plane(3)=OmAxis(3)*coeff+cos_om;
    10101013   
    1011     M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1];
     1014    M1=[cos(PlaneAngle(1)) -sin(PlaneAngle(1)) 0;sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1];
    10121015    M=M1;
    10131016    if checkM2
    1014         M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))];
    1015         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)
     1017        M2=[1 0 0;0 cos(PlaneAngle(2)) -sin(PlaneAngle(2));0 sin(PlaneAngle(2)) cos(PlaneAngle(2))];
     1018        M=M1*M2;% first rotate in the x,y plane with angle PlaneAngle(1), then slant around the new x axis0 with angle PlaneAngle(2)
    10161019    end
    10171020    norm_plane=M*[0 0 1]';
     
    10381041end
    10391042InterpMesh=min(DX,DY);%mesh used for interpolation in a slanted plane
    1040 % if strcmp(ObjectData.Type,'plane_z')
    1041 %     InterpMesh=10*InterpMesh;%TODO: temporary, to shorten computation
    1042 % end
    10431043
    10441044%% extrema along each axis
     
    11261126% end
    11271127icell_grid=[];% field cell index which defines the grid
    1128 if ~strcmp(ObjectData.ProjMode,'projection')&& ~strcmp(ObjectData.Type,'plane_z')% TODO:rationalize
    1129     %% define the new coordinates in case of interpolation on a imposed grid
    1130     if ~testYMin
    1131         errormsg='min Y value not defined for the projection grid';return
    1132     end
    1133     if ~testYMax
    1134         errormsg='max Y value not defined for the projection grid';return
    1135     end
    1136     if ~testXMin
    1137         errormsg='min X value not defined for the projection grid';return
    1138     end
    1139     if ~testXMax
    1140         errormsg='max X value not defined for the projection grid';return
    1141     end
    1142 else
     1128if strcmp(ObjectData.ProjMode,'projection')
    11431129    %% case of a grid requested by the input field
    11441130    for icell=1:numel(CellInfo)% TODO: recalculate coordinates here to get the bounds in the rotated coordinates
     
    11621148        end
    11631149    end
    1164     if ~isempty(find(check_grid))% if a grid is requested by the input field
     1150    if ~isempty(find(check_grid,1))% if a grid is requested by the input field
    11651151        if isempty(icell_grid)%  if the grid is not given by cell #icell_grid
    11661152            if ~isfield(FieldData,'XMax')
     
    11691155        end
    11701156    end
    1171 end
    1172 if ~isempty(find(check_grid))||~strcmp(ObjectData.ProjMode,'projection')%no existing gridded data used
     1157else
     1158    %% define the new coordinates in case of interpolation on a imposed grid
     1159    if ~testYMin
     1160        errormsg='min Y value not defined for the projection grid';return% %%%%%%%%%%%%%%%%        A REVOIR
     1161
     1162    end
     1163    if ~testYMax
     1164        errormsg='max Y value not defined for the projection grid';return
     1165    end
     1166    if ~testXMin
     1167        errormsg='min X value not defined for the projection grid';return
     1168    end
     1169    if ~testXMax
     1170        errormsg='max X value not defined for the projection grid';return
     1171    end
     1172end
     1173if ~isempty(find(check_grid,1))||~strcmp(ObjectData.ProjMode,'projection')%no existing gridded data used
    11731174    if isempty(icell_grid)||~strcmp(ObjectData.ProjMode,'projection')%no existing gridded data used
    11741175        AYName='coord_y';
    11751176        AXName='coord_x';
    1176         if strcmp(ObjectData.ProjMode,'projection')||strcmp(ObjectData.Type,'plane_z')
     1177        if strcmp(ObjectData.ProjMode,'projection')||strcmp(ObjectData.Type,'plane')
    11771178            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
    11781179            ProjData.coord_x=[FieldData.XMin FieldData.XMax];
     
    14721473            end
    14731474           
    1474             if isequal(ProjMode{icell},'projection')% && (~testangle || test90y || test90x)
     1475            if strcmp(ProjMode{icell},'projection')% && (~testangle || test90y || test90x)
    14751476                if  NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax% no range restriction
    14761477                    ListVarName=[ListVarName FieldData.ListVarName(VarIndex)];
     
    16531654                        if isfield(FieldData,'VarAttribute')&&length(FieldData.VarAttribute)>=ivar
    16541655                            VarAttribute{length(ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
    1655                         end;
     1656                        end
    16561657                        ProjData.(FFName)=isnan(ProjData.(VarName));%detact NaN (points outside the interpolation range)
    16571658                        ProjData.(VarName)(ProjData.(FFName))=0; %set to 0 the NaN data
     
    16611662                    VarDimName=[VarDimName {DimCell}];
    16621663                    VarAttribute{numel(ListVarName)}.Role='errorflag';
    1663                 elseif ~testangle
     1664                elseif ~testangle % 3Dcase without change of angle
    16641665                    % unstructured z coordinate
    16651666                    test_sup=(Coord{1}>=ObjectData.Coord(1,3));
     
    16831684                    % determine the boundaries of the projected field,
    16841685                    % first find the 8 summits of the initial volume in the
    1685                     PlaneAngle=ObjectData.Angle*pi/180;
     1686                    PlaneAngle=[0 0 0];% default
     1687                    PlaneAngle(1:numel(ObjectData.Angle))=ObjectData.Angle*pi/180;
    16861688                    % new coordinates
    16871689                    Coord{1}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(1)});%initial z coordinates
    16881690                    Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});%initial y coordinates
    16891691                    Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});%initial x coordinates
    1690                     summit=zeros(3,8);% initialize summit coordinates
    1691                     summit(1,1:4)=[Coord{3}(1) Coord{3}(end) Coord{3}(1) Coord{3}(end)];%square
    1692                     summit(2,1:4)=[Coord{2}(1) Coord{2}(1) Coord{2}(end) Coord{2}(end)];% square at z= Coord{1}(1)
    1693                     summit(1:2,5:8)=summit(1:2,1:4);
    1694                     summit(3,:)=[Coord{1}(1)*ones(1,4) Coord{1}(end)*ones(1,4)];
    1695                     %Mrot_inv=rodrigues(-PlaneAngle);
    1696                     newsummit=zeros(3,8);% initialize the rotated summit coordinates
     1692%                     summit=zeros(3,8);% initialize summit coordinates
     1693%                     summit(1,1:4)=[Coord{3}(1) Coord{3}(end) Coord{3}(1) Coord{3}(end)];%square
     1694%                     summit(2,1:4)=[Coord{2}(1) Coord{2}(1) Coord{2}(end) Coord{2}(end)];% square at z= Coord{1}(1)
     1695%                     summit(1:2,5:8)=summit(1:2,1:4);
     1696%                     summit(3,:)=[Coord{1}(1)*ones(1,4) Coord{1}(end)*ones(1,4)];
     1697%                     %Mrot_inv=rodrigues(-PlaneAngle);
     1698%                     newsummit=zeros(3,8);% initialize the rotated summit coordinates
    16971699                    ObjectData.Coord=ObjectData.Coord';% set ObjectData.Coord as a vertical vector
    16981700                    if size(ObjectData.Coord,1)<3
     
    17001702                    end
    17011703                   
    1702                     M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1];
    1703                     M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))];
    1704                     M=M2*M1;
     1704%                     M1=[cos(PlaneAngle(1)) sin(PlaneAngle(1)) 0;-sin(PlaneAngle(1)) cos(PlaneAngle(1)) 0;0 0 1];
     1705%                     M2=[1 0 0;0 cos(PlaneAngle(2)) sin(PlaneAngle(2));0 -sin(PlaneAngle(2)) cos(PlaneAngle(2))];
     1706%                     M=M1*M2;
    17051707                    M_inv=inv(M);
    17061708                   
    1707                     for isummit=1:8% TODO: introduce a function for rotation of n points (to use also for scattered data)
    1708                         newsummit(:,isummit)=M*(summit(:,isummit)-(ObjectData.Coord));
    1709                     end
    1710                     coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:));% set of coordinqtes in the projection plane
    1711                     coord_y_proj=min(newsummit(2,:)):InterpMesh: max(newsummit(2,:));
    1712                     coord_z_proj=-width:width;
     1709%                     for isummit=1:8% TODO: introduce a function for rotation of n points (to use also for scattered data)
     1710%                         newsummit(:,isummit)=M*(summit(:,isummit)-(ObjectData.Coord));
     1711%                     end
     1712%                     coord_x_proj=min(newsummit(1,:)):InterpMesh: max(newsummit(1,:));% set of coordinqtes in the projection plane
     1713%                     coord_y_proj=min(newsummit(2,:)):InterpMesh: max(newsummit(2,:));
     1714%                     coord_z_proj=-width:width;
     1715                    coord_x_proj=ObjectData.RangeX(1):InterpMesh:ObjectData.RangeX(2);% set of coordinqtes in the projection plane
     1716                     coord_y_proj=ObjectData.RangeY(1):InterpMesh:ObjectData.RangeY(2);
     1717                    coord_z_proj=-floor(ObjectData.RangeInterp/InterpMesh):InterpMesh:floor(ObjectData.RangeInterp/InterpMesh);
    17131718                    %Mrot=rodrigues(PlaneAngle);% inverse rotation matrix
    1714                     Origin=M_inv*[coord_x_proj(1);coord_y_proj(1);coord_z_proj(1)]+ObjectData.Coord;
     1719                    %Origin=M_inv*[coord_x_proj(1);coord_y_proj(1);coord_z_proj(1)]+ObjectData.Coord;
    17151720                    npx=numel(coord_x_proj);
    17161721                    npy=numel(coord_y_proj);
     
    17221727                    iX=[coord_x_proj(end)-coord_x_proj(1);0;0]/(npx-1);
    17231728                    iY=[0;coord_y_proj(end)-coord_y_proj(1);0]/(npy-1);
     1729                    if npz==1
     1730                        iZ=[0;0;0];
     1731                    else
    17241732                    iZ=[0;0;coord_z_proj(end)-coord_z_proj(1)]/(npz-1);
     1733                    end
    17251734                    %                     iX(1:2)=[cosphi -sinphi;sinphi cosphi]*iX(1:2);
    17261735                    %                     iY(1:2)=[-cosphi -sinphi;sinphi cosphi]*iY(1:2);
     
    17311740                   
    17321741                    [Grid_x,Grid_y,Grid_z]=meshgrid(0:npx-1,0:npy-1,0:npz-1);
     1742                    %[Grid_x,Grid_y,Grid_z]=meshgrid(0:npz-1,0:npy-1,0:npx-1);
    17331743                    if ismatrix(Grid_x)% add a singleton in case of a single z value
    17341744                        Grid_x=shiftdim(Grid_x,-1);
     
    17361746                        Grid_z=shiftdim(Grid_z,-1);
    17371747                    end
    1738                     XI=Origin(1)+ix(1)*Grid_x+iy(1)*Grid_y+iz(1)*Grid_z;
    1739                     YI=Origin(2)+ix(2)*Grid_x+iy(2)*Grid_y+iz(2)*Grid_z;
    1740                     ZI=Origin(3)+ix(3)*Grid_x+iy(3)*Grid_y+iz(3)*Grid_z;
    1741                     [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1});% mesh in the initial coordinates
     1748                    XI=ObjectData.Coord(1)+ix(1)*Grid_x+iy(1)*Grid_y+iz(1)*Grid_z;
     1749                    YI=ObjectData.Coord(2)+ix(2)*Grid_x+iy(2)*Grid_y+iz(2)*Grid_z;
     1750                    ZI=ObjectData.Coord(3)+ix(3)*Grid_x+iy(3)*Grid_y+iz(3)*Grid_z;
     1751                  %  [X,Y,Z]=meshgrid(Coord{3},Coord{2},Coord{1});% grid of initial coordinates
    17421752                    for ivar=VarIndex
    17431753                        VarName=FieldData.ListVarName{ivar};
     
    17451755                        VarDimName=[VarDimName {{'coord_y','coord_x'}}];
    17461756                        VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
    1747                         FieldData.(VarName)=permute(FieldData.(VarName),[2 3 1]);
     1757                        FieldData.(VarName)=permute(FieldData.(VarName),[2 3 1]); %coordinate permutation needed to use interp3
     1758                        indexnan=isnan(FieldData.(VarName));
     1759                        FieldData.(VarName)(indexnan)=0;%set to zero the undefined values
    17481760                        ProjData.coord_x=coord_x_proj;
    17491761                        ProjData.coord_y=coord_y_proj;
    1750                         ProjData.(VarName)=interp3(X,Y,Z,double(FieldData.(VarName)),XI,YI,ZI,'*linear');
    1751                         ProjData.(VarName)=nanmean(ProjData.(VarName),3);
     1762                        ProjData.(VarName)=interp3(Coord{3},Coord{2},Coord{1},double(FieldData.(VarName)),ZI,YI,XI,'*linear');
     1763                        ProjData.(VarName)=nanmean(ProjData.(VarName),1);
    17521764                        ProjData.(VarName)=squeeze(ProjData.(VarName));
    17531765                    end
Note: See TracChangeset for help on using the changeset viewer.