function read_in_data(infile)
% reads input data for mesh generator

global  numlps
global  meltyp
global  ndiv
global  join
global  xcor
global  ycor
global  rn
global  numel
global  lnr
global  lnp
global  np
global  xlc ylc nplc  % loop corner info
%     open input data file

fid=0;
if (infile)
    [fid,message] = fopen(infile,'r');
    if(fid==-1)
        disp(message);
    end
else
     filename = input('Mesh generation file name: ','s');
    [fid,message] = fopen(filename,'r');
    if(fid==-1)
        disp(message);
    end
end

%    READING AND INITIALIZATION OF DATA

meshsize = fscanf(fid,'%d',2);

numlps = meshsize(1);
meltyp = meshsize(2);
%
%
% FIX MELTYPE AT 3 FOR THE CURRENT IMPLEMENTATION
%
%
meltyp = 3;

if(numlps<=0)
    error('numlps must be greater than zero')
end

if(meltyp<=0)
    error('meltyp incorrect value')
end

maxsnp = 8;
maxside = 4;
maxdim = 2;

rn = zeros(1,maxsnp);
xcor = zeros(numlps,maxsnp);
ycor = zeros(numlps,maxsnp);
ndiv = zeros(numlps,maxside);
join = zeros(numlps,maxside,maxdim);

for i = 1:numlps

    ndiv_vals = fscanf(fid,'%d',2);

    ndiv(i,1) = ndiv_vals(1);
    ndiv(i,2) = ndiv_vals(2);

    join_vals = fscanf(fid,'%d',8);

    join(i,1,1) = join_vals(1);
    join(i,1,2) = join_vals(2);
    join(i,2,1) = join_vals(3);
    join(i,2,2) = join_vals(4);
    join(i,3,1) = join_vals(5);
    join(i,3,2) = join_vals(6);
    join(i,4,1) = join_vals(7);
    join(i,4,2) = join_vals(8);

    for j=1:8
        cor_vals = fscanf(fid,'%lf',2);
        xcor(i,j) = cor_vals(1);
        ycor(i,j) = cor_vals(2);
    end

end

numdiv2 = 0;
numbig  = 0;

for i=1:numlps
    ndiv(i,3)=ndiv(i,1);
    ndiv(i,4)=ndiv(i,2);
    numdiv2 = numdiv2+ndiv(i,1)*ndiv(i,2);
    numbig=max([numbig ndiv(i,1) ndiv(i,2)]);
end

if(meltyp==3)
    numel = 2*numdiv2;
end
if(meltyp==4)
    numel = numdiv2;
end
if(meltyp==6)
    numel = 2*numdiv2;
end
if(meltyp==9)
    numel = numdiv2;
end

lnr = zeros(numlps,maxside,numbig);
lnp = zeros(numlps,maxside,numbig);

np = zeros(numel,meltyp);

xlc=xcor(:,[1,3,5,7]); % x coords of loop corners
ylc=ycor(:,[1,3,5,7]); % y coords of loop corners
nplc=zeros(numlps,4);  % global node numbers of loop corners

status = fclose(fid);
