Ignore:
Timestamp:
Aug 17, 2012, 11:47:16 PM (9 years ago)
Author:
sommeria
Message:

various bugs corrected. get_field now used in a passive way from uvmat: variable names are transferred from get_field to uvmat.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/calc_field_interp.m

    r516 r517  
    11
    2 %'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them
     2%'calc_field': defines fields (velocity, vort, div...) from civ data and calculate them
    33%---------------------------------------------------------------------
    44% [DataOut,VarAttribute,errormsg]=calc_field_interp(Coord_tps,NbSites,SubRange,FieldVar,Operation,Coord_interp)
     
    1515% Coord_interp: coordiantes of sites on which the fields need to be calculated
    1616
    17 function [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,FieldVar,Operation,XI,YI)
     17function [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,Operation,XI,YI)
    1818
    1919VarVal=[];
     
    2121VarAttribute={};
    2222errormsg='';
    23 check_u=0;
    24 check_v=0;
    25 for ilist=1:length(Operation)
    26     switch Operation{ilist}
    27         case {'U'}
    28            check_u=1;
    29         case {'V'}
    30             check_v=1;
    31           case {'vec(U,V)','norm(U,V)'} 
    32              check_u=1;
    33              check_v=1;
     23InputVarList={};
     24if ischar(Operation),Operation={Operation};end
     25for ilist=1:numel(Operation)
     26    r=regexp(Operation{ilist},'(?<Operator>(^vec|^norm))\((?<UName>.+),(?<VName>.+)\)$','names');
     27    if isempty(r) % the operation is the variable
     28        if isempty(find(strcmp(Operation{ilist},InputVarList)));
     29            InputVarList=[InputVarList Operation{ilist}];
     30        end
     31        Operator{ilist}='';
     32    else
     33        UName{ilist}=r.UName;
     34        VName{ilist}=r.VName;
     35        if isempty(find(strcmp(r.UName,InputVarList)));
     36            InputVarList=[InputVarList UName{ilist}];
     37        end
     38        if isempty(find(strcmp(r.VName,InputVarList)));
     39            InputVarList=[InputVarList VName{ilist}];
     40        end
     41        Operator{ilist}=r.Operator;
    3442    end
    3543end
    36 if check_u
    37     F_u = TriScatteredInterp(Coord,FieldVar(:,1),'linear');
     44%create interpolator for linear interpolation
     45if exist('XI','var')
     46    for ilist=1:numel(InputVarList)
     47        F.(InputVarList{ilist})=TriScatteredInterp(Coord,Data.(InputVarList{ilist}),'linear');
     48    end
    3849end
    39 if check_v
    40     F_v = TriScatteredInterp(Coord,FieldVar(:,2),'linear');
    41 end
    42 for ilist=1:length(Operation)
     50for ilist=1:numel(Operation)
    4351    nbvar=numel(ListVarName);
    44     switch Operation{ilist}
    45         case 'vec(U,V)'
    46             VarVal{nbvar+1}=F_u(XI,YI);
    47             VarVal{nbvar+2}=F_v(XI,YI);
    48             ListVarName{nbvar+1}='U';
    49             ListVarName{nbvar+2}='V';
     52    switch Operator{ilist}
     53        case 'vec'
     54            if exist('XI','var')
     55                VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
     56                VarVal{nbvar+2}=F.(VName{ilist})(XI,YI);
     57            else
     58                VarVal{nbvar+1}=Data.(UName{ilist});
     59                VarVal{nbvar+2}=Data.(VName{ilist});
     60            end
     61            ListVarName{nbvar+1}=UName{ilist};
     62            ListVarName{nbvar+2}=VName{ilist};
    5063            VarAttribute{nbvar+1}.Role='vector_x';
    5164            VarAttribute{nbvar+2}.Role='vector_y';
    52         case 'U'
    53             VarVal{nbvar+1}=F_u(XI,YI);
    54             ListVarName{nbvar+1}='U';
     65            %         case 'U'
     66            %             VarVal{nbvar+1}=F_u(XI,YI);
     67            %             ListVarName{nbvar+1}='U';
     68            %             VarAttribute{nbvar+1}.Role='scalar';
     69            %         case 'V'
     70            %             VarVal{nbvar+1}=F_v(XI,YI);
     71            %             ListVarName{nbvar+1}='V';
     72            %             VarAttribute{nbvar+1}.Role='scalar';
     73        case 'norm'
     74            if exist('XI','var')
     75                U2=F.(UName{ilist})(XI,YI).*F.(UName{ilist})(XI,YI);
     76                V2=F.(VName{ilist})(XI,YI).*F.(VName{ilist})(XI,YI);
     77            else
     78                U2=Data.(UName{ilist}).*Data.(UName{ilist});
     79                V2=Data.(VName{ilist}).*Data.(VName{ilist});
     80            end
     81            VarVal{nbvar+1}=sqrt(U2+V2);
     82            ListVarName{nbvar+1}='norm';
    5583            VarAttribute{nbvar+1}.Role='scalar';
    56         case 'V'
    57             VarVal{nbvar+1}=F_v(XI,YI);
    58             ListVarName{nbvar+1}='V';
    59             VarAttribute{nbvar+1}.Role='scalar';
    60         case 'norm(U,V)'
    61             VarVal{nbvar+1}=sqrt(F_u(XI,YI).*F_u(XI,YI)+F_v(XI,YI).*F_v(XI,YI));
    62             ListVarName{nbvar+1}='norm(U,V)';
    63             VarAttribute{nbvar+1}.Role='scalar';
     84        otherwise
     85            if ~isempty(Operation{ilist})
     86                if exist('XI','var')
     87                    VarVal{nbvar+1}=F.(Operation{ilist})(XI,YI);
     88                else
     89                    VarVal{nbvar+1}= Data.(Operation{ilist});
     90                end
     91                ListVarName{nbvar+1}=Operation{ilist};
     92                VarAttribute{nbvar+1}.Role='scalar';
     93            end
    6494    end
    6595end
     96% put an error flag to indicate NaN data
     97if exist('XI','var')
    6698nbvar=numel(ListVarName);
    6799ListVarName{nbvar+1}='FF';
    68100VarVal{nbvar+1}=isnan(VarVal{nbvar});
    69101VarAttribute{nbvar+1}.Role='errorflag';
     102end
    70103
    71104% Attr_FF.Role='errorflag';
Note: See TracChangeset for help on using the changeset viewer.