/* code integrates 1d Schroedinger eqn. using a shooting method. */
/* adapted to cope with NUM identical potential wells */
/* general well width W, depth V, well wall thickness A*/
/* setting NUM to 1 and A to zero gives single well of width 2*W*/
/* setting NMUM to 1 and A non-zero gives chemical bond potential */
/* setting NUM large yields 1d crystal - see band structure */
/* potential is set to zero outside the well - until it encounters a hardwall *//* at LARGE_X */
 
#include <math.h>
#include <stdio.h>

/* function prototypes */
double find_energy(double ue, double le);
double wavefunction_at_large_x(double e);
void integrate(int i, double e);
double potential(double x);

/* parameters of the potential */

#define V -1.0
#define W 10.0
#define A 0.0
#define NUM 1
#define POT_LEN (NUM*(W+A))
#define LARGE_X (POT_LEN+W)
 
/* largest number of steps allowed */

#define MAX 100000
/* xstep and wavefunction accuracy */

#define DX 0.01
#define RESOL 0.001

/* find all solutions up to MAX_ENERGY assuming the average gap between */
/* sucessive allowed energies is ENERGY_GAP; */

#define MAX_ENERGY 0.0
#define MIN_ENERGY V
#define ENERGY_GAP 0.0005

/* PARITY may be 0 or 1 depending on whether the wavefunction is odd or even */

#define PARITY 0

/* wavefunction */
double psi[MAX],dpsi[MAX];

/* number of steps */
int N;

int main(void){
FILE *fp;
char s[25];
int i;

double energy,norm;
int outerloopcount=0;

/* search intervals for energy and resolution */
double upper_energy=MIN_ENERGY,lower_energy=MIN_ENERGY;

/* sign of wavefunction at large x using upper and lower energies */
double upper_energy_sign,lower_energy_sign;

N=(int)LARGE_X/DX;
if(N>MAX){(void)printf("Increase size of arrays - MAX variable\n");exit(1);}

lower_energy_sign=wavefunction_at_large_x(lower_energy);

/* loop incrementing the upper_energy until an interval is found which */
/* contains a possible allowed energy */

while(upper_energy<MAX_ENERGY){
outerloopcount++;

/*compute wavefunction at large x from upper energy*/
upper_energy_sign=wavefunction_at_large_x(upper_energy);

if(upper_energy_sign*lower_energy_sign<0.0){
(void)printf("Energy interval %lg--%lg contains an allowed energy\n",lower_energy,upper_energy);
(void)printf("Number of iterations required to locate interval %d\n",outerloopcount);
fflush(stdout);

/* have found an interval containing an odd number of allowed*/
/* energies - hopefully just one ! */

energy=find_energy(upper_energy,lower_energy);
lower_energy=upper_energy+ENERGY_GAP;
lower_energy_sign=upper_energy_sign;
(void)printf("Energy is %lg\n",energy);

/* normalize wavefunction */
norm=0.0;
for(i=0;i<N;i++)
norm+=psi[i]*psi[i]*DX;

for(i=0;i<N;i++)
psi[i]/=sqrt(norm);

/* write out wavefunction to file labelled by energy*/
sprintf(s,"E=%lgwave",energy);
fp=fopen(s,"w");
for(i=0;i<N;i++){
fprintf(fp,"%lg %lg\n",i*DX,psi[i]);}

fclose(fp);

fflush(stdout);
outerloopcount=0;
}

upper_energy+=(ENERGY_GAP);
}

(void)printf("Finished\n");
return(0);
}

double find_energy(double ue, double le){
/* uses bisection to locate true energy */
double e,sl,su,s;
int count=0;

sl=wavefunction_at_large_x(le);
su=wavefunction_at_large_x(ue);

do{
count++;

e=(le+ue)*0.5;
s=wavefunction_at_large_x(e);
if((s*sl)>0.0){
sl=s;
le=e;
}
else{
su=s;
ue=e;
}
}while(fabs(s)>RESOL);

(void)printf("Number of iterations to locate energy %d\n",count);
(void)fflush(stdout);
return(e);
}

double wavefunction_at_large_x(double e){
/* integrates out to large x and returns value of psi there */

int i;
psi[0]=(1-PARITY);
dpsi[0]=PARITY;
for(i=0;i<N-1;i++){
integrate(i,e);
}

return(psi[N-1]);
}

void integrate(int i, double e){

/* calculates one update of pis and dpsi using simple Euler scheme */
psi[i+1]=psi[i]+DX*dpsi[i];
dpsi[i+1]=dpsi[i]+DX*(potential(i*DX)-e)*psi[i];

return;
}

double potential(double x){
/* calculates potential at point x */
int idum;
double xdum;

if(x>LARGE_X) return(10.0/RESOL);
if(x>POT_LEN) return(0);

/* first find which well you are in */
idum=(int)(x/(A+W));
/* next find position within this well */
xdum=x-idum*(A+W);

if(xdum<A/2) return(0);
if(xdum>(W+A/2)) return(0);
return(V);

}