zsample.c

/* file: zsample.c						*/

#include <stdio.h>
#include "zeros.c"

/*  Sample program. 

    Uses DFZERO + zeros.c  interval zero finder. 
    
    Solve zeros for a boundary value problem:

    The sample program calculates the energies and states for a particle in the
    following potential: output for the states is in Mathematica format.
    
        V(x)
    
     ~~~^		^~~~~~  inf
	|		|
	|		|
	|	 _______|       Vo
	|	 
	|________
    		 a       --->x


    						output file: zsample.out    
*/

FILE *fp;

double VV=0.;
double MM, AA,HH;



double gke(double *e) {		     /* function whose zeros we to find	*/
 double k,E,ka, kk_sin_ka, cos_ka;	
 E=*e;
 k=sqrt(2*MM*E)/HH;
 ka=sqrt(2*MM*(VV-E))/HH;
 kk_sin_ka=(ka/k)*sin(k*AA);
 cos_ka=cos(k*AA);
 E=exp(2*ka*AA)*(kk_sin_ka+cos_ka)+(kk_sin_ka-cos_ka);
 return(E);
}


int SampleQP() {    /* A SAMPLE PROGRAM FOR DFZERO                       */
float f; 
double a,b,re,ae,k,ka,E,A,B,ex; 
int i,j=0;
zeroLink *Z;
ddfortran* Fun;
Fun=&gke;
re=.00001; ae=.00001;
HH=1.;
MM=1.;
if (VV==0.) VV=4;
printf("System is h/2pi=%g   m=%g   V=%g . Half length of the box, a= ?\n",HH,MM,VV);
scanf("%f",&f); 
AA=1.*f;
a=0.0001;
b=VV-.0001;
puts("Calling zero finder\n");
Z=getPointsInterval(Fun, a,b, re,ae);
if (Z==0) {
  puts("NO E<V energy found for this well. Try another.\n");
  return(SampleQP());
}
fprintf(fp,"System is h/2pi=%g   m=%g   V=%g   a=%g\n",HH,MM,VV,AA);
fprintf(fp,"Found(energies):\n");
puts("Found(energies):\n");
if (Z) for(;Z->next;Z=Z->next) {
	printf("(%22.20g,%d)  \n",Z->val,Z->isZero);
	fprintf(fp,"(%22.20g,%d)  \n",Z->val,Z->isZero);
}
if (Z) {
	printf("(%22.20g,%d)  ",Z->val,Z->isZero);
	fprintf(fp,"(%22.20g,%d)  ",Z->val,Z->isZero);
}
puts("\nThe wavefunctions are:\n");
fprintf(fp,"\nThe wavefunctions are:\n");

Z=firstLink(Z);
if (Z) for(;Z->next;Z=Z->next,j++) 
  { LABEL:
   E=Z->val; 
   k=sqrt(2*MM*E)/HH;
   ka=sqrt(2*MM*(VV-E))/HH;
   ex=exp(ka*AA);
   A=1000000/ka*k*cos(k*AA)/ex/(1.+ex*ex);
   B=-1.00000*k/ka*cos(k*AA)*ex*ex/(ex+1./ex);
   printf("psi%d[x_]:=Sin[%22.20f x] /; x<=%g \n",j,k,AA);
   printf("psi%d[x_]:=%22.20f*10^-6 Exp[%22.20f x]+ %22.20f Exp[%22.20f x]  /; x>%g \n",
    j,A,ka,B,-ka,AA);
   fprintf(fp,"psi%d[x_]:=Sin[%22.20f x] /; x<=%g \n",j,k,AA);
   fprintf(fp,"psi%d[x_]:=%22.20f*10^-6 Exp[%22.20f x]+ %22.20f Exp[%22.20f x]  /; x>%g \n",
    j,A,ka,B,-ka,AA); 
  } 
 if (Z) goto LABEL ;
 printf("\n(*where*) a=%g;\n",AA);
 fprintf(fp,"\n(*where*) a=%g;\n\n",AA);
 fprintf(fp,"Plot[psi0[x],{x,0,2*a}]\n\n");
 puts("Plot[psi0[x],{x,0,2*a}]\n\n");
 return(1);
}
 


int main() {
 float f;
 fp=fopen("zsample.dat","w");  /* external fp */
 
 printf("V= ?\n");
 scanf("%f",&f);
 VV=1.*f;
 SampleQP();
 
 return(fclose(fp)); 
}