Changeset 1121 for trunk/src/proj_field.m

Timestamp:

Dec 20, 2022, 2:00:18 PM (21 months ago)

Author:

sommeria

Message:

merge_proj_volume added

File:

• trunk/src/proj_field.m (modified) (16 diffs)

trunk/src/proj_field.m

@@ -1072 +1072 @@
 end
 
-%% initiate Matlab  structure for physical fieldopen
+%% initiate Matlab  structure for physical field
 [ProjData,errormsg]=proj_heading(FieldData,ObjectData);
 if ~isempty(errormsg)
@@ -1345 +1345 @@
             % two cases of projection for scattered coordinates
             switch ProjMode{icell}
-                case 'projection'
+                case 'projection'% ptoject scattered data without interpolation
                     nbvar=0;
                     %nbvar=numel(ProjData.ListVarName);
@@ -1366 +1366 @@
                         end
                     end
-                case 'interp_lin'%interpolate data on a regular grid
+                case 'interp_lin'%interpolate scattered data on a regular grid
                     if isfield(CellInfo{icell},'VarIndex_errorflag')
                         VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};
@@ -1382 +1382 @@
                     end
                     % interpolate and calculate field on the grid
-                    
+           
                     [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp([coord_X coord_Y],FieldData,CellInfo{icell}.FieldName,XI,YI);
                     
@@ -1500 +1500 @@
             end
             
-            if strcmp(ProjMode{icell},'projection')% && (~testangle || test90y || test90x)
+            if strcmp(ProjMode{icell},'projection')% project gridded data without interpolation
                 if  NbDim==2 && ~testXMin && ~testXMax && ~testYMin && ~testYMax% no range restriction
                     ListVarName=[ListVarName FieldData.ListVarName(VarIndex)];
@@ -1606 +1606 @@
                     end
                 end
-            else       % case with interpolation on a grid
+            else       % project gridded data with interpolation on a grid
                 if NbDim==2 %2D case
                     if isequal(ProjMode{icell},'interp_tps')
@@ -1778 +1778 @@
 
 %-----------------------------------------------------------------
-ProjData=FieldData;%default output
+ProjData=[];%default output
 
 %% axis origin
@@ -1809 +1809 @@
 DY=0; %default 
 DZ=0;
-if isfield(ObjectData,'DX')&~isempty(ObjectData.DX)
+ProjMode=ObjectData.ProjMode;
+if isfield(ObjectData,'DX')&&~isempty(ObjectData.DX)
      DX=abs(ObjectData.DX);%mesh of interpolation points 
 end
-if isfield(ObjectData,'DY')&~isempty(ObjectData.DY)
+if isfield(ObjectData,'DY')&&~isempty(ObjectData.DY)
      DY=abs(ObjectData.DY);%mesh of interpolation points 
 end
-if isfield(ObjectData,'DZ')&~isempty(ObjectData.DZ)
+if isfield(ObjectData,'DZ')&&~isempty(ObjectData.DZ)
      DZ=abs(ObjectData.DZ);%mesh of interpolation points 
 end
@@ -1872 +1873 @@
 %-----------------------------------------------------------------
 idimvar=0;
+
+[CellInfo,NbDimArray,errormsg]=find_field_cells(FieldData);
 % LOOP ON GROUPS OF VARIABLES SHARING THE SAME DIMENSIONS
 % CellVarIndex=cells of variable index arrays
@@ -1878 +1881 @@
 nbcoord=0;%number of added coordinate variables brought by projection
 nbvar=0;
-for icell=1:length(CellVarIndex)
-    NbDim=NbDimVec(icell);
+for icell=1:length(CellInfo)
+    NbDim=NbDimArray(icell);
     if NbDim<3
         continue
     end
-    VarIndex=CellVarIndex{icell};%  indices of the selected variables in the list FieldData.ListVarName
-    VarType=VarTypeCell{icell};
+    VarIndex=CellInfo{icell}.VarIndex;%  indices of the selected variables in the list FieldData.ListVarName
+    VarType=CellInfo{icell}.VarType;
     ivar_X=VarType.coord_x;
-    ivar_Y=VarType.coord_y;
+    ivar_Y=VarType.coorCeld_y;
     ivar_Z=VarType.coord_z;
     ivar_U=VarType.vector_x;
@@ -1897 +1900 @@
     ivar_F=VarType.warnflag;
     ivar_FF=VarType.errorflag;
-    check_unstructured_coord=~isempty(ivar_X) && ~isempty(ivar_Y);
+    %check_unstructured_coord=~isempty(ivar_X) && ~isempty(ivar_Y);
     DimCell=FieldData.VarDimName{VarIndex(1)};
     if ischar(DimCell)
@@ -1904 +1907 @@
 
 %% case of input fields with unstructured coordinates
-    if check_unstructured_coord
+    if strcmp(CellInfo{icell}.CoordType,'scattered')
         XName=FieldData.ListVarName{ivar_X};
         YName=FieldData.ListVarName{ivar_Y};
-        eval(['coord_x=FieldData.' XName ';'])
-        eval(['coord_y=FieldData.' YName ';'])
+        coord_x=FieldData.(XName);
+        coord_y=FieldData.(YName);
         if length(ivar_Z)==1
             ZName=FieldData.ListVarName{ivar_Z};
-            eval(['coord_z=FieldData.' ZName ';'])
-        end
-
+            coord_z=FieldData.(ZName);
+        end
         % translate  initial coordinates
         coord_x=coord_x-ObjectData.Coord(1,1);
@@ -1988 +1990 @@
                 VarName=FieldData.ListVarName{ivar};
                 if ivar==ivar_X %x coordinate
-                    eval(['ProjData.' VarName '=coord_X;'])
+                    ProjData.(VarName)=coord_X;
                 elseif ivar==ivar_Y % y coordinate
-                    eval(['ProjData.' VarName '=coord_Y;'])
+                    ProjData.(VarName)=coord_Y;
                 elseif isempty(ivar_Z) || ivar~=ivar_Z % other variables (except Z coordinate wyhich is not reproduced)
-                    eval(['ProjData.' VarName '=FieldData.' VarName ';'])
+                    ProjData.(VarName)=FieldData.(VarName);
                 end
                 if isempty(ivar_Z) || ivar~=ivar_Z 
@@ -2285 +2287 @@
                     test_interp_tps=0;
                 end
-                eval(['ProjData.' AYName '=[coord_y_proj(1) coord_y_proj(end)];']) %record the new (projected ) y coordinates
-                eval(['ProjData.' AXName '=[coord_x_proj(1) coord_x_proj(end)];']) %record the new (projected ) x coordinates
+                ProjData.(AYName)=[coord_y_proj(1) coord_y_proj(end)]; %record the new (projected ) y coordinates
+                ProjData.(AXName)=[coord_x_proj(1) coord_x_proj(end)]; %record the new (projected ) x coordinates
                 for ivar=VarIndex
                     VarName=FieldData.ListVarName{ivar};
@@ -2354 +2356 @@
     if testangle
         if isempty(ivar_V)
-            msgbox_uvmat('ERROR','v velocity component missing in proj_field.m')
+            errormsg='v velocity component missing in proj_field.m';
             return
         end
@@ -2372 +2374 @@
     end
 end
+if isempty(ProjData.ListVarName)
+    errormsg='input field is not 3D, no volume projection';
+    ProjData=[];
+end
 
 %------------------------------------------------------------------------