/* code integrates 1d Schroedinger eqn. using a shooting method. */
/* adapted to double well */

#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);


/* xstep and wavefunction accuracy */

#define dx 0.001
#define RESOL 0.01

/* parameters of potential */

#define V0 -10.0
#define A 1.0
#define B 0.0

/* 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 V0
#define ENERGY_GAP 1.0

/* point where wavefunction should be small */
#define XLONG (A+B)

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

/* wavefunction */
double psi[100000],dpsi[100000];
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;

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

N=(int)(XLONG/dx);
if(N>=100000){printf("too many steps for arrays\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;
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(2*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*0.1);
}

(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];
}

double potential(double x){
/* calculates potential at point x */

/* harmonic oscillator */
/*return(0.25*x*x); */

if(x<=B) return(0);
if((x>B) && (x<(A+B))) return(V0);
if(x>=(A+B)) return(10000.0);

}