%
% erik_plot
%
%   Main program for making the overview plots
%   Data should already be loaded using routine erik_read_cluster
%
%   panel 1 - spinplane angle
%   panel 2 - perpendicular component in spinplane
%   panel 3 - perpendicular component out of spinplane
%   panel 4 - parallel component in spinplane
%   panel 5 - spacecraft potential
%   panel 6 - angle between Sun and B-field in spinplane
%

zoom = input('Zoomed in (1 = yes, 0 = no)?: ');
sc_no = input('Which S/C (1-4)?: ');
filter = input('LP filter freq (0 = no filter)?: ');
disp('Working...');

% Get the Cluster data (already loaded from disk)
command = sprintf('B_ds = B_fgm_%1.0f;\0',sc_no); eval(command);
command = sprintf('E_ds = E_vect_%1.0f;\0',sc_no); eval(command);
command = sprintf('EDI_ds = EDI_vect_%1.0f;\0',sc_no); eval(command);
command = sprintf('CIS_ds = CIS_v_%1.0f;\0',sc_no); eval(command);
command = sprintf('N_e = Ne_%1.0f;\0',sc_no); eval(command);
command = sprintf('SC_p = SCp_%1.0f;\0',sc_no); eval(command);
command = sprintf('ph = ph_vect_%1.0f;\0',sc_no); eval(command);
command = sprintf('E_12 = E_12_%1.0f;\0',sc_no); eval(command);
command = sprintf('E_34 = E_34_%1.0f;\0',sc_no); eval(command);
command = sprintf('LAT = LAT_vect_%1.0f;\0',sc_no); eval(command);
command = sprintf('R = R_vect_%1.0f;\0',sc_no); eval(command);

% Extract Invariant latitude and altitude
LAT_start = LAT(1, 2);
LAT_end = LAT(length(LAT), 2);
R_start = R(1, 2);
R_end = R(length(R), 2);

% Get the angles
theta_xy = erik_angle(B_ds, 'xy');
theta_x = erik_angle(B_ds, 'x');
theta_x = [E_ds(:,1) interp1(theta_x(:,1), theta_x(:,2), E_ds(:,1), ...
'linear', 'extrap')];

% Get date and time
start_time = E_ds(1,1);
date = datestr(datenum(fromepoch(start_time)));
date = date(1:11);

% Get EFW and EDI in FAC
E_pp = erik_Epp(B_ds, E_ds);
E_pp_EDI = erik_Epp(B_ds, EDI_ds);

% Calculate E-field from CIS data (E=-VxB)
E_pp_CIS = av_cross(CIS_ds, B_ds);
E_pp_CIS(:, 2) = E_pp_CIS(:, 2) * (-1e-3);
E_pp_CIS(:, 3) = E_pp_CIS(:, 3) * (-1e-3);
E_pp_CIS(:, 4) = E_pp_CIS(:, 4) * (-1e-3);

% Get CIS in FAC
E_pp_CIS = erik_Epp(B_ds, E_pp_CIS);

