function [E_pp] = erik_Epp(B_ds, E_ds)
%
% function [E_pp] = erik_Epp(B_ds, E_ds)
%
% Input:  B and E-fields in despun coordinates
% Output: E-field in parallel and perpendicular coord, 
%   if EFW (E_ds miss the z component)
%   [E_pp] = [t  E_perp1  E_perp2  E_perp2_para=0  E_para_Ez=0  E_para_perp2=0]
%   if CIS or EDI (all components present)
%   [E_pp] = [t  E_perp1  E_perp2 E_para_perp2=0]
%
%   All references are to Parallel Electric Fields in Space
%   thesis paper by Erik Bergman, 2003
%
%   See Table 4.1 for definitions
%   E_perp1			if E_z = 0
%   E_perp2  		if E_z = 0
%   E_perp2_para=0  	if E_para = 0
%   E_para_Ez=0  		if E_z = 0
%   E_para_perp2=0	if E_perp2 = 0
%


%See Eq. (4.1)
B_perp1 = av_norm(av_cross(B_ds, [0 0 1]));
%See Eq. (4.2)
B_perp2 = av_norm(av_cross(B_perp1, B_ds));
%See Eq. (4.3)
B_para = av_norm(B_ds);

% Interpolate the B data so we have the same number of points as the E data
B_perp1_h = [E_ds(:,1) interp1(B_perp1(:,1), B_perp1(:,[2 3 4]), E_ds(:,1), ...
'linear', 'extrap')];
B_perp2_h = [E_ds(:,1) interp1(B_perp2(:,1), B_perp2(:,[2 3 4]), E_ds(:,1), ...
'linear', 'extrap')];
B_para_h = [E_ds(:,1) interp1(B_para(:,1), B_para(:,[2 3 4]), E_ds(:,1), ...
'linear', 'extrap')];

% If E_ds miss z-comp then from EFW
if length(E_ds(1, :)) == 3
  % Assuming E_z = 0
  E_ds(:, 4) = 0;
  E_perp1 = av_dot(B_perp1_h, E_ds);
  E_perp2 = av_dot(B_perp2_h, E_ds);
  E_para = av_dot(B_para_h, E_ds);

  E_pp = E_perp1;
  E_pp(:, 3) = E_perp2(:, 2);
  E_pp(:, 5) = E_para(:, 2);

  % Assuming E_para = 0, see Eq. (4.6)
  E_ds(:, 4) = -(E_ds(:, 2) .* B_para_h(:, 2) + E_ds(:, 3) .* ...
B_para_h(:, 3)) ./ B_para_h(:, 4);
  E_perp2 = av_dot(B_perp2_h, E_ds);
  E_pp(:, 4) = E_perp2(:, 2);

  % Assuming E_perp2 = 0, see Eq. (4.10)
  E_ds(:, 4) = -(E_ds(:, 2) .* B_perp2_h(:, 2) + E_ds(:, 3) .* ...
B_perp2_h(:, 3)) ./ B_perp2_h(:, 4);
  E_para = av_dot(B_para_h, E_ds);
  E_pp(:, 6) = E_para(:, 2);

% Else from CIS or EDI
else
  E_perp1 = av_dot(B_perp1_h, E_ds);
  E_perp2 = av_dot(B_perp2_h, E_ds);
  
  % Assuming E_perp2 = 0, see Eq. (4.10)
  E_ds(:, 4) = -(E_ds(:, 2) .* B_perp2_h(:, 2) + E_ds(:, 3) .* ...
B_perp2_h(:, 3)) ./ B_perp2_h(:, 4);
  E_para = av_dot(B_para_h, E_ds);
  
  E_pp = E_perp1;
  E_pp(:, 3) = E_perp2(:, 2);
  E_pp(:, 4) = E_para(:, 2);
end
