Commit d7f7a6dc authored by Jürgen Matzka's avatar Jürgen Matzka
Browse files

Version von obs251

parent b2fd85aa
......@@ -12,10 +12,10 @@
clear all
close all
station = ['TTB1l'];
station = ['WNG0g'];
%station = ['BFO1s';'VNA1s']
starttime = '2018-03-28';
stoptime = '2018-04-02'; % nicht für Varlab
starttime = '2018-03-10';
stoptime = '2018-04-05'; % nicht für Varlab
% -------------------------------------------------------------------------
......
......@@ -10,14 +10,16 @@
clear all
close all
station = ['VNA1s']
%station = ['SMA0g']
%station = ['SHE0g';'TTB1l']
%station = ['VNA1s']
%station = ['NGK0g']
%station = ['TDC1d']
%station = ['SHE0g']
%station = ['WNG0g']
%station = ['NGK0g';'NGK1g';'NGK2g';'NGK6g';'WNG0g';'WNG1g';'WNG2g';'WNG3g';'VNA1s';'SHE0g';'VSS0g';'TTB0g';'TTB1g';'VSS1g';'PUL0g'; 'TEL0g'];
starttime = '2018-03-29';
stoptime = '2018-04-02'; % nicht für Varlab
%station = ['TTB1l']
%station = ['BFO1s']
%station = ['NGK0g';'NGK1g';'NGK2g';'NGK6g';'WNG0g';'WNG1g';'WNG2g';'WNG3g';'VNA1s';'SHE0g';'VSS0g';'TTB0g';'TTB1l';'VSS1g';'SMA0g';'PUL0g';'TEL0g'];
starttime = '2017-05-08';
stoptime = '2017-05-2'; % nicht für Varlab
% -------------------------------------------------------------------------
......@@ -30,19 +32,7 @@ ppmsec2cdf('starttime',starttime,'stoptime',stoptime,'station',station);
% -------------------------------------------------------------------------
% cdf2def
% -------------------------------------------------------------------------
%
produce_IAGA_2002_min = 'P'; % kann auch 'n' sein, warum?
produce_IAGA_2002_sec = '0'; % kann auch 'n' sein, warum?
output_resolution = '0';
calculate_IAGA_F = 'n';
output_format = 'c';
cdf2def('starttime',starttime,'stoptime',stoptime,...
'station',station,...
'produce_IAGA_2002_min',produce_IAGA_2002_min,...
'produce_IAGA_2002_sec',produce_IAGA_2002_min,...
'calculate_IAGA_F', calculate_IAGA_F,...
'output_format', output_format,...
'output_resolution', output_resolution);
% - produce_IAGA_2002_min/sec '0' (default) no output
% 'P' provisional data
% 'D' definitive data
......@@ -62,9 +52,6 @@ cdf2def('starttime',starttime,'stoptime',stoptime,...
% 'r' means raw, only raw-file 1 second
% and IAGA V
%
% -------------------------------------------------------------------------
% varlab
......@@ -75,6 +62,7 @@ varlab_01('parfile','parameters_win_JM.txt','date',starttime,...
% % -------------------------------------------------------------------------
% % Start programs
% % -------------------------------------------------------------------------
......
......@@ -10,14 +10,16 @@
clear all
close all
station = ['SMA0g']
%station = ['SHE0g';'TTB1l']
%station = ['VNA1s']
%station = ['NGK0g']
%station = ['TDC1d']
station = ['SHE0g']
%station = ['WNG0g']
%station = ['NGK0g';'NGK1g';'NGK2g';'NGK6g';'WNG0g';'WNG1g';'WNG2g';'WNG3g';'VNA1s';'SHE0g';'VSS0g';'TTB0g';'TTB1g';'VSS1g';'PUL0g'; 'TEL0g'];
starttime = '2018-04-02';
stoptime = '2018-04-02'; % nicht für Varlab
%station = ['TTB1l']
%station = ['BFO1s']
%station = ['NGK0g';'NGK1g';'NGK2g';'NGK6g';'WNG0g';'WNG1g';'WNG2g';'WNG3g';'VNA1s';'SHE0g';'VSS0g';'TTB0g';'TTB1l';'VSS1g';'SMA0g';'PUL0g';'TEL0g'];
starttime = '2017-05-08';
stoptime = '2017-05-11'; % nicht für Varlab
% -------------------------------------------------------------------------
......@@ -30,19 +32,7 @@ ppmsec2cdf('starttime',starttime,'stoptime',stoptime,'station',station);
% -------------------------------------------------------------------------
% cdf2def
% -------------------------------------------------------------------------
%
produce_IAGA_2002_min = 'P'; % kann auch 'n' sein, warum?
produce_IAGA_2002_sec = '0'; % kann auch 'n' sein, warum?
output_resolution = '0';
calculate_IAGA_F = 'n';
output_format = 'c';
cdf2def('starttime',starttime,'stoptime',stoptime,...
'station',station,...
'produce_IAGA_2002_min',produce_IAGA_2002_min,...
'produce_IAGA_2002_sec',produce_IAGA_2002_min,...
'calculate_IAGA_F', calculate_IAGA_F,...
'output_format', output_format,...
'output_resolution', output_resolution);
% - produce_IAGA_2002_min/sec '0' (default) no output
% 'P' provisional data
% 'D' definitive data
......@@ -62,9 +52,6 @@ cdf2def('starttime',starttime,'stoptime',stoptime,...
% 'r' means raw, only raw-file 1 second
% and IAGA V
%
% -------------------------------------------------------------------------
% varlab
......@@ -75,6 +62,7 @@ varlab_01('parfile','parameters_win_JM.txt','date',starttime,...
% % -------------------------------------------------------------------------
% % Start programs
% % -------------------------------------------------------------------------
......
......@@ -21,8 +21,8 @@
% attachment: same file
clear all
startdate = jd2000(2018, 03, 29); % HERE START DATE FOR non-online
enddate = jd2000(2018, 04, 02); % HERE END DATE FOR non-online
startdate = jd2000(2018, 04, 26); % HERE START DATE FOR non-online
enddate = jd2000(2018, 04, 26); % HERE END DATE FOR non-online
station = 'VNA1';
%station = 'WNG1';
......
This diff is collapsed.
% datalogger code 'a' maybe for tabulated values
% datalogger code 'b' maybe for photographic paper
% datalogger code 'c' for PCD
% datalogger code 'd' for ADAM module
% datalogger code 'g' for GFZ 0.5 Hz
% datalogger code 's' for seismic
% datalogger code 'l' for LEMI
% datalogger code 'f' for Flare
%
% 2013-04-15 jrgm:
% make function calcualte_baseline.m fit for Matlab2013 cdfread:
% %informcdf.Variables{i,3},1 contains the length of the array
%
% make function calcualte_ime_error_absolutes.m fit for Matlab2013 cdfread:
% %informcdf.Variables{i,3},1 contains the length of the array
clear all
close all
figure
%station = 'VSS0g'; % c or d or g or s or l for datalogger
%station = 'WNG1g'; % c or d or g or s or l for datalogger
%station = 'NGK0g'; % c or d or g or s or l for datalogger
% station = 'TDC1d'; % c or d or g or s or l for datalogger
station = 'VNA1s'; % c or d or g or s or l for datalogger
% station = 'BFO1s'; % c or d or g or s or l for datalogger
%station = 'SHE0g'; % c or d or g or s or l for datalogger
%station = 'TTB0g'; % c or d or g or s or l for datalogger
%station = 'TTB1l'; % c or d or g or s or l for datalogger
%station = 'VSS0g'; % c or d or g or s or l for datalogger
tempsource = 'c'; % c for cdf
temp_correct = 'y'; % apply temp-correction function to variations 'y'
apply_time_lag = 'y'; %reads from raw1(K61)
find_time_error = 'n';
timeshift_array = [-300:1:300]';
number_DI = 1000; % read in all DI-measurements for this year
number_DI = 5; % read in the last x measurements in this year
number_DI = 0; % use date range from below
year = 2018;
start_mmdd = '0325';
end___mmdd = '0325';
ppmsource = 'c'; % 's' for spreadsheet, 'c' for cdf, here always 'c'
% If 's' is to be used, then change this in calculate_baselines.m
return_array = []; %returnes baseliens etc from function
opsys = computer;
% if strcmp(opsys(1:2), 'GL') == 1 %removed because XLS-files can only be
% processed on Windows systems
% path1 = '/nfs/magobs/';
% bs = '/'; % slash
% else
if strcmp(opsys(1:2), 'PC') == 1
path1 = 'Y:';
bs = '\'; % backslash
else
diplay('operating system not recognised');
end
path_DI = [path1 bs 'BPDATA' bs 'DI_' station(1:3) ...
num2str(year, '%4.4d') bs];
if strcmp(station(1:3), 'MSN')
path_DI = [path1 bs 'BPDATA' bs 'DI_' 'NGK' ...
num2str(year, '%4.4d') bs];
end %MSN
path_cdf = [path1 bs 'laptop_data' bs station(1:4) bs 'cdf' bs];
if strcmp(pwd, 'O:\jmat\0programme_win')
path_DI = ['O:' bs 'jmat' bs '0DI_measurements' bs 'DI_' ...
station(1:3) ...
num2str(year, '%4.4d') bs];
path_cdf = ['O:' bs 'jmat' bs station(1:4) bs 'cdf' bs];
end
%month = 1; removed 27/2-2012
%day = 7;
% % START Identify and convert Achims GFZ-DI-schemes
% % with name e.g. ng0_140102.xls
% %
% %
% if strcmp(station(1:3), 'NGK') || strcmp(station(1:3), 'WNG') ...
% || strcmp(station(1:3), 'VRE') || strcmp(station(1:3), 'SHE') ...
% || strcmp(station(1:3), 'PAG') || strcmp(station(1:3), 'SUR') ...
% || strcmp(station(1:3), 'ODE') || strcmp(station(1:3), 'KMH') ...
% || strcmp(station(1:3), 'YAK') || strcmp(station(1:3), 'MGD') ...
% || strcmp(station(1:3), 'ABG') || strcmp(station(1:3), 'HYD') ...
% || strcmp(station(1:3), 'PNL') || strcmp(station(1:3), 'SMA') ...
% || strcmp(station(1:3), 'LOC') || strcmp(station(1:3), 'PET')
% filename_criterion = [path_DI lower(station(1:2)) station(4) ...
% '_' num2str(year-2000, '%2.2d') '*.xls'];
% filename_length = 14;
%
%
% % find all files in directory that satisfy criterion
% B = dir(filename_criterion);
% C = struct2cell(B);
%
% %remove filernames with wrong length, this is because I add comment to file
% %name as indicator for bad measuremmtns and this changes length of filename
% for i = 1:size(C, 2)
% length_here = size(char(C{1,i}), 2);
% if length_here == filename_length
% D(i) = 1;
% else
% D(i) = 0;
% end
% end
% C = C(:,D == 1);
% filenames = char(C{1,:});
% clear D
%
% % I skip this part here: remove filenames that have '0' or '1' at 'end-11'
%
% % reduce number of DI-measurements to selected date range or number
% if number_DI ~= 1000
% if number_DI == 0
% for i = 1:size(C, 2)
% characters_here = char(C{1,i});
% if str2num(characters_here(end-7:end-4)) >= ...
% str2num(start_mmdd) &&...
% str2num(characters_here(end-7:end-4)) <= ...
% str2num(end___mmdd)
% D(i) = 1;
% else
% D(i) = 0;
% end
%
% end
% C = C(:,D == 1);
% filenames = char(C{1,:});
% clear D
% end
% else
% filenames1 = filenames;
% filenames = filenames1(end - number_DI+1:end,:);
% end
%
% % call function to convert
% for i = 1:size(filenames, 1)
% DI_file = [filenames(i,:)];
% if exist([path_DI DI_file], 'file')
% % convert
% [return_string] = ...
% calculate_baseline_convert_GFZ_form(station, path_DI, DI_file);
% % move GFZ-files to directory 'converted_GFZ'
% if exist([path_DI 'converted_GFZ'],'dir') == 0
% mkdir([path_DI 'converted_GFZ'])
% end
% for j = 1:size(filenames, 1)
% % incode here moving of filenames to 'conver...
% % folder!!!!!!!!!!!!!
% end
% end
% end %i
% clear B C D filenames filenames1 filename
%
% end
% %
% %
% %
% % END Identify and convert Achims GFZ-DI-schemes
%
%
%read in filenames for DI-files from directory, remove unrelevant filenames
%from list
if strcmp(station(1:3), 'TDC') || strcmp(station(1:3), 'VNA') ...
|| strcmp(station(1:3), 'BFO') || strcmp(station(1:3), 'MSN')...
|| strcmp(station(1:3), 'NGK') || strcmp(station(1:3), 'WNG') ...
|| strcmp(station(1:3), 'TTB') || strcmp(station(1:3), 'VRE')...
|| strcmp(station(1:3), 'SHE') || strcmp(station(1:3), 'VSS');
filename_criterion = [path_DI 'DI_' station(1:3) '_'...
num2str(year, '%4.4d') '*.xls'];
filename_length = 23;
% elseif strcmp(station(1:3), 'WNG')
%
% filename_criterion = [path_DI 'DI_' station(1:3) '0_'...
% num2str(year-2000, '%2.2d') '*.xls'];
% filename_length = 22;
else
filename_criterion = [path_DI '*.xls'];
filename_length = 12;
end
B = dir(filename_criterion);
C = struct2cell(B);
for i = 1:size(C, 2) %remove filernames with wrong length
length_here = size(char(C{1,i}), 2);
if length_here == filename_length
D(i) = 1;
else
D(i) = 0;
end
end
C = C(:,D == 1);
filenames = char(C{1,:});
clear D
for i = 1:size(C, 2) %remove filenames that have '0' or '1' at 'end-11'
characters_here = char(C{1,i});
if strcmp(characters_here(end-11), '0') || ...
strcmp(characters_here(end-11), '1')
D(i) = 1;
else
D(i) = 0;
end
end
C = C(:,D == 1);
filenames = char(C{1,:});
clear D
% reduce number of DI-measurements
if number_DI ~= 1000
if number_DI == 0
for i = 1:size(C, 2) %remove filenames that have '0' or '1' at 'end-11'
characters_here = char(C{1,i});
if str2num(characters_here(end-11:end-8)) >= ...
str2num(start_mmdd) &&...
str2num(characters_here(end-11:end-8)) <= ...
str2num(end___mmdd)
D(i) = 1;
else
D(i) = 0;
end
end
C = C(:,D == 1);
filenames = char(C{1,:});
clear D
end
else
filenames1 = filenames;
filenames = filenames1(end - number_DI+1:end,:);
end
%
%
% call function for list of filenames
%
%
for i = 1:size(filenames, 1)
DI_file = [path_DI filenames(i,:)];
cdf_file = [path_cdf station(1:4) '_' num2str(year, '%4.4d')...
filenames(i,end-11:end-8) '.cdf'];
if exist(DI_file, 'file')
if strcmp(find_time_error, 'y')
[return_string] = ...
calculate_time_error_absolutes(DI_file, cdf_file, station,...
ppmsource, tempsource, temp_correct, timeshift_array);
else
if strcmp(station(1:3), 'VRE')
[return_string] = ...
calculate_baselines_VRE(DI_file, cdf_file, station,...
ppmsource, tempsource, temp_correct, apply_time_lag);
return_string = char(return_string{1});
else
[return_string] = ...
calculate_baselines(DI_file, cdf_file, station,...
ppmsource, tempsource, temp_correct, apply_time_lag);
return_string = char(return_string{1});
end
end
end
return_array = [return_array; return_string]
end %i
save('recent-DI', 'return_array')
%% compare IAF-files with absolute measurements
clear all
station = 'NAQ6d';
station = 'NGK0g';
%calculate IAF files for timeinterval:
timeinterval = [2014; 1; 17]; %start year, start month, number of months
timeinterval = [2018; 1; 3]; %start year, start month, number of months
% option to apply time correction
time_correction_str = 'n';
......
function [p]=TTB_rotate(date,PATH)
% Read data
TTB1=read_data(PATH,'TTB1',date);
TTB0=read_data(PATH,'TTB0',date);
% Calculate F from rotated values
%F=@(V,p) sqrt((V.HNvar.*Sc(V.T1,p(3)).*cosd(p(1))).^2+(V.HEvar.*Sc(V.T1,p(4))-V.HNvar.*Sc(V.T1,p(3)).*sind(p(1))).^2+(V.Zvar.*Sc(V.T1,p(5))).^2)+p(2);
%F=@(V,p) sqrt((V.HNvar.*cosd(p(1))).^2+(V.HEvar-V.HNvar.*sind(p(1))).^2+(V.Zvar.^2))+p(2);
%F=@(V,p) sqrt(((V.HNvar-V.HEvar./sind(p(1))).*cosd(p(1))).^2+(V.HEvar).^2+(V.Zvar.^2));
F=@(V,p) sqrt( ((V.HNvar-V.HEvar.*sind(p(1)))./cosd(p(1))).^2 + (V.HEvar).^2 + (V.Zvar.^2) ) +p(2);
% Function for misfit value
O=@(TTB1,p,n) nansum(abs(F(TTB1,p)-TTB1.Fsc).^n)^(1/n);
% Fit
n=2; % Least-squares
p=fminsearch(@(p) O(TTB1,p,n),[0 0]);
O(TTB1,p,n)
% Plot
clf
subplot(4,1,1)
plot(TTB1.time,TTB1.Fsc-TTB1.Fcalc,'.')
title('Original \Delta F')
subplot(4,1,2)
plot(TTB1.time,TTB1.Fsc-F(TTB1,p),'.')
title('Rotated \Delta F')
subplot(4,1,3)
plot(TTB1.time,-TTB1.T1)
% Now temperature
% Temperature scaling
Sc=@(V,T,p) V.*(p(1).*T.^2+ones);
%Sc=@(V,T,p) V+p(1).*T;
%Sc=@(V,T,p) V;
% Calculate F from rotated values
r=p(1:2);
F=@(V,p) sqrt( ((Sc(V.HNvar,V.T1,p(1))-Sc(V.HEvar,V.T1,p(1)).*sind(p(1)))./cosd(p(1))).^2 + (Sc(V.HEvar,V.T1,p(1))).^2 + (Sc(V.Zvar,V.T1,p(1)).^2) );
%F=@(V,p) sqrt((Sc(V.HNvar,V.T1,p(1)).*cosd(r(1))).^2+(Sc(V.HEvar,V.T1,p(1))-Sc(V.HNvar,V.T1,p(2)).*sind(r(1))).^2+(Sc(V.Zvar,V.T1,p(3))).^2)+r(2);
%F=@(V,p) sqrt((Sc(V.HNvar,V.T1,p(1)).*cosd(r(1))).^2+(Sc(V.HEvar,V.T1,p(1))-Sc(V.HNvar,V.T1,p(1)).*sind(r(1))).^2+(Sc(V.Zvar,V.T1,p(1))).^2)+r(2);
% Function for misfit value
O=@(TTB1,p,n) nansum(abs(F(TTB1,p)-TTB1.Fsc).^n)^(1/n);
% Fit
n=2; % Least-squares
p=fminsearch(@(p) O(TTB1,p,n),[0 0 0]);
O(TTB1,p,n)
% Plot
subplot(4,1,4)
plot(TTB1.time,TTB1.Fsc-F(TTB1,p),'.')
title('Rotated \Delta F')
figure
plot(TTB1.time,TTB0.HNvar-(Sc(TTB1.HNvar,TTB1.T1,p(1)).*cosd(r(1))));
p=[r p];
end
\ No newline at end of file
......@@ -55,7 +55,7 @@ function [ return_string ] = ...
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% TO DO: Zeitfehler des Variometers bercksichtigen
% TO DO: Zeitfehler des Variometers ber?cksichtigen
%
%
......@@ -596,7 +596,11 @@ if exist(DI_file, 'file') == 2;
if strcmp(station(1:3), 'TTB')
ppmsource = 's'; % 's' for spreadsheet, 'c' for cdf
end
if strcmp(station(1:3), 'SMA')...
&& timeDI(1) > jd2000(2017, 5, 11)...
&& timeDI(1) < jd2000(2017, 5, 12)
ppmsource = 's'; % 's' for spreadsheet, 'c' for cdf
end
......@@ -697,6 +701,28 @@ if exist(DI_file, 'file') == 2;
end
end
if strcmp(station(1:3), 'SMA')
mire_nr = raw1{8,6};
pillar_nr = 1;
deltaI = 0.0000 + 0.0000;
deltaD = 0.0000 + 0.0000;
if timeDI(1) >= jd2000(2017, 05, 01)
deltappm =11;% Sockeldiff. Apr. 2018
end
if mire_nr == 1
mire_azim = 328.64259 % in gon % 295.77833 deg
elseif mire_nr == 2
mire_azim = 315.4611 % in gon % 283.9150 deg
elseif mire_nr == 3
mire_azim = 328.2722 % in gon % 295.4450 deg
elseif mire_nr == 4
mire_azim = 328.9944 % in gon % 296.0950 deg
else
error('This mire doesnt exist')
end
end
if strcmp(station(1:3), 'NAQ') && timeDI(1) < jd2000(2010, 08, 21, 12)
mire_nr = raw1{8,6};
......@@ -885,7 +911,7 @@ if exist(DI_file, 'file') == 2;
% mire_azim = 292.5622; % GPS Azimuth 263.306/0.9 = 292.5622
% mire_azim = 292.5511; % SunAZ1 263.296/0.9 = 292.5511
mire_azim = 291.1300; % SunAZ2 262.017/0.9 = 291.1300
% HIER MUSS EVENTUELL NOCH GPS-Azimut 263.17? fr
% HIER MUSS EVENTUELL NOCH GPS-Azimut 263.17? f?r
% 2013 extra gemacht werden
elseif pillar_nr == 3 && mire_nr == 1 % position 2015 and mire M-L9
% mire_azim = 284.1011; % GPS Azimuty 255.691 deg = 255.691/0.9 = 284.1011
......@@ -915,8 +941,8 @@ if exist(DI_file, 'file') == 2;
deltaD = 0.0000 + 0.0000;
deltappm = 0.0 + 257.02; % instrument_diff = 0
mire_azim = 24.0883; %preliminary
% Azimut aus altem Spreadsheet: 201 40' 49''