program bands
	integer b
	real e(125)
	real vec(3)
	real hamilr(125,125),hamili(125,125)
	real m
	real zr(125,125),zi(125,125),fv1(125,125),fv2(125,125),fm1(125,125)

	m=0.0

	open(24,file='gaas.dat',status='new')

	do 1 b=0,40
		call kgenerator(m,vec)
		write(*,*) vec
		call hamiltonian(vec,hamilr,hamili)
		call ch(125,125,hamilr,hamili,e,1,zr,zi,fv1,fv2,fm1,ier)
		write(*,*) ier
		e=e-8.8010
		write(24,20) m,(e(l),l=1,15)
		m=m+1.0
      1 continue
     20 format(16f10.4)
	end

	subroutine hamiltonian(vec,hamilr,hamili)
	real vec(3)
	complex hamil(125,125)
	real hamilr(125,125),hamili(125,125)
	integer i,j
	complex element

	do 2 i=-62,62
		do 3 j=-62,62
			call el(i,j,vec,element)
	      		hamil(i+63,j+63)=element
		3 continue
      2 continue
	
	hamilr=real(hamil)
	hamili=aimag(hamil)
	
	return
	end

! kgenerator gives the starting kvector of the different symmetry points to be calculated
	subroutine kgenerator(m,vec)
	real vec(3)
	real m
	real a

	if (m.LE.10) then
		vec=[0.5-m/20,0.5-m/20,0.5-m/20]
	elseif ((m.GT.10).AND.(m.LE.20)) then
   		vec=[(m-10)/10,0.,0.]
	elseif ((m.GT.20).AND.(m.LE.25)) then			
   		vec=[1.,(m-20)/10,0.]
	elseif ((m.GT.25).AND.(m.LE.30)) then		
   		vec=[1-(m-25)/20,0.5+(m-25)/20,0.]
	else				
   		vec=[.75-(m-30)/(40/3),.75-(m-30)/(40/3),0.]
	endif

	return
	end

! el gives the hamiltonian matrix element given an index i,j and the symmetry vec
	subroutine el(i,j,vec,element)
	integer i,j
	real vec(3),out(3)
	complex vout
	complex element
	real pi,hbar,a,mo

	pi=3.141592654
	hbar=(6.626e-34)/(2*pi)
	a=5.64e-10
	mo=9.109e-31

	if (i.Eq.j) then
		call kval(i,out)
		element=(hbar*hbar)/(2*mo)*dot_product(out+vec,out+vec)*(4*pi*pi/(a*a))/1.602e-19
		write(*,*) hbar
		write(*,*) element
		
	else
		call vg(i-j,vout)
		element=vout
		
	endif

	return
	end
	
! vg gives the pseudopotential for a given RPL given by its index m
	subroutine vg(m,vout)
	real Vs,Va
	complex vout
	integer m,a
	real pi
	real out(3)
	real s(3)

	pi=3.141592654

	call kval(m,out)
	a=dot_product(out,out)

	s=[0.125,0.125,0.125]

	if (a.EQ.3) then
   		Vs=-0.23
		Va=0.07
	elseif (a.EQ.4) then
		Vs=0.0
		Va=0.05
	elseif (a.EQ.8) then
		Vs=0.01
		Va=0.0
	elseif (a.EQ.11) then
		Vs=0.06
		Va=0.01
	else
		Vs=0.0
		Va=0.0
	endif

	vout=cmplx(13.6*Vs*cos(2*pi*dot_product(out,s)),-13.6*Va*sin(2*pi*dot_product(out,s)))

	return
	end

! kval gives unique RPL's based on an index
! kval generates RPL's in such a way that kval(i)-kval(j)=kval(i-j), which is handy
! The idea for this subroutine was taken from http://vcsel.micro.uiuc.edu/~adanner/pseudopotential.htm
	subroutine kval(m,out)
	integer m
	real lim,n
	real vecb1(3),vecb2(3),vecb3(3)
	real hvec,kvec,lvec
	real out(3)

	lim=5
	n=(lim*lim*lim)-1

	vecb1(1)=-1
	vecb1(2)=1
	vecb1(3)=1

	vecb2(1)=1
	vecb2(2)=-1
	vecb2(3)=1

	vecb3(1)=1
	vecb3(2)=1
	vecb3(3)=-1
	

	hvec=floor((m+n/2)/(lim*lim))-floor(lim/2)
	kvec=floor(mod(m+n/2,lim*lim)/lim)-floor(lim/2)
	lvec=mod(m+n/2,lim)-floor(lim/2)
	
	out=hvec*vecb1+kvec*vecb2+lvec*vecb3



	return
	end