% Read in the first hdf file to obtain grid and other attributes 
    sd_id = hdfsd( 'start', [filenames(1,:)], 'rdonly');
    alpha = hdfsd('readattr', sd_id, 5);
    beta  = hdfsd('readattr', sd_id, 6);
    R     = hdfsd('readattr', sd_id, 7);
   
    % Read in the X grid
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Grid X'));
    X = hdfsd('readdata',sds_id,zeros(1,2),[],[181; 27]);
    
    % Read in the Y grid
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Grid Y'));
    Y = hdfsd('readdata',sds_id,zeros(1,2),[],[181 27]);
   
    % Calculate Radius
    radius = sqrt(X.^2 + Y.^2);
    
    % Calculate co-latitude
    theta = atan2(y,x)
    % Note: the following syntax might not work... May need to explicitly loop over theta for periodic point:
    theta[theta < 0] = theta[theta < 0] + 2*pi

    hdfsd('end',sd_id);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Now allocate & read all the data stored in a list of MIX output files:     
    f_gridded = zeros(2,length(filenames),27,181);
    x_gridded = f_gridded; 
    y_gridded = f_gridded;
   
    for i=1:length(filenames)
    
    sd_id = hdfsd( 'start', [filenames(i,:)], 'rdonly');
    stime(i) = datenum(double(hdfsd('readattr', sd_id, 2)));
    
    % Read in energy: keV (kilo - electron volt)
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Average energy North [keV]'));
    ENERGY = hdfsd('readdata',sds_id,zeros(1,2),[],[181 27]);
 
    % Read in flux
    sds_id = hdfsd( 'select', sd_id,hdfsd('nametoindex',sd_id,'Number flux North [1/cm^2 s]'));
    FLUX = hdfsd('readdata',sds_id,zeros(1,2),[],[181 27]);
    
    f_gridded(1,i,:,:) = ENERGY;
    f_gridded(2,i,:,:) = FLUX;

    %x_gridded(1,i,:,:,1) = X;
    %x_gridded(2,i,:,:,1) = X;
    %y_gridded(1,i,:,:,1) = Y;
    %y_gridded(2,i,:,:,1) = Y;
    
    %Close the file buffer when you are finished.
    hdfsd('end',sd_id);

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