% Additional plots if zoom selected
if zoom == 1
  figure(2);
  clf;
  ph = erik_fixphase(ph, E_pp(:, 1));
  h2(1) = subplot(4, 1, 1); plot( ph(:, 2), E_pp(:, 5));
  set(gca ,'XTick', 0:45:360);
  curr_axis = axis;
  axis([curr_axis(1) curr_axis(2) -7 7]);
  ylabel('mV/m');
  strTitle = sprintf('phase vs E_p_a_r_a on S/C-%1.0f', sc_no);
  title(strTitle);
  theta = theta_x(:, 2) + ph(:, 2);
  h2(2) = subplot(4, 1, 2); plot(theta, E_pp(:, 5));
  set(gca ,'XTick', 0:45:360);
  curr_axis = axis;
  axis([curr_axis(1) curr_axis(2) -7 7]);
  ylabel('mV/m');
  strTitle = sprintf('Angle between B-field and spinning x-axis ...
vs E_p_a_r_a on S/C-%1.0f', sc_no);
  title(strTitle);


  h2(3) = subplot(4, 1, 3); plot(E_12(:,1), E_12(:, 2));
  curr_axis = axis;
  axis([curr_axis(1) curr_axis(2) -10 10]);
  ylabel('Volts');
  strTitle = sprintf('S/C-%1.0f, raw probe data',sc_no);
  title(strTitle);
  h2(4) = subplot(4, 1, 4); plot(E_34(:,1), E_34(:, 2));
  curr_axis = axis;
  axis([curr_axis(1) curr_axis(2) -10 10]);
  ylabel('Volts');
  add_timeaxis(h2([3 4]), date);
  legend(h2(3), 'E_1_2', 0);
  legend(h2(4), 'E_3_4', 0);

end

% Filter the data to get rid of wake effects
if filter ~= 0
  E_pp = av_tsfilt(E_pp, 0, filter, [], 3);
end

% Main plots
figure(1);
clf;
h(1) = subplot(6, 1, 1); plot(theta_xy(:,1), theta_xy(:, 2));
ylabel('[Degrees]');
grid;
strTitle = sprintf('S/C-%1.0f, ILAT=[%0.05g, %0.05g], R=[%0.05g, ...
%0.05g], filter=%0.5g\0',...
sc_no, LAT_start, LAT_end, R_start, R_end, filter);
if zoom == 1
  strTitle = [strTitle ', ZOOMED']
end
title(strTitle);

h(2) = subplot(6, 1, 2);
hold on;
plot(E_pp(:, 1), E_pp(:, 2));
plot(E_pp_EDI(:, 1), E_pp_EDI(:, 2), '+r');
plot(E_pp_CIS(:, 1), E_pp_CIS(:, 2), '+g');
hold off;
curr_axis = axis;
axis([curr_axis(1) curr_axis(2) -20 20]);
ylabel('E_p_e_r_p_1 [mV/m]');
grid;

h(3) = subplot(6, 1, 3);
hold on;
plot(E_pp(:, 1), E_pp(:, 4));
plot(E_pp_EDI(:, 1), E_pp_EDI(:, 3), '+r');
plot(E_pp_CIS(:, 1), E_pp_CIS(:, 3), '+g');
hold off;
curr_axis = axis;
axis([curr_axis(1) curr_axis(2) -20 20]);
ylabel('E_p_e_r_p_2 [mV/m]');
grid;

h(4) = subplot(6, 1, 4);
hold on;
plot(E_pp(:, 1), E_pp(:, 5));
plot(E_pp_EDI(:, 1), E_pp_EDI(:, 4), '+r');
plot(E_pp_CIS(:, 1), E_pp_CIS(:, 4), '+g');
hold off;
curr_axis = axis;
axis([curr_axis(1) curr_axis(2) -3 3]);
ylabel('E_p_a_r_a [mV/m]');
grid;

h(5) = subplot(6, 1, 5); plot(SC_p(:, 1), SC_p(:, 2));
ylabel('[Volts]');
grid;

h(6) = subplot(6, 1, 6); plot(theta_x(:, 1), theta_x(:, 2));
ylabel('[Degrees]');
grid;

add_timeaxis(h, date);
legend(h(1), 'Angle between spin plane and B-field', 0);
legend(h(2), 'EFW', 'EDI', 'CIS', 0);
legend(h(3), 'EFW if E_p_a_r_a = 0', 'EDI', 'CIS', 0);
legend(h(4), 'EFW if E_z = 0', 'EDI if E_p_e_r_p_2 = 0', ...
'CIS if E_p_e_r_p_2 = 0', 0);
legend(h(5), 'S/C potential', 0);
legend(h(6), 'Angle between Sun and B-field', 0);

clear command theta_xy theta_x theta B_ds E_ds N_e E_pp E_pp_EDI filter ...
start_time date ph E_12 E_34 SC_p LAT R;
