%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% erik_read_cluster.m
%
% Matlab script that reads Cluster EFW/STAFF/FGM/Ephemeris/CIS/PEACE etc
% data using the isGetDataLite interface and do various preparations of
% that data, e.g. de-spin.
%
% Created by Jan-Erik Wahlund, IRF-Uppsala, 2001 (modified 2002).
% Modified by Michiko Morooka, IRF-Uppsala, 2002.
% Modified by Erik Bergman, IRF-Uppsala, HT 2002 - just so it fits me like a jumpsuit size small
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Open dbh server and assign a database handle (db) identifier.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%DATABASE = input('Give database as string (e.g. panurge:16) :');
DATABASE = 'disco:10';
db = Mat_DbOpen(DATABASE);

% Ask necessary information to extract data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t_start = input('Start time? Format [YYYY MM DD hh mm ss]: ');

% Since fft, inv-fft etc is done, the edges are subject to
% ringing effects which should be avoided. About 100 seconds of
% data before (and after) the selected time interval are affected.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t_start_staff = t_start + [0 0 0 0 -3 0];
t_start_fgm = t_start + [0 0 0 0 -3 0];

Dt = input('Duration in seconds?: ');
SC_status = input('Which S/C? Format [sc1 sc2 sc3 sc4] (yes=1,no=0): ');

reply = input('Normal mode data (25 Hz)? (yes=1,no=0): ');

if reply
f_E_str = '10Hz';
f_B_str = '0-10Hz';
else
reply = input('Full burst mode data (8 kHz)? (yes=1,no=0): ');
if reply
f_E_str = '8kHz';
f_B_str = '0-180Hz';
else
f_E_str = '180Hz';
f_B_str = '0-180Hz';
end
end

%INST_status = ...
% input('Read data? Format [EFW STAFF FGM CIS PEACE EDI] (yes=1,no=0): ');

INST_status = [1 0 1 1 0 1];

% Now data extraction starts.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if INST_status(1) == 1
% EFW electric field data.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for sc_no = 1:4
if SC_status(sc_no) == 1

notes = sprintf('Reading EFW on S/C-%1.0f\0',sc_no); disp(notes);
sc_name = sprintf('%1.0f\0',sc_no);

if ((sc_no == 1) & (toepoch(t_start) >= toepoch([2002 02 01 00 00 00])) )

