Changeset 404 for trunk/src/calc_field.m


Ignore:
Timestamp:
Apr 30, 2012, 6:46:45 PM (12 years ago)
Author:
sommeria
Message:

various bugs corrected. nc2struct_toolbox suppressed (correspond to obsolete versions of Matlab)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/calc_field.m

    r399 r404  
     1
     2
    13%'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them
    24%---------------------------------------------------------------------
    … …  
    3537% the scalar is calculated from other fields, as explicited below
    3638
    37 list_field={'velocity';...%image correlation corresponding to a vel vector
     39%% list of field options implemented
     40FieldOptions={'velocity';...%image correlation corresponding to a vel vector
    3841    'ima_cor';...%image correlation corresponding to a vel vector
    3942    'norm_vel';...%norm of the velocity
    … …  
    4851errormsg=[]; %default error message
    4952if ~exist('FieldList','var')
    50     DataOut=list_field;% gives the list of possible fields in the absence of input
     53    DataOut=FieldOptions;% gives the list of possible field inputs in the absence of input
     54    return
     55end
     56
     57%% check input
     58if ~exist('DataIn','var')
     59    DataIn=[];
     60end
     61if ischar(FieldList)
     62    FieldList={FieldList};%convert a string input to a cell with one string element
     63end
     64check_grid=0;
     65check_der=0;
     66check_calc=ones(size(FieldList));
     67for ilist=1:length(FieldList)
     68    switch FieldList{ilist}
     69        case {'u','v'}
     70            check_grid=1;% needs a regular grid
     71        case{'vort','div','strain'}% needs spatial derivatives spatial derivatives
     72            check_der=1;
     73        case {'velocity','norm_vel'};
     74        otherwise
     75            check_calc(ilist)=0;
     76    end
     77end
     78FieldList=FieldList(check_calc==1);
     79if isempty(FieldList)
     80    DataOut=DataIn;
     81    return
     82end
     83if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
     84    nbcoord=3;
    5185else
    52     if ~exist('DataIn','var')
    53         DataIn=[];
    54     end
    55     if ischar(FieldList)
    56         FieldList={FieldList};%convert a string input to a cell with one string element
    57     end
    58     if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
    59         nbcoord=3;
    60     else
    61         nbcoord=2;
    62     end
    63     ListVarName={};
    64     ValueList={};
    65     RoleList={};
    66     units_cell={};
    67    
    68     %% interpolation with new civ data
    69     if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps')&& exist('Coord_interp','var')
    70         DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
    71         for ilist=1:numel(DataOut.ListGlobalAttribute)
    72             DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
    73         end
    74         DataOut.ListVarName={'coord_y','coord_x','FF'};
    75         DataOut.VarDimName{1}='coord_y';
    76         DataOut.VarDimName{2}='coord_x';
    77         XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates
    78         YMax=max(max(DataIn.SubRange(2,:,:)));
    79         XMin=min(min(DataIn.SubRange(1,:,:)));
    80         YMin=min(min(DataIn.SubRange(2,:,:)));
    81         check_der=0;
    82         check_val=0;
    83         nb_sites=size(Coord_interp,1);
    84         nb_coord=size(Coord_interp,2);
     86    nbcoord=2;
     87end
     88ListVarName={};
     89ValueList={};
     90RoleList={};
     91units_cell={};
     92
     93%% interpolation with new civ data
     94if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var') || check_grid ||check_der)
     95    %create a default grid if needed
     96    if  ~exist('Coord_interp','var')
     97        coord_x=DataIn.XMin:DataIn.Mesh:DataIn.XMax;
     98        coord_y=DataIn.XMin:DataIn.Mesh:DataIn.YMax;
     99        [XI,YI]=meshgrid(coord_x,coord_y);
     100        XI=reshape(XI,[],1);
     101        YI=reshape(YI,[],1);
     102        Coord_interp=[XI YI];
     103    end
     104    DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; %reproduce global attribute
     105    for ilist=1:numel(DataOut.ListGlobalAttribute)
     106        DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
     107    end
     108    DataOut.ListVarName={'coord_y','coord_x','FF'};
     109    DataOut.VarDimName{1}='coord_y';
     110    DataOut.VarDimName{2}='coord_x';
     111    XMax=max(max(DataIn.SubRange(1,:,:)));% extrema of the coordinates
     112    YMax=max(max(DataIn.SubRange(2,:,:)));
     113    XMin=min(min(DataIn.SubRange(1,:,:)));
     114    YMin=min(min(DataIn.SubRange(2,:,:)));
     115%     check_der=0;
     116%     check_val=0;
     117    nb_sites=size(Coord_interp,1);
     118    nb_coord=size(Coord_interp,2);
     119    %initialise output
     120    for ilist=1:length(FieldList)
     121        switch FieldList{ilist}
     122            case 'velocity'
     123                DataOut.U=zeros(nb_sites,1);
     124                DataOut.V=zeros(nb_sites,1);
     125            otherwise
     126  %          case{'vort','div','strain'}% case of spatial derivatives
     127                DataOut.(FieldList{ilist})=zeros(nb_sites,1);
     128%                 otherwise % case of a scalar
     129%                     check_val=1;
     130%                     DataOut.(FieldList{ilist})=zeros(size(Coord_interp,1));
     131        end
     132    end
     133    nbval=zeros(nb_sites,1);
     134    NbSubDomain=size(DataIn.SubRange,3);
     135    %DataIn.Coord_tps=DataIn.Coord_tps(1:end-3,:,:);% suppress the 3 zeros used to fit with the dimensions of variables
     136    for isub=1:NbSubDomain
     137        nbvec_sub=DataIn.NbSites(isub);
     138        check_range=(Coord_interp >=ones(nb_sites,1)*DataIn.SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*DataIn.SubRange(:,2,isub)');
     139        ind_sel=find(sum(check_range,2)==nb_coord);
     140        %rho smoothing parameter
     141        %                 epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites
     142        %                 ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
     143        nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
     144        if check_val
     145            EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
     146        end
     147        if check_der
     148            [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
     149        end
    85150        for ilist=1:length(FieldList)
    86151            switch FieldList{ilist}
    87152                case 'velocity'
    88                     check_val=1;
    89                     DataOut.U=zeros(nb_sites,1);
    90                     DataOut.V=zeros(nb_sites,1);
    91                 case{'vort','div','strain'}% case of spatial derivatives
    92                     check_der=1;
    93                     DataOut.(FieldList{ilist})=zeros(nb_sites,1);
    94 %                 otherwise % case of a scalar
    95 %                     check_val=1;
    96 %                     DataOut.(FieldList{ilist})=zeros(size(Coord_interp,1));
     153                    ListFields={'U', 'V'};
     154                    VarAttributes{1}.Role='vector_x';
     155                    VarAttributes{2}.Role='vector_y';
     156                    DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
     157                    DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
     158                case 'u'
     159                    ListFields={'U'};
     160                    VarAttributes{1}.Role='scalar';
     161                    DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
     162                case 'v'
     163                    ListFields={'V'};
     164                    VarAttributes{1}.Role='scalar';
     165                    DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
     166                case 'norm_vel'
     167                    ListFields={'norm_vel'};
     168                    VarAttributes{1}.Role='scalar';
     169                    V=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
     170                    V=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
     171                    DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V);
     172                case 'vort'
     173                    ListFields={'vort'};
     174                    VarAttributes{1}.Role='scalar';
     175                    DataOut.vort(ind_sel)=DataOut.vort(ind_sel)+EMDY *DataIn.U_tps(1:nbvec_sub+3,isub)-EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
     176                case 'div'
     177                    ListFields={'div'};
     178                    VarAttributes{1}.Role='scalar';
     179                    DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDY *DataIn.V_tps(1:nbvec_sub+3,isub);
     180                case 'strain'
     181                    ListFields={'strain'};
     182                    VarAttributes{1}.Role='scalar';
     183                    DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
    97184            end
    98185        end
    99         nbval=zeros(nb_sites,1);
    100         NbSubDomain=size(DataIn.SubRange,3);
    101         %DataIn.Coord_tps=DataIn.Coord_tps(1:end-3,:,:);% suppress the 3 zeros used to fit with the dimensions of variables
    102         for isub=1:NbSubDomain
    103             nbvec_sub=DataIn.NbSites(isub);
    104             check_range=(Coord_interp >=ones(nb_sites,1)*DataIn.SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*DataIn.SubRange(:,2,isub)');
    105             ind_sel=find(sum(check_range,2)==nb_coord);
    106             %rho smoothing parameter
    107             %                 epoints = Coord_interp(ind_sel) ;% coordinates of interpolation sites
    108             %                 ctrs=DataIn.Coord_tps(1:nbvec_sub,:,isub);%(=initial points) ctrs
    109             nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
    110             if check_val
    111                 EM = tps_eval(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
    112             end
    113             if check_der
    114                 [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),DataIn.Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
    115             end
    116             for ilist=1:length(FieldList)
    117                 switch FieldList{ilist}
    118                     case 'velocity'
    119                         ListFields={'U', 'V'};
    120                         VarAttributes{1}.Role='vector_x';
    121                         VarAttributes{2}.Role='vector_y';
    122                         DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    123                         DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
    124                     case 'u'
    125                         ListFields={'U'};
    126                         VarAttributes{1}.Role='scalar';
    127                         DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *DataIn.U_tps(1:nbvec_sub+3,isub);
    128                     case 'v'
    129                         ListFields={'V'};
    130                         VarAttributes{1}.Role='scalar';
    131                         DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *DataIn.V_tps(1:nbvec_sub+3,isub);
    132                     case 'vort'
    133                         ListFields={'vort'};
    134                         VarAttributes{1}.Role='scalar';
    135                         DataOut.vort(ind_sel)=DataOut.vort(ind_sel)+EMDY *DataIn.U_tps(1:nbvec_sub+3,isub)-EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
    136                     case 'div'
    137                         ListFields={'div'};
    138                         VarAttributes{1}.Role='scalar';
    139                         DataOut.div(ind_sel)=DataOut.div(ind_sel)+EMDX*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDY *DataIn.V_tps(1:nbvec_sub+3,isub);
    140                     case 'strain'
    141                         ListFields={'strain'};
    142                         VarAttributes{1}.Role='scalar';
    143                         DataOut.strain(ind_sel)=DataOut.strain(ind_sel)+EMDY*DataIn.U_tps(1:nbvec_sub+3,isub)+EMDX *DataIn.V_tps(1:nbvec_sub+3,isub);
    144                 end
    145             end
    146             DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
    147 %            DataOut.FF=reshape(DataOut.FF,numel(yI),numel(xI));
    148             nbval(nbval==0)=1;
    149 %             switch FieldList{1}
    150 %                 case {'velocity','u','v'}
    151 %                     DataOut.U=reshape(DataOut.U./nbval,numel(yI),numel(xI));
    152 %                     DataOut.V=reshape(DataOut.V./nbval,numel(yI),numel(xI));
    153 %                 case 'vort'
    154 %                     DataOut.vort=reshape(DataOut.vort,numel(yI),numel(xI));
    155 %                 case 'div'
    156 %                     DataOut.div=reshape(DataOut.div,numel(yI),numel(xI));
    157 %                 case 'strain'
    158 %                     DataOut.strain=reshape(DataOut.strain,numel(yI),numel(xI));
    159 %             end
    160             DataOut.ListVarName=[DataOut.ListVarName ListFields];
    161             for ilist=3:numel(DataOut.ListVarName)
    162                 DataOut.VarDimName{ilist}={'coord_y','coord_x'};
    163             end
    164             DataOut.VarAttribute={[],[]};
    165             DataOut.VarAttribute{3}.Role='errorflag';
    166             DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];
    167         end
    168     else
    169        
    170         %% civx data
    171         DataOut=DataIn;
    172         for ilist=1:length(FieldList)
    173             if ~isempty(FieldList{ilist})
    174                 [VarName,Value,Role,units]=feval(FieldList{ilist},DataIn);%calculate field with appropriate function named FieldList{ilist}
    175                 ListVarName=[ListVarName VarName];
    176                 ValueList=[ValueList Value];
    177                 RoleList=[RoleList Role];
    178                 units_cell=[units_cell units];
    179             end
    180         end
    181         %erase previous data (except coordinates)
    182         for ivar=nbcoord+1:length(DataOut.ListVarName)
    183             VarName=DataOut.ListVarName{ivar};
    184             DataOut=rmfield(DataOut,VarName);
    185         end
    186         DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
    187         if isfield(DataOut,'VarDimName')
    188             DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
    189         else
    190             errormsg='element .VarDimName missing in input data';
    191             return
    192         end
    193         DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
    194         %append new data
    195         DataOut.ListVarName=[DataOut.ListVarName ListVarName];
    196         for ivar=1:length(ListVarName)
    197             DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
    198             DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
    199             DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
    200             DataOut.(ListVarName{ivar})=ValueList{ivar};
    201         end
    202     end
    203 end
     186        DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
     187        nbval(nbval==0)=1;
     188        DataOut.ListVarName=[DataOut.ListVarName ListFields];
     189        for ilist=3:numel(DataOut.ListVarName)
     190            DataOut.VarDimName{ilist}={'coord_y','coord_x'};
     191        end
     192        DataOut.VarAttribute={[],[]};
     193        DataOut.VarAttribute{3}.Role='errorflag';
     194        DataOut.VarAttribute=[DataOut.VarAttribute VarAttributes];
     195    end
     196else
     197
     198    %% civx data
     199    DataOut=DataIn;
     200    for ilist=1:length(FieldList)
     201        if ~isempty(FieldList{ilist})
     202            [VarName,Value,Role,units]=feval(FieldList{ilist},DataIn);%calculate field with appropriate function named FieldList{ilist}
     203            ListVarName=[ListVarName VarName];
     204            ValueList=[ValueList Value];
     205            RoleList=[RoleList Role];
     206            units_cell=[units_cell units];
     207        end
     208    end
     209    %erase previous data (except coordinates)
     210    for ivar=nbcoord+1:length(DataOut.ListVarName)
     211        VarName=DataOut.ListVarName{ivar};
     212        DataOut=rmfield(DataOut,VarName);
     213    end
     214    DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
     215    if isfield(DataOut,'VarDimName')
     216        DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
     217    else
     218        errormsg='element .VarDimName missing in input data';
     219        return
     220    end
     221    DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
     222    %append new data
     223    DataOut.ListVarName=[DataOut.ListVarName ListVarName];
     224    for ivar=1:length(ListVarName)
     225        DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
     226        DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
     227        DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
     228        DataOut.(ListVarName{ivar})=ValueList{ivar};
     229    end
     230end
     231
    204232
    205233
Note: See TracChangeset for help on using the changeset viewer.