function [Pinfo] = loadfile(filename)

pathname = pwd;

[fpath, fname, extfile] = fileparts(fullfile(pathname, filename));

rda_filename = fullfile(pathname, filename);
fid = fopen(rda_filename, 'r');


rda.sparname = fname;
coordfilespath = fullfile(fpath, 'lcmodel_outputs', fname);
rda.coordpath = coordfilespath;
rda.spectrapath = fpath;
rdafilepath = fullfile(fpath, fname);
rda.sparfilepath = rdafilepath;

exam1 = strsplit(rda.spectrapath, 'exam_1');
if ~isempty(exam1)
    rda.mainpath = exam1{1};
else
    rda.mainpath = rda.spectrapath;
end

patientno = strsplit(rda.mainpath, filesep);
if length(patientno) > 1
    rda.name = patientno{end-1};
else
    rda.name = fname;
end

rdaheadinfo = textread([pathname,'/',filename], '%s');

head_start_text = '>>> Begin of header <<<';
head_end_text   = '>>> End of header <<<';

tline = fgets(fid);

while (isempty(strfind(tline , head_end_text)))

    tline = fgets(fid);

    if ( isempty(strfind (tline , head_start_text)) + isempty(strfind (tline , head_end_text )) == 2)


        % Store this data in the appropriate format

        occurence_of_colon = findstr(':',tline);
        variable = tline(1:occurence_of_colon-1) ;
        value    = tline(occurence_of_colon+1 : length(tline)) ;

        switch variable
            case { 'PatientID' , 'PatientName' , 'StudyDescription' , 'PatientBirthDate' , 'StudyDate' , 'StudyTime' , 'PatientAge' , 'SeriesDate' , ...
                    'SeriesTime' , 'SeriesDescription' , 'ProtocolName' , 'PatientPosition' , 'ModelName' , 'StationName' , 'InstitutionName' , ...
                    'DeviceSerialNumber', 'InstanceDate' , 'InstanceTime' , 'InstanceComments' , 'SequenceName' , 'SequenceDescription' , 'Nucleus' ,...
                    'TransmitCoil' }
                eval(['rda.' , variable , ' = value; ']);
            case { 'PatientSex' }
                % Sex converter! (int to M,F,U)
                switch value
                    case 0
                        rda.sex = 'Unknown';
                    case 1
                        rda.sex = 'Male';
                    case 2

                        rda.sex = 'Female';
                end

            case {  'SeriesNumber' , 'InstanceNumber' , 'AcquisitionNumber' , 'NumOfPhaseEncodingSteps' , 'NumberOfRows' , 'NumberOfColumns' , 'VectorSize' }
                %Integers
                eval(['rda.' , variable , ' = str2num(value); ']);
            case { 'PatientWeight' , 'TR' , 'TE' , 'TM' , 'DwellTime' , 'NumberOfAverages' , 'MRFrequency' , 'MagneticFieldStrength' , 'FlipAngle' , ...
                    'SliceThickness' ,  'FoVHeight' , 'FoVWidth' , 'PercentOfRectFoV' , 'PixelSpacingRow' , 'PixelSpacingCol'}
                value = strrep(value, ',', '.');
                %Floats
                eval(['rda.' , variable , ' = str2num(value); ']);
            case {'SoftwareVersion[0]' }
                rda.software_version = value;
            case {'CSIMatrixSize[0]' }
                rda.CSIMatrix_Size(1) = str2num(value);
            case {'CSIMatrixSize[1]' }
                rda.CSIMatrix_Size(2) = str2num(value);
            case {'CSIMatrixSize[2]' }
                rda.CSIMatrix_Size(3) = str2num(value);
            case {'PositionVector[0]' }
                value = strrep(value, ',', '.');
                out.PositionVector(1) = str2num(value);
            case {'PositionVector[1]' }
                value = strrep(value, ',', '.');
                out.PositionVector(2) = str2num(value);
            case {'PositionVector[2]' }
                value = strrep(value, ',', '.');
                out.PositionVector(3) = str2num(value);
            case {'RowVector[0]' }
                value = strrep(value, ',', '.');
                out.RowVector(1) = str2num(value);
            case {'RowVector[1]' }
                value = strrep(value, ',', '.');
                out.RowVector(2) = str2num(value);
            case {'RowVector[2]' }
                value = strrep(value, ',', '.');
                out.RowVector(3) = str2num(value);
            case {'ColumnVector[0]' }
                value = strrep(value, ',', '.');
                out.ColumnVector(1) = str2num(value);
            case {'ColumnVector[1]' }
                value = strrep(value, ',', '.');
                out.ColumnVector(2) = str2num(value);
            case {'ColumnVector[2]' }
                value = strrep(value, ',', '.');
                out.ColumnVector(3) = str2num(value);

            otherwise
                % We don't know what this variable is.  Report this just to keep things clear
                %disp(['Unrecognised variable ' , variable ]);
        end

    else
        % Don't bother storing this bit of the output
    end

