source: trunk/src/calc_field_tps.m @ 516

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

improvement of calc-field and combination of two fields

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            DataOut.U=zeros(nb_sites,1);
51            DataOut.V=zeros(nb_sites,1);
52            VarAttribute{1}.Role='vector_x';
53            VarAttribute{2}.Role='vector_y';
54        case {'U','V','norm'}
55            check_grid=1;
56            DataOut.(OperationType)=zeros(nb_sites,1);
57            VarAttribute{1}.Role='scalar';
58        case {'curl','div','strain'}
59            check_der=1;
60            DataOut.(OperationType)=zeros(nb_sites,1);
61            VarAttribute{1}.Role='scalar';
62    end
63end
64Attr_FF.Role='errorflag';
65VarAttribute=[VarAttribute {Attr_FF}];
66
67%% loop on subdomains
68for isub=1:NbSubDomain
69    nbvec_sub=NbSites(isub);
70    check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
71    ind_sel=find(sum(check_range,2)==nb_coord);
72    nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
73    if check_grid
74        EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
75    end
76    if check_der
77        [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
78    end
79    ListVar={};
80    for ilist=1:length(Operation)
81        var_count=numel(ListVar);
82        switch Operation{ilist}
83            case 'vec(U,V)'
84                ListVar=[ListVar {'U', 'V'}];
85                VarAttribute{var_count+1}.Role='vector_x';
86                VarAttribute{var_count+2}.Role='vector_y';
87                DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
88                DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
89            case 'U'
90                ListVar=[ListVar {'U'}];
91                VarAttribute{var_count+1}.Role='scalar';
92                DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
93            case 'V'
94                ListVar=[ListVar {'V'}];
95                VarAttribute{var_count+1}.Role='scalar';
96                DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
97            case 'norm(U,V)'
98                ListVar=[ListVar {'norm'}];
99                VarAttribute{var_count+1}.Role='scalar';
100                U=DataOut.U(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,1);
101                V=DataOut.V(ind_sel)+EM *FieldVar(1:nbvec_sub+3,isub,2);
102                DataOut.norm(ind_sel)=sqrt(U.*U+V.*V);
103            case 'curl(U,V)'
104                ListVar=[ListVar {'curl'}];
105                VarAttribute{var_count+1}.Role='scalar';
106                DataOut.curl(ind_sel)=DataOut.curl(ind_sel)-EMDY *FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
107            case 'div(U,V)'
108                ListVar=[ListVar {'div'}];
109                VarAttribute{var_count+1}.Role='scalar';
110                DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*FieldVar(1:nbvec_sub+3,isub,1)+EMDY *FieldVar(1:nbvec_sub+3,isub,2);
111            case 'strain(U,V)'
112                ListVar=[ListVar {'strain'}];
113                VarAttribute{var_count+1}.Role='scalar';
114                DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*FieldVar(1:nbvec_sub+3,isub,1)+EMDX *FieldVar(1:nbvec_sub+3,isub,2);
115        end
116    end
117end
118DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
119nbval(nbval==0)=1;% to avoid division by zero for averaging
120ListFieldOut=fieldnames(DataOut);
121for ifield=1:numel(ListFieldOut)
122    DataOut.(ListFieldOut{ifield})=DataOut.(ListFieldOut{ifield})./nbval;
123    DataOut.(ListFieldOut{ifield})=reshape(DataOut.(ListFieldOut{ifield}),npy,npx);
124end
125
126
127
128
Note: See TracBrowser for help on using the repository browser.