Changeset 1094


Ignore:
Timestamp:
Mar 25, 2021, 9:37:04 AM (3 years ago)
Author:
sommeria
Message:

phys_polar debugged

Location:
trunk/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/plot_field.m

    r1093 r1094  
    768768            test_interp_X=0; %default, regularly meshed X coordinate
    769769            test_interp_Y=0; %default, regularly meshed Y coordinate
    770             if isfield(Data,'VarAttribute')
    771                 if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end)},'units')
    772                     x_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end)}.units;
    773                 end
    774                 if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end-1) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)},'units')
    775                     y_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)}.units;
    776                 end
    777             end
     770%             if isfield(Data,'VarAttribute')
     771%                 if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end)},'units')
     772%                     x_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end)}.units;
     773%                 end
     774%                 if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end-1) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)},'units')
     775%                     y_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)}.units;
     776%                 end
     777%             end
    778778            if numel(Coord_y)>2
    779779                DCoord_y=diff(Coord_y);
    … …  
    790790                DCoord_x_min=min(DCoord_x);
    791791                DCoord_x_max=max(DCoord_x);
    792                 if sign(DCoord_x_min)~=sign(DCoord_x_max);% =1 for increasing values, 0 otherwise
     792                if sign(DCoord_x_min)~=sign(DCoord_x_max)% =1 for increasing values, 0 otherwise
    793793                    errormsg=['errror in plot_field.m: non monotonic dimension variable ' Data.ListVarName{VarRole.coord(2)} ];
    794794                    return
    … …  
    821821    end
    822822    %define coordinates as CoordUnits, if not defined as attribute for each variable
     823%     if isfield(Data,'VarAttribute')&& numel(Data.VarAttribute)>=1 && isfield(Data.VarAttribute{1},'unit')
     824%         y_units=Data.VarAttribute{1}.unit;
     825%     end
    823826    if isfield(Data,'CoordUnit')
    824827        if isempty(x_units)
    … …  
    827830        if isempty(y_units)
    828831            y_units=Data.CoordUnit;
     832        end
     833    elseif isfield(Data,'VarAttribute')
     834        if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end)},'units')
     835            x_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end)}.units;
     836        end
     837        if numel(Data.VarAttribute)>=CellInfo{icell}.CoordIndex(end-1) && isfield(Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)},'units')
     838            y_units=Data.VarAttribute{CellInfo{icell}.CoordIndex(end-1)}.units;
    829839        end
    830840    end

  • trunk/src/proj_field.m

    r1093 r1094  
    11251125ProjMode=num2cell(blanks(numel(CellInfo)));
    11261126ProjMode=regexprep(ProjMode,' ',ObjectData.ProjMode);
    1127 %ProjMode=cell(size(CellInfo));
    1128 % for icell=1:numel(CellInfo)
    1129 %     ProjMode{icell}=ObjectData.ProjMode;% projection mode of the plane object
    1130 % end
    11311127icell_grid=[];% field cell index which defines the grid
    11321128if strcmp(ObjectData.ProjMode,'projection')
    … …  
    11921188        [XI,YI]=meshgrid(coord_x_proj,coord_y_proj);%grid in the new coordinates
    11931189        ProjData.VarDimName={AYName,AXName};
    1194         %         XI=ObjectData.Coord(1,1)+(X)*cos(PlaneAngle(3))-YI*sin(PlaneAngle(3));%corresponding coordinates in the original system
    1195         %         YI=ObjectData.Coord(1,2)+(X)*sin(PlaneAngle(3))+YI*cos(PlaneAngle(3));
    11961190    else% we use the existing grid from field cell #icell_grid
    11971191        NbDim=NbDimArray(icell_grid);
    … …  
    12041198        ProjData.(AXName)=FieldData.(AXName); % new (projected ) y coordinates
    12051199    end
    1206     ProjData.ListVarName={AYName,AXName};
    1207    
     1200    ProjData.ListVarName={AYName,AXName};   
    12081201    ProjData.VarAttribute{1}.Role='coord_y';
    12091202    ProjData.VarAttribute{2}.Role='coord_x';
     1203            YAttribute=FieldData.VarAttribute{CellInfo{icell_grid}.CoordIndex(NbDim-1)};
     1204        XAttribute=FieldData.VarAttribute{CellInfo{icell_grid}.CoordIndex(NbDim)};
     1205    if ~testangle
     1206        if isfield(YAttribute,'units')
     1207            ProjData.VarAttribute{1}.units=YAttribute.units;
     1208        end
     1209         if isfield(XAttribute,'units')
     1210            ProjData.VarAttribute{2}.units=XAttribute.units;
     1211         end
     1212    end
    12101213end
    12111214
    … …  
    12441247            coord_x=FieldData.(CellInfo{icell}.XName);% initial x coordinates
    12451248            coord_y=FieldData.(CellInfo{icell}.YName);% initial y coordinates
    1246      
     1249           
    12471250            if check3D
    12481251                coord_z=FieldData.(CellInfo{icell}.ZName);
    … …  
    13481351                        VarName_FF=FieldData.ListVarName{CellInfo{icell}.VarIndex_errorflag};
    13491352                        indsel=find(FieldData.(VarName_FF)==0);
     1353                        if isempty(indsel)
     1354                            errormsg='bad projection plane: no data found in the projection domain';
     1355                            return
     1356                        end
    13501357                        coord_X=coord_X(indsel);
    13511358                        coord_Y=coord_Y(indsel);
    … …  
    16251632                    end
    16261633                    [X,Y]=meshgrid(Coord{2},Coord{1});%initial coordinates
    1627                     %name of error flag variable
    1628 %                     FFName='FF';%default name (if not already used)
    1629 %                     if isfield(ProjData,'FF')
    1630 %                         ind=1;
    1631 %                         while isfield(ProjData,['FF_' num2str(ind)])
    1632 %                             ind=ind+1;
    1633 %                         end
    1634 %                         FFName=['FF_' num2str(ind)];% append an index to the name of error flag, FF_1,FF_2...
    1635 %                     end
     1634
    16361635                    % project all variables in the cell
    16371636                    for ivar=VarIndex
    … …  
    16561655                            VarAttribute{length(ListVarName)+nbcoord}=FieldData.VarAttribute{ivar};
    16571656                        end
    1658 %                         ProjData.(FFName)=isnan(ProjData.(VarName));%detact NaN (points outside the interpolation range)
    1659 %                         ProjData.(VarName)(ProjData.(FFName))=0; %set to 0 the NaN data
    1660                     end
    1661                     %update list of variables with error flag
    1662 %                     ListVarName=[ListVarName FFName];
    1663 %                     VarDimName=[VarDimName {DimCell}];
    1664 %                     VarAttribute{numel(ListVarName)}.Role='errorflag';
     1657                    end
    16651658                elseif ~testangle % 3Dcase without change of angle
    16661659                    % unstructured z coordinate
    … …  
    16911684                    Coord{2}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(2)});%initial y coordinates
    16921685                    Coord{3}=FieldData.(FieldData.ListVarName{CellInfo{icell}.CoordIndex(3)});%initial x coordinates
    1693 
    16941686                    coord_x_proj=ObjectData.RangeX(1):InterpMesh:ObjectData.RangeX(2);% set of coordinates in the projection plane
    16951687                    coord_y_proj=ObjectData.RangeY(1):InterpMesh:ObjectData.RangeY(2);
    … …  
    17001692                    YI_proj=M(2,1)*XI+M(2,2)*YI+ObjectData.Coord(1,2);
    17011693                    ZI_proj=M(3,1)*XI+M(3,2)*YI+ObjectData.Coord(1,3);
    1702                    
    1703    
    17041694                    for ivar=VarIndex
    17051695                        VarName=FieldData.ListVarName{ivar};
    … …  
    17631753    end
    17641754end
    1765 % %prepare substraction in case of two input fields
    1766 % SubData.ListVarName={};
    1767 % SubData.VarDimName={};
    1768 % SubData.VarAttribute={};
    1769 % check_remove=zeros(size(ProjData.ListVarName));
    1770 % for iproj=1:numel(ProjData.VarAttribute)
    1771 %     if isfield(ProjData.VarAttribute{iproj},'CheckSub')&&isequal(ProjData.VarAttribute{iproj}.CheckSub,1)
    1772 %         VarName=ProjData.ListVarName{iproj};
    1773 %         SubData.ListVarName=[SubData.ListVarName {VarName}];
    1774 %         SubData.VarDimName=[SubData.VarDimName ProjData.VarDimName{iproj}];
    1775 %         SubData.VarAttribute=[SubData.VarAttribute ProjData.VarAttribute{iproj}];
    1776 %         SubData.(VarName)=ProjData.(VarName);
    1777 %         check_remove(iproj)=1;
    1778 %     end
    1779 % end
    1780 % if ~isempty(find(check_remove))
    1781 %     ind_remove=find(check_remove);
    1782 %     ProjData.ListVarName(ind_remove)=[];
    1783 %     ProjData.VarDimName(ind_remove)=[];
    1784 %     ProjData.VarAttribute(ind_remove)=[];
    1785 %     ProjData=sub_field(ProjData,[],SubData);
    1786 % end
    17871755
    17881756%-----------------------------------------------------------------

  • trunk/src/series/merge_proj.m

    r1093 r1094  
    253253end
    254254OutputPath=fullfile(Param.OutputPath,Param.Experiment,Param.Device);
     255
    255256for index=1:NbField
    256257        update_waitbar(WaitbarHandle,index/NbField)

  • trunk/src/series/sliding_average.m

    r1093 r1094  
    2525%             .RUN =0 for GUI input, =1 for function activation
    2626%             .RunMode='local','background', 'cluster': type of function  use
    27 %             
    28 %    .IndexRange: set the file or frame indices on which the action must be performed
    29 %    .FieldTransform: .TransformName: name of the selected transform function
     27%             900
     28%    .IndexRange: set the file or frame indices on which the action must be performseriesed
     29%    .FieldTransform: .TransformName: name of the select39ed transform function
    3030%                     .TransformPath:   path  of the selected transform function
    3131%    .InputFields: sub structure describing the input fields withfields
    … …  
    3434%              .FieldName_1: name of the second field in case of two input series
    3535%              .VelType_1: velocity type of the second field in case of two input series
    36 %              .Coord_y: name of y coordinate variable
     36%             uvmat .Coord_y: name of y coordinate variable
    3737%              .Coord_x: name of x coordinate variable
    3838%    .ProjObject: %sub structure describing a projection object (read from ancillary GUI set_object)
    3939%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    40 
    4140%=======================================================================
    4241% Copyright 2008-2021, LEGI UMR 5519 / CNRS UGA G-INP, Grenoble, France
    … …  
    6261if isstruct(Param) && isequal(Param.Action.RUN,0)
    6362    ParamOut.AllowInputSort='off';% allow alphabetic sorting of the list of input file SubDir (options 'off'/'on', 'off' by default)
    64     ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to max (options 'off'/'on', 'off' by default)
     63    ParamOut.WholeIndexRange='off';% prescribes the file index ranges from min to mseriesax (options 'off'/'on', 'off' by default)
    6564    ParamOut.NbSlice=1; %nbre of slices ('off' by default)
    6665    ParamOut.VelType='off';% menu for selecting the velocity type (options 'off'/'one'/'two',  'off' by default)
    6766    ParamOut.FieldName='one';% menu for selecting the field (s) in the input file(options 'off'/'one'/'two', 'off' by default)
    6867    ParamOut.FieldTransform = 'on';%can use a transform function
    69     ParamOut.ProjObject='off';%can use projection object(option 'off'/'on',
     68    ParamOut.ProjObject='off';%can use projection object39(option 'off'/'on',
    7069    ParamOut.Mask='off';%can use mask option   (option 'off'/'on', 'off' by default)
    7170    ParamOut.OutputDirExt='.tfilter';%set the output dir extension
    … …  
    122121        return
    123122    end
    124     [FileInfo{iview},MovieObject{iview}]=get_file_info(filecell{iview,1});
     123    [FileInfo{iview},MovieObject{iview}]=get_file_info(filecell{iview,1});900
    125124    FileType{iview}=FileInfo{iview}.FileType;
    126125    CheckImage{iview}=~isempty(find(strcmp(FileType{iview},ImageTypeOptions)));% =1 for images
    … …  
    138137if size(time,1)>1
    139138    diff_time=max(max(diff(time)));
    140     if diff_time>0
     139    if diff_time>0series
    141140        msgbox_uvmat('WARNING',['times of series differ by (max) ' num2str(diff_time)])
    142141    end   
    … …  
    177176
    178177%% Set field names and velocity types
    179 InputFields{1}=[];%default (case of images)
     178InputFields{1}=[];%default (case of images)series
    180179if isfield(Param,'InputFields')
    181180    InputFields{1}=Param.InputFields;
    … …  
    187186%% initialisation
    188187T=24.2; %main wave period
     188t0=3; % time for motion start (torus at its maximum x)
    189189NbPeriod=2; %number of periods for the sliding average
    190190omega=2*pi/T;
     191amplitude=2.5; %oscillation amplitude
     192Lscale=15;%diameter of the torus, length scale for normalisation
     193Uscale=amplitude*omega;
    191194
    192195DataOut.ListGlobalAttribute= {'Conventions','Time'};
    193196DataOut.Conventions='uvmat';
    194 DataOut.ListVarName={'coord_y','coord_x','Umean','Vmean','Ucos','Vcos','DUDXsin','DUDYsin','DVDXsin','DVDYsin','Ustokes','Vstokes'};
     197DataOut.ListVarName={'coord_y','coord_x','Umean','Vmean','Ucos','Vcos','DUDXsin','DUDXcos','DUDYsin','DVDXsin','DVDXcos'...
     198    ,'DVDYsin','Ustokes','Vstokes'};
    195199DataOut.VarDimName={'coord_y','coord_x',{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},...
    196     {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};
     200    {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'},...
     201    {'coord_y','coord_x'},{'coord_y','coord_x'},{'coord_y','coord_x'}};
    197202
    198203%%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
    … …  
    207212Time_end=Data.Time;
    208213dt=(Time_end-Time_1)/(NbField-1); %time interval
    209 NpTime=round(NbPeriod*T/dt+1)
    210 
     214NpTime=round(NbPeriod*T/dt+1);
     215
     216OutputPath=fullfile(Param.OutputPath,Param.Experiment,Param.Device);
     217RootFileOut=RootFile{1};
     218NomTypeOut='_1';
    211219%%%%%%%%%%%%%%%% loop on field indices %%%%%%%%%%%%%%%%
    212220disp('loop for filtering started')
    … …  
    219227    [Field,tild,errormsg] = read_field(filecell{1,index},FileType{iview},InputFields{iview},frame_index{iview}(index));
    220228   
    221     %%%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
     229    %%%%%%%%%%% MAIN RUNNING OPERATIONS  %%%%%%%%%%%%
    222230    if index==1 %first field
    223         DataOut.coord_x=Field.coord_x;
    224         DataOut.coord_y=Field.coord_y;
     231        DataOut.coord_x=Field.coord_x/Lscale;
     232        DataOut.coord_y=Field.coord_y/Lscale;
    225233        npy=numel(DataOut.coord_y);
    226234        npx=numel(DataOut.coord_x);
    … …  
    229237        Ucos=zeros(NpTime,npy,npx);
    230238        Vcos=zeros(NpTime,npy,npx);
     239        DUDXcos=zeros(NpTime,npy,npx);
    231240        DUDXsin=zeros(NpTime,npy,npx);
    232241        DUDYsin=zeros(NpTime,npy,npx);
     242        DVDXcos=zeros(NpTime,npy,npx);
    233243        DVDXsin=zeros(NpTime,npy,npx);
    234244        DVDYsin=zeros(NpTime,npy,npx);
    235245    end
    236       Time(index)=Field.Time;
     246    Time(index)=Field.Time-t0;%time from the start of the motion
    237247    Umean=circshift(Umean,[-1 0 0]); %shift U by ishift along the first index
    238         Vmean=circshift(Vmean,[-1 0 0]); %shift U by ishift along the first index
    239         Ucos=circshift(Ucos,[-1 0 0]); %shift U by ishift along the first index
    240         Vcos=circshift(Vcos,[-1 0 0]); %shift U by ishift along the first index
    241         DUDXsin=circshift(DUDXsin,[-1 0 0]);
    242         DUDYsin=circshift(DUDYsin,[-1 0 0]);
    243         DVDXsin=circshift(DVDXsin,[-1 0 0]);
    244         DVDYsin=circshift(DVDYsin,[-1 0 0]);
    245         Umean(end,:,:)=Field.U;
    246         Vmean(end,:,:)=Field.V;
    247         Ucos(end,:,:)=Field.U*cos(omega*Time(index));
    248         Vcos(end,:,:)=Field.V*cos(omega*Time(index));
    249         DUDXsin(end,:,:)=Field.DUDX*sin(omega*Time(index));
    250         DUDYsin(end,:,:)=Field.DUDY*sin(omega*Time(index));
    251         DVDXsin(end,:,:)=Field.DVDX*sin(omega*Time(index));
    252         DVDYsin(end,:,:)=Field.DVDY*sin(omega*Time(index));
    253         DataOut.Time=Time(index)-(NpTime-1)*dt/2;
    254         DataOut.Umean=squeeze(nanmean(Umean,1));
    255         DataOut.Vmean=squeeze(nanmean(Vmean,1));
    256         DataOut.Ucos=2*squeeze(nanmean(Ucos,1));
    257         DataOut.Vcos=2*squeeze(nanmean(Vcos,1));
    258         DataOut.DUDXsin=2*squeeze(nanmean(DUDXsin,1));
    259         DataOut.DUDYsin=2*squeeze(nanmean(DUDYsin,1));
    260         DataOut.DVDXsin=2*squeeze(nanmean(DVDXsin,1));
    261         DataOut.DVDYsin=2*squeeze(nanmean(DVDYsin,1));
    262         DataOut.Ustokes=(1/omega)*(DataOut.Ucos.*DataOut.DUDXsin+DataOut.Vcos.*DataOut.DUDYsin);
    263         DataOut.Vstokes=(1/omega)*(DataOut.Ucos.*DataOut.DVDXsin+DataOut.Vcos.*DataOut.DVDYsin);
     248    Vmean=circshift(Vmean,[-1 0 0]); %shift U by ishift along the first index
     249    Ucos=circshift(Ucos,[-1 0 0]); %shift U by ishift along the first index
     250    Vcos=circshift(Vcos,[-1 0 0]); %shift U by ishift along the first index
     251    DUDXcos=circshift(DUDXcos,[-1 0 0]);
     252    DUDXsin=circshift(DUDXsin,[-1 0 0]);
     253    DUDYsin=circshift(DUDYsin,[-1 0 0]);       
     254    DVDXcos=circshift(DVDXcos,[-1 0 0]);
     255    DVDXsin=circshift(DVDXsin,[-1 0 0]);
     256    DVDYsin=circshift(DVDYsin,[-1 0 0]);       
     257    Umean(end,:,:)=Field.U;
     258    Vmean(end,:,:)=Field.V;
     259    Ucos(end,:,:)=Field.U*cos(omega*Time(index));
     260    Vcos(end,:,:)=Field.V*cos(omega*Time(index));
     261    DUDXcos(end,:,:)=Field.DUDX*cos(omega*Time(index));
     262    DUDXsin(end,:,:)=Field.DUDX*sin(omega*Time(index));
     263    DUDYsin(end,:,:)=Field.DUDY*sin(omega*Time(index));% ParamOut=[];%default output
     264
     265    DVDXcos(end,:,:)=Field.DVDX*cos(omega*Time(index));
     266    DVDXsin(end,:,:)=Field.DVDX*sin(omega*Time(index));
     267    DVDYsin(end,:,:)=Field.DVDY*sin(omega*Time(index));
     268    DataOut.Time=(Time(index)-(NpTime-1)*dt/2)/T;%time inperiods from the beginning of the oscillation (torus at max abscissa)
     269    DataOut.Umean=(1/Uscale)*squeeze(nanmean(Umean,1));
     270    DataOut.Vmean=(1/Uscale)*squeeze(nanmean(Vmean,1));
     271    DataOut.Ucos=2*squeeze(nanmean(Ucos,1));
     272    DataOut.Vcos=2*squeeze(nanmean(Vcos,1));
     273    DataOut.DUDXcos=2*squeeze(nanmean(DUDXcos,1));
     274    DataOut.DUDXsin=2*squeeze(nanmean(DUDXsin,1));
     275    DataOut.DUDYsin=2*squeeze(nanmean(DUDYsin,1));
     276    DataOut.DVDXcos=2*squeeze(nanmean(DVDXcos,1));
     277    DataOut.DVDXsin=2*squeeze(nanmean(DVDXsin,1));
     278    DataOut.DVDYsin=2*squeeze(nanmean(DVDYsin,1));
     279    DataOut.Ustokes=(1/omega)*(1/Uscale)*(DataOut.Ucos.*DataOut.DUDXsin+DataOut.Vcos.*DataOut.DUDYsin);
     280    DataOut.Vstokes=(1/omega)*(1/Uscale)*(DataOut.Ucos.*DataOut.DVDXsin+DataOut.Vcos.*DataOut.DVDYsin);
     281
    264282    % writing the result file as netcdf file
    265     if index-round(NpTime/2)>=1
    266         OutputFile=fullfile_uvmat(RootPath{1},OutputDir,RootFile{1},FileExtOut,NomType{1},index-round(NpTime/2));
    267         %case of netcdf input file , determine global attributes
    268         errormsg=struct2nc(OutputFile,DataOut); %save result file
    269         if isempty(errormsg)
    270             disp([OutputFile ' written']);
    271         else
    272             disp(['error in writting result file: ' errormsg])
    273         end
    274     end
    275 end
    276 %%%%%%%%%%%%%%%% end loop on field indices %%%%%%%%%%%%%%%%
    277 
     283    i1=i1_series{1}(index);
     284    OutputFile=fullfile_uvmat(OutputPath,OutputDir,RootFileOut,'.nc',NomTypeOut,i1);
     285    errormsg=struct2nc(OutputFile, DataOut);
     286    if isempty(errormsg)
     287        disp([OutputFile ' written'])
     288    else
     289        disp(errormsg)
     290    end
     291end
     292   

  • trunk/src/transform_field/phys_polar.m

    r1093 r1094  
    5252    dlg_title = 'set the parameters for the polar coordinates';
    5353    num_lines= 2;
    54     def     = { '[0 0]';'0';'0';'+'};
     54    def     = { '[0 0]';'';'0';'+'};
    5555    if isfield(XmlData,'TransformInput')
    5656        if isfield(XmlData.TransformInput,'PolarCentre')
    … …  
    8787DataCell{2}=[];%default
    8888checkpixel(1)=0;
    89 if isfield(DataCell{1},'CoorUnit')&& strcmp(DataCell{1}.CoorUnit,'px')
     89if isfield(DataCell{1},'CoordUnit')&& strcmp(DataCell{1}.CoordUnit,'pixel')
    9090    checkpixel(1)=1;
    9191end
    … …  
    101101if nargin==4% case of two input fields
    102102    checkpixel(2)=0;
    103 if isfield(DataCell{2},'CoorUnit')&& strcmp(DataCell{2}.CoorUnit,'px')
     103if isfield(DataCell{2},'CoordUnit')&& strcmp(DataCell{2}.CoordUnit,'pixel')
    104104    checkpixel(2)=1;
    105105end
    … …  
    124124        end
    125125    end
    126     if isfield(XmlData.TransformInput,'PolarReferenceRadius') && isnumeric(XmlData.TransformInput.PolarReferenceRadius)
     126    if isfield(XmlData.TransformInput,'PolarReferenceRadius') && ~isempty(XmlData.TransformInput.PolarReferenceRadius)
    127127        radius_offset=XmlData.TransformInput.PolarReferenceRadius;
    128128    end
    129129    if radius_offset > 0
    130130        angle_scale=radius_offset; %the azimuth is rescale in terms of the length along the reference radius
    131         check_degree=0; %the output has the same unit asthe input
     131        check_degree=0; %the output has the same unit as the input
    132132    else
    133133        angle_scale=180/pi; %polar angle in degrees
    … …  
    144144
    145145nbvar=0;%counter for the number of output variables
    146 nbcoord=0;%counter for the number of variables for radial coordiantes (case of multiple field inputs)
     146nbcoord=0;%counter for the number of variablecheck_degrees for radial coordiantes (case of multiple field inputs)
    147147nbgrid=0;%counter for the number of gridded fields (all linearly interpolated on the same output polar grid)
    148148nbscattered=0;%counter of scattered fields
    … …  
    157157        return
    158158    end
    159     % end
    160159    %transform of X,Y coordinates for vector fields
    161160    if isfield(DataCell{ifield},'ZIndex')&& ~isempty(DataCell{ifield}.ZIndex)
    … …  
    180179                Data.VarAttribute{nbvar-1}.Role='coord_x';
    181180                check_unit=1;
    182                 if isfield(DataCell{ifield},'CoordUnit')
    183                     Data=rmfield(Data,'CoordUnit');
    184                     Data.VarAttribute{nbvar-1}.unit=DataCell{ifield}.CoordUnit;
    185                 elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
    186                     Data.VarAttribute{nbvar-1}.unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
     181                %unit of output field
     182                if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
     183                    radius_unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
     184                elseif isfield(DataCell{ifield},'CoordUnit')
     185                    radius_unit=DataCell{ifield}.CoordUnit;
    187186                else
    188                     check_unit=0;
    189                 end
     187                    radius_unit='';
     188                end
     189                Data.VarAttribute{nbvar-1}.units=radius_unit;
     190                if check_degree
     191                     Data.VarAttribute{nbvar}.units='degree';
     192                else %case of a reference radius
     193                    Data.VarAttribute{nbvar}.units=radius_unit;
     194                    Data.CoordUnit=radius_unit;
     195                end
     196%                 if isfield(DataCell{ifield},'CoordUnit')
     197%                     Data=rmfield(Data,'CoordUnit');
     198%                     Data.VarAttribute{nbvar-1}.unit=DataCell{ifield}.CoordUnit;
     199%                 elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
     200%                     Data.VarAttribute{nbvar-1}.unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
     201%                 else
     202%                     check_unit=0;
     203%                 end
    190204                Data.VarAttribute{nbvar}.Role='coord_y';
    191                 if check_degree
    192                 Data.VarAttribute{nbvar}.unit='degree';
    193                 elseif check_unit
    194                     Data.VarAttribute{nbvar}.unit=Data.VarAttribute{nbvar-1}.unit;
    195                 end
     205%                 if check_degree
     206%                 Data.VarAttribute{nbvar}.units='degree';
     207%                 elseif check_unit
     208%                     Data.VarAttribute{nbvar}.units=Data.VarAttribute{nbvar-1}.units;
     209%                 end
    196210 
    197211                %transform u,v into polar coordinates
    … …  
    250264                if nbgrid==0% no gridded data yet, introduce the coordinate variables common to all gridded data
    251265                    nbcoord=nbcoord+1;%add new radial coordinates for the first gridded field
    252                     radius_name = rename_indexing(radius_name,Data.ListVarName);
    253                     theta_name = rename_indexing(theta_name,Data.ListVarName);
    254                     Data.ListVarName = [Data.ListVarName {radius_name} {theta_name}];
     266                    radius_name = rename_indexing(radius_name,Data.ListVarName);% add an index to the name, or increment an existing index,
     267                    theta_name = rename_indexing(theta_name,Data.ListVarName);% if the proposed Name already exists in the list
     268                    Data.ListVarName = [Data.ListVarName {radius_name} {theta_name}];%add polar coordinates to the list of variables
    255269                    Data.VarDimName=[Data.VarDimName {radius_name} {theta_name}];
    256270                    nbvar=nbvar+2;
    257271                    if check_reverse
    258                                             Data.VarAttribute{nbvar-1}.Role='coord_y';
    259                     Data.VarAttribute{nbvar}.Role='coord_x';
     272                        Data.VarAttribute{nbvar-1}.Role='coord_y';
     273                        Data.VarAttribute{nbvar}.Role='coord_x';
    260274                    else
    261                     Data.VarAttribute{nbvar-1}.Role='coord_x';
    262                     Data.VarAttribute{nbvar}.Role='coord_y';
     275                        Data.VarAttribute{nbvar-1}.Role='coord_x';
     276                        Data.VarAttribute{nbvar}.Role='coord_y';
    263277                    end
    264278                    check_unit=1;
    265                     if isfield(DataCell{ifield},'CoordUnit')
    266                         Data.VarAttribute{nbvar-1}.unit=DataCell{ifield}.CoordUnit;
    267                     elseif isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
    268                         Data.VarAttribute{nbvar-1}.unit=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
     279
     280                    if isfield(XmlData,'GeometryCalib')&& isfield(XmlData.GeometryCalib,'CoordUnit')
     281                        Data.VarAttribute{nbvar-1}.units=XmlData.GeometryCalib.CoordUnit;% states that the output is in unit defined by GeometryCalib, then erased all projection objects with different units
     282                    elseif isfield(DataCell{ifield},'CoordUnit')
     283                        Data.VarAttribute{nbvar-1}.units=DataCell{ifield}.CoordUnit;%radius in coord units
    269284                    else
    270285                        check_unit=0;
    271286                    end
    272287                    if check_degree
    273                         Data.VarAttribute{nbvar}.unit='degree';
     288                        Data.VarAttribute{nbvar}.units='degree';%angle in degree
    274289                    elseif check_unit
    275                         Data.VarAttribute{nbvar}.unit=Data.VarAttribute{nbvar-1}.unit;
     290                        Data.VarAttribute{nbvar}.units=Data.VarAttribute{nbvar-1}.units;% angle in coord unit (normalised by reference radiuss)
    276291                    end
    277292                end
    … …  
    282297                    FieldName{nbgrid}=DataCell{ifield}.ListVarName{CellInfo{icell}.VarIndex_scalar};
    283298                    A{nbgrid}=DataCell{ifield}.(FieldName{nbgrid});
    284 %                     Data.ListVarName=[Data.ListVarName {FieldName{nbgrid}}];
    285 %                     Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}}];
    286299                    nbpoint(nbgrid)=numel(A{nbgrid});
    287300                    check_scalar(nbgrid)=1;
    … …  
    296309                    A{nbgrid+1}=DataCell{ifield}.(FieldName{nbgrid+1});
    297310                    A{nbgrid+2}=DataCell{ifield}.(FieldName{nbgrid+2});
    298                    % Data.ListVarName=[Data.ListVarName {'U_r','U_theta'}];
     311                    % Data.ListVarName=[Data.ListVarName {'U_r','U_theta'}];
    299312                    %Data.VarDimName=[Data.VarDimName {{theta_name,radius_name}} {{theta_name,radius_name}}];
    300313                    Data.VarAttribute{nbvar+1}.Role='vector_x';

  • trunk/src/uvmat.m

    r1093 r1094  
    15381538XmlData.LIFCalib.Ray1Coord=LineData{1}.Coord;
    15391539XmlData.LIFCalib.Ray2Coord=LineData{2}.Coord;
     1540if numel(LineData)<3
     1541    msgbox_uvmat('ERROR','draw a reference line of direct laser illumination (without dye absorbsion)');
     1542    return
     1543end
    15401544XmlData.LIFCalib.RefLineCoord=LineData{3}.Coord;
    15411545
    … …  
    15451549y=linspace(UvData.Field.Coord_y(1),UvData.Field.Coord_y(2),nby)-nby/2;
    15461550[X,Y]=meshgrid(x,y);
    1547 coeff_quad=0.15*4/(nbx*nbx);% image luminosity reduced by 10% at the edge
    1548 UvData.Field.A=double(UvData.Field.A).*(1+coeff_quad*(X.*X+Y.*Y));
     1551%coeff_quad=0.15*4/(nbx*nbx);% image luminosity reduced by 10% at the edge
     1552%UvData.Field.A=double(UvData.Field.A).*(1+coeff_quad*(X.*X+Y.*Y));
    15491553
    15501554%% display the current image in polar axes with origin at the  illumination source
    … …  
    15641568    y_ref=y_ref-y0;
    15651569    [theta_ref,r_ref] = cart2pol(x_ref,y_ref);%theta_ref  and r_ref are the polar coordinates of the points on the line
    1566     theta_ref=theta_ref*180/pi;% theta_ref in radians
     1570    theta_ref=theta_ref*180/pi;% theta_ref in degrees
    15671571    figure(10)
    15681572    plot(theta_ref,r_ref)
    … …  
    15701574    ylabel('radius from light source')
    15711575    title('ref line in polar coordinates')
    1572     azimuth_ima=linspace(DataPol.Coord_y(1),DataPol.Coord_y(2),size(DataPol.A,1));%array of angular index on the transformed image
     1576    azimuth_ima=linspace(DataPol.Coord_y(1),DataPol.Coord_y(2),size(DataPol.A,1))-360;%array of angular indices on the transformed image
    15731577    dist_source = interp1(theta_ref,r_ref,azimuth_ima);% get the polar position of the reference line
    15741578    dist_source_pixel=round(size(DataPol.A,2)*(dist_source-DataPol.Coord_x(1))/(DataPol.Coord_x(2)-DataPol.Coord_x(1)));
    … …  
    37193723%% choose and read a second field FileName_1 if defined
    37203724ParamOut_1=[];
    3721 if numel(UvData.FileInfo)>1
     3725if ~isempty(FileName_1)
     3726    if numel(UvData.FileInfo)==1
     3727        UvData.FileInfo{2}=UvData.FileInfo{1};
     3728    end
    37223729    VelType_1=[];%default
    37233730    FieldName_1=[];
    … …  
    37363743        NomType_1=get(handles.NomType,'String');
    37373744    end
    3738     if strcmp(UvData.FileInfo{2}.FieldType,'image')
     3745    if strcmp(UvData.FileInfo{2}.FileType,'image')
    37393746        FieldName_1='image';
    37403747        frame_index_1=1;%default
    … …  
    37613768    end
    37623769    switch UvData.FileInfo{2}.FileType
    3763         case {'civx','civdata','netcdf','pivdata_fluidimage'};
     3770        case {'civx','civdata','netcdf','pivdata_fluidimage'}
    37643771            list_fields=get(handles.FieldName_1,'String');% list menu fields
    37653772            FieldName_1= list_fields{get(handles.FieldName_1,'Value')}; % selected field
    … …  
    38493856%% update the display menu for the second velocity type (second menuline)
    38503857test_veltype_1=0;
    3851 if isempty(FileName_1)
    3852 elseif ~test_keepdata_1
    3853     if strcmp(UvData.FileInfo{2}.FieldType,'civdata')&& ~strcmp(FieldName_1,'get_field...')
     3858if ~isempty(FileName_1)
     3859
     3860    if strcmp(UvData.FileInfo{2}.FileType,'civdata')&& ~strcmp(FieldName_1,'get_field...')
    38543861        test_veltype_1=1;
    38553862        set(handles.VelType_1,'Visible','on')
Note: See TracChangeset for help on using the changeset viewer.