end


%
% So now we should have got to the point after the header text
%
% Siemens documentation suggests that the data should be in a double complex format (8bytes for real, and 8 for imaginary?)
%

bytes_per_point = 16;
complex_data = fread(fid , rda.CSIMatrix_Size(1) * rda.CSIMatrix_Size(1) *rda.CSIMatrix_Size(1) *rda.VectorSize * 2 , 'double');
%fread(fid , 1, 'double');  %This was a check to confirm that we had read all the data (it passed!)
fclose(fid);

% Now convert this data into something meaningful

%Reshape so that we can get the real and imaginary separated
hmm = reshape(complex_data,  2 , rda.VectorSize , rda.CSIMatrix_Size(1) ,  rda.CSIMatrix_Size(2) ,  rda.CSIMatrix_Size(3) );

%Combine the real and imaginary into the complex matrix
fids = complex(hmm(1,:,:,:,:),hmm(2,:,:,:,:));
for rows=1:rda.CSIMatrix_Size(1)
    for cols=1:rda.CSIMatrix_Size(2)
        fids(:,:,rows,cols)=fids(:,:,rows,cols)';
    end
end

specs = fftshift(fft(fids,[],2),2);



% make calculations for the output mrs structure
sz = size(specs);
dwelltime = rda.DwellTime/1000000;
spectralwidth=1/dwelltime;
txfrq = rda.MRFrequency*1000000;
dims.t = 1;
dims.subSpecs = 0;
Bo = rda.MagneticFieldStrength;
rawAverages = rda.NumberOfAverages;
averages = 1;
subspecs =1;
rawSubspecs = 'na';
date = rda.StudyDate;
seq = rda.SequenceDescription;
TE = rda.TE;
TR = rda.TR;
pointsToLeftShift = 'N/A';

%Calculate t and ppm arrays using the calculated parameters:
f=[(-spectralwidth/2)+(spectralwidth/(2*sz(1))):spectralwidth/(sz(1)):(spectralwidth/2)-(spectralwidth/(2*sz(1)))];
ppm=-f/(Bo*42.577);
ppm=ppm+4.6082;

t=[0:dwelltime:(sz(1)-1)*dwelltime];

%FILLING IN DATA STRUCTURE
out.fids=fids;
out.specs=specs;
out.sz=sz;
out.ppm=ppm;
out.t=t;
out.spectralwidth=spectralwidth;
out.dwelltime=dwelltime;
out.txfrq=txfrq;
out.date=date;
out.dims=dims;
out.Bo=Bo;
out.averages=averages;
out.rawAverages=rawAverages;
out.subspecs=subspecs;
out.rawSubspecs=rawSubspecs;
out.seq=seq;
out.te=TE;
out.tr=TR;
out.pointsToLeftshift=pointsToLeftShift;
out.Header=rda;
out.slithickness= rda.SliceThickness;
out.nrow= rda.CSIMatrix_Size(1);
out.ncol= rda.CSIMatrix_Size(2);
out.nsli= rda.CSIMatrix_Size(3);
out.FOV = rda.CSIMatrix_Size(1)*rda.PixelSpacingRow;

