Ignore:
Timestamp:
Jun 24, 2016, 8:49:08 PM (7 years ago)
Author:
sommeria
Message:

various

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/src/script_readlvm.m

    r950 r959  
    11%% get the input file
    2 fileinput=uigetfile_uvmat('pick an input file','/fsnet/project/coriolis/2016/16MILESTONE/Data');
     2project='/fsnet/project/coriolis/2015/15MINI_MEDDY/PROBES';
     3if ~exist(project,'dir')
     4    project='U:\project\coriolis\2015\15MINI_MEDDY\PROBES';%windows
     5end
     6fileinput=uigetfile_uvmat('pick an input file',project);
    37[Path,Name,Ext]=fileparts(fileinput);
    48
    … …  
    711if strcmp(Ext,'.lvm')
    812    disp(['reading ' fileinput '...'])
    9 Data=read_lvm(fileinput)
     13    Data=read_lvm(fileinput)
    1014elseif strcmp(Ext,'.nc')
    1115    disp(['reading ' fileinput '...'])
    … …  
    1519end
    1620
    17 %% save as netcdf file if it has been opened as .lvm
     21%% save netcdf file as .lvm
    1822if strcmp(Ext,'.lvm')
    1923    OutputFile=fullfile(Path,[Name '.nc']);
    … …  
    2731end
    2832
    29 %% use calibration data stored in an xml file
    30 ProbeDoc=[];
    31 XmlFile=fullfile(Path,[Name '.xml']);
     33%% use calibration stored in a specified xml file; always use the same file
     34ProbeDoc=[];
     35%XmlFile=fullfile(Path,[Name '.xml']);
     36XmlFile=fullfile(Path,'calib.xml'); 
    3237if exist(XmlFile,'file')
    33     ProbeDoc=xml2struct(XmlFile);
    34 end
     38    ProbeDoc=xml2struct(XmlFile)
     39else
     40    disp('no calibration file .xml detected')
     41end
     42
    3543% if isfield(Data,'Position'), C2,C4,C6
    36 %     Min=1;
    37 %     Data.Position=Data.Position+PositionMin;
     44%     Min=1; Data.Position=Data.Position+PositionMin;
    3845% end
    39 % a5=-2.134088,b5=1010.1611,
    40 % Data.C2=Data.C2;
    41 % Data.C5=a5*Data.C5+b5;
    42 % Data.C6=Data.C2;
     46% a5=-2.134088,b5=1010.1611, Data.C2=Data.C2; Data.C5=a5*Data.C5+b5; Data.C6=Data.C2;
    4347% ylabelstring='density drho (kg/m3)';
    4448
    45 %% treqnsform conductivity probe signals: calibration + filter at 10 Hz
    46 ylabelstring='conductivity signal (volts)';
    47 clist=0;% counter of conductivity probes
     49
     50%% transform temperature probe signals
     51if isfield(ProbeDoc,'T5')&& ~isempty(ProbeDoc.T5)  % if temperature calibration exists; see calibT.m
     52   
     53    Data.T5=exp((Data.T5 - ProbeDoc.T5.a)./ProbeDoc.T5.b);% transform volt signal into temperature
     54    Data.T5=filter(ones(1,60)/60,1,Data.T5); % filter the signal to 4 Hz
     55    figure(6); set(6,'name','temperature'); plot(Data.Time,Data.T5);
     56    %plot(Data.Time,Data.T5,Data.Time,20+Data.Position/100)
     57    title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', temperature'])
     58    xlabel('Time(s)'); ylabel('temperature (degree C)'); %ylim([20 37])
     59    grid on
     60end
     61
     62%% check camera signal
     63ind_start=find(Data.Trig_Cam>3.5,1,'first')
     64disp(['camera starts at time ' num2str(Data.Time(ind_start))])
     65%% transform and filter conductivity probe signals into [temperature-corrected] density
     66ylabelstring='conductivity signal (volts)'; clist=0;% counter of conductivity probes
    4867for ilist=1:numel(Data.ListVarName)
    4968    if isequal(Data.ListVarName{ilist}(1),'C');% if the var name begins by 'C'
    … …  
    5170        CName{clist}=Data.ListVarName{ilist};
    5271        if isfield(ProbeDoc,CName{clist})&& ~isempty(ProbeDoc.(CName{clist}))
    53             a=ProbeDoc.(CName{clist}).a;
    54             b=ProbeDoc.(CName{clist}).b;
    55             Data.(CName{clist})=a*Data.(CName{clist})+b;% transform volt signal into density
     72       
     73            a=ProbeDoc.(CName{clist}).a; b=ProbeDoc.(CName{clist}).b;
     74            % Data.(CName{clist})=a*Data.(CName{clist})+b;% volts STRAIGHT into density (if const T)
     75           
     76            %% BUT now need to modify conductivity, density due to temperature effect
     77            %% first put into conductivity - using just C5 conductivity calibration for now (22/7/15)
     78            ac1 = 0.5766e4; bc1 = 2.9129e4; %% this was found later, w/ different gain. Use the one below instead
     79            %% ac1 = 0.4845e4; bc1 = 2.064e4; %% this was w/ the gain when we did the experiments; but it doesn't work?
     80            Data.(CName{clist})=ac1*Data.(CName{clist})+bc1; %% voltage translated into conductivity via calibration
     81             
     82            %% read in temperature...
     83            refT = 23; T = Data.T5;
     84            Hewitt_fit = [2.2794885e-11 -6.2634979e-9 1.5439826e-7 7.8601061e-5 2.1179818e-2];
     85            bfit = reshape(polyval(Hewitt_fit,T(:)),size(T)); bref = reshape(polyval(Hewitt_fit,refT(:)),size(refT));
     86            sigTsig18 = (1+bfit.*(T-18)); sigREFsig18 = (1+bref.*(refT-18));
     87            modfactor = sigTsig18./sigREFsig18;
     88            modfactor = 1./modfactor;  %% now in sigREF/sigT, which (* sigT) below to get sigREF, which can convert to rho
     89            Data.(CName{clist})=Data.(CName{clist}).*modfactor %% = temp corrected conductivity (= as if at reference temp)
     90           
     91            %% now we need to put the modified conductivity back into a (modified) voltage, then estimate density directly.
     92            ac2 = 1.7343e-4; bc2 = -5.0048;
     93            Data.(CName{clist})=ac2*Data.(CName{clist})+bc2; % temp-corrected voltage
     94            Data.(CName{clist})=a*Data.(CName{clist})+b; % now finally voltage into density
    5695            Data.(CName{clist})=filter(ones(1,20)/20,1,Data.(CName{clist})); % filter the signal to 10 Hz
    5796            ylabelstring='density drho (g/cm3)';
     97                 
    5898        end
    5999    end
    60100end
    61101
    62 %% plot conductivity probe signals
    63 figure(1)
    64 set(1,'name','conductivity')
     102%% plot conductivity signals
     103figure(1); set(1,'name','conductivity')
     104bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz
    65105plot_string='plot(';
    66106for clist=1:numel(CName)
     107    Data.(CName{clist})=filter(ones(1,bandwidth2)/bandwidth2,1,Data.(CName{clist}));%low pass filter
    67108    plot_string=[plot_string 'Data.Time,Data.' CName{clist} ','];
    68109end
    … …  
    74115xlabel('Time(s)')
    75116ylabel(ylabelstring)
     117        ylim([1 1.02])
    76118grid on
    77119
    78120if isfield(Data,'Position')
    79     %% plot  motor position
     121    %% plot motor position
    80122    figure(2)
    81123    set(2,'name','position')
    … …  
    88130    hold on
    89131    plot(Data.Time,Data.Speed)
    90    
    91     %%plot first profile
    92     figure(4)
    93     index1=find(Data.Speed<0,1); %start of descent
    94     Speed=Data.Speed(index1:end);
    95     index2=index1-1+find(Speed>0,1); %end of descent
    96     plot(Data.Position(index1:index2),Data.C1(index1:index2))
    97    
     132       
    98133    %% plot conductivity probe profiles (limited to downward motion, Data.Speed<0)
    99134    if ~isempty(ProbeDoc)
    … …  
    119154        xlabel('Z(cm)')
    120155        ylabel(ylabelstring)
     156        ylim([1 1.02])
    121157        grid on
    122158    end
    123     %     plot(Z,Data.C2(Data.Speed<0),Z,Data.C3(Data.Speed<0),Z,Data.C4(Data.Speed<0),Z,Data.C5(Data.Speed<0))
    124     %     legend({'C2';'C3';'C4';'C5'})
    125     %     title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', profiles conductivity probes'])
    126     %     xlabel('Z(cm)')
    127     %     %ylabel('signal (Volt)')
    128     %     ylabel(ylabelstring)
    129     %     grid on
    130    
    131    
    132     %     set(3,'name','profiles')
    133     %     PositionMin=1; %position of probes at the bottom (Z=1 cm)
    134     %     Data.Position=Data.Position-min(Data.Position)+PositionMin;
    135     %     Z=Data.Position(Data.Speed<0);
    136     %     plot(Z,Data.C2(Data.Speed<0),Z,Data.C5(Data.Speed<0),Z,Data.C6(Data.Speed<0))
    137     %     legend({'C2';'C5';'C6'})
    138     %     title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', profiles conductivity probes'])
    139     %     xlabel('Z(cm)')
    140     %     %ylabel('signal (Volt)')
    141     %     ylabel(ylabelstring)
    142     %     grid on
    143    
    144 end
    145 
    146 
    147 
    148 
    149 %% plot temperature signals
    150 % figure(4)
    151 % set(4,'name','temperature')
    152 % plot(Data.Time,Data.T2,Data.Time,Data.T5,Data.Time,Data.T6)
    153 % legend({'T2';'T5';'T6'})
    154 % title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', signal conductivity probes'])
    155 % xlabel('Time(s)')
    156 % ylabel('signal (Volt)')
    157 % grid on
     159       
     160end
     161
     162%%%%
     163figure(4)
     164bandwidth1=480; % corresponds to 2 cm (1 cm/s with 240 pts/s)
     165bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz
     166C5_filter_low=filter(ones(1,bandwidth2)/bandwidth2,1,Data.C5);%low pass filter
     167C5_filter=filter(ones(1,bandwidth1)/bandwidth1,1,C5_filter_low);%low pass filter
     168C5_filter=C5_filter_low-C5_filter;% high pass filter
     169C5_filter(Data.Speed>-0.1)=NaN;
     170plot(Data.Time,C5_filter)
     171ylim([-0.001 0.001])
     172hold on
     173plot(Data.Time,Data.Position/100000)
     174grid on
     175hold off
     176title('C5 filtered')
     177
     178%%%%
     179figure(5)
     180bandwidth1=480; % corresponds to 2 cm (1 cm/s with 240 pts/s)
     181bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz
     182C3_filter_low=filter(ones(1,bandwidth2)/bandwidth2,1,Data.C3);%low pass filter
     183C3_filter=filter(ones(1,bandwidth1)/bandwidth1,1,C3_filter_low);%low pass filter
     184C3_filter=C3_filter_low-C3_filter;% high pass filter
     185C3_filter(Data.Speed>-0.1)=NaN;
     186plot(Data.Time,C3_filter)
     187ylim([-0.001 0.001])
     188hold on
     189plot(Data.Time,Data.Position/100000)
     190grid on
     191hold off
     192title('C3 filtered')
    158193
    159194%% plot velocity (ADV) signals
    … …  
    178213
    179214
     215
Note: See TracChangeset for help on using the changeset viewer.