Changeset 866 for trunk/src/proj_field.m


Ignore:
Timestamp:
Feb 8, 2015, 7:26:02 AM (10 years ago)
Author:
sommeria
Message:

stereo_civ_updated

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/proj_field.m

    r864 r866  
    120120end
    121121
     122% %% take the difference of projected input fields if relevant
     123% [CellInfo,NbDim,errormsg]=find_field_cells(ProjData);
     124% ind_remove=zeros(size(ProjData.ListVarName));
     125% ivar=[];
     126% ivar_1=[];
     127% for icell=1:numel(CellInfo)
     128%     if ~isempty(CellInfo{icell})
     129%         % if two scalar are in the same cell
     130%         if isfield(CellInfo{icell},'VarIndex_scalar') && numel(CellInfo{icell}.VarIndex_scalar)==2 && ProjData.VarAttribute{CellInfo{icell}.VarIndex_scalar(2)}.CheckSub;
     131%             ivar=[ivar CellInfo{icell}.VarIndex_scalar(1)];
     132%             ivar_1=[ivar_1 CellInfo{icell}.VarIndex_scalar(2)];
     133%         end
     134%         % if two vector u components are in the same cell
     135%         if isfield(CellInfo{icell},'VarIndex_vector_x') && numel(CellInfo{icell}.VarIndex_vector_x)==2 && ProjData.VarAttribute{CellInfo{icell}.VarIndex_vector_x(2)}.CheckSub;
     136%             ivar=[ivar CellInfo{icell}.VarIndex_vector_x(1)];
     137%             ivar_1=[ivar_1 CellInfo{icell}.VarIndex_vector_x(2)];
     138%         end
     139%          % if two vector v components are in the same cell
     140%         if isfield(CellInfo{icell},'VarIndex_vector_y') && numel(CellInfo{icell}.VarIndex_vector_y)==2 && ProjData.VarAttribute{CellInfo{icell}.VarIndex_vector_y(2)}.CheckSub;
     141%             ivar=[ivar CellInfo{icell}.VarIndex_vector_y(1)];
     142%             ivar_1=[ivar_1 CellInfo{icell}.VarIndex_vector_y(2)];
     143%         end
     144%     end
     145% end
     146% % subtract fields if relevant
     147% for imod=1:numel(ivar)
     148%         VarName=ProjData.ListVarName{ivar(imod)};
     149%         VarName_1=ProjData.ListVarName{ivar_1(imod)};
     150%         ProjData.(VarName)=double(ProjData.(VarName))-double(ProjData.(VarName_1));
     151%         ind_remove(ivar_1(imod))=1;
     152% end
     153% ProjData.ListVarName(find(ind_remove))=[];
     154% ProjData.VarDimName(find(ind_remove))=[];
     155% ProjData.VarAttribute(find(ind_remove))=[];
     156
    122157%-----------------------------------------------------------------
    123158%project on a set of points
    … …  
    387422            XName=FieldData.ListVarName{ivar_X};
    388423            YName=FieldData.ListVarName{ivar_Y};
    389             coord_x=FieldData.(XName);
    390             coord_y=FieldData.(YName);
     424            coord_x=FieldData.(FieldData.ListVarName{ivar_X});
     425            coord_y=FieldData.(FieldData.ListVarName{ivar_Y});
    391426            % image or 2D matrix
    392427        case 'grid' %case of structured coordinates
    … …  
    584619                costheta=cos(theta(isegment));
    585620                sintheta=sin(theta(isegment));
    586 %                 XIsegment=LineCoord(isegment,1)+DX_edge*costheta:DX*costheta:LineCoord(isegment+1,1));
    587 %                 YIsegment=(LineCoord(isegment,2)+DX_edge*sintheta:DX*sintheta:LineCoord(isegment+1,2));
     621                %                 XIsegment=LineCoord(isegment,1)+DX_edge*costheta:DX*costheta:LineCoord(isegment+1,1));
     622                %                 YIsegment=(LineCoord(isegment,2)+DX_edge*sintheta:DX*sintheta:LineCoord(isegment+1,2));
    588623                NbInterval=floor((dlength(isegment)-DX_edge)/DX);
    589624                LastX=DX_edge+DX*NbInterval;
    … …  
    618653ProjData.ListVarName={};
    619654ProjData.VarDimName={};
     655check_abscissa=0;
    620656for icell=1:length(CellInfo)
    621657    % list of variable types to be projected
    … …  
    631667    %% identify vector components
    632668    testU=isfield(CellInfo{icell},'VarIndex_vector_x') &&isfield(CellInfo{icell},'VarIndex_vector_y') ;% test for vectors
    633 %     if testU
    634 %         UName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x};
    635 %         VName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y};
    636 %         vector_x=FieldData.(UName);
    637 %         vector_y=FieldData.(VName);
    638 %     end
     669    %     if testU
     670    %         UName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_x};
     671    %         VName=FieldData.ListVarName{CellInfo{icell}.VarIndex_vector_y};
     672    %         vector_x=FieldData.(UName);
     673    %         vector_y=FieldData.(VName);
     674    %     end
    639675    %identify error flag
    640676    errorflag=0; %default, no error flag
    … …  
    654690        %case of unstructured coordinates
    655691        case 'scattered'
    656             XName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
    657             YName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
    658             coord_x=FieldData.(XName);
    659             coord_y=FieldData.(YName);
    660             ProjData.ListVarName=[ProjData.ListVarName {XName}];
    661             ProjData.VarDimName=[ProjData.VarDimName {XName}];
    662             nbvar=numel(ProjData.ListVarName);
    663             ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
     692%             XName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
     693%             YName= FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)};
     694            coord_x=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)});
     695            coord_y=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)});
     696           
    664697            if isequal(ProjMode,'projection')
    665698                if width==0
    … …  
    667700                    return
    668701                end
     702                ProjData.ListVarName=[ProjData.ListVarName {FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}}];
     703                ProjData.VarDimName=[ProjData.VarDimName {FieldData.ListVarName{CellInfo{icell}.CoordIndex(end-1)}}];
     704                nbvar=numel(ProjData.ListVarName);
     705                ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
    669706                % select the (non false) input data located in the band of projection
    670707                flagsel=(errorflag==0) & ((coord_y -yinf(1))*(xinf(2)-xinf(1))>(coord_x-xinf(1))*(yinf(2)-yinf(1))) ...
    … …  
    678715                Xproj=(coord_x-ObjectData.Coord(1,1))*costheta + (coord_y-ObjectData.Coord(1,2))*sintheta; %projection on the line
    679716                [Xproj,indsort]=sort(Xproj);% sort points by increasing absissa along the projection line
    680                 ProjData.(XName)=Xproj;
     717                ProjData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)})=Xproj;
    681718                for ivar=1:numel(VarIndex)
    682719                    ProjData.(VarName{ivar})=FieldData.(VarName{ivar})(flagsel);% restrict vrtibles to the projection band
    683720                    ProjData.(VarName{ivar})=ProjData.(VarName{ivar})(indsort);% sort by absissa
    684721                    ProjData.ListVarName=[ProjData.ListVarName VarName{ivar}];
    685                     ProjData.VarDimName=[ProjData.VarDimName {XName}];
     722                    ProjData.VarDimName=[ProjData.VarDimName {FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)}}];
    686723                    ProjData.VarAttribute{nbvar+ivar}=FieldData.VarAttribute{VarIndex(ivar)};%reproduce var attribute
    687724                    if isfield(ProjData.VarAttribute{nbvar+ivar},'Role')
    … …  
    695732                end
    696733            elseif isequal(ProjMode,'interp_lin')  %filtering %linear interpolation:
     734                if ~check_abscissa
     735                    XName=FieldData.ListVarName{CellInfo{icell}.CoordIndex(end)};
     736                    ProjData.ListVarName=[ProjData.ListVarName {XName}];
     737                    ProjData.VarDimName=[ProjData.VarDimName {XName}];
     738                    nbvar=numel(ProjData.ListVarName);
     739                    ProjData.VarAttribute{nbvar}.long_name='abscissa along line';
     740                    check_abscissa=1; % define abcissa only once
     741                end
    697742                if ~isequal(errorflag,0)
    698743                    VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};
    … …  
    704749                        FieldData.(VarName)=FieldData.(VarName)(indsel);
    705750                    end
    706                 end                   
     751                end
    707752                [ProjVar,ListFieldProj,VarAttribute,errormsg]=calc_field_interp([coord_x coord_y],FieldData,CellInfo{icell}.FieldName,XI,YI);
    708753                ProjData.X=Xproj;
     754                nbvar=numel(ProjData.ListVarName);
    709755                ProjData.ListVarName=[ProjData.ListVarName ListFieldProj];
    710756                ProjData.VarAttribute=[ProjData.VarAttribute VarAttribute];
    … …  
    713759                    if isfield(VarAttribute{ivar},'Role')
    714760                        if  strcmp(VarAttribute{ivar}.Role,'vector_x');
    715                             ivar_U=ivar;
     761                            ivar_U=ivar+nbvar;
    716762                        elseif strcmp(VarAttribute{ivar}.Role,'vector_y');
    717                             ivar_V=ivar;
     763                            ivar_V=ivar+nbvar;
    718764                        end
    719765                    end
    … …  
    741787                    if isfield(VarAttribute{ivar},'Role')
    742788                        if  strcmp(VarAttribute{ivar}.Role,'vector_x');
    743                             ivar_U=ivar;
     789                            ivar_U=ivar+nbvar;
    744790                        elseif strcmp(VarAttribute{ivar}.Role,'vector_y');
    745                             ivar_V=ivar;
     791                            ivar_V=ivar+nbvar;
    746792                        end
    747793                    end
    … …  
    861907                end
    862908            end
    863     end 
     909    end
    864910    if ~isempty(ivar_U) && ~isempty(ivar_V)
    865         vector_x =ProjData.(ProjData.ListVarName{nbvar+ivar_U});
    866          ProjData.(ProjData.ListVarName{nbvar+ivar_U}) =cos(theta)*vector_x+sin(theta)*ProjData.(ProjData.ListVarName{nbvar+ivar_V});
    867           ProjData.(ProjData.ListVarName{nbvar+ivar_V}) =-sin(theta)*vector_x+cos(theta)*ProjData.(ProjData.ListVarName{nbvar+ivar_V});
    868     end
    869        
     911        vector_x =ProjData.(ProjData.ListVarName{ivar_U});
     912        ProjData.(ProjData.ListVarName{ivar_U}) =cos(theta)*vector_x+sin(theta)*ProjData.(ProjData.ListVarName{ivar_V});
     913        ProjData.(ProjData.ListVarName{ivar_V}) =-sin(theta)*vector_x+cos(theta)*ProjData.(ProjData.ListVarName{ivar_V});
     914    end
    870915end
    871916
    … …  
    12301275                   
    12311276                    if isfield(CellInfo{icell},'CheckSub') && CellInfo{icell}.CheckSub && ~isempty(vector_x_proj)
    1232                         ProjData.(ListVarName{vector_x_proj})=ProjData.(ListVarName{vector_x_proj})-VarVal{1};
    1233                         ProjData.(ListVarName{vector_y_proj})=ProjData.(ListVarName{vector_y_proj})-VarVal{2};
     1277                        ProjData.(FieldData.ListVarName{vector_x_proj})=ProjData.(FieldData.ListVarName{vector_x_proj})-VarVal{1};
     1278                        ProjData.(FieldData.ListVarName{vector_y_proj})=ProjData.(FieldData.ListVarName{vector_y_proj})-VarVal{2};
     1279                        ListVarName={};% no new variable
     1280                        VarAttribute={};
    12341281                    else
    12351282                        VarDimName=cell(size(ListVarName));
    … …  
    12431290                        end
    12441291                    end
     1292                    vector_x_proj=CellInfo{icell}.VarIndex_vector_x; %preserve for next cell
     1293                    vector_y_proj=CellInfo{icell}.VarIndex_vector_y; %preserve for next cell
    12451294            end
    12461295           
    … …  
    14911540                            ListVarName=[ListVarName VarName];
    14921541                            VarAttribute{length(ListVarName)}=FieldData.VarAttribute{ivar}; %reproduce the variable attributes
    1493                             eval(['ProjData.' VarName '=squeeze(FieldData.' VarName '(iz,:,:));'])% select the z index iz
     1542                            ProjData.(VarName)=squeeze(FieldData.(VarName)(iz,:,:));% select the z index iz
    14941543                            %TODO : do a vertical average for a thick plane
    14951544                            if test_interp(2) || test_interp(3)
    1496                                 eval(['ProjData.' VarName '=interp2(Coord{3},Coord{2},ProjData.' VarName ',Coord_x,Coord_y'');'])
     1545                                ProjData.(VarName)=interp2(Coord{3},Coord{2},ProjData.(VarName),Coord_x,Coord_y);
    14971546                            end
    14981547                        end
    … …  
    15351584                if ~isempty(ivar_W)
    15361585                    WName=FieldData.ListVarName{ivar_W};
    1537                     eval(['ProjData.' VName '=ProjData.' VName '+ ProjData.' WName '*sin(Theta);'])%
    1538                     eval(['ProjData.' WName '=NormVec_X*ProjData.' UName '+ NormVec_Y*ProjData.' VName '+ NormVec_Z* ProjData.' WName ';']);
     1586                    ProjData.(VName)=ProjData.(VName)+ ProjData.(WName)*sin(Theta);%
     1587                    ProjData.(WName)=NormVec_X*ProjData.(UName)+ NormVec_Y*ProjData.(VName)+ NormVec_Z* ProjData.(WName);
    15391588                end
    15401589            end
Note: See TracChangeset for help on using the changeset viewer.