disp('read p34');
[t_E1234, E1234] = isGetDataLite( db, t_start, Dt, ...
'Cluster', sc_name, 'efw', 'E', 'p34', f_E_str, 'hx');
dummy = zeros(size(E1234'));
tmp1 = cat(1, dummy, dummy);
E1234 = cat(1, tmp1, E1234');
else
disp('read p1234');
[t_E1234, E1234] = isGetDataLite( db, t_start, Dt, ...
'Cluster', sc_name, 'efw', 'E', 'p1234', f_E_str, 'hx');
end

E1234 = real( double(E1234) );
E_efw_vect = [ t_E1234'; detrend(E1234(3,:),0); detrend(E1234(2,:),0) ]';

[t_ph, ph] = isGetDataLite( db, t_start, Dt, ...
'Cluster', sc_name, 'ephemeris', 'phase', ' ', ' ', ' ');
ph = real( double(ph) );
ph_vect = [ t_ph ph ];
E_vect = EFW_despin( E_efw_vect, ph_vect );

[t_E12, E_12] = isGetDataLite( db, t_start, Dt, ...
'Cluster', sc_name, 'efw', 'E', 'p12', f_E_str, 'hx');
[t_E34, E_34] = isGetDataLite( db, t_start, Dt, ...
'Cluster', sc_name, 'efw', 'E', 'p34', f_E_str, 'hx');
E_12 = [ t_E12 E_12 ];
E_34 = [ t_E34 E_34 ];

if ((sc_no == 1) & (toepoch(t_start) >= toepoch([2002 02 01 00 00 00])) )
Ne = SC_pot1( db, t_start, Dt, sc_no );
else
Ne = erik_sc_potential( db, t_start, Dt, sc_no );
SCp = Ne(:, [1 3]);
Ne = Ne(:, [1 2]);
end

command = sprintf('E_vect_%1.0f = E_vect;\0',sc_no); eval(command);
command = sprintf('Ne_%1.0f = Ne;\0', sc_no); eval(command);
command = sprintf('SCp_%1.0f = SCp;\0', sc_no); eval(command);
command = sprintf('E_12_%1.0f = E_12;\0', sc_no); eval(command);
command = sprintf('E_34_%1.0f = E_34;\0', sc_no); eval(command);

clear E1234 t_E1234 t_ph ph E_efw_vect ph_vect Ne SCp E_vect t_E12 t_E34 E_12 E_34
end
end

disp('E-field data read & despun & offset-corrected!');

clear command notes sc_name sc_no

end % EFW_status end

if INST_status(2) == 1
% STAFF SC data.
%%%%%%%%%%%%%%%%
% Have you heard this one:
% Calibrated SC data must be read by a potens of 2 samples!!!
% Sigh...

disp('No of E-points');
% Some n_samples is needed. Don't care from which S/C.
if SC_status(1) == 1, n_samples = length( E_vect_1(:,1) )
elseif SC_status(2) == 1, n_samples = length( E_vect_2(:,1) )
elseif SC_status(3) == 1, n_samples = length( E_vect_3(:,1) )
elseif SC_status(4) == 1, n_samples = length( E_vect_4(:,1) )
end
pot_read = ceil( log10( n_samples ) / log10(2) );
disp('No of B-points');
n_read = 2^pot_read

Dt_staff = Dt * n_read/n_samples + 360;

%
% Big STAFF reading problem! As soon as there is a data gap reading
% stops due to the fft-inverse-fft-procedure for STAFF SC calibrated
% data. Therefore we have too check end of data and loop the reading
% procedure assuming there are short data gaps with less than 1 min
% duration. This means the inverse-fft-filter effects are visable near
% those data gaps. Sigh... you can't have it all.
%

% NOTE!!!
% x is spin axis for WEC ! STAFF has y as z and vice versa !

% Despin STAFF data. Uses WEC coordinates!
% Same phase data as EFW.

for sc_no = 1:4
if SC_status(sc_no) == 1

notes = sprintf('Reading STAFF on S/C-%1.0f\0',sc_no); disp(notes);
sc_name = sprintf('%1.0f\0',sc_no);

[t_Bsc, Bsc] = isGetDataLite( db, t_start_staff, Dt_staff, ...
'Cluster', sc_name, 'staff', 'B_SC', 'Bx_By_Bz', f_B_str, ' ');
Bsc = real(double(Bsc));
B_spin = Bsc(1,:);

B_staff_vect = [ t_Bsc'; Bsc(3,:); Bsc(2,:) ]';
[t_ph, ph] = isGetDataLite( db, t_start, Dt_staff, ...
'Cluster', sc_name, 'ephemeris', 'phase', ' ', ' ', ' ');
ph = real( double(ph) );
ph_vect = [ t_ph ph ];

B_vect = EFW_despin( B_staff_vect, ph_vect );
B_vect = [ B_vect B_spin' ];

% Read enough data? If not, start looping...
n_staff = length(t_Bsc);
while n_samples > n_staff
% One minute later we try read again.

t_new_start = fromepoch(t_Bsc(length(t_Bsc)) + 60);
Dt_read = t_Bsc(length(t_Bsc)) - t_Bsc(1);
Dt_new = Dt_staff - Dt_read - 60;

[t_Bsc, Bsc] = isGetDataLite( db, t_new_start, Dt_new, ...
'Cluster', sc_name, 'staff', 'B_SC', 'Bx_By_Bz', f_B_str, ' ');
Bsc = real(double(Bsc));
B_spin = Bsc(1,:);

B_staff_vect = [ t_Bsc'; Bsc(3,:); Bsc(2,:) ]';
[t_ph, ph] = isGetDataLite( db, t_new_start, Dt_new, ...
'Cluster', sc_name, 'ephemeris', 'phase', ' ', ' ', ' ');
ph = real( double(ph) );
ph_vect = [ t_ph ph ];

B_tmp = EFW_despin( B_staff_vect, ph_vect );
B_tmp = [ B_tmp B_spin' ];

B_vect = cat( 1, B_vect, B_tmp );

n_staff = n_staff + length(t_Bsc);
end

command = sprintf('B_vect_%1.0f = B_vect;\0', sc_no); eval(command);

end
end

clear t_Bsc Bsc t_ph B_staff_vect ph_vect ph ...
Dt_staff n_read pot_read n_samples ...
B_spin Dt_new Dt_read n_staff t_new_start B_tmp ...
B_vect

disp('Calibrated B-field data read & despun!');

clear STAFF_status f_B_str f_E_str
clear command notes sc_name sc_no

end % STAFF status end

if INST_status(3) == 1
% FGM data.
%%%%%%%%%%%

disp('FGM B-field in DSC coordinates!');

for sc_no = 1:4
if SC_status(sc_no) == 1

sc_name = sprintf('C%1.0f\0',sc_no);
sc_name2 = sprintf('B_xyz_gse__C%1.0f_PP_FGM\0',sc_no);

[t_sc, B_fgm] = isGetDataLite( db, t_start, Dt, ...
'CSDS_PP', sc_name, 'FGM', sc_name2, ' ', ' ', ' ');

if isempty(t_sc)
B_fgm = NaN .* zeros(3,1);
t_sc = NaN;
end
B_fgm = real(double(B_fgm));
t_sc = real(double(t_sc));
B_fgm = [t_sc B_fgm'];
B_gse = B_fgm;

% Note: Latitude and Longitude of Spacecraft from SCDS on Web.
% Using tables & routines compiled by Andris Vaivads.
if not(isnan(B_fgm(1,1))), B_fgm =c_gse2dsc( B_fgm, [B_fgm(1,1) sc_no] );, end

command = sprintf('B_fgm_%1.0f = B_fgm;\0', sc_no); eval(command);
command = sprintf('B_gse_%1.0f = B_gse;\0', sc_no); eval(command);

end
end

clear t_sc command sc_name sc_name2 sc_no B_fgm B_gse

end % FGM status end

% Ephemeris.
%%%%%%%%%%%%
disp('Reading Ephemeris');

for sc_no = 1:4
if SC_status(sc_no) == 1
sc_name = sprintf('%1.0f\0',sc_no);

[t_R, R] = isGetDataLite( db, t_start, Dt, ...
'Cluster', sc_name, 'ephemeris', 'position', ' ', ' ', ' ');
R = real( double(R) );
R_tot = sqrt( R(1,:).^2 + R(2,:).^2 + R(3,:).^2 );
R_tot = R_tot ./ 6370;
R_vect = [ t_R R_tot' ];

[t_V, V] = isGetDataLite( db, t_start, Dt, ...
'Cluster', sc_name, 'ephemeris', 'velocity', ' ', ' ', ' ');
V = real( double(V) );
V_tot = sqrt( V(1,:).^2 + V(2,:).^2 + V(3,:).^2 );
V_vect = [ t_V V_tot' ];

%[t_MLT, MLT] = isGetDataLite( db, t_start, Dt, ...
% 'Cluster', sc_name, 'ephemeris', 'mlt', ' ', ' ', ' ');
%MLT = real( double(MLT) );
%MLT_vect = [ t_MLT MLT ];

[t_L, L] = isGetDataLite( db, t_start, Dt, ...
'Cluster', sc_name, 'ephemeris', 'l_shell', ' ', ' ', ' ');
L = real( double(L) );
L_vect = [ t_L L ];

[t_LAT, LAT] = isGetDataLite( db, t_start, Dt, ...
'Cluster', sc_name, 'ephemeris', 'inv_lat', ' ', ' ', ' ');
LAT = real( double(LAT) );
LAT_vect = [ t_LAT LAT ];

[t_ph, ph] = isGetDataLite( db, t_start, Dt, ...
'Cluster', sc_name, 'ephemeris', 'phase', ' ', ' ', ' ');
ph = real( double(ph) );
ph_vect = [ t_ph ph ];

command = sprintf('V_vect_%1.0f = V_vect;\0',sc_no); eval(command);
command = sprintf('L_vect_%1.0f = L_vect;\0',sc_no); eval(command);
command = sprintf('R_vect_%1.0f = R_vect;\0',sc_no); eval(command);
%command = sprintf('MLT_vect_%1.0f = MLT_vect;\0',sc_no); eval(command);
command = sprintf('LAT_vect_%1.0f = LAT_vect;\0',sc_no); eval(command);
command = sprintf('ph_vect_%1.0f = ph_vect;\0',sc_no); eval(command);

end
end

clear t_R R R_tot t_V V V_tot t_MLT MLT t_L L ...
t_LAT LAT ph_vect ...
V_vect L_vect R_vect MLT_vect LAT_vect ...
sc_name sc_no command

if INST_status(4) == 1
% CIS data.
%%%%%%%%%%%

disp('read CIS-PP data ...');

for sc_no = 1:4
if SC_status(sc_no) == 1

sc_name = sprintf('C%1.0f\0',sc_no);
sc_name2 = sprintf('V_HIA_xyz_gse__C%1.0f_PP_CIS\0',sc_no);

[t_sc, CIS_v] = isGetDataLite( db, t_start, Dt, ...
'CSDS_PP', sc_name, 'CIS', sc_name2, ' ', ' ', ' ');

if isempty(t_sc)
CIS_v = NaN .* zeros(3,1);
t_sc = NaN;
end

CIS_v = real(double(CIS_v));
t_sc = real(double(t_sc));
CIS_v = [t_sc CIS_v'];

if not(isnan(CIS_v(1,1))), CIS_v =c_gse2dsc( CIS_v, [CIS_v(1,1) sc_no] );, end

command = sprintf('CIS_v_%1.0f = CIS_v;\0',sc_no); eval(command);

end
end

clear sc_name sc_no command t_sc sc_name2 CIS_v

end % CIS status end

if INST_status(5) == 1
% PEACE data.
%%%%%%%%%%%%%%

disp('read PEACE-PP data ...');

for sc_no = 1:4
if SC_status(sc_no) == 1
sc_name = sprintf('C%1.0f\0',sc_no);
sc_name2 = sprintf('V_e_xyz_gse__C%1.0f_PP_PEA\0',sc_no);

[t_sc, PEACE_v] = isGetDataLite( db, t_start, Dt, ...
'CSDS_PP', sc_name, 'PEACE', sc_name2, ' ', ' ', ' ');

if isempty(t_sc)
PEACE_v = NaN .* zeros(3,1);
t_sc = NaN;
end

PEACE_v = real(double(PEACE_v));
t_sc = real(double(t_sc));
PEACE_v = [t_sc PEACE_v'];

command = sprintf('PEACE_v_%1.0f = PEACE_v;\0',sc_no); eval(command);
end
end

clear sc_name sc_no command t_sc sc_name2 PEACE_v

end % PEACE status end

if INST_status(6) == 1
% EDI E_perp data.
%%%%%%%%%%%%%%

disp('read EDI-PP data ...');

for sc_no = 1:4
if SC_status(sc_no) == 1
sc_name = sprintf('C%1.0f\0',sc_no);
sc_name2 = sprintf('E_xyz_gse__C%1.0f_PP_EDI\0',sc_no);

[t_sc, EDI_vect] = isGetDataLite( db, t_start, Dt, ...
'CSDS_PP', sc_name, 'EDI', sc_name2, ' ', ' ', ' ');

if isempty(t_sc)
EDI_vect = NaN .* zeros(3,1);
t_sc = NaN;
end

EDI_vect = real(double(EDI_vect));
t_sc = real(double(t_sc));
EDI_vect = [t_sc EDI_vect'];

% Note: Latitude and Longitude of Spacecraft from SCDS on Web.
% Using tables & routines compiled by Andris Vaivads.
if not(isnan(EDI_vect(1,1))), EDI_vect =c_gse2dsc( EDI_vect, [EDI_vect(1,1) sc_no] );, end

command = sprintf('EDI_vect_%1.0f = EDI_vect;\0',sc_no); eval(command);
end
end

clear sc_name sc_no command t_sc sc_name2 EDI_vect

end % EDI status end

clear reply INST_status % SC_status

% Do your favourite analysis!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('Ready for analysis!');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%