#include<iostream>
#include<cmath>
#include<fstream>
double U_K(double K,double a,double u);
using namespace std;
int main()
{
cout<<"Enter set of possible lattice vector's i.e., 'N' & 'U'\t:\t"<<endl;
int N ; double u ;
cin>>N>>u;
double a = 1;
double K[2*N+1];
int j=0;
//assigning value to 'K' i.e., lattice vectors (-N to N)
for(int i = -N; i <=N ; i++)
{
K[j] = (2*3.14*i)/a;
j++;
}
int r = 2*N+1, c = 2*N+1;
double mat[r][c],k=0;
ofstream ofs("data_cmp_U-2_N_3.txt");
//varying wave-vector i.e., "k"
for(double k = (-3.14)/a ; k < (3.14/a) ; k += 0.1)
{
for(int i=0; i < r ; i++)
{
for(int j=0; j < c ; j++)
{
if (i == j) mat[i][j] = (0.5)*pow((k - K[i]),2);
else mat[i][j] = U_K(K[j] - K[i],a,u);
}
}
double ftr;
//digonalizing matrix to find energy eigen-values
for(int i=0;i<r;i++)
{
ftr =0;
for(int j=0 ; j<c ;j++)
{
if(i>j)
{
ftr = mat[i][j]/mat[j][j];
for( int k=0; k<c;k++)
{
mat[i][k] -= ftr*mat[j][k];
}
}
}
}
for(int i=r-1;i>=0;i--)
{
ftr =0;
for(int j=c-1 ; j>=0 ;j--)
{
if(i<j)
{
ftr = mat[i][j]/mat[j][j];
for( int k=0; k<c;k++)
{
mat[i][k] -= ftr*mat[j][k];
}
}
}
}
ofs<<k<<"\t\t";
for(int i =0 ; i< r;i++)
{
for(int j =0 ; j < c ; j++)
{
if(i == j) ofs<<mat[i][j]<<"\t\t";
}
}
ofs<<endl;
}
return 0;
}
double U_K(double K,double a,double u)
{
return u*exp(-K)/a;
}
I2luY2x1ZGU8aW9zdHJlYW0+CiNpbmNsdWRlPGNtYXRoPgojaW5jbHVkZTxmc3RyZWFtPgoKZG91YmxlIFVfSyhkb3VibGUgSyxkb3VibGUgYSxkb3VibGUgdSk7Cgp1c2luZyBuYW1lc3BhY2Ugc3RkOwppbnQgbWFpbigpCnsKCWNvdXQ8PCJFbnRlciBzZXQgb2YgcG9zc2libGUgbGF0dGljZSB2ZWN0b3IncyBpLmUuLCAnTicgJiAnVSdcdDpcdCI8PGVuZGw7CglpbnQgTiA7IGRvdWJsZSB1IDsKCWNpbj4+Tj4+dTsKCWRvdWJsZSBhID0gMTsKCWRvdWJsZSBLWzIqTisxXTsKCQoJaW50IGo9MDsKCS8vYXNzaWduaW5nIHZhbHVlIHRvICdLJyBpLmUuLCBsYXR0aWNlIHZlY3RvcnMgKC1OIHRvIE4pCglmb3IoaW50IGkgPSAtTjsgaSA8PU4gOyBpKyspCgl7CgkJS1tqXSA9ICgyKjMuMTQqaSkvYTsKCQlqKys7Cgl9CgkKCWludCByID0gMipOKzEsIGMgPSAyKk4rMTsKCWRvdWJsZSBtYXRbcl1bY10saz0wOwoJCglvZnN0cmVhbSBvZnMoImRhdGFfY21wX1UtMl9OXzMudHh0Iik7CgkvL3Zhcnlpbmcgd2F2ZS12ZWN0b3IgaS5lLiwgImsiCglmb3IoZG91YmxlIGsgPSAoLTMuMTQpL2EgOyBrIDwgKDMuMTQvYSkgOyBrICs9IDAuMSkKCXsKCQlmb3IoaW50IGk9MDsgaSA8IHIgOyBpKyspCgkJewoJCQlmb3IoaW50IGo9MDsgaiA8IGMgOyBqKyspCgkJCXsKCQkJCWlmIChpID09IGopIG1hdFtpXVtqXSA9ICgwLjUpKnBvdygoayAtIEtbaV0pLDIpOwoJCQkJZWxzZSBtYXRbaV1bal0gPSBVX0soS1tqXSAtIEtbaV0sYSx1KTsKCQkJfQoJCX0KCQkKCQlkb3VibGUgZnRyOwoJCS8vZGlnb25hbGl6aW5nIG1hdHJpeCB0byBmaW5kIGVuZXJneSBlaWdlbi12YWx1ZXMKCQlmb3IoaW50IGk9MDtpPHI7aSsrKQoJCXsKCQkJZnRyID0wOwoJCQlmb3IoaW50IGo9MCA7IGo8YyA7aisrKQoJCQl7CgkJCQlpZihpPmopCgkJCQl7CgkJCQkJZnRyID0gbWF0W2ldW2pdL21hdFtqXVtqXTsKCQkJCQlmb3IoIGludCBrPTA7IGs8YztrKyspCgkJCQkJewoJCQkJCQltYXRbaV1ba10gLT0gZnRyKm1hdFtqXVtrXTsKCQkJCQl9CgkJCQl9CgkJCQkKCQkJfQoJCX0KCQkKCQlmb3IoaW50IGk9ci0xO2k+PTA7aS0tKQoJCXsJCgkJCWZ0ciA9MDsKCQkJZm9yKGludCBqPWMtMSA7IGo+PTAgO2otLSkKCQkJewoJCQkJaWYoaTxqKSAKCQkJCXsKCQkJCQlmdHIgPSBtYXRbaV1bal0vbWF0W2pdW2pdOwoJCQkJCWZvciggaW50IGs9MDsgazxjO2srKykKCQkJCQl7CgkJCQkJCW1hdFtpXVtrXSAtPSBmdHIqbWF0W2pdW2tdOwoJCQkJCX0KCQkJCX0KCQkJfQoKCQl9CgkJb2ZzPDxrPDwiXHRcdCI7CgkJZm9yKGludCBpID0wIDsgaTwgcjtpKyspCgkJewoJCQlmb3IoaW50IGogPTAgOyBqIDwgYyA7IGorKykKCQkJewkKCQkJCWlmKGkgPT0gaikgb2ZzPDxtYXRbaV1bal08PCJcdFx0IjsKCQkJfQkKCQl9CgkJb2ZzPDxlbmRsOwkJCgl9CgkKCQoJcmV0dXJuIDA7Cn0KCmRvdWJsZSBVX0soZG91YmxlIEssZG91YmxlIGEsZG91YmxlIHUpCnsKCXJldHVybiB1KmV4cCgtSykvYTsKfQ==