source: trunk/src/calc_field_tps.m @ 660

Last change on this file since 660 was 651, checked in by sommeria, 11 years ago

various bugs corrected. Introduction of campaign lists in Open menu

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