Changeset 397 for trunk/src/proj_field.m


Ignore:
Timestamp:
Apr 26, 2012, 8:59:09 AM (12 years ago)
Author:
sommeria
Message:

civ_matlab and patch improved, changes in the management of interpolation (still in progress).
adapatation to movies (use of VideoReader?)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/proj_field.m

    r388 r397  
    11%'proj_field': projects the field on a projection object
    22%--------------------------------------------------------------------------
    3 %  function [ProjData,errormsg]=proj_field(FieldData,ObjectData)
     3%  function [ProjData,errormsg]=proj_field(FieldData,ObjectData,FieldName)
    44%
    55% OUTPUT:
    … …  
    2020%    .ProjMode=mode of projection ;
    2121%    .CoordUnit: 'px', 'cm' units for the coordinates defining the object
    22 %    .Phi  angle of rotation (=0 by default)
     22%    .Angle (  angles of rotation (=[0 0 0] by default)
    2323%    .ProjAngle=angle of projection;
    2424%    .DX,.DY,.DZ=increments along each coordinate
    … …  
    8282function [ProjData,errormsg]=proj_field(FieldData,ObjectData,FieldName)
    8383errormsg='';%default
     84if ~exist('FieldName','var')
     85    FieldName='';
     86end
    8487%% case of no projection (object is used only as graph display)
    8588if isfield(ObjectData,'ProjMode') && (isequal(ObjectData.ProjMode,'none')||isequal(ObjectData.ProjMode,'mask_inside')||isequal(ObjectData.ProjMode,'mask_outside'))
     89    if ~isempty(FieldName)
     90        [ProjData,errormsg]=calc_field(FieldName,FieldData);
     91    else
    8692    ProjData=[];
     93    end
    8794    return
    8895end
    … …  
    123130%%%%%%%%%%
    124131
    125 %% apply projection depending on the object style
     132%% apply projection depending on the object type
    126133switch ObjectData.Type
    127134    case 'points'
    … …  
    136143        end
    137144    case 'plane'
    138             [ProjData,errormsg] = proj_plane(FieldData,ObjectData);
     145            [ProjData,errormsg] = proj_plane(FieldData,ObjectData,FieldName);
    139146    case 'volume'
    140147        [ProjData,errormsg] = proj_volume(FieldData,ObjectData);
    … …  
    892899if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end;
    893900
    894 %axis origin
     901%% axis origin
    895902if isempty(ObjectData.Coord)
    896903    ObjectData.Coord(1,1)=0;%origin of the plane set to [0 0] by default
    … …  
    899906end
    900907
    901 %rotation angles
     908%% rotation angles
    902909PlaneAngle=[0 0 0];
    903910norm_plane=[0 0 1];
    … …  
    921928testangle=~isequal(PlaneAngle,[0 0 0]);% && ~test90y && ~test90x;%=1 for slanted plane
    922929
    923 %mesh sizes DX and DY
     930%% mesh sizes DX and DY
    924931DX=0;
    925932DY=0; %default
    … …  
    936943end
    937944
    938 %extrema along each axis
     945%% extrema along each axis
    939946testXMin=0;
    940947testXMax=0;
    … …  
    958965end
    959966
    960 % initiate Matlab  structure for physical field
     967%% initiate Matlab  structure for physical field
    961968[ProjData,errormsg]=proj_heading(FieldData,ObjectData);
    962969ProjData.NbDim=2;
    … …  
    968975    ProjData.Mesh=FieldData.Mesh;
    969976end
    970 
    971977error=0;%default
    972978flux=0;
    … …  
    10171023   
    10181024    %% case of input fields with unstructured coordinates
     1025    coord_z=0;%default
    10191026    if testX
    10201027        XName=FieldData.ListVarName{ivar_X};
    10211028        YName=FieldData.ListVarName{ivar_Y};
    1022         eval(['coord_x=FieldData.' XName ';'])
    1023         eval(['coord_y=FieldData.' YName ';'])
     1029        coord_x=FieldData.(XName);
     1030        coord_y=FieldData.(YName);
    10241031        if length(ivar_Z)==1
    10251032            ZName=FieldData.ListVarName{ivar_Z};
    1026             eval(['coord_z=FieldData.' ZName ';'])
     1033            coord_z=FieldData.(ZName);
    10271034        end
    10281035       
    … …  
    10511058       
    10521059        %rotate coordinates if needed:
     1060        Psi=PlaneAngle(1);
     1061        Theta=PlaneAngle(2);
     1062        Phi=PlaneAngle(3);
    10531063        if testangle && ~test90y && ~test90x;%=1 for slanted plane
    1054             %             coord_X=coord_x;
    1055             %             coord_Y=coord_y;
    1056             %             if ~isequal(Theta,0)
    1057             %                 coord_Y=coord_Y *cos(Theta);
    1058             %             end
    1059             %         else
    10601064            coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi));
    10611065            coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta);
    1062             %         end
    1063             %         if ~isempty(ivar_Z)
    10641066            coord_Y=coord_Y+coord_z *sin(Theta);
    1065             %         end
    1066             %         if testangle
    10671067            coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER
     1068
    10681069            coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi));
    10691070        else
    … …  
    11631164            FF=zeros(1,length(coord_y_proj)*length(coord_x_proj));
    11641165            testFF=0;
     1166            FieldName
    11651167            if ~isempty(FieldName)
    1166                 FieldData=calc_field(FieldName,FieldData,XI,YI);
    1167             end
    1168             for ivar=VarIndex
    1169                 VarName=FieldData.ListVarName{ivar};
    1170                 if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1)
    1171                     ivar_new=ivar_new+1;
    1172                     ProjData.ListVarName=[ProjData.ListVarName {VarName}];
     1168                ProjData=calc_field(FieldName,FieldData,[XI YI]);
     1169            else
     1170                for ivar=VarIndex
     1171                    VarName=FieldData.ListVarName{ivar};
     1172                    if ~( ivar==ivar_X || ivar==ivar_Y || ivar==ivar_Z || ivar==ivar_F || ivar==ivar_FF || test_anc(ivar)==1)
     1173                        ivar_new=ivar_new+1;
     1174                        ProjData.ListVarName=[ProjData.ListVarName {VarName}];
     1175                        ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
     1176                        if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
     1177                            ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar};
     1178                        end
     1179                        if  ~isequal(ivar_FF,0)
     1180                            FieldData.(VarName)=FieldData.(VarName)(indsel);
     1181                        end
     1182                        %                     if isfield(FieldData,[VarName '_tps'])
     1183                        %                         [XI,YI]=meshgrid(coord_x_proj,coord_y_proj');
     1184                        %                         XI=reshape(XI,[],1);
     1185                        %                         YI=reshape(YI,[],1);
     1186                        %
     1187                        ProjData.(VarName)=griddata_uvmat(double(coord_X),double(coord_Y),double(FieldData.(VarName)),coord_x_proj,coord_y_proj',rho);
     1188                        varline=reshape(ProjData.(VarName),1,length(coord_y_proj)*length(coord_x_proj));
     1189                        FFlag= isnan(varline); %detect undefined values NaN
     1190                        indnan=find(FFlag);
     1191                        if~isempty(indnan)
     1192                            varline(indnan)=zeros(size(indnan));
     1193                            ProjData.(VarName)=reshape(varline,length(coord_y_proj),length(coord_x_proj));
     1194                            FF(indnan)=ones(size(indnan));
     1195                            testFF=1;
     1196                        end
     1197                        if ivar==ivar_U
     1198                            ivar_U=ivar_new;
     1199                        end
     1200                        if ivar==ivar_V
     1201                            ivar_V=ivar_new;
     1202                        end
     1203                        if ivar==ivar_W
     1204                            ivar_W=ivar_new;
     1205                        end
     1206                    end
     1207                end
     1208                if testFF
     1209                    ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
     1210                    ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
    11731211                    ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1174                     if isfield(FieldData,'VarAttribute') && length(FieldData.VarAttribute) >=ivar
    1175                         ProjData.VarAttribute{ivar_new+nbcoord}=FieldData.VarAttribute{ivar};
    1176                     end
    1177                     if  ~isequal(ivar_FF,0)
    1178                         FieldData.(VarName)=FieldData.(VarName)(indsel);
    1179                     end
    1180 %                     if isfield(FieldData,[VarName '_tps'])
    1181 %                         [XI,YI]=meshgrid(coord_x_proj,coord_y_proj');
    1182 %                         XI=reshape(XI,[],1);
    1183 %                         YI=reshape(YI,[],1);
    1184 %                         
    1185                    if ~isempty(FieldName)
    1186                         ProjData.(VarName)=griddata_uvmat(double(coord_X),double(coord_Y),double(FieldData.(VarName)),coord_x_proj,coord_y_proj',rho);
    1187                     end
    1188                     varline=reshape(ProjData.(VarName),1,length(coord_y_proj)*length(coord_x_proj));
    1189                     FFlag= isnan(varline); %detect undefined values NaN
    1190                     indnan=find(FFlag);
    1191                     if~isempty(indnan)
    1192                         varline(indnan)=zeros(size(indnan));
    1193                         ProjData.(VarName)=reshape(varline,length(coord_y_proj),length(coord_x_proj));
    1194                         FF(indnan)=ones(size(indnan));
    1195                         testFF=1;
    1196                     end
    1197                     if ivar==ivar_U
    1198                         ivar_U=ivar_new;
    1199                     end
    1200                     if ivar==ivar_V
    1201                         ivar_V=ivar_new;
    1202                     end
    1203                     if ivar==ivar_W
    1204                         ivar_W=ivar_new;
    1205                     end
    1206                 end
    1207             end
    1208             if testFF
    1209                 ProjData.FF=reshape(FF,length(coord_y_proj),length(coord_x_proj));
    1210                 ProjData.ListVarName=[ProjData.ListVarName {'FF'}];
    1211                 ProjData.VarDimName=[ProjData.VarDimName {DimCell}];
    1212                 ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
     1212                    ProjData.VarAttribute{ivar_new+1+nbcoord}.Role='errorflag';
     1213                end
    12131214            end
    12141215        end
Note: See TracChangeset for help on using the changeset viewer.