program kotrafo;
uses dos,crt;
const grad=3;
var c: array[1..grad] of real;
    pxyz: array[1..100,1..grad] of real;
    de: array[1..grad] of extended;
    b0,b,bh,bp,vec:array[1..grad,1..grad] of extended;
    stelle:array [1..grad] of integer;
    abc:word;
    vektor:array[1..grad] of extended;
    anz,anz1,danz,i,itest,it2,panz:integer;
    f1:array[-1..1,-1..1,-1..1] of word;
    cf,f2:array[-3..3,-3..3,-3..3] of word;
    xh,yh,zh,r:real;
    x,y,z:array [1..100] of real;

procedure eingabe; {++++++++++++++++++++++++++++++++++++++++++++++++++}
begin
writeln('Das Programm erzeugt ein reziprokes Gitter und die dazugeh”rige');
writeln('erste Brillouinzone.');
writeln('Geben Sie die Koordinaten der Gittergrundvektoren der primitiven');
writeln('Elementarzelle ein!');
writeln('Koordinaten des Gittervektors a : ');
write('ax : ');readln(b0[1,1]);
write('ay : ');readln(b0[2,1]);
write('az : ');readln(b0[3,1]);
writeln;
writeln('Koordinaten des Gittervektors b : ');
write('bx : ');readln(b0[1,2]);
write('by : ');readln(b0[2,2]);
write('bz : ');readln(b0[3,2]);
writeln;
writeln('Koordinaten des Gittervektors c : ');
write('cx : ');readln(b0[1,3]);
write('cy : ');readln(b0[2,3]);
write('cz : ');readln(b0[3,3]);
writeln;
end;

procedure reziprok;  {+++++++++++++++++++++++++++++++++++++++++++++++++}
var sp:extended;
    i1,i2:integer;
    strx,stry,strz:string;
begin
sp:=b0[1,1]*(b0[2,2]*b0[3,3]-b0[3,2]*b0[2,3]);
sp:=sp+b0[2,1]*(b0[3,2]*b0[1,3]-b0[1,2]*b0[3,3]);
sp:=sp+b0[3,1]*(b0[1,2]*b0[2,3]-b0[2,2]*b0[1,3]);
b[1,1]:=(b0[2,2]*b0[3,3]-b0[3,2]*b0[2,3])/sp;
b[2,1]:=(b0[3,2]*b0[1,3]-b0[1,2]*b0[3,3])/sp;
b[3,1]:=(b0[1,2]*b0[2,3]-b0[2,2]*b0[1,3])/sp;
b[1,2]:=(b0[2,3]*b0[3,1]-b0[3,3]*b0[2,1])/sp;
b[2,2]:=(b0[3,3]*b0[1,1]-b0[1,3]*b0[3,1])/sp;
b[3,2]:=(b0[1,3]*b0[2,1]-b0[2,3]*b0[1,1])/sp;
b[1,3]:=(b0[2,1]*b0[3,2]-b0[3,1]*b0[2,2])/sp;
b[2,3]:=(b0[3,1]*b0[1,2]-b0[1,1]*b0[3,2])/sp;
b[3,3]:=(b0[1,1]*b0[2,2]-b0[2,1]*b0[1,2])/sp;
writeln('Es ergeben sich folgende reziproke Gittervektoren :');
writeln;
str(b[1,1]:7:4,strx);
str(b[2,1]:7:4,stry);
str(b[3,1]:7:4,strz);
writeln(' a" : (',strx,',',stry,',',strz,')');
writeln;
str(b[1,2]:7:4,strx);
str(b[2,2]:7:4,stry);
str(b[3,2]:7:4,strz);
writeln(' b" : (',strx,',',stry,',',strz,')');
writeln;
str(b[1,3]:7:4,strx);
str(b[2,3]:7:4,stry);
str(b[3,3]:7:4,strz);
writeln(' c" : (',strx,',',stry,',',strz,')');
writeln;
writeln('In den reziproken Gittervektoren ist der Faktor 2*pi nicht enthalten.');
end;

procedure makevektor;    {+++++++++++++++++++++++++++++++++++++++++++++++++}
var i1:integer;
begin
for i1:=1 to grad do begin
vektor[i1]:=b[i1,abc];
end;
end;

