   % 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);