finderinfo=find(ismember(rdaheadinfo, 'MidSlabPosition[1,1]:')==1);
out.voxoffap = str2num(rdaheadinfo{finderinfo+1});
finderinfo=find(ismember(rdaheadinfo, 'MidSlabPosition[0,0]:')==1);
out.voxofflr = str2num(rdaheadinfo{finderinfo+1});
finderinfo=find(ismember(rdaheadinfo, 'MidSlabPosition[2,2]:')==1);
out.voxoffcc = str2num(rdaheadinfo{finderinfo+1});


finderinfo=find(ismember(rdaheadinfo, 'RowVector[1]:')==1);
finderinfo= strrep(rdaheadinfo{finderinfo+1}, ',', '.');
ap_row= str2num(finderinfo);

finderinfo=find(ismember(rdaheadinfo, 'RowVector[0]:')==1);
finderinfo= strrep(rdaheadinfo{finderinfo+1}, ',', '.');
lr_row = str2num(finderinfo);

finderinfo=find(ismember(rdaheadinfo, 'RowVector[2]:')==1);
finderinfo= strrep(rdaheadinfo{finderinfo+1}, ',', '.');
cc_row = str2num(finderinfo);

finderinfo=find(ismember(rdaheadinfo, 'ColumnVector[1]:')==1);
finderinfo= strrep(rdaheadinfo{finderinfo+1}, ',', '.');
ap_col= str2num(finderinfo);

finderinfo=find(ismember(rdaheadinfo, 'ColumnVector[0]:')==1);
finderinfo= strrep(rdaheadinfo{finderinfo+1}, ',', '.');
lr_col = str2num(finderinfo);

finderinfo=find(ismember(rdaheadinfo, 'ColumnVector[2]:')==1);
finderinfo= strrep(rdaheadinfo{finderinfo+1}, ',', '.');
cc_col = str2num(finderinfo);

r = [lr_row ap_row cc_row]; % RowVector (Readout / Right)
c = [lr_col ap_col cc_col]; % ColumnVector (Phase / Down)
n = cross(r, c);

% Calculate Angles in Radians
% Angle LR
% Formula: atan2(c_z, n_z)
out.voxanglr = atan2(c(3), n(3));

% Angle AP
% Formula: atan2(-r_z, sqrt(c_z^2 + n_z^2))
out.voxangap = atan2(-r(3), sqrt(c(3)^2 + n(3)^2));

% Angle CC
% Formula: atan2(r_y, r_x)
out.voxangcc = atan2(r(2), r(1));


finderinfo=find(ismember(rdaheadinfo, 'SlabThickness[2]:')==1);
out.apVOI= str2num(rdaheadinfo{finderinfo+1});
finderinfo=find(ismember(rdaheadinfo, 'SlabThickness[1]:')==1);
out.lrVOI = str2num(rdaheadinfo{finderinfo+1});
finderinfo=find(ismember(rdaheadinfo, 'SlabThickness[0]:')==1);
out.ccVOI = str2num(rdaheadinfo{finderinfo+1});


Pinfo=out;
Pinfo.ext='rda';

% Auto-detect MRI file
if exist(fullfile(pwd, 'sub1_Bet.nii'), 'file')
    Pinfo.MRI = fullfile(pwd, 'sub1_Bet.nii');
elseif exist(fullfile(pwd, 'sub1_Bet.nii.gz'), 'file')
    Pinfo.MRI = fullfile(pwd, 'sub1_Bet.nii.gz');
else
    warning('MRI file not found. Please set Pinfo.MRI manually.');
end


%FILLING IN THE FLAGS
out.flags.writtentostruct=1;
out.flags.gotparams=1;
out.flags.leftshifted=0;
out.flags.filtered=0;
out.flags.zeropadded=0;
out.flags.freqcorrected=0;
out.flags.phasecorrected=0;
out.flags.averaged=1;
out.flags.addedrcvrs=1;
out.flags.subtracted=0;
out.flags.writtentotext=0;
out.flags.downsampled=0;
if out.dims.subSpecs==0
    out.flags.isFourSteps=0;
else
    out.flags.isFourSteps=(out.sz(out.dims.subSpecs)==4);
end



end