procedure makesum;    {++++++++++++++++++++++++++++++++++++++++++++++++++++}
var i1:integer;
begin
for i1:=1 to grad do begin
vec[i1,3]:=vec[i1,1]+vec[i1,2];
end;
end;

function sqrvec(xx:word):real; {+++++++++++++++++++++++++++++++++++++++++++}
var i1:integer;
    rr:real;
begin
rr:=0;
for i1:=1 to grad do begin
rr:=rr+sqr(vec[i1,xx]);
end;
sqrvec:=rr;
end;


procedure reduktion;   {++++++++++++++++++++++++++++++++++++++++++++++++++}
var i1:integer;
    strx,stry,strz:string;
begin
anz1:=0;
repeat
anz:=0;

repeat
danz:=0;
abc:=1;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=2;makevektor;
for i1:=1 to grad do vec[i1,2]:=vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,1]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

repeat
danz:=0;
abc:=1;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=2;makevektor;
for i1:=1 to grad do vec[i1,2]:=-vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,1]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

repeat
danz:=0;
abc:=1;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=3;makevektor;
for i1:=1 to grad do vec[i1,2]:=vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,1]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

repeat
danz:=0;
abc:=1;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=3;makevektor;
for i1:=1 to grad do vec[i1,2]:=-vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,1]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

{.......................................................}

repeat
danz:=0;
abc:=2;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=1;makevektor;
for i1:=1 to grad do vec[i1,2]:=vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,2]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

repeat
danz:=0;
abc:=2;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=1;makevektor;
for i1:=1 to grad do vec[i1,2]:=-vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,2]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

repeat
danz:=0;
abc:=2;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=3;makevektor;
for i1:=1 to grad do vec[i1,2]:=vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,2]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

repeat
danz:=0;
abc:=2;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=3;makevektor;
for i1:=1 to grad do vec[i1,2]:=-vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,2]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

{.......................................................}

repeat
danz:=0;
abc:=3;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=2;makevektor;
for i1:=1 to grad do vec[i1,2]:=vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,3]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

repeat
danz:=0;
abc:=3;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=2;makevektor;
for i1:=1 to grad do vec[i1,2]:=-vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,3]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

repeat
danz:=0;
abc:=3;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=1;makevektor;
for i1:=1 to grad do vec[i1,2]:=vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,3]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

repeat
danz:=0;
abc:=3;makevektor;
for i1:=1 to grad do vec[i1,1]:=vektor[i1];
abc:=1;makevektor;
for i1:=1 to grad do vec[i1,2]:=-vektor[i1];
makesum;
if sqrvec(3)<sqrvec(1) then begin
for i1:=1 to grad do b[i1,3]:=vec[i1,3];inc(anz);inc(danz);end;
until danz=0;

{.......................................................}

anz1:=anz1+anz;
until anz=0;

                           
if anz1>0 then begin
writeln('Es ergeben sich folgende reduzierte reziproke Gittervektoren :');
writeln;
str(b[1,1]:7:4,strx);
str(b[2,1]:7:4,stry);
str(b[3,1]:7:4,strz);
writeln(' a" : (',strx,',',stry,',',strz,')');
writeln;
str(b[1,2]:7:4,strx);
str(b[2,2]:7:4,stry);
str(b[3,2]:7:4,strz);
writeln(' b" : (',strx,',',stry,',',strz,')');
writeln;
str(b[1,3]:7:4,strx);
str(b[2,3]:7:4,stry);
str(b[3,3]:7:4,strz);
writeln(' c" : (',strx,',',stry,',',strz,')');
writeln;
end;
end;

function test(j1,j2,j3,j4,j5,j6:integer):integer;{+++++++++++++++++++++++++}
var rr:real;
    cc:integer;
begin
rr:=(j1*b[1,1]+j2*b[1,2]+j3*b[1,3])*(j4*b[1,1]+j5*b[1,2]+j6*b[1,3])+
    (j1*b[2,1]+j2*b[2,2]+j3*b[2,3])*(j4*b[2,1]+j5*b[2,2]+j6*b[2,3])+
    (j1*b[3,1]+j2*b[3,2]+j3*b[3,3])*(j4*b[3,1]+j5*b[3,2]+j6*b[3,3]);
rr:=rr/(sqr(j4*b[1,1]+j5*b[1,2]+j6*b[1,3])+
    sqr(j4*b[2,1]+j5*b[2,2]+j6*b[2,3])+sqr(j4*b[3,1]+j5*b[3,2]+j6*b[3,3]));

