source: trunk/src/script_readlvm.m @ 1134

Last change on this file since 1134 was 1123, checked in by sommeria, 14 months ago

various improvements

File size: 8.6 KB
Line 
1%% get the input file
2project='/fsnet/project/coriolis/2016/16CREST';
3% if ~exist(project,'dir')
4%     project='U:\project\coriolis\2015\15MINI_MEDDY\PROBES';%windows
5% end
6fileinput=uigetfile_uvmat('pick an input file',project);
7[Path,Name,Ext]=fileparts(fileinput);
8
9%% read the input file
10FileType='';
11if strcmp(Ext,'.lvm')
12    disp(['reading ' fileinput '...'])
13    Data=read_lvm(fileinput)
14elseif strcmp(Ext,'.nc')
15    disp(['reading ' fileinput '...'])
16    Data=nc2struct(fileinput)
17else
18    disp('invalid input file extension')
19end
20
21%% save netcdf file as .lvm
22if strcmp(Ext,'.lvm')
23    OutputFile=fullfile(Path,[Name '.nc']);
24    errormsg=struct2nc(OutputFile,Data);% copy the data in a netcdf file
25    if isempty(errormsg)
26        disp([OutputFile ' written'])
27    else
28        disp(errormsg)
29    end
30    [success,msg] = fileattrib(OutputFile,'+w','g')% allow writing access for the group of users
31end
32
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'); 
37if exist(XmlFile,'file')
38    ProbeDoc=xml2struct(XmlFile)
39else
40    disp('no calibration file .xml detected')
41end
42
43% if isfield(Data,'Position'), C2,C4,C6
44%     Min=1; Data.Position=Data.Position+PositionMin;
45% end
46% a5=-2.134088,b5=1010.1611, Data.C2=Data.C2; Data.C5=a5*Data.C5+b5; Data.C6=Data.C2;
47% ylabelstring='density drho (kg/m3)';
48
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
67for ilist=1:numel(Data.ListVarName)
68    if isequal(Data.ListVarName{ilist}(1),'C');% if the var name begins by 'C'
69        clist=clist+1;
70        CName{clist}=Data.ListVarName{ilist};
71        if isfield(ProbeDoc,CName{clist})&& ~isempty(ProbeDoc.(CName{clist}))
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
95            Data.(CName{clist})=filter(ones(1,20)/20,1,Data.(CName{clist})); % filter the signal to 10 Hz
96            ylabelstring='density drho (g/cm3)';
97                 
98        end
99    end
100end
101
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
105plot_string='plot(';
106for clist=1:numel(CName)
107    Data.(CName{clist})=filter(ones(1,bandwidth2)/bandwidth2,1,Data.(CName{clist}));%low pass filter
108    plot_string=[plot_string 'Data.Time,Data.' CName{clist} ','];
109end
110plot_string(end)=')';
111eval(plot_string)
112legend(CName')
113htitle=title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', conductivity probes']);
114set(htitle,'Interpreter','none')% desable tex interpreter
115xlabel('Time(s)')
116ylabel(ylabelstring)
117        ylim([1 1.02])
118grid on
119
120if isfield(Data,'Position')
121    %% plot motor position
122    figure(2)
123    set(2,'name','position')
124    plot(Data.Time,Data.Position)
125    htitle=title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', probe position ']);
126    set(htitle,'Interpreter','none')% desable tex interpreter
127    xlabel('Time(s)')
128    ylabel('Z (cm)')
129    grid on
130    hold on
131    plot(Data.Time,Data.Speed)
132       
133    %% plot conductivity probe profiles (limited to downward motion, Data.Speed<0)
134    if ~isempty(ProbeDoc)
135        figure(3)
136        set(3,'name','profiles')
137        Zmotor=Data.Position(Data.Speed<0);
138        Z=Zmotor*ones(1,numel(CName));%motor position transformed in a matrix with a columnfor each probe
139        Zmax=max(Zmotor);
140        plot_string='plot(';
141        for clist=1:numel(CName)
142            if isfield(ProbeDoc,CName{clist})&& ~isempty(ProbeDoc.(CName{clist})) && size(ProbeDoc.(CName{clist}).Position,2)>=2 % if at least two positions are defined to indicate that the probe moves
143               Zprobe=Zmotor-Zmax+ProbeDoc.(CName{clist}).Position(2,3);%upper position of the probe
144               Zprobe(Zprobe<ProbeDoc.(CName{clist}).Position(1,3))=ProbeDoc.(CName{clist}).Position(1,3);
145                Z(:,clist)=Zprobe;% add to z the first z position of the chosen probe (given in the xml file)
146                plot_string=[plot_string 'Z(:,' num2str(clist) '),Data.' CName{clist} '(Data.Speed<0),'];
147            end
148        end
149        plot_string(end)=')';
150        eval(plot_string)
151        legend(CName')
152        htitle=title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', conductivity probes']);
153        set(htitle,'Interpreter','none')% desable tex interpreter
154        xlabel('Z(cm)')
155        ylabel(ylabelstring)
156        ylim([1 1.02])
157        grid on
158    end
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')
193
194%% plot velocity (ADV) signals
195% figure(5)
196% set(5,'name','velocity')
197% plot(Data.Time,Data.ADV_X,Data.Time,Data.ADV_Y,Data.Time,Data.ADV_Z1,Data.Time,Data.ADV_Z2)
198% legend({'ADV_X';'ADV_Y';'ADV_Z1';'ADV_Z2'})
199% title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', signal velocity'])
200% xlabel('Time(s)')
201% ylabel('signal (Volt)')
202% grid on
203
204%% plot velocity interface signals
205% figure(6)
206% set(6,'name','interface')
207% plot(Data.Time,Data.I1,Data.Time,Data.I2,Data.Time,Data.I3,Data.Time,Data.I4)
208% legend({'I1';'I2';'I3';'I4'})
209% title([Data.Experiment ', '  Data.FileName ', Time=' Data.DateTime ', signal interface (ultrasound)'])
210% xlabel('Time(s)')
211% ylabel('signal (Volt)')
212% grid on
213
214
215
Note: See TracBrowser for help on using the repository browser.