Changeset 397 for trunk/src/proj_field.m
- Timestamp:
- Apr 26, 2012, 8:59:09 AM (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/src/proj_field.m
r388 r397 1 1 %'proj_field': projects the field on a projection object 2 2 %-------------------------------------------------------------------------- 3 % function [ProjData,errormsg]=proj_field(FieldData,ObjectData )3 % function [ProjData,errormsg]=proj_field(FieldData,ObjectData,FieldName) 4 4 % 5 5 % OUTPUT: … … 20 20 % .ProjMode=mode of projection ; 21 21 % .CoordUnit: 'px', 'cm' units for the coordinates defining the object 22 % . Phi angle of rotation (=0by default)22 % .Angle ( angles of rotation (=[0 0 0] by default) 23 23 % .ProjAngle=angle of projection; 24 24 % .DX,.DY,.DZ=increments along each coordinate … … 82 82 function [ProjData,errormsg]=proj_field(FieldData,ObjectData,FieldName) 83 83 errormsg='';%default 84 if ~exist('FieldName','var') 85 FieldName=''; 86 end 84 87 %% case of no projection (object is used only as graph display) 85 88 if 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 86 92 ProjData=[]; 93 end 87 94 return 88 95 end … … 123 130 %%%%%%%%%% 124 131 125 %% apply projection depending on the object style132 %% apply projection depending on the object type 126 133 switch ObjectData.Type 127 134 case 'points' … … 136 143 end 137 144 case 'plane' 138 [ProjData,errormsg] = proj_plane(FieldData,ObjectData );145 [ProjData,errormsg] = proj_plane(FieldData,ObjectData,FieldName); 139 146 case 'volume' 140 147 [ProjData,errormsg] = proj_volume(FieldData,ObjectData); … … 892 899 if isfield(ObjectData,'ProjMode'),ProjMode=ObjectData.ProjMode; end; 893 900 894 % axis origin901 %% axis origin 895 902 if isempty(ObjectData.Coord) 896 903 ObjectData.Coord(1,1)=0;%origin of the plane set to [0 0] by default … … 899 906 end 900 907 901 % rotation angles908 %% rotation angles 902 909 PlaneAngle=[0 0 0]; 903 910 norm_plane=[0 0 1]; … … 921 928 testangle=~isequal(PlaneAngle,[0 0 0]);% && ~test90y && ~test90x;%=1 for slanted plane 922 929 923 % mesh sizes DX and DY930 %% mesh sizes DX and DY 924 931 DX=0; 925 932 DY=0; %default … … 936 943 end 937 944 938 % extrema along each axis945 %% extrema along each axis 939 946 testXMin=0; 940 947 testXMax=0; … … 958 965 end 959 966 960 % initiate Matlab structure for physical field967 %% initiate Matlab structure for physical field 961 968 [ProjData,errormsg]=proj_heading(FieldData,ObjectData); 962 969 ProjData.NbDim=2; … … 968 975 ProjData.Mesh=FieldData.Mesh; 969 976 end 970 971 977 error=0;%default 972 978 flux=0; … … 1017 1023 1018 1024 %% case of input fields with unstructured coordinates 1025 coord_z=0;%default 1019 1026 if testX 1020 1027 XName=FieldData.ListVarName{ivar_X}; 1021 1028 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); 1024 1031 if length(ivar_Z)==1 1025 1032 ZName=FieldData.ListVarName{ivar_Z}; 1026 eval(['coord_z=FieldData.' ZName ';'])1033 coord_z=FieldData.(ZName); 1027 1034 end 1028 1035 … … 1051 1058 1052 1059 %rotate coordinates if needed: 1060 Psi=PlaneAngle(1); 1061 Theta=PlaneAngle(2); 1062 Phi=PlaneAngle(3); 1053 1063 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 % end1059 % else1060 1064 coord_X=(coord_x *cos(Phi) + coord_y* sin(Phi)); 1061 1065 coord_Y=(-coord_x *sin(Phi) + coord_y *cos(Phi))*cos(Theta); 1062 % end1063 % if ~isempty(ivar_Z)1064 1066 coord_Y=coord_Y+coord_z *sin(Theta); 1065 % end1066 % if testangle1067 1067 coord_X=(coord_X *cos(Psi) - coord_Y* sin(Psi));%A VERIFIER 1068 1068 1069 coord_Y=(coord_X *sin(Psi) + coord_Y* cos(Psi)); 1069 1070 else … … 1163 1164 FF=zeros(1,length(coord_y_proj)*length(coord_x_proj)); 1164 1165 testFF=0; 1166 FieldName 1165 1167 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'}]; 1173 1211 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 1213 1214 end 1214 1215 end
Note: See TracChangeset
for help on using the changeset viewer.