if rr<1 then cc:=0 else cc:=1;
test:=cc;
end;


procedure f1set;   {++++++++++++++++++++++++++++++++++++++++++++++++++++++}
var i1,i2,i3,i4,i5,i6,testwert:integer;
begin
for i1:=-1 to 1 do begin
for i2:=-1 to 1 do begin
for i3:=-1 to 1 do begin
f1[i1,i2,i3]:=0;end;end;end;

for i1:=-1 to 1 do begin
for i2:=-1 to 1 do begin
for i3:=-1 to 1 do begin
if ((i1<>0)or(i2<>0)or(i3<>0)) then begin
testwert:=0;
for i4:=-1 to 1 do begin
for i5:=-1 to 1 do begin
for i6:=-1 to 1 do begin
if (((i4<>0)or(i5<>0)or(i6<>0))and((i1<>i4)or(i2<>i5)or(i3<>i6))) then begin
testwert:=testwert+test(i1,i2,i3,i4,i5,i6);end;end;end;end;
if testwert=0 then begin f1[i1,i2,i3]:=1;inc(itest);end;end;end;end;end;
end;

procedure f2set;     {+++++++++++++++++++++++++++++++++++++++++++++++++++}
var i1,i2,i3,i4,i5,i6,testwert:integer;
begin
for i1:=-3 to 3 do begin
for i2:=-3 to 3 do begin
for i3:=-3 to 3 do begin
f2[i1,i2,i3]:=0;
if((abs(i1)<>abs(i2))and(abs(i1)<>abs(i3))and(abs(i2)<>abs(i3))
and(i1<>0)and(i2<>0)and(i3<>0)) then begin f2[i1,i2,i3]:=1;inc(it2);end;
end;end;end;
end;

procedure matrixtrafo(j1:integer); {++++++++++++++++++++++++++++++++}
var bb:array[1..grad] of real;
    dede:real;
    j2,jmax:integer;
begin
if (j1<grad) then begin
jmax:=j1;
for j2:=j1+1 to grad do begin
if (abs(b[jmax,j1])<abs(b[j2,j1])) then jmax:=j2;end;
for j2:=1 to grad do begin
bb[j2]:=b[j1,j2];b[j1,j2]:=b[jmax,j2];b[jmax,j2]:=bb[j2];end;
dede:=de[j1];de[j1]:=de[jmax];de[jmax]:=dede;
end;
end;

procedure koeff;{+++++++++++++++++++++++++++++++++++++++++++++}
var detb:real;
    faktor:extended;
    ka,ki,i1,i2:integer;
begin
for i1:=1 to grad do begin
matrixtrafo(i1);
i2:=0;
repeat
inc(i2);
stelle[i1]:=i2;
until b[i1,i2]<>0;
 if i1<grad then begin
  for ka:=i1+1 to grad do begin
  faktor:=b[ka,i2]/b[i1,i2];
   for ki:=i2 to grad do begin
   b[ka,ki]:=b[ka,ki]-faktor*b[i1,ki];
   end;
  b[ka,i2]:=0;
  de[ka]:=de[ka]-faktor*de[i1];
  end;
 end;
end;
end;

procedure koeff2;   {+++++++++++++++++++++++++++++++++++++++++++++}
var i1,i2:integer;
begin
for i1:=grad downto 2 do begin
  c[stelle[i1]]:=de[i1]/b[i1,stelle[i1]];
    for i2:=i1-1 downto 1 do begin
    de[i2]:=de[i2]-c[stelle[i1]]*b[i2,stelle[i1]];
    end;
end;
c[stelle[1]]:=de[1]/b[1,stelle[1]];
end;

procedure feldtausch;  {+++++++++++++++++++++++++++++++++++++++++++++++++}
var i1,i2:integer;
begin
for i1:=1 to 3 do begin
for i2:=1 to 3 do begin
bh[i1,i2]:=b[i1,i2];
end;end;
end;

