source: trunk/src/calc_field_tps.m @ 534

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

various bugs corrected

File size: 5.5 KB
Line 
1
2%'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them
3%---------------------------------------------------------------------
4% [DataOut,VarAttribute,errormsg]=calc_field_tps(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 [DataOut,VarAttribute,errormsg]=calc_field_tps(Coord_tps,NbSites,SubRange,FieldVar,Operation,Coord_interp)
18
19%list of defined scalars to display in menus (in addition to 'ima_cor').
20% a type is associated to each scalar:
21%              'discrete': related to the individual velocity vectors, not interpolated by patch
22%              'vel': calculated from velocity components, continuous field (interpolated with velocity)
23%              'der': needs spatial derivatives
24%              'var': the scalar name corresponds to a field name in the netcdf files
25% a specific variable name for civ1 and civ2 fields are also associated, if
26% the scalar is calculated from other fields, as explicited below
27errormsg='';
28   
29%% nbre of subdomains
30if ndims(Coord_interp)==3
31    nb_coord=size(Coord_interp,3);
32    npx=size(Coord_interp,2);
33    npy=size(Coord_interp,1);
34    nb_sites=npx*npy;
35    Coord_interp=reshape(Coord_interp,nb_sites,nb_coord);
36else
37    nb_coord=size(Coord_interp,2);
38    nb_sites=size(Coord_interp,1);
39end
40NbSubDomain=size(Coord_tps,3);
41nbval=zeros(nb_sites,1);
42
43%% list of operations
44check_grid=0;
45check_der=0;
46for ilist=1:length(Operation)
47    OperationType=regexprep(Operation{ilist},'(.+','');
48    switch OperationType
49        case 'vec'
50            check_grid=1;
51            DataOut.U=zeros(nb_sites,1);
52            DataOut.V=zeros(nb_sites,1);
53            VarAttribute{1}.Role='vector_x';
54            VarAttribute{2}.Role='vector_y';
55        case {'U','V','norm'}
56            check_grid=1;
57            DataOut.(OperationType)=zeros(nb_sites,1);
58            VarAttribute{1}.Role='scalar';
59        case {'curl','div','strain'}
60            check_der=1;
61            DataOut.(OperationType)=zeros(nb_sites,1);
62            VarAttribute{1}.Role='scalar';
63    end
64end
65Attr_FF.Role='errorflag';
66VarAttribute=[VarAttribute {Attr_FF}];
67
68%% loop on subdomains
69for isub=1:NbSubDomain
70    nbvec_sub=NbSites(isub);
71    check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
72    ind_sel=find(sum(check_range,2)==nb_coord);
73    nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
74    if check_grid
75        EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
76    end
77    if check_der
78        [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
79    end
80    ListVar={};
81    for ilist=1:length(Operation)
82        var_count=numel(ListVar);
83        switch Operation{ilist}
84            case 'vec(U,V)'
85                ListVar=[ListVar {'U', 'V'}];
86                VarAttribute{var_count+1}.Role='vector_x';
87                VarAttribute{var_count+2}.Role='vector_y';
88                DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
89                DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
90            case 'U'
91                ListVar=[ListVar {'U'}];
92                VarAttribute{var_count+1}.Role='scalar';
93                DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
94            case 'V'
95                ListVar=[ListVar {'V'}];
96                VarAttribute{var_count+1}.Role='scalar';
97                DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
98            case 'norm(U,V)'
99                ListVar=[ListVar {'norm'}];
100                VarAttribute{var_count+1}.Role='scalar';
101                U=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
102                V=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
103                DataOut.norm(ind_sel)=sqrt(U.*U+V.*V);
104            case 'curl(U,V)'
105                ListVar=[ListVar {'curl'}];
106                VarAttribute{var_count+1}.Role='scalar';
107                DataOut.curl(ind_sel)=DataOut.curl(ind_sel)-EMDY *FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
108            case 'div(U,V)'
109                ListVar=[ListVar {'div'}];
110                VarAttribute{var_count+1}.Role='scalar';
111                DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*FieldVar(1:nbvec_sub+3,isub,1)+EMDY *FieldVar(1:nbvec_sub+3,isub,2);
112            case 'strain(U,V)'
113                ListVar=[ListVar {'strain'}];
114                VarAttribute{var_count+1}.Role='scalar';
115                DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
116        end
117    end
118end
119DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
120nbval(nbval==0)=1;% to avoid division by zero for averaging
121ListFieldOut=fieldnames(DataOut);
122for ifield=1:numel(ListFieldOut)
123    DataOut.(ListFieldOut{ifield})=DataOut.(ListFieldOut{ifield})./nbval;
124    DataOut.(ListFieldOut{ifield})=reshape(DataOut.(ListFieldOut{ifield}),npy,npx);
125end
126
127
128
129
Note: See TracBrowser for help on using the repository browser.