source: trunk/src/calc_field.m @ 534

%'calc_field': defines fields (velocity, vort, div...) from civx data and calculate them
%---------------------------------------------------------------------
% [DataOut,errormsg]=calc_field(FieldList,DataIn,Coord_interp)
%
% OUTPUT:
% Scal: matlab vector representing the scalar values (length nbvec defined by var_read)
%      if no input, Scal=list of programmed scalar names (to put in menus)
%      if only the field name is put as input, vec_A=type of scalar, which can be:
%                   'discrete': related to the individual velocity vectors, not interpolated by patch
%                   'vel': scalar calculated solely from velocity components
%                   'der': needs spatial derivatives
%                   'var': the scalar name directly corresponds to a field name in the netcdf files
% error: error flag
%      error = 0; OK
%      error = 1; the prescribed scalar cannot be read or calculated from available fields
%
% INPUT:
% FieldList: cell array of strings representing the name(s) of the field(s) to calculate
% DataIn: structure representing the field, as defined in check_field_srtructure.m
% Coord_interp(:,nb_coord) optional set of coordinates to interpolate the field (use with thin plate shell)
%
% FUNCTION related
% varname_generator.m: determines the field names to read in the netcdf
% file, depending on the scalar
function [FieldOptions,ColorList]=calc_field
%function [DataOut,errormsg]=calc_field(FieldList,DataIn,Coord_interp)

%list of defined scalars to display in menus (in addition to 'ima_cor').
% a type is associated to each scalar:
%              'discrete': related to the individual velocity vectors, not interpolated by patch
%              'vel': calculated from velocity components, continuous field (interpolated with velocity)
%              'der': needs spatial derivatives
%              'var': the scalar name corresponds to a field name in the netcdf files
% a specific variable name for civ1 and civ2 fields are also associated, if
% the scalar is calculated from other fields, as explicited below

%% list of field options implemented
FieldOptions={'vec(U,V)';...%image correlation corresponding to a vel vector
    'C';...%image correlation corresponding to a vel vector
    'norm(U,V)';...%norm of the velocity
    'curl(U,V)';...%vorticity
    'div(U,V)';...%divergence
    'strain(U,V)';...%rate of strain
    'U';... %u velocity component
    'V';... %v velocity component
    'w';... %w velocity component
    'w_normal';... %w velocity component normal to the plane
    'error'}; %error associated to a vector (for stereo or patch)
ColorList={'C';...%image correlation corresponding to a vel vector
    'norm(U,V)';...%norm of the velocity
    'U';... %u velocity component
    'V';... %v velocity component
    }
% errormsg=[]; %default error message
% if ~exist('FieldList','var')
%     DataOut=FieldOptions;% gives the list of possible field inputs in the absence of input
%     return
% end
return
%% check input
if ~exist('DataIn','var')
    DataIn=[];
end
if ischar(FieldList)
    FieldList={FieldList};%convert a string input to a cell with one string element
end
check_grid=0;
check_der=0;
check_calc=ones(size(FieldList));
for ilist=1:length(FieldList)
    switch FieldList{ilist}
        case {'u','v','velocity','norm_vel','ima_cor'}
            check_grid=1;% needs a regular grid
        case{'vort','div','strain'}% needs spatial derivatives spatial derivatives
            check_der=1;
%         case {'velocity','norm_vel','ima_cor'};
        otherwise
            check_calc(ilist)=0;
    end
end
FieldList=FieldList(check_calc==1);
if isfield(DataIn,'Z')&& isequal(size(DataIn.Z),size(DataIn.X))
    nbcoord=3;
else
    nbcoord=2;
end
ListVarName={};
ValueList={};
RoleList={};
units_cell={};

%% reproduce global attributes
DataOut.ListGlobalAttribute=DataIn.ListGlobalAttribute; 
for ilist=1:numel(DataOut.ListGlobalAttribute)
    DataOut.(DataOut.ListGlobalAttribute{ilist})=DataIn.(DataIn.ListGlobalAttribute{ilist});