procedure vektor2(j:integer);  {++++++++++++++++++++++++++++++++++++++++++}
begin
xh:=0;yh:=0;zh:=0;
if j=-3 then begin xh:=-bh[1,3];yh:=-bh[2,3];zh:=-bh[3,3];end;
if j=-2 then begin xh:=-bh[1,2];yh:=-bh[2,2];zh:=-bh[3,2];end;
if j=-1 then begin xh:=-bh[1,1];yh:=-bh[2,1];zh:=-bh[3,1];end;
if j=1 then begin xh:=bh[1,1];yh:=bh[2,1];zh:=bh[3,1];end;
if j=2 then begin xh:=bh[1,2];yh:=bh[2,2];zh:=bh[3,2];end;
if j=3 then begin xh:=bh[1,3];yh:=bh[2,3];zh:=bh[3,3];end;
end;

function ecke(j1,j2,j3:integer):integer;{++++++++++++++++++++++++++++}
var j4,j5,j6:integer;
begin
if abs(j1)=1 then j4:=j1;
if abs(j2)=1 then j4:=j2;
if abs(j3)=1 then j4:=j3;
if abs(j1)=2 then j5:=j1;
if abs(j2)=2 then j5:=j2;
if abs(j3)=2 then j5:=j3;
if abs(j1)=3 then j6:=j1;
if abs(j2)=3 then j6:=j2;
if abs(j3)=3 then j6:=j3;
j5:=round(j5/2);j6:=round(j6/3);
ecke:=f1[j4,j5,j6];
end;

function kante(j1,j2:integer):integer;{++++++++++++++++++++++++++++++++++}
var j4,j5,j6:integer;
begin
j4:=0;j5:=0;j6:=0;
if abs(j1)=1 then j4:=j1;
if abs(j2)=1 then j4:=j2;
if abs(j1)=2 then j5:=j1;
if abs(j2)=2 then j5:=j2;
if abs(j1)=3 then j6:=j1;
if abs(j2)=3 then j6:=j2;
j5:=round(j5/2);j6:=round(j6/3);
kante:=f1[j4,j5,j6];
end;

procedure makeb1punkte;  {++++++++++++++++++++++++++++++++++++++++++++++++}
var i1,i2,i3,ih:integer;
begin
for i1:=-3 to 3 do begin
for i2:=-3 to 3 do begin
for i3:=-3 to 3 do begin
cf[i1,i2,i3]:=0;end;end;end;
for i1:=-3 to 3 do begin
for i2:=-3 to 3 do begin
for i3:=-3 to 3 do begin
if f2[i1,i2,i3]<>0 then begin
vektor2(i1);
bp[1,1]:=xh;bp[1,2]:=yh;bp[1,3]:=zh;
vektor2(i2);
bp[2,1]:=xh;bp[2,2]:=yh;bp[2,3]:=zh;
vektor2(i3);
bp[3,1]:=xh;bp[3,2]:=yh;bp[3,3]:=zh;

{if ((i1=1)and(i2=-2)and(i3=3)) then begin
readln;end;}
if ecke(i1,i2,i3)=1 then begin
for ih:=1 to 3 do b[1,ih]:=bp[1,ih];
if kante(i1,i2)=1 then begin
for ih:=1 to 3 do b[2,ih]:=bp[2,ih]+bp[1,ih];end else begin
for ih:=1 to 3 do b[2,ih]:=bp[2,ih];end;
for ih:=1 to 3 do b[3,ih]:=bp[3,ih]+bp[2,ih]+bp[1,ih];
for ih:=1 to 3 do de[ih]:=sqr(b[ih,1])+sqr(b[ih,2])+sqr(b[ih,3]);
koeff;
koeff2;
inc(panz);
for ih:=1 to 3 do pxyz[panz,ih]:=c[ih];
end else begin
for ih:=1 to 3 do b[1,ih]:=bp[1,ih];
if kante(i1,i2)=1 then begin
for ih:=1 to 3 do b[2,ih]:=bp[2,ih]+bp[1,ih];end else begin
for ih:=1 to 3 do b[2,ih]:=bp[2,ih];end;
for ih:=1 to 3 do b[3,ih]:=bp[3,ih];
for ih:=1 to 3 do de[ih]:=sqr(b[ih,1])+sqr(b[ih,2])+sqr(b[ih,3]);
koeff;
koeff2;
inc(panz);
for ih:=1 to 3 do pxyz[panz,ih]:=c[ih];

