function [k,f,imf,vg,imk,ll,vl]=gme

imf=[];vg=[];imk=[];ll=[];vl=[];
%%%%%%%%%%%%%%%%%%%%%%%%% User defined parameters
code='pwe2d';                                                   % 2D Plane Wave Expansion method; Uncomment if chosen
%code='gme';                                                    % Guide mode expansion method (3D); Uncomment if chosen

nrun=2;                                                         % Nmber of runs
jrun=[0 1];                                                     % Vector with j values for each run
%sy='kz';                                                       % Change vertical symmetry in runs; Uncomment if chosen
sy='xy';                                                        % Change horizontal symmetry in runs; Uncomment if chosen

ofilename={'pwe2dag_j0.out','pwe2dag_j1.out'};                  % Output file names (number of filenames=nrun)

folderb='/Users/alejandrogiacomotti/DisqueD/matlab/gme/';        % Choose folder for binaries (gmeag.e and pwe2dag.e)
folder='/Users/alejandrogiacomotti/DisqueD/matlab/gme/2DCylinders/';    % Choose folder

nbands=10;
%%%%%%%%%%%%%%%%%%%%%%%%% Define GME parameters

jlattice=1;                                                                                 %jlattice
jbasis=1;       alength1=0.38;      alength2=0.;        alength3=0.;        alength4=0.;    %jbasis,alength1,alength2,alength3,alength4
eps1=9.0;       eps2=1.0;                                                                  %eps1,eps2=dielectric constants (must be REAL and POSITIVE)
jparxy=0;       jparkz=-1;                                                                  %jparxy=0 H/TE (even),=1 E/TM (odd), =-1 both; jparkz=-1,0,1
akzmin=0.;      akzmax=0.;          nkz=0;                                                  %akzmin,akzmax,nkz=values of kz (used only if jparxy=-1)
nk=20;                                                                                      %nk=number of k-point subdivisions along each line in BZ 
npw=97;                                                                                    %npw=number of plane waves in expansion
jeigenv=0;      jkeigenv=1;         jstate=5;                                               %jeigenv=0 eigenvalues only, =1 eigenvectors;jkeigenv,jstate
factor=1.;                                                                                  %factor=multiplication factor for bands
xmin=-3.;       xmax=3.;            ymin=-3.;           ymax=3.;            ngrid=50;       %xmin,xmax,ymin,ymax,ngrid=for plot of field              
jfreephotons=0;                                                                             %jfreephotons=1 free photons (average epsilon)
jloop=0;        varmin=0.;          varmax=1.;        nvar=0;                               %jloop=0 no loop, =1 loop; varmin,varmax,nvar=values of loop variable
path=1231;                                                                                  %path (added by AG) *

% values of jlattice and jbasis:
%  jlattice=1 --> square lattice:
%    jbasis=1 circular rods of eps1 and radius=alength in material with eps2
%    jbasis=2 square rods of eps1 and side=alength (non rotated) in eps2
%    jbasis=3 square rods of eps1 and side=alength1 (45 deg rotated) in eps2
%    jbasis=99 1D lattice of stripes of length alength1 and eps1 in eps2
%  jlattice=2 --> triangular lattice:
%    jbasis=1 circular rods of eps1 and radius=alength in material with eps2
%    jbasis=2 graphite structure of circular rods of eps1 in eps2
%  jlattice=4 --> supercell for W1 waveguide, period along x=aret, along y=aret*(alength2+alength3)
%    jbasis=0 no defect: folded bands of 2D triangular lattice (my routine)
%    jbasis=1  linear defect: alength2=1 --> W1 waveguide (Agio routine)
%    jbasis=11 linear defect: alength2=1 --> W1 waveguide (my routine)
%    jbasis=2 no defect: folded bands of 2D triangular lattice (Agio routine)
% 
% possible values of npw consistent with symmetry: 
% square --> 1 9 13 25 29 37 45 49 57 61 69 81 89 97      145      293
% triang --> 1 7 13 19 31 37 43 55 61 73 85 91 97 109 121 127 139 151 163 169 187 199 211 223
% W1 waveguide --> depends on supercell period along GM!
%
% * Variable path is given by a sequence of the following special points
%
% 1D slabs (jlattice=1, jbasis=99): Gamma=1; X=2; Y=3; M=4;
% Paths in 1D: 12,13 (only for GME)
%
% Square lattice (jlattice=1): Gamma=1; X=2; M=3;
% Paths in square 2D lattice: 12 (only for GME), 21, 13, 213, 2132, 1231
%
% Triangular lattice (jlattice=2): Gamma=1; M=2; K=3;
% Paths in triangular 2D lattice: 21, 13, 213, 2132, 1231, 1321 (only for GME)
%
% Paths in W1 (jlattice=4): 12

%%%%%%%%%%%%%%%%%%%%%%%%% Run GME
cd(folder)                                                      % Change to folder
system(['cp ' folderb code 'ag.e ' folder code 'ag.e'])         % Copy binary to current folder  
ifilename=[code 'ag.dat'];                                      % Input file name
for n=1:nrun
    if (nrun>1 & sy=='kz') jparkz=jrun(n);end                   % Sets j value for each run: run over different kz symmetry
    if (nrun>1 & sy=='xy') jparxy=jrun(n);end                   % Sets j value for each run: run over different xy symmetry
    wdatafile(ifilename);    
    system(['./' code 'ag.e'])                                  % Run binary
    system(['cp ' code 'ag.out ' ofilename{n}])                     
