Index: trunk/src/script_readlvm.m
===================================================================
--- trunk/src/script_readlvm.m	(revision 950)
+++ trunk/src/script_readlvm.m	(revision 959)
@@ -1,4 +1,8 @@
 %% get the input file
-fileinput=uigetfile_uvmat('pick an input file','/fsnet/project/coriolis/2016/16MILESTONE/Data');
+project='/fsnet/project/coriolis/2015/15MINI_MEDDY/PROBES';
+if ~exist(project,'dir')
+    project='U:\project\coriolis\2015\15MINI_MEDDY\PROBES';%windows
+end
+fileinput=uigetfile_uvmat('pick an input file',project);
 [Path,Name,Ext]=fileparts(fileinput);
 
@@ -7,5 +11,5 @@
 if strcmp(Ext,'.lvm')
     disp(['reading ' fileinput '...'])
-Data=read_lvm(fileinput)
+    Data=read_lvm(fileinput)
 elseif strcmp(Ext,'.nc')
     disp(['reading ' fileinput '...'])
@@ -15,5 +19,5 @@
 end
 
-%% save as netcdf file if it has been opened as .lvm
+%% save netcdf file as .lvm
 if strcmp(Ext,'.lvm')
     OutputFile=fullfile(Path,[Name '.nc']);
@@ -27,23 +31,38 @@
 end
 
-%% use calibration data stored in an xml file
-ProbeDoc=[];
-XmlFile=fullfile(Path,[Name '.xml']);
+%% use calibration stored in a specified xml file; always use the same file
+ProbeDoc=[]; 
+%XmlFile=fullfile(Path,[Name '.xml']);
+XmlFile=fullfile(Path,'calib.xml');  
 if exist(XmlFile,'file')
-    ProbeDoc=xml2struct(XmlFile);
-end
+    ProbeDoc=xml2struct(XmlFile)
+else
+    disp('no calibration file .xml detected')
+end
+
 % if isfield(Data,'Position'), C2,C4,C6
-%     Min=1;
-%     Data.Position=Data.Position+PositionMin;
+%     Min=1; Data.Position=Data.Position+PositionMin;
 % end
-% a5=-2.134088,b5=1010.1611,
-% Data.C2=Data.C2;
-% Data.C5=a5*Data.C5+b5;
-% Data.C6=Data.C2;
+% a5=-2.134088,b5=1010.1611, Data.C2=Data.C2; Data.C5=a5*Data.C5+b5; Data.C6=Data.C2;
 % ylabelstring='density drho (kg/m3)';
 
-%% treqnsform conductivity probe signals: calibration + filter at 10 Hz
-ylabelstring='conductivity signal (volts)';
-clist=0;% counter of conductivity probes
+
+%% transform temperature probe signals 
+if isfield(ProbeDoc,'T5')&& ~isempty(ProbeDoc.T5)  % if temperature calibration exists; see calibT.m
+    
+    Data.T5=exp((Data.T5 - ProbeDoc.T5.a)./ProbeDoc.T5.b);% transform volt signal into temperature 
+    Data.T5=filter(ones(1,60)/60,1,Data.T5); % filter the signal to 4 Hz 
+    figure(6); set(6,'name','temperature'); plot(Data.Time,Data.T5);
+    %plot(Data.Time,Data.T5,Data.Time,20+Data.Position/100)
+    title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', temperature'])
+    xlabel('Time(s)'); ylabel('temperature (degree C)'); %ylim([20 37])
+    grid on
+end
+
+%% check camera signal
+ind_start=find(Data.Trig_Cam>3.5,1,'first')
+disp(['camera starts at time ' num2str(Data.Time(ind_start))])
+%% transform and filter conductivity probe signals into [temperature-corrected] density
+ylabelstring='conductivity signal (volts)'; clist=0;% counter of conductivity probes
 for ilist=1:numel(Data.ListVarName)
     if isequal(Data.ListVarName{ilist}(1),'C');% if the var name begins by 'C'
@@ -51,18 +70,40 @@
         CName{clist}=Data.ListVarName{ilist};
         if isfield(ProbeDoc,CName{clist})&& ~isempty(ProbeDoc.(CName{clist}))
-            a=ProbeDoc.(CName{clist}).a;
-            b=ProbeDoc.(CName{clist}).b;
-            Data.(CName{clist})=a*Data.(CName{clist})+b;% transform volt signal into density
+        
+            a=ProbeDoc.(CName{clist}).a; b=ProbeDoc.(CName{clist}).b;
+            % Data.(CName{clist})=a*Data.(CName{clist})+b;% volts STRAIGHT into density (if const T)
+            
+            %% BUT now need to modify conductivity, density due to temperature effect
+            %% first put into conductivity - using just C5 conductivity calibration for now (22/7/15)
+            ac1 = 0.5766e4; bc1 = 2.9129e4; %% this was found later, w/ different gain. Use the one below instead
+            %% ac1 = 0.4845e4; bc1 = 2.064e4; %% this was w/ the gain when we did the experiments; but it doesn't work?
+            Data.(CName{clist})=ac1*Data.(CName{clist})+bc1; %% voltage translated into conductivity via calibration
+             
+            %% read in temperature... 
+            refT = 23; T = Data.T5;
+            Hewitt_fit = [2.2794885e-11 -6.2634979e-9 1.5439826e-7 7.8601061e-5 2.1179818e-2]; 
+            bfit = reshape(polyval(Hewitt_fit,T(:)),size(T)); bref = reshape(polyval(Hewitt_fit,refT(:)),size(refT));
+            sigTsig18 = (1+bfit.*(T-18)); sigREFsig18 = (1+bref.*(refT-18)); 
+            modfactor = sigTsig18./sigREFsig18;
+            modfactor = 1./modfactor;  %% now in sigREF/sigT, which (* sigT) below to get sigREF, which can convert to rho
+            Data.(CName{clist})=Data.(CName{clist}).*modfactor %% = temp corrected conductivity (= as if at reference temp)
+            
+            %% now we need to put the modified conductivity back into a (modified) voltage, then estimate density directly. 
+            ac2 = 1.7343e-4; bc2 = -5.0048;
+            Data.(CName{clist})=ac2*Data.(CName{clist})+bc2; % temp-corrected voltage
+            Data.(CName{clist})=a*Data.(CName{clist})+b; % now finally voltage into density
             Data.(CName{clist})=filter(ones(1,20)/20,1,Data.(CName{clist})); % filter the signal to 10 Hz
             ylabelstring='density drho (g/cm3)';
