Changeset 364 for trunk/src/civ_matlab.m


Ignore:
Timestamp:
Jan 10, 2012, 5:06:30 PM (12 years ago)
Author:
sommeria
Message:

civ2, fix2, patch2 introduced under Matlab (no interpolation and rotation done yet)
small corrections in plot_field and read_civdata

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/civ_matlab.m

    r363 r364  
    5353if isfield (Param,'Civ1')
    5454    check_civ1=1;% test for further use of civ1 results
     55    % %% prepare images
     56    par_civ1=Param.Civ1;
     57    if isfield(par_civ1,'reverse_pair')
     58        if par_civ1.reverse_pair
     59            if ischar(par_civ1.ImageB)
     60                temp=par_civ1.ImageA;
     61                par_civ1.ImageA=imread(par_civ1.ImageB);
     62            end
     63            if ischar(temp)
     64                par_civ1.ImageB=imread(temp);
     65            end
     66        end
     67    else
     68        if ischar(par_civ1.ImageA)
     69            par_civ1.ImageA=imread(par_civ1.ImageA);
     70        end
     71        if ischar(par_civ1.ImageB)
     72            par_civ1.ImageB=imread(par_civ1.ImageB);
     73        end
     74    end
     75   
    5576    % caluclate velocity data (y and v in indices, reverse to y component)
    56     [xtable ytable utable vtable ctable F result_conv errormsg] = civ (Param.Civ1);
    57 
    58    % to try the reverse_pair method, uncomment below
    59 %     [xtable1 ytable1 utable1 vtable1 ctable1 F1 result_conv1 errormsg1] = civ (Param.Civ1);
    60 %     Param.Civ1.reverse_pair=1;
    61 %     [xtable2 ytable2 utable2 vtable2 ctable2 F2 result_conv2 errormsg2] = civ (Param.Civ1);
    62 %     xtable=[xtable1; xtable2];
    63 %     ytable=[ytable1; ytable2];
    64 %     utable=[utable1; -utable2];
    65 %     vtable=[vtable1; -vtable2];
    66 %     ctable=[ctable1; ctable2];
    67 %     F=[F1; F2];
    68 %     result_conv=[result_conv1; result_conv2];
    69 %     errormsg=[errormsg1; errormsg2];
     77    [xtable ytable utable vtable ctable F result_conv errormsg] = civ (par_civ1);
    7078   
    71    
     79    % to try the reverse_pair method, uncomment below
     80    %     [xtable1 ytable1 utable1 vtable1 ctable1 F1 result_conv1 errormsg1] = civ (Param.Civ1);
     81    %     Param.Civ1.reverse_pair=1;
     82    %     [xtable2 ytable2 utable2 vtable2 ctable2 F2 result_conv2 errormsg2] = civ (Param.Civ1);
     83    %     xtable=[xtable1; xtable2];
     84    %     ytable=[ytable1; ytable2];
     85    %     utable=[utable1; -utable2];
     86    %     vtable=[vtable1; -vtable2];
     87    %     ctable=[ctable1; ctable2];
     88    %     F=[F1; F2];
     89    %     result_conv=[result_conv1; result_conv2];
     90    %     errormsg=[errormsg1; errormsg2];
    7291    if ~isempty(errormsg)
    7392        return
    … …  
    8099    end
    81100    Data.ListGlobalAttribute=[Data.ListGlobalAttribute Civ1_param];% {'Civ1_Time','Civ1_Dt'}];
    82     Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_C','Civ1_F'};%  cell array containing the names of the fields to record
     101    Data.ListVarName={'Civ1_X','Civ1_Y','Civ1_U','Civ1_V','Civ1_F','Civ1_C'};%  cell array containing the names of the fields to record
    83102    Data.VarDimName={'NbVec1','NbVec1','NbVec1','NbVec1','NbVec1','NbVec1'};
    84103    Data.VarAttribute{1}.Role='coord_x';
    … …  
    93112    Data.Civ1_C=reshape(ctable,[],1);
    94113    Data.Civ1_F=reshape(F,[],1);
    95     Data.CivStage=1;
    96    
     114    Data.CivStage=1; 
    97115else
    98116    if exist('ncfile','var')
    … …  
    111129    else
    112130        Data=nc2struct(CivFile);%read civ1 and fix1 data in the existing netcdf file
    113 %         if isfield(Param,'Fix1')
    114 %             Data=nc2struct(CivFile,ListVarCiv1);%read civ1 data in the existing netcdf file
    115 %         else
    116 %             Data=nc2struct(CivFile,ListVarFix1);%read civ1 and fix1 data in the existing netcdf file
    117 %         end
    118131    end
    119132end
    … …  
    181194    par_civ2=Param.Civ2;
    182195    if ~isfield (Param,'Civ1') || ~strcmp(Param.Civ1.ImageA,par_civ2.ImageA)
    183         par_civ2.ImageA=imread(par_civ2.ImageA);%read first image if not already done for civ1
     196        par_civ2.ImageA=imread(Param.Civ2.ImageA);%read first image if not already done for civ1
    184197    else
    185         par_civ2.ImageA=Param.Civ1.ImageA;
     198        par_civ2.ImageA=par_civ1.ImageA;
    186199    end
    187200    if ~isfield (Param,'Civ1') || ~strcmp(Param.Civ1.ImageB,par_civ2.ImageB)
    188         par_civ2.ImageB=imread(par_civ2.ImageB);%read second image if not already done for civ1
     201        par_civ2.ImageB=imread(Param.Civ2.ImageB);%read second image if not already done for civ1
    189202         else
    190         par_civ2.ImageB=Param.Civ1.ImageB;
     203        par_civ2.ImageB=par_civ1.ImageB;
    191204    end
    192205    ibx2=ceil(par_civ2.Bx/2);
    … …  
    269282    Data.Civ2_Time=str2double(par_civ2.Time);
    270283    Data.Civ2_Dt=str2double(par_civ2.Dt);
    271     Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_C','Civ2_F'}];%  cell array containing the names of the fields to record
     284    nbvar=numel(Data.ListVarName);
     285    Data.ListVarName=[Data.ListVarName {'Civ2_X','Civ2_Y','Civ2_U','Civ2_V','Civ2_F','Civ2_C'}];%  cell array containing the names of the fields to record
    272286    Data.VarDimName=[Data.VarDimName {'NbVec2','NbVec2','NbVec2','NbVec2','NbVec2','NbVec2'}];
    273     Data.VarAttribute{1}.Role='coord_x';
    274     Data.VarAttribute{2}.Role='coord_y';
    275     Data.VarAttribute{3}.Role='vector_x';
    276     Data.VarAttribute{4}.Role='vector_y';
    277     Data.VarAttribute{5}.Role='warnflag';
     287    Data.VarAttribute{nbvar+1}.Role='coord_x';
     288    Data.VarAttribute{nbvar+2}.Role='coord_y';
     289    Data.VarAttribute{nbvar+3}.Role='vector_x';
     290    Data.VarAttribute{nbvar+4}.Role='vector_y';
     291    Data.VarAttribute{nbvar+5}.Role='warnflag';
    278292    Data.Civ2_X=reshape(xtable,[],1);
    279293    Data.Civ2_Y=reshape(size(par_civ2.ImageA,1)-ytable+1,[],1);
    … …  
    290304    for ilist=1:length(ListFixParam)
    291305        ParamName=ListFixParam{ilist};
    292         ListName=['Fix1_' ParamName];
     306        ListName=['Fix2_' ParamName];
    293307        eval(['Data.ListGlobalAttribute=[Data.ListGlobalAttribute ''' ParamName '''];'])
    294308        eval(['Data.' ListName '=Param.Fix2.' ParamName ';'])
    … …  
    308322        Data.VarAttribute{nbvar}.Role='errorflag';   
    309323        Data.Civ2_FF=fix(Param.Fix2,Data.Civ2_F,Data.Civ2_C,Data.Civ2_U,Data.Civ2_V);
    310         Data.CivStage=5;   
     324        Data.CivStage=Data.CivStage+1;   
    311325    end
    312326   
    … …  
    316330if isfield (Param,'Patch2')
    317331    Data.ListGlobalAttribute=[Data.ListGlobalAttribute {'Patch2_Rho','Patch2_Threshold','Patch2_SubDomain'}];
    318     Data.Patch2_Rho=str2double(Param.Patch2.Rho);
    319     Data.Patch2_Threshold=str2double(Param.Patch2.Threshold);
    320     Data.Patch2_SubDomain=str2double(Param.Patch2.SubDomain);
     332    Data.Patch2_Rho=Param.Patch2.SmoothingParam;
     333    Data.Patch2_Threshold=Param.Patch2.MaxDiff;
     334    Data.Patch2_SubDomain=Param.Patch2.SubdomainSize;
    321335    Data.ListVarName=[Data.ListVarName {'Civ2_U_Diff','Civ2_V_Diff','Civ2_X_SubRange','Civ2_Y_SubRange','Civ2_NbSites','Civ2_X_tps','Civ2_Y_tps','Civ2_U_tps','Civ2_V_tps'}];
    322336    Data.VarDimName=[Data.VarDimName {'NbVec2','NbVec2',{'NbSubDomain2','Two'},{'NbSubDomain2','Two'},'NbSubDomain2',...
    … …  
    337351      Data.Civ2_V_Diff(ind_good)=Data.Civ2_V(ind_good)-Vres;
    338352      Data.Civ2_FF(ind_good)=FFres;
    339       Data.CivStage=6;                             
    340 end   
     353      Data.CivStage=Data.CivStage+1;                             
     354end 
    341355
    342356%% write result in a netcdf file if requested
    … …  
    379393shifty=-par_civ.Shifty;% sign minus because image j index increases when y decreases
    380394if isfield(par_civ,'Grid')
    381     if ischar(par_civ.Grid)
     395    if ischar(par_civ.Grid)%read the drid file if the input is a file name
    382396        par_civ.Grid=dlmread(par_civ.Grid);
    383397        par_civ.Grid(1,:)=[];%the first line must be removed (heading in the grid file)
    … …  
    422436check_MaxIma=isfield(par_civ,'MaxIma') && ~isempty(par_civ.MaxIma);
    423437
    424 %% prepare images
    425 if isfield(par_civ,'reverse_pair')
    426     if par_civ.reverse_pair
    427         if ischar(par_civ.ImageB)
    428             temp=par_civ.ImageA;
    429             par_civ.ImageA=imread(par_civ.ImageB);
    430         end
    431         if ischar(temp)
    432             par_civ.ImageB=imread(temp);
    433         end
    434     end
    435 else
    436     if ischar(par_civ.ImageA)
    437         par_civ.ImageA=imread(par_civ.ImageA);
    438     end
    439     if ischar(par_civ.ImageB)
    440         par_civ.ImageB=imread(par_civ.ImageB);
    441     end
    442 end
     438% %% prepare images
     439% if isfield(par_civ,'reverse_pair')
     440%     if par_civ.reverse_pair
     441%         if ischar(par_civ.ImageB)
     442%             temp=par_civ.ImageA;
     443%             par_civ.ImageA=imread(par_civ.ImageB);
     444%         end
     445%         if ischar(temp)
     446%             par_civ.ImageB=imread(temp);
     447%         end
     448%     end
     449% else
     450%     if ischar(par_civ.ImageA)
     451%         par_civ.ImageA=imread(par_civ.ImageA);
     452%     end
     453%     if ischar(par_civ.ImageB)
     454%         par_civ.ImageB=imread(par_civ.ImageB);
     455%     end
     456% end
    443457
    444458[npy_ima npx_ima]=size(par_civ.ImageA);
    … …  
    660674    end
    661675end
    662 return
    663 
    664 
    665 FF=double(FF);
    666676
    667677
Note: See TracChangeset for help on using the changeset viewer.