end;

%%%%%%%%%%%%%%%%%%%%%%%%% Read data    
if (code(1:3)=='gme') ll=textread([folder 'gmeag_ll.out']);end;
vl=textread([folder code 'ag_vl.out']);

if (code(1:3)=='gme')
for n=1:nrun
fid = fopen([folder ofilename{n}]);
ycell=textscan(fid,'%n %n %n %n %n','headerLines',7);
    for j=1:nbands
        k(n,:,j)=cell2mat(ycell(1));
        f(n,:,j)=cell2mat(ycell(2));
        imf(n,:,j)=cell2mat(ycell(3));
        vg(n,:,j)=cell2mat(ycell(4));
        imk(n,:,j)=cell2mat(ycell(5));
        textscan(fid,'%*[^\n]',1);
        ycell=textscan(fid,'%n %n %n %n %n');
    end; 
    fclose(fid);
end
else
for n=1:nrun
fid = fopen([folder ofilename{n}]);
ycell=textscan(fid,'%n %n','headerLines',6);
    for j=1:nbands
        k(n,:,j)=cell2mat(ycell(1));
        f(n,:,j)=cell2mat(ycell(2));
        textscan(fid,'%*[^\n]',1);
        ycell=textscan(fid,'%n %n %n %n %n');
    end; 
    fclose(fid);
end    
end

%%%%%%%%%%%%%%%%%%%%%%%%% Plot data  
% Comment unwanted lines

% Band diagrams: two runs superimposed 
figure(1);
%plot(squeeze(k(1,:,1)),squeeze(f(1,:,:)),'ko',squeeze(k(2,:,1)),squeeze(f(2,:,:)),'k.') % point plot
plot(squeeze(k(1,:,1)),squeeze(f(1,:,:)),'r',squeeze(k(2,:,1)),squeeze(f(2,:,:)),'b') % line plot
hold on 
%plot(ll(:,1),ll(:,2),'k--')                                    % Light line: comment if unwanted         
plot(vl(1:2,1),vl(1:2,2),'k--')
plot(vl(3:4,1),vl(3:4,2),'k--')
hold off
ylim([0 0.7])  
set(gca,'XTick',[])
text(0,0,'\Gamma','VerticalAlignment','Top','HorizontalAlignment','Center','fontsize',14);
text(vl(1,1),0,'X','VerticalAlignment','Top','HorizontalAlignment','Center','fontsize',14);
text(vl(3,1),0,'M','VerticalAlignment','Top','HorizontalAlignment','Center','fontsize',14);
text(1,0,'\Gamma','VerticalAlignment','Top','HorizontalAlignment','Center','fontsize',14);
set(gca,'fontsize',14')
ylabel('\omegaa/(2\pic)');
%%%%%%%%%%%%%%%%%%%%%%%%% End of main routine

function wdatafile(file)
fid = fopen(file, 'w');
file
if (code(1:3)=='gme') 
fprintf(fid, '%d\n', jlattice);
fprintf(fid, '%d %.2f %.2f %.2f %.2f\n',jbasis,alength1,alength2,alength3,alength4);
fprintf(fid, '%.2f %.2f\n', eps1_clad1,eps2_clad1);
fprintf(fid, '%.2f %.2f\n', eps1_core,eps2_core);
fprintf(fid, '%.2f %.2f\n', eps1_clad3,eps2_clad3);
fprintf(fid, '%.3f\n', d);
fprintf(fid, '%d %d\n', npw,nalpha);
fprintf(fid, '%d %d\n', jparxy,jparkz);
fprintf(fid, '%d\n', nk);
fprintf(fid, '%d %d %d\n', jeigenv,jlossmin,jlossband);
fprintf(fid, '%d %.1e %d\n', nom,error,niter);
fprintf(fid, '%.2f %d\n', factor,jshift);
fprintf(fid, '%d\n', jfreephotons);
fprintf(fid, '%d %.2f %.2f %d\n', jloop,varmin,varmax,nvar);
fprintf(fid, '%d %d %d\n', jfields,jkeigenv,jstate);
fprintf(fid, '%.1f %.1f %d %.1f %.1f %d %.1f %.1f %d\n',xmin,xmax,nxgrid,ymin,ymax,nygrid,zmin,zmax,nzgrid);
fprintf(fid, '%d\n', path);
else
fprintf(fid, '%d\n', jlattice);
fprintf(fid, '%d %.2f %.2f %.2f %.2f\n',jbasis,alength1,alength2,alength3,alength4);
fprintf(fid, '%.2f %.2f\n', eps1,eps2);
fprintf(fid, '%d %d\n', jparxy,jparkz);
fprintf(fid, '%.2f %.2f %d\n', akzmin,akzmax,nkz);
fprintf(fid, '%d\n', nk);
fprintf(fid, '%d\n', npw);
fprintf(fid, '%d %d %d\n', jeigenv,jkeigenv,jstate);
fprintf(fid, '%.2f\n', factor);
fprintf(fid, '%.1f %.1f %.1f %.1f %d\n',xmin,xmax,ymin,ymax,ngrid);
fprintf(fid, '%d\n', jfreephotons);
fprintf(fid, '%d %.2f %.2f %d\n', jloop,varmin,varmax,nvar);
fprintf(fid, '%d\n', path);
end
fclose(fid)
end
end