+                  
         end
     end
 end
 
-%% plot conductivity probe signals
-figure(1)
-set(1,'name','conductivity')
+%% plot conductivity signals
+figure(1); set(1,'name','conductivity')
+bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz
 plot_string='plot(';
 for clist=1:numel(CName)
+    Data.(CName{clist})=filter(ones(1,bandwidth2)/bandwidth2,1,Data.(CName{clist}));%low pass filter
     plot_string=[plot_string 'Data.Time,Data.' CName{clist} ','];
 end
@@ -74,8 +115,9 @@
 xlabel('Time(s)')
 ylabel(ylabelstring)
+        ylim([1 1.02])
 grid on
 
 if isfield(Data,'Position')
-    %% plot  motor position
+    %% plot motor position
     figure(2)
     set(2,'name','position')
@@ -88,12 +130,5 @@
     hold on
     plot(Data.Time,Data.Speed)
-    
-    %%plot first profile
-    figure(4)
-    index1=find(Data.Speed<0,1); %start of descent
-    Speed=Data.Speed(index1:end);
-    index2=index1-1+find(Speed>0,1); %end of descent
-    plot(Data.Position(index1:index2),Data.C1(index1:index2))
-    
+       
     %% plot conductivity probe profiles (limited to downward motion, Data.Speed<0)
     if ~isempty(ProbeDoc)
@@ -119,41 +154,41 @@
         xlabel('Z(cm)')
         ylabel(ylabelstring)
+        ylim([1 1.02])
         grid on
     end
-    %     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))
-    %     legend({'C2';'C3';'C4';'C5'})
-    %     title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', profiles conductivity probes'])
-    %     xlabel('Z(cm)')
-    %     %ylabel('signal (Volt)')
-    %     ylabel(ylabelstring)
-    %     grid on
-    
-    
-    %     set(3,'name','profiles')
-    %     PositionMin=1; %position of probes at the bottom (Z=1 cm)
-    %     Data.Position=Data.Position-min(Data.Position)+PositionMin;
-    %     Z=Data.Position(Data.Speed<0);
-    %     plot(Z,Data.C2(Data.Speed<0),Z,Data.C5(Data.Speed<0),Z,Data.C6(Data.Speed<0))
-    %     legend({'C2';'C5';'C6'})
-    %     title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', profiles conductivity probes'])
-    %     xlabel('Z(cm)')
-    %     %ylabel('signal (Volt)')
-    %     ylabel(ylabelstring)
-    %     grid on
-    
-end
-
-
-
-
-%% plot temperature signals
-% figure(4)
-% set(4,'name','temperature')
-% plot(Data.Time,Data.T2,Data.Time,Data.T5,Data.Time,Data.T6)
-% legend({'T2';'T5';'T6'})
-% title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', signal conductivity probes'])
-% xlabel('Time(s)')
-% ylabel('signal (Volt)')
-% grid on
+        
+end
+
+%%%%
+figure(4)
+bandwidth1=480; % corresponds to 2 cm (1 cm/s with 240 pts/s)
+bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz
+C5_filter_low=filter(ones(1,bandwidth2)/bandwidth2,1,Data.C5);%low pass filter
+C5_filter=filter(ones(1,bandwidth1)/bandwidth1,1,C5_filter_low);%low pass filter
+C5_filter=C5_filter_low-C5_filter;% high pass filter
+C5_filter(Data.Speed>-0.1)=NaN;
+plot(Data.Time,C5_filter)
+ylim([-0.001 0.001])
+hold on
+plot(Data.Time,Data.Position/100000)
+grid on
+hold off
+title('C5 filtered')
+
+%%%%
+figure(5)
+bandwidth1=480; % corresponds to 2 cm (1 cm/s with 240 pts/s)
+bandwidth2=60; % corresponds to 0.25 cm (1 cm/s with 240 pts/s), removes 4 Hertz
+C3_filter_low=filter(ones(1,bandwidth2)/bandwidth2,1,Data.C3);%low pass filter
+C3_filter=filter(ones(1,bandwidth1)/bandwidth1,1,C3_filter_low);%low pass filter
+C3_filter=C3_filter_low-C3_filter;% high pass filter
+C3_filter(Data.Speed>-0.1)=NaN;
+plot(Data.Time,C3_filter)
+ylim([-0.001 0.001])
+hold on
+plot(Data.Time,Data.Position/100000)
+grid on
+hold off
+title('C3 filtered')
 
 %% plot velocity (ADV) signals
@@ -178,2 +213,3 @@
 
 
+
