source: trunk/src/calc_field_interp.m @ 533

Last change on this file since 533 was 533, checked in by sommeria, 12 years ago

bugs repaired in proj_field

File size: 4.2 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')
102    nbvar=numel(ListVarName);
103    ListVarName{nbvar+1}='FF';
104    VarVal{nbvar+1}=isnan(VarVal{nbvar});
105    VarAttribute{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.