/* Solve arbitrary potential problem in 1-D    */
/* by inverting scattering matrix h[][]        */
/* of quantum boundary conditions              */
/* scattering solutions only.                  */
/* The number of lattice points in potential   */
/* region is (MAX-2)/2, the array determines   */
/* the coefficients R,T, A_i, B_i              */
/* for i=1,2,...,(MAX-2)/2                     */
/* dx should be chosen accordingly             */
/* the lattice points are at matching boundary */
/* points                                      */
#include<stdlib.h>
#include<stdio.h>
#include<complex.h>
#include<math.h>
#define PI 3.1415926
#define MAX 100   /* maximum array size */
#define N 99    /* 1 less than MAX; passed to inversion algorithm */

_Complex double h[MAX][MAX], hinv[MAX][MAX];
_Complex double R,T,A,B;
double dx,x,E,dE;

double V(double x); /* potential function used to get k[] */
double k[MAX/2-1]; /* array to hold momenta */
double X[MAX/2-1]; /* array to hold boundary match points */

int i,j,l,m;
char c;
void cinverse(_Complex double (*A)[N+1],_Complex double (*Ainv)[N+1]);

double modulus(_Complex  double c);


main(){

_Complex double (*h_ptr)[MAX], (*h_invptr)[N+1];
 
 
 dx=0.1; /* axis increments */
 dE=0.05;
 for(l=1;l<=500;l++){
E=(double)l*dE;
 /* set the incoming wave energy */

k[0]=sqrt(2.0*E);
 for(m=1;m<=MAX/2-1;m++){  /* creation of k-vectors */
   k[m]=csqrt(2.0*E+2.0*V((double)m*dx-dx/2.0)) ;
   X[m]=(double)m*dx;}     /* set up array of matching points */
X[0]=0.0;


/* create the matrix to be inverted */

  for (i = 0; i < MAX; i++)
    {
      for(j = 0;j < MAX; j++){
          hinv[i][j]=0.0;
          h[i][j]=0.0;
      }
    }

  for(m=0;m<MAX-1;m++){  /* create diagonal matrix entries */
x=(double)m*dx;
if((m/2)*2==m)
h[m][m]=cexp(-I*k[m/2]*X[m/2]);
else
  h[m][m]=-k[(m+1)/2]*cexp(I*k[(m+1)/2]*X[m/2]);}
  /* fix the first two */
h[0][0]=-1.0;
h[1][1]=k[1];
  for(m=0;m<MAX-1;m++){  /* create off diagonal upper matrix entries */
x=(double)m*dx;
if((m/2)*2==m)
h[m][m+1]=-cexp(I*k[m/2+1]*X[m/2]);
else
  h[m][m+1]=k[(m+1)/2]*cexp(-I*k[(m+1)/2]*X[(m-1)/2]);}
  /* fix first two entries */
h[0][1]=1.0;
h[1][2]=-k[1];
  for(m=1;m<MAX;m++){  /* create off diagonal lower matrix entries */
x=(double)m*dx;
if((m/2)*2==m)
h[m][m-1]=cexp(I*k[m/2]*X[m/2]);
else
  h[m][m-1]=-k[(m-1)/2]*cexp(-I*k[(m-1)/2]*X[(m-1)/2]);}
  /* fix first two entries */
h[1][0]=k[0];
h[1][2]=-k[1];
  for(m=1;m<MAX-2;m++){  /* create far off diagonal upper matrix entries */
if((m/2)*2==m)
  h[m][m+2]=conj(h[m][m+1]);}
  for(m=3;m<MAX;m++){  /* create far off diagonal lower matrix entries */
if((m/2)*2!=m)
  h[m][m-2]=-conj(h[m][m-1]);}
  /* set first one */
h[0][2]=1.0;
/* fix the T portion */
h[MAX-1][MAX-1]=-k[0]*cexp(I*k[0]*X[MAX/2-1]);
h[MAX-2][MAX-1]=-cexp(I*k[0]*X[MAX/2-1]);



/* initialize the first addresses for pointers */ 
h_ptr=&h[0];
h_invptr=&hinv[0];

cinverse(h_ptr,h_invptr);  /* invert the matrix */

 /* compute the transmission coefficient */
R=hinv[0][0]+k[0]*hinv[0][1];
T=hinv[MAX-1][0]+k[0]*hinv[MAX-1][1];

printf("%g\t%g\n",E,modulus(T));


 }/* end loop over incident energies*/

 return;      /* end program */

}



void cinverse(_Complex  double (*A)[N+1], _Complex  double (*Ainv)[N+1]){

  /* Program to compute the inverse of a complex matrix     */
  /* (*A)[] is a pointer to matrix to be inverted           */
  /* (*Ainv)[] is a pointer to array that holds the inverse */
  /* program will hang if A[0][0]=0.0                       */
  /* this clobbers the original matrix                      */

int i,j,k,l,m,col,row;
_Complex  double factor, sum;
_Complex  double Aux[N+1][N+1],one[N+1][N+1],Aux2[N+1][N+1];

/* all matrices are (N+1)X(N+1) */
/* a[][] is to be inverted, one[][] is an auxiliary unit matrix */ 


for(i=0;i<=N;i++){
for(j=0;j<=N;j++){
Aux[i][j]=0.0;Aux2[i][j]=0.0;
if(i==j)one[i][j]=1.0;
else one[i][j]=0.0;}}



/* First we perform the row reductions 
reducing the lower triangle to zeros */
for(k=0;k<=N;k++){
int m=k;
int u=0; 
do{u+1;}while(*(*(A+k)+k+u)==0.0);
m=k+u;
for(l=m+1;l<=N;l++){
if(*(*(A+k)+m)!=0.0){
factor=*(*(A+k)+l) / *(*(A+k)+m);
for(j=0;j<=N;j++){
*(*(A+j)+l)=*(*(A+j)+l)-factor * *(*(A+j)+m);
one[j][l]=one[j][l]-factor*one[j][m];}}}
for(j=0;j<N;j++){
Aux[j][k]=*(*(A+j)+k);
*(*(A+j)+k)=*(*(A+j)+m);
*(*(A+j)+m)=Aux[j][k];
/* apply to unit matrix */
Aux2[j][k]=one[j][k];
one[j][k]=one[j][m];
one[j][m]=Aux2[j][k];}
}

/* Now reduce the upper triangle 
to zeros by row reduction */

for(k=N;k>=0;k--){
for(l=k-1;l>=0;l--){
factor=*(*(A+k)+l)/ *(*(A+k)+k);
for( j=0;j<=N;j++){
*(*(A+j)+l)=*(*(A+j)+l)-factor * *(*(A+j)+k);
one[j][l]=one[j][l]-factor*one[j][k];}}}

/* Now put together the inverse */

for(col=0;col<=N;col++){
for(row=0;row<=N;row++){
one[col][row]=one[col][row]/ *(*(A+row)+row);
*(*(Ainv+col)+row)=one[col][row];}}

return ;
}

double V(double x)
{
    if(x<1.0 || x>3.0)
	return(0.0);
    else
	return(4.0);
}



double modulus(_Complex  double c)
{
double re,im;
 re=creal(c);
 im=cimag(c);
return(re*re+im*im);
}