if kante(i2,i3)=1 then begin
for ih:=1 to 3 do b[1,ih]:=bp[1,ih];
if kante(i1,i2)=1 then begin
for ih:=1 to 3 do b[2,ih]:=bp[2,ih]+bp[1,ih];end else begin
for ih:=1 to 3 do b[2,ih]:=bp[2,ih];end;
for ih:=1 to 3 do b[3,ih]:=bp[3,ih]+bp[2,ih];
for ih:=1 to 3 do de[ih]:=sqr(b[ih,1])+sqr(b[ih,2])+sqr(b[ih,3]);
koeff;
koeff2;
if (sqr(c[1])+sqr(c[2])+sqr(c[3])<
sqr(pxyz[panz,1])+sqr(pxyz[panz,2])+sqr(pxyz[panz,3])) then begin
for ih:=1 to 3 do pxyz[panz,ih]:=c[ih];end;end;

if kante(i1,i3)=1 then begin
for ih:=1 to 3 do b[1,ih]:=bp[1,ih];
if kante(i1,i2)=1 then begin
for ih:=1 to 3 do b[2,ih]:=bp[2,ih]+bp[1,ih];end else begin
for ih:=1 to 3 do b[2,ih]:=bp[2,ih];end;
for ih:=1 to 3 do b[3,ih]:=bp[3,ih]+bp[1,ih];
for ih:=1 to 3 do de[ih]:=sqr(b[ih,1])+sqr(b[ih,2])+sqr(b[ih,3]);
koeff;
koeff2;
if (sqr(c[1])+sqr(c[2])+sqr(c[3])<
sqr(pxyz[panz,1])+sqr(pxyz[panz,2])+sqr(pxyz[panz,3])) then begin
for ih:=1 to 3 do pxyz[panz,ih]:=c[ih];end;end;
end;
cf[i1,i2,i3]:=panz;
end;end;end;end;
end;

procedure select;   {+++++++++++++++++++++++++++++++++++++++++++++++++}
var i1,i2,i3:integer;
    xx,yy,zz:real;
begin
anz:=0;
for i1:=1 to panz do begin
xx:=pxyz[i1,1];yy:=pxyz[i1,2];zz:=pxyz[i1,3];
if anz=0 then begin
inc(anz);
x[anz]:=xx;y[anz]:=yy;z[anz]:=zz;end else begin
i3:=0;
for i2:=1 to anz do begin
if (sqr(xx-x[i2])+sqr(yy-y[i2])+sqr(zz-z[i2])<0.001) then inc(i3);
end;
if i3=0 then begin
inc(anz);
x[anz]:=xx;y[anz]:=yy;z[anz]:=zz;end;
end;end;
end;

procedure ausgabe;
var tex:text;
    name,bbz,str1,str2,str3:string;
    rbz,fbz:real;
    i1,i2,i3,i4:integer;
begin {
for i1:=-3 to 3 do begin
for i2:=-3 to 3 do begin
for i3:=-3 to 3 do begin
if cf[i1,i2,i3]<>0 then begin
i4:=cf[i1,i2,i3];
str(pxyz[i4,1]:7:4,str1);
str(pxyz[i4,2]:7:4,str2);
str(pxyz[i4,3]:7:4,str3);
writeln(i1,' ',i2,' ',i3,'    ',str1,'  ',str2,'  ',str3);
end;end;end;end;
readln;}
write('Radius der B-Zonen-Ecken : ');readln(rbz);
write('B-Zone um Faktor pi gestaucht (j/n) ');readln(bbz);
if bbz='n' then fbz:=pi else fbz:=1;
write('Dateiname zur Speicherung : ');readln(name);
assign(tex,name);rewrite(tex);
for i:=1 to anz do begin
writeln(tex,7);writeln(tex,rbz);
writeln(tex,x[i]*fbz);writeln(tex,y[i]*fbz);writeln(tex,z[i]*fbz);
end;
close(tex);
end;



begin    {++++++++++ Hauptprogramm +++++++++++++++++++++++++++++++++++}
clrscr;
itest:=0;it2:=0;panz:=0;
eingabe;
reziprok;
reduktion;
f1set;
f2set;
feldtausch;
makeb1punkte;
select;
writeln('Die erste Brillouinzone ist ein K”rper mit ',anz,' Ecken und ',itest,' Fl„chen.');
writeln;
{readln;}
ausgabe;
writeln;
writeln('Programmende.');
readln;
end.