read_extender(filename)

function f = read_extender(filename)
%READ_EXTENDER(filename) Reads EXTENDER output files
%   READ_EXTENDER Reads the extended_mesh, extended_field, and points.out
%   files generated by the EXTENDER code and returns them as structures.
%
%   Note that EXTENDER outputs the magnetic field (H) but this is converted
%   to the magnetic induction (B) through multiplication by mu0.
%
%   Example:
%       extended_mesh=read_extender('extended_mesh.test');
%
%   Written by:     S.Lazerson ([e-mail])
%   Version:        1.0
%   Date:           5/11/11

mu0=4*pi*1e-7;
% First we check to see what type of file it
if strcmp('extended_mesh',filename(1:13))
    f.datatype='EXTENDER_MESH';
elseif strcmp('extended_field',filename(1:14))
    f.datatype='EXTENDER_FIELD';
elseif strcmp('.out',filename(length(filename)-3:length(filename)))
    f.datatype='EXTENDER_POINTS';
else
    f=-1;
    disp(['ERROR: Unknown filetype ' strtrim(filename)]);
    return
end
fid=fopen(filename,'r');
% Now we read the file based on datatype
switch f.datatype
    case 'EXTENDER_MESH'
        type=fgetl(fid);
        if ~strcmp('ConstrainedGeometry',strtrim(type))
            disp(['ERROR: ' strtrim(filename) ' is of type ' strtrim(type)]);
            disp('       Expected ConstrainedGeometry');
            f=-2;
        end
        temp=fscanf(fid,'%d',[1 2]);
        f.npoints=temp(1);
        f.ncells=temp(2);
        temp=fscanf(fid,'%f',[1 3]);
        f.aminor=temp(1);
        f.delta_phi=temp(2);
        f.zminor=temp(3);
        temp=fscanf(fid,'%f',[1 3]);
        f.rmin=temp(1);
        f.rmax=temp(2);
        f.zmax=temp(3);
        f.zmin=-temp(3);
        f.nfp=fscanf(fid,'%d',[1 1]);
        temp=fscanf(fid,'%d',[1 3]);
        f.nr=temp(1);
        f.nphi=temp(2);
        f.nz=temp(3);
        temp=fscanf(fid,'%d',[1 2]);
        f.napoints=temp(1);
        f.nacells=temp(2);
        temp=fscanf(fid,'%d',[1 2]);
        f.bpoints=temp(1);
        f.icells=temp(2);
        temp=fscanf(fid,'%d',[2 f.npoints]);
        f.point_activity=temp(1,:);
        f.point_address=temp(2,:);
        temp=fscanf(fid,'%d',[2 f.ncells]);
        f.cell_activity=temp(1,:);
        f.cell_address=temp(2,:);
        f.apoints=fscanf(fid,'%d',[1 f.napoints]);
        f.acells=fscanf(fid,'%d',[1 f.nacells]);
        % Create the grid
        f.raxis=f.rmin+f.aminor.*(0:1/(f.nr-1):1);
        f.zaxis=-f.zmax+f.zminor.*(0:1/(f.nz-1):1);
        f.phi=f.delta_phi.*(0:1/(f.nphi-1):1);
    case 'EXTENDER_FIELD'
        type=fgetl(fid);
        if ~strcmp('ConstrainedVField',strtrim(type))
            disp(['ERROR: ' strtrim(filename) ' is of type ' strtrim(type)]);
            disp('       Expected ConstrainedVField');
            f=-2;
        end
        f.napoints=fscanf(fid,'%d',[1 1]);
        f.nfp=fscanf(fid,'%d',[1 1]);
        f.nr=fscanf(fid,'%d',[1 1]);
        f.nphi=fscanf(fid,'%d',[1 1]);
        f.nz=fscanf(fid,'%d',[1 1]);
        f.rmin=fscanf(fid,'%f',[1 1]);
        f.rmax=fscanf(fid,'%f',[1 1]);
        f.zmax=fscanf(fid,'%f\n',[1 1]);
        f.parityline=strtrim(fgetl(fid));
        f.br=fscanf(fid,'%f',[1 f.napoints]).*mu0;
        f.bphi=fscanf(fid,'%f',[1 f.napoints]).*mu0;
        f.bz=fscanf(fid,'%f',[1 f.napoints]).*mu0;
    case 'EXTENDER_POINTS';
        f.index=[];
        f.r=[];
        f.phi=[];
        f.z=[];
        f.br=[];
        f.bphi=[];
        f.bz=[];
        f.modb=[];
        f.inside=[];
        while ~feof(fid)
            line=fgetl(fid);
            if isempty(strfind(line,'#'))
                temp=sscanf(line,'%d %f %f %f %f %f %f %f %d');
                f.index  = [ f.index  temp(1)];
                f.r      = [ f.r      temp(2)];
                f.phi    = [ f.phi    temp(3)];
                f.z      = [ f.z      temp(4)];
                f.br     = [ f.br     temp(5)];
                f.bphi   = [ f.bphi   temp(6)];
                f.bz     = [ f.bz     temp(7)];
                f.modb   = [ f.modb   temp(8)];
                f.inside = [ f.inside temp(9)];
            end
        end
        % Convert from H to B
        f.br=mu0.*f.br;
        f.bphi=mu0.*f.bphi;
        f.bz=mu0.*f.bz;
        % Calculate Cartesian components
        f.npoints=length(f.br);
        f.x=f.r.*cos(f.phi);
        f.y=f.r.*sin(f.phi);
        f.bx=f.br.*cos(f.phi)-f.bphi.*sin(f.phi);
        f.by=f.bphi.*cos(f.phi)+f.br.*sin(f.phi);
end
fclose(fid);
return
end