source: trunk/src/calc_field_interp.m @ 521

Last change on this file since 521 was 521, checked in by sommeria, 9 years ago

various bugs corrected

File size: 4.0 KB
Line 
1
2%'calc_field': defines fields (velocity, vort, div...) from civ data and calculate them
3%---------------------------------------------------------------------
4% [DataOut,VarAttribute,errormsg]=calc_field_interp(Coord_tps,NbSites,SubRange,FieldVar,Operation,Coord_interp)
5%
6% OUTPUT:
7% DataOut: structure representing the output fields
8%
9% INPUT:
10% Coord_tps:
11% NbSites
12% SubRange
13% FieldVar
14% Operation: cell array representing the list of operations (eg div, rot..)
15% Coord_interp: coordiantes of sites on which the fields need to be calculated
16
17function [VarVal,ListVarName,VarAttribute,errormsg]=calc_field_interp(Coord,Data,Operation,XI,YI)
18
19VarVal=[];
20ListVarName={};
21VarAttribute={};
22errormsg='';
23InputVarList={};
24if ischar(Operation),Operation={Operation};end
25check_skipped=zeros(size(Operation));
26Operator=cell(size(Operation));
27for ilist=1:numel(Operation)
28    r=regexp(Operation{ilist},'(?<Operator>(^vec|^norm))\((?<UName>.+),(?<VName>.+)\)$','names');
29    if isempty(r) % the operation is the variable
30        if ~isfield(Data,Operation{ilist})
31            check_skipped(ilist)=1;
32        else
33            if isempty(find(strcmp(Operation{ilist},InputVarList)));
34                InputVarList=[InputVarList Operation{ilist}];% the variable is added to the list if it is not already in the list
35            end
36            Operator{ilist}='';
37        end
38    else
39        if ~isfield(Data,r.UName)||~isfield(Data,r.VName)
40            check_skipped(ilist)=1;
41        else
42            UName{ilist}=r.UName;
43            VName{ilist}=r.VName;
44            if isempty(find(strcmp(r.UName,InputVarList)));
45                InputVarList=[InputVarList UName{ilist}]; %the variable is added to the list if it is not already in the list
46            end
47            if isempty(find(strcmp(r.VName,InputVarList)));
48                InputVarList=[InputVarList VName{ilist}]; %the variable is added to the list if it is not already in the list
49            end
50            Operator{ilist}=r.Operator;
51        end
52    end
53end
54%create interpolator for linear interpolation
55if exist('XI','var')
56    for ilist=1:numel(InputVarList)
57        F.(InputVarList{ilist})=TriScatteredInterp(Coord,Data.(InputVarList{ilist}),'linear');
58    end
59end
60for ilist=1:numel(Operation)
61    if ~check_skipped(ilist)
62    nbvar=numel(ListVarName);
63    switch Operator{ilist}
64        case 'vec'
65            if exist('XI','var')
66                VarVal{nbvar+1}=F.(UName{ilist})(XI,YI);
67                VarVal{nbvar+2}=F.(VName{ilist})(XI,YI);
68            else
69                VarVal{nbvar+1}=Data.(UName{ilist});
70                VarVal{nbvar+2}=Data.(VName{ilist});
71            end
72            ListVarName{nbvar+1}=UName{ilist};
73            ListVarName{nbvar+2}=VName{ilist};
74            VarAttribute{nbvar+1}.Role='vector_x';
75            VarAttribute{nbvar+2}.Role='vector_y';
76        case 'norm'
77            if exist('XI','var')
78                U2=F.(UName{ilist})(XI,YI).*F.(UName{ilist})(XI,YI);
79                V2=F.(VName{ilist})(XI,YI).*F.(VName{ilist})(XI,YI);
80            else
81                U2=Data.(UName{ilist}).*Data.(UName{ilist});
82                V2=Data.(VName{ilist}).*Data.(VName{ilist});
83            end
84            VarVal{nbvar+1}=sqrt(U2+V2);
85            ListVarName{nbvar+1}='norm';
86            VarAttribute{nbvar+1}.Role='scalar';
87        otherwise
88            if ~isempty(Operation{ilist})
89                if exist('XI','var')
90                    VarVal{nbvar+1}=F.(Operation{ilist})(XI,YI);
91                else
92                    VarVal{nbvar+1}= Data.(Operation{ilist});
93                end
94                ListVarName{nbvar+1}=Operation{ilist};
95                VarAttribute{nbvar+1}.Role='scalar';
96            end
97    end
98    end
99end
100% put an error flag to indicate NaN data
101if exist('XI','var')
102nbvar=numel(ListVarName);
103ListVarName{nbvar+1}='FF';
104VarVal{nbvar+1}=isnan(VarVal{nbvar});
105VarAttribute{nbvar+1}.Role='errorflag';
106end
107
108% Attr_FF.Role='errorflag';
109% VarAttribute=[VarAttribute {Attr_FF}];
110
111
112
113
114
Note: See TracBrowser for help on using the repository browser.