%brillion zones

m=1
c=1
a=1

x=input("Enter the brillion zone number ")

brillion=fopen('zone.dat','a+')

for k=(x*pi)/a:(2*pi)/500:(x*pi)/a;
  theta=k*(a/2);
  w=2*sqrt(c/m)*abs(sin(theta));	
  fprintf(brillion,"%f\t%f\n",theta,w)
endfor

fclose (brillion);