end

%% interpolation with new civ data
[CellVarIndex,NbDimVec,VarTypeCell,errormsg]=find_field_cells(DataIn);
icell_tps=[];
for icell=1:numel(CellVarIndex)
    VarType=VarTypeCell{icell};
    if NbDimVec(icell)>=2 && ~isempty(VarType.subrange_tps)      
        icell_tps=[icell_tps icell];
    end
end

%if isfield(DataIn,'SubRange') && isfield(DataIn,'Coord_tps') && (exist('Coord_interp','var') || check_grid ||check_der)
if ~isempty(icell_tps)
    %create a default grid if needed
    if  ~exist('Coord_interp','var')
        for ind=1:numel(icell_tps)
            SubRange=DataIn.(DataIn.ListVarName{VarType{icell_tps(ind)}.subrange_tps});
            XMax(ind)=max(max(SubRange(1,:,:)));% extrema of the coordinates
            YMax(ind)=max(max(SubRange(2,:,:)));
            XMin(ind)=min(min(SubRange(1,:,:)));
            YMin(ind)=min(min(SubRange(2,:,:))); 
        end
        XMax=max(XMax);
        YMax=max(YMax);
        XMin=min(XMin);
        YMin=min(YMin);
        if ~isfield(DataIn,'Mesh')
            DataIn.Mesh=sqrt(2*(XMax-XMin)*(YMax-YMin)/numel(DataIn.Coord_tps));
            % adjust the mesh to a value 1, 2 , 5 *10^n
            ord=10^(floor(log10(DataIn.Mesh)));%order of magnitude
            if DataIn.Mesh/ord>=5
                DataIn.Mesh=5*ord;
            elseif DataIn.Mesh/ord>=2
                DataIn.Mesh=2*ord;
            else
                DataIn.Mesh=ord;
            end
        end
        coord_x=XMin:DataIn.Mesh:XMax;% increase the recommanded mesh to  visualisation
        coord_y=YMin:DataIn.Mesh:YMax;
        DataOut.coord_x=[XMin XMax];
        DataOut.coord_y=[YMin YMax];
        [XI,YI]=meshgrid(coord_x,coord_y);
        Coord_interp=cat(3,XI,YI);%[XI YI];
    end
    npx=size(Coord_interp,2);
    npy=size(Coord_interp,1);
    Coord_interp=reshape(Coord_interp,npx*npy,size(Coord_interp,3));

    %initialise output
    nb_sites=size(Coord_interp,1);
    nb_coord=size(Coord_interp,2);
    nbval=zeros(nb_sites,1);
    NbSubDomain=size(DataIn.SubRange,3);
    DataOut.ListVarName={'coord_y','coord_x'};
    DataOut.VarDimName={{'coord_y'},{'coord_x'}};
    DataOut.VarAttribute{1}=[];
    DataOut.VarAttribute{2}=[];
    for ilist=1:length(FieldList)
        switch FieldList{ilist}
            case 'velocity'
                DataOut.U=zeros(nb_sites,1);
                DataOut.V=zeros(nb_sites,1);
            otherwise
                DataOut.(FieldList{ilist})=zeros(nb_sites,1);
        end
    end
    
    % interpolate data in each subdomain
    for icell=icell_tps
        for isub=1:NbSubDomain
            nbvec_sub=DataIn.(DataIn.ListVarName{VarType{icell}.nbsites_tps});
            SubRange=DataIn.(DataIn.ListVarName{VarType{icell}.subrange_tps});
            check_range=(Coord_interp >=ones(nb_sites,1)*SubRange(:,1,isub)' & Coord_interp<=ones(nb_sites,1)*SubRange(:,2,isub)');
            ind_sel=find(sum(check_range,2)==nb_coord);
            nbval(ind_sel)=nbval(ind_sel)+1;% records the number of values for eacn interpolation point (in case of subdomain overlap)
            Coord_tps=DataIn.(DataIn.ListVarName{VarType{icell}.coord_tps});
            if check_grid
                EM = tps_eval(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the velocity from tps 'sources'
            end
            if check_der
                [EMDX,EMDY] = tps_eval_dxy(Coord_interp(ind_sel,:),Coord_tps(1:nbvec_sub,:,isub));%kernels for calculating the spatial derivatives from tps 'sources'
            end
            ListVar={};
            U_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_x_tps});
            V_tps=DataIn.(DataIn.ListVarName{VarType{icell}.vector_y_tps});
            for ilist=1:length(FieldList)
                var_count=numel(ListVar);
                switch FieldList{ilist}
                    case 'velocity'
                        ListVar=[ListVar {'U', 'V'}];
                        VarAttribute{var_count+1}.Role='vector_x';
                        VarAttribute{var_count+2}.Role='vector_y';
                        DataOut.U(ind_sel)=DataOut.U(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
                        DataOut.V(ind_sel)=DataOut.V(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
                    case 'u'
                        ListVar=[ListVar {'u'}];
                        VarAttribute{var_count+1}.Role='scalar';
                        DataOut.u(ind_sel)=DataOut.u(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
                    case 'v'
                        ListVar=[ListVar {'v'}];
                        VarAttribute{var_count+1}.Role='scalar';
                        DataOut.v(ind_sel)=DataOut.v(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
                    case 'norm_vel'
                        ListVar=[ListVar {'norm_vel'}];
                        VarAttribute{var_count+1}.Role='scalar';
                        U=DataOut.U(ind_sel)+EM *U_tps(1:nbvec_sub+3,isub);
                        V=DataOut.V(ind_sel)+EM *V_tps(1:nbvec_sub+3,isub);
                        DataOut.norm_vel(ind_sel)=sqrt(U.*U+V.*V);
                    case 'vort'
                        ListVar=[ListVar {'vort'}];
                        VarAttribute{var_count+1}.Role='scalar';
                        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);
                    case 'div'
                        ListVar=[ListVar {'div'}];
                        VarAttribute{var_count+1}.Role='scalar';
                        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);
                    case 'strain'
                        ListVar=[ListVar {'strain'}];
                        VarAttribute{var_count+1}.Role='scalar';
                        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);
                end
            end
        end
        DataOut.FF=nbval==0; %put errorflag to 1 for points outside the interpolation rang
        nbval(nbval==0)=1;% to avoid division by zero for averaging
        if isempty(find(strcmp('FF',DataOut.ListVarName),1))% if FF is not already listed
            DataOut.ListVarName=[DataOut.ListVarName {'FF'}];
            DataOut.VarDimName=[DataOut.VarDimName {{'coord_y','coord_x'}}];
            DataOut.VarAttribute{length(DataOut.ListVarName)}.Role='errorflag';
        end
        DataOut.ListVarName=[DataOut.ListVarName ListVar];
        for ifield=1:numel(ListVar)
            VarDimName{ifield}={'coord_y','coord_x'};
            DataOut.(ListVar{ifield})=DataOut.(ListVar{ifield})./nbval;
            DataOut.(ListVar{ifield})=reshape(DataOut.(ListVar{ifield}),npy,npx);
        end
        DataOut.FF=reshape(DataOut.FF,npy,npx);
        DataOut.VarDimName=[DataOut.VarDimName VarDimName];
        DataOut.VarAttribute=[DataOut.VarAttribute VarAttribute];
    end
else

    %% civx data
    DataOut=DataIn;
    for ilist=1:length(FieldList)
        if ~isempty(FieldList{ilist})
            [VarName,Value,Role,units]=feval(FieldList{ilist},DataIn);%calculate field with appropriate function named FieldList{ilist}
            ListVarName=[ListVarName VarName];
            ValueList=[ValueList Value];
            RoleList=[RoleList Role];
            units_cell=[units_cell units];
        end
    end
    %erase previous data (except coordinates)
    for ivar=nbcoord+1:length(DataOut.ListVarName)
        VarName=DataOut.ListVarName{ivar};
        DataOut=rmfield(DataOut,VarName);
    end
    DataOut.ListVarName=DataOut.ListVarName(1:nbcoord);
    if isfield(DataOut,'VarDimName')
        DataOut.VarDimName=DataOut.VarDimName(1:nbcoord);
    else
        errormsg='element .VarDimName missing in input data';
        return
    end
    DataOut.VarAttribute=DataOut.VarAttribute(1:nbcoord);
    %append new data
    DataOut.ListVarName=[DataOut.ListVarName ListVarName];
    for ivar=1:length(ListVarName)
        DataOut.VarDimName{nbcoord+ivar}=DataOut.VarDimName{1};
        DataOut.VarAttribute{nbcoord+ivar}.Role=RoleList{ivar};
        DataOut.VarAttribute{nbcoord+ivar}.units=units_cell{ivar};
        DataOut.(ListVarName{ivar})=ValueList{ivar};
    end
end



%%%%%%%%%%%%% velocity fieldn%%%%%%%%%%%%%%%%%%%%
function [VarName,ValCell,Role,units_cell]=velocity(DataIn)
VarName={};
ValCell={};
Role={};
units_cell={};
if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
    units=[DataIn.CoordUnit '/' DataIn.TimeUnit];
else
    units='pixel';
end
if isfield(DataIn,'U')
    VarName=[VarName {'U'}];
    ValCell=[ValCell {DataIn.U}];
    Role=[Role {'vector_x'}];
    units_cell=[units_cell {units}];
end
if isfield(DataIn,'V')
    VarName=[VarName {'V'}];
    ValCell=[ValCell {DataIn.V}];
    Role=[Role {'vector_y'}];
    units_cell=[units_cell {units}];
end
if isfield(DataIn,'W')
    VarName=[VarName {'W'}];
    ValCell=[ValCell {DataIn.W}];
    Role=[Role {'vector_z'}];
    units_cell=[units_cell {units}];
end
if isfield(DataIn,'F')
    VarName=[VarName {'F'}];
    ValCell=[ValCell {DataIn.F}];
    Role=[Role {'warnflag'}];
    units_cell=[units_cell {[]}];
end
if isfield(DataIn,'FF')
    VarName=[VarName,{'FF'}];
    ValCell=[ValCell {DataIn.FF}];
    Role=[Role {'errorflag'}];
    units_cell=[units_cell {[]}];
end

%%%%%%%%%%%%% ima cor%%%%%%%%%%%%%%%%%%%%
function [VarName,ValCell,Role,units]=ima_cor(DataIn)
VarName={};
ValCell={};
Role={};
units={};
if isfield(DataIn,'C')
    VarName{1}='C';
    ValCell{1}=DataIn.C;
    Role={'ancillary'};
    units={[]};
end

%%%%%%%%%%%%% norm_vec %%%%%%%%%%%%%%%%%%%%
function [VarName,ValCell,Role,units]=norm_vel(DataIn)
VarName={};
ValCell={};
Role={};
units={};
if isfield(DataIn,'U') && isfield(DataIn,'V')
    VarName{1}='norm_vel';
    ValCell{1}=DataIn.U.*DataIn.U+ DataIn.V.*DataIn.V;
    if isfield(DataIn,'W') && isequal(size(DataIn.W),size(DataIn.U))
        ValCell{1}=ValCell{1}+DataIn.W.*DataIn.W;
    end
    ValCell{1}=sqrt(ValCell{1});
    Role{1}='scalar';
    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
    else
        units={'pixel'};
    end
end



%%%%%%%%%%%%% vorticity%%%%%%%%%%%%%%%%%%%%
function [VarName,ValCell,Role,units]=vort(DataIn)
VarName={};
ValCell={};
Role={};
units={};
if isfield(DataIn,'DjUi')
    VarName{1}='vort';
    ValCell{1}=DataIn.DjUi(:,1,2)-DataIn.DjUi(:,2,1);  %vorticity
    siz=size(ValCell{1});
    ValCell{1}=reshape(ValCell{1},siz(1),1);
    Role{1}='scalar';
    if isfield(DataIn,'TimeUnit')
        units={[DataIn.TimeUnit '-1']};
    else
        units={[]};
    end
end

%%%%%%%%%%%%% divergence%%%%%%%%%%%%%%%%%%%%
function [VarName,ValCell,Role,units]=div(DataIn)
VarName={};
ValCell={};
Role={};
units={};
if isfield(DataIn,'DjUi')
    VarName{1}='div';
    ValCell{1}=DataIn.DjUi(:,1,1)+DataIn.DjUi(:,2,2); %DUDX+DVDY
    siz=size(ValCell{1});
    ValCell{1}=reshape(ValCell{1},siz(1),1);
    Role{1}='scalar';
    if isfield(DataIn,'TimeUnit')
        units={[DataIn.TimeUnit '-1']};
    else
        units={[]};
    end
end

%%%%%%%%%%%%% strain %%%%%%%%%%%%%%%%%%%%
function [VarName,ValCell,Role,units]=strain(DataIn)
VarName={};
ValCell={};
Role={};
units={};
if isfield(DataIn,'DjUi')
    VarName{1}='strain';
    ValCell{1}=DataIn.DjUi(:,1,2)+DataIn.DjUi(:,2,1);%DVDX+DUDY
    siz=size(ValCell{1});
    ValCell{1}=reshape(ValCell{1},siz(1),1);
    Role{1}='scalar';
    if isfield(DataIn,'TimeUnit')
        units={[DataIn.TimeUnit '-1']};
    else
        units={[]};
    end
end

%%%%%%%%%%%%% u %%%%%%%%%%%%%%%%%%%%
function [VarName,ValCell,Role,units]=u(DataIn)
VarName={};
ValCell={};
Role={};
units={};
if isfield(DataIn,'U')
    VarName{1}='U';
    ValCell{1}=DataIn.U;
    Role{1}='scalar';
    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
    else
        units={'pixel'};
    end
end

%%%%%%%%%%%%% v %%%%%%%%%%%%%%%%%%%%
function [VarName,ValCell,Role,units]=v(DataIn)
VarName={};
ValCell={};
Role={};
units={};
if isfield(DataIn,'V')
    VarName{1}='V';
    ValCell{1}=DataIn.V;
    Role{1}='scalar';
    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
    else
        units={'pixel'};
    end
end

%%%%%%%%%%%%% w %%%%%%%%%%%%%%%%%%%%
function [VarName,ValCell,Role,units]=w(DataIn)
VarName={};
ValCell={};
Role={};
units={};
if isfield(DataIn,'W')
    VarName{1}='W';
    ValCell{1}=DataIn.W;
    Role{1}='scalar';%will remain unchanged by projection
    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
    else
        units={'pixel'};
    end
end

%%%%%%%%%%%%% w_normal %%%%%%%%%%%%%%%%%%%%
function [VarName,ValCell,Role,units]=w_normal(DataIn)
VarName={};
ValCell={};
Role={};
units={};
if isfield(DataIn,'W')
    VarName{1}='W';
    ValCell{1}=DataIn.W;
    Role{1}='vector_z';%will behave like a vector component  by projection
    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
    else
        units={'pixel'};
    end
end

%%%%%%%%%%%%% error %%%%%%%%%%%%%%%%%%%%
function [VarName,ValCell,Role,units]=error(DataIn)
VarName={};
ValCell={};
Role={};
units={};
if isfield(DataIn,'E')
    VarName{1}='E';
    ValCell{1}=DataIn.E;
    Role{1}='ancillary'; %TODO CHECK units in actual fields
    if isfield(DataIn,'CoordUnit') && isfield(DataIn,'TimeUnit')
        units={[DataIn.CoordUnit '/' DataIn.TimeUnit]};
    else
        units={'pixel'};
    end
end