source: trunk/src/calc_field_interp.m @ 579

Last change on this file since 579 was 579, checked in by sommeria, 8 years ago

bugs corrected for reading vort div

File size: 5.8 KB
RevLine 
[575]1%'calc_field_interp': defines fields (velocity, vort, div...) from civ data and calculate them
2% for projection with linear interpolation
[515]3%---------------------------------------------------------------------
[579]4% [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,FieldName,XI,YI)
[515]5%
6% OUTPUT:
[575]7% VarVal: array giving the values of the calculated field
[579]8% ListVarName: corresponding list of variable names
9% VarAttribute: corresponding list of variable attributes, each term #ilist is of the form VarAttribute{ilist}.tag=value
[515]10%
11% INPUT:
[579]12% Coord(nbpoints,2): matrix of x,y coordinates of theinput data points
13% Data: inputfield structure
14% FieldName: string representing the field to calculate, or cell array of fields (as displayed in uvmat/FieldName)
15% XI, YI: set of x and y coordiantes where the fields need to be linearly interpolated
[515]16
[579]17function [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,FieldName,XI,YI)
[515]18
[579]19VarVal={};
[515]20ListVarName={};
21VarAttribute={};
22errormsg='';
[517]23InputVarList={};
[579]24if ischar(FieldName),FieldName={FieldName};end
25check_skipped=zeros(size(FieldName));% default, =1 to mark the variables which can be calculated
26Operator=cell(size(FieldName));
27for ilist=1:numel(FieldName)
28    Operator{ilist}='';%default empty operator (vec, norm,...)
29    r=regexp(FieldName{ilist},'(?<Operator>(^vec|^norm|^curl|^div|^strain))\((?<UName>.+),(?<VName>.+)\)$','names');% analyse the field name
30    if isempty(r) % the field name is a variable itself
31        if ~isfield(Data,FieldName{ilist})
32            check_skipped(ilist)=1; %variable not found
[521]33        else
[579]34            if isempty(find(strcmp(FieldName{ilist},InputVarList), 1));
35                InputVarList=[InputVarList FieldName{ilist}];% the variable is added to the list of input variables if it is not already in the list
[521]36            end
[517]37        end
38    else
[579]39        if ~isfield(Data,r.UName)||~isfield(Data,r.VName)%needed input variable not found
[521]40            check_skipped(ilist)=1;
[579]41        elseif strcmp(r.Operator,'curl')||strcmp(r.Operator,'div')||strcmp(r.Operator,'strain')
42            Operator{ilist}=r.Operator;
43            switch r.Operator
44                case 'curl'% case of CivX data format
45                    UName{ilist}='vort';
46                    Data.vort=Data.DjUi(:,1,2)-Data.DjUi(:,2,1);
47                case 'div'
48                    UName{ilist}='div';
49                    Data.div=Data.DjUi(:,1,1)+Data.DjUi(:,2,2);
50                case 'strain'
51                    UName{ilist}='strain';
52                    Data.strain=Data.DjUi(:,1,2)+Data.DjUi(:,2,1);
53            end
54            InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
[521]55        else
56            UName{ilist}=r.UName;
57            VName{ilist}=r.VName;
58            if isempty(find(strcmp(r.UName,InputVarList)));
59                InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
60            end
[579]61            if isempty(find(strcmp(r.VName,InputVarList), 1));
[521]62                InputVarList=[InputVarList VName{ilist}]; %the variable is added to the list if it is not already in the list
63            end
64            Operator{ilist}=r.Operator;
[517]65        end
[515]66    end
67end
[517]68%create interpolator for linear interpolation
69if exist('XI','var')
70    for ilist=1:numel(InputVarList)
71        F.(InputVarList{ilist})=TriScatteredInterp(Coord,Data.(InputVarList{ilist}),'linear');
72    end
[515]73end
[579]74for ilist=1:numel(FieldName)
[521]75    if ~check_skipped(ilist)
[533]76        nbvar=numel(ListVarName);
77        switch Operator{ilist}
78            case 'vec'
[517]79                if exist('XI','var')
[533]80                    VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
81                    VarVal{nbvar+2}=F.(VName{ilist})(XI,YI);
[517]82                else
[533]83                    VarVal{nbvar+1}=Data.(UName{ilist});
84                    VarVal{nbvar+2}=Data.(VName{ilist});
[517]85                end
[533]86                ListVarName{nbvar+1}=UName{ilist};
87                ListVarName{nbvar+2}=VName{ilist};
88                VarAttribute{nbvar+1}.Role='vector_x';
89                VarAttribute{nbvar+2}.Role='vector_y';
90            case 'norm'
91                if exist('XI','var')
92                    U2=F.(UName{ilist})(XI,YI).*F.(UName{ilist})(XI,YI);
93                    V2=F.(VName{ilist})(XI,YI).*F.(VName{ilist})(XI,YI);
94                else
95                    U2=Data.(UName{ilist}).*Data.(UName{ilist});
96                    V2=Data.(VName{ilist}).*Data.(VName{ilist});
97                end
98                VarVal{nbvar+1}=sqrt(U2+V2);
99                ListVarName{nbvar+1}='norm';
[517]100                VarAttribute{nbvar+1}.Role='scalar';
[579]101            case {'curl','div','strain'}
102                if exist('XI','var')
103                    VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
104                else
105                    VarVal{nbvar+1}=Data.(UName{ilist});
106                end
107                ListVarName{nbvar+1}=UName{ilist};
108                VarAttribute{nbvar+1}.Role='scalar';
[533]109            otherwise
[579]110                if ~isempty(FieldName{ilist})
[533]111                    if exist('XI','var')
[579]112                        VarVal{nbvar+1}=F.(FieldName{ilist})(XI,YI);
[533]113                    else
[579]114                        VarVal{nbvar+1}= Data.(FieldName{ilist});
[533]115                    end
[579]116                    ListVarName{nbvar+1}=FieldName{ilist};
[533]117                    VarAttribute{nbvar+1}.Role='scalar';
118                end
119        end
[515]120    end
121end
[517]122% put an error flag to indicate NaN data
[579]123if exist('XI','var')&&~isempty(VarVal)
[533]124    nbvar=numel(ListVarName);
125    ListVarName{nbvar+1}='FF';
126    VarVal{nbvar+1}=isnan(VarVal{nbvar});
127    VarAttribute{nbvar+1}.Role='errorflag';
[517]128end
[515]129
130% Attr_FF.Role='errorflag';
131% VarAttribute=[VarAttribute {Attr_FF}];
132
133
134
135
136
Note: See TracBrowser for help on using the repository browser.