function JouleHeating = JouleHeating_MIX(xi,yi,pot,ped)

% This function estimates the Joule heating using the ionospheric output
% (MIX) from the LFM-MIX/CMIT simulation. The Joule Heating in the
% ionosphere is defined as
%             JouleHeating = SigmaP * E^2
% where SigmaP is the Pedersen conductance in the ionosphere, E is the
% electric field defined as E = -grad(potential)
%
% The inputs are:
%            (xi,yi) : the MIX grid
%             pot    : the ionospheric potential from MIX [V]
%             ped    : the Pedersen conductance from MIX  [mho]
% The output is :
%             JouleHeating : Joule Heating in the ionosphere(on the same
%                            grid (xi,yi) ) [mW/m^2]
%
% Example:
%         filename='C:\temp\oldgc1_mix_2000-01-01T02-00-00Z.hdf'
%         xi = double(hdfread(filename, '/Grid X', 'Index', {[1  1],[1  1],[181   27]}));
%         yi = double(hdfread(filename, '/Grid Y', 'Index', {[1  1],[1  1],[181   27]}));
%         pot = double(hdfread(filename, '/Potential North [V]', 'Index', {[1  1],[1  1],[181   27]}));
%         ped = double(hdfread(filename, '/Pedersen conductance North [S]', 'Index', {[1  1],[1  1],[181   27]}));
%         JH = JouleHeating_MIX(xi,yi,pot,ped);

dtheta=1.5;                       % Define a regular grid (X,Y) for interpolation
theta1=pi*(-50:dtheta:50)/180;    % On this regular grid, grid points are evenly spaced
theta2=pi*(-50:dtheta:50)/180;    % The difference between grid points is dtheta, default value is 1.5 degrees,
[T1 T2]=meshgrid(theta1,theta2);  % which is close to the MIX resolution (dtheta should not be too small)

Rion=1.02;                        % Relative radius of the Earth ionosphere
RE=6.35e6;                        % Radius of the earth [m]
X=Rion.*sin(T1);                  % Cartesien grid X
Y=Rion.*sin(T2);                  % Cartesien grid Y

POT=griddata(xi,yi,pot,X,Y);      % Potential on the interpolated grid (X,Y)
PED=griddata(xi,yi,ped,X,Y);      % Pedersen conductance on the interpolated grid (X,Y)
warning off;

[Ex Ey]=gradient(-POT);           % E = - grad(potential)
Ex=Ex/(RE*Rion*dtheta*pi/180);    % Ex -- [V/m^2]
Ey=Ey/(RE*Rion*dtheta*pi/180);    % Ey -- [V/m^2]

Joule=(Ex.^2+Ey.^2).*PED;         % Joule heating = sigmaP * E^2 [W/m^2]

JouleHeating=griddata(X,Y,Joule*1e3,xi,yi); % Interpolation back to the original grid
                                            % Note the 1e3 factor here
                                            % makes the JouleHeating have the unit [mW/m^2]
                                            % Now JouleHeating variable is
                                            % on the grid (xi,yi)

[Nx Ny]=size(JouleHeating);                 % Set the NaN points to be zero for plotting
for i=1:Nx
    for j=1:Ny
        if (isnan(JouleHeating(i,j))==1)
            JouleHeating(i,j)=0;
        end
    end
end

end
