DECLARE SUB Pager ()
DECLARE SUB Picture1 ()
   PI = 4! * ATN(1!)
30 REM                                                                  
40 REM BARRIER PROBLEM                                                  
50 REM G E COPELAND                                                     
60 REM DEPT OF PHYSICS                                                  
70 REM OLD DOMINION UNIVERSITY                                          
80 REM NORFOLK VA 23508                                                 
90 REM                                                                  
100 REM                                                                  
110 DIM x(370), Y(370), ENERG(370)
120 REM                                                                  
130 PRINT CHR$(12);
140 PRINT "The PERIODIC  POTENTIAL BARRIER PROBLEM"
150 PRINT
160 PRINT "We will examine the QUANTUM MECHANICAL character of a POTENTIAL"
170 PRINT "FUNCTION that is periodic. It consists of repeating units"
180 PRINT "that are barriers of width b separated by regions of zero"
190 PRINT "potential of width a. The barriers are of heigth V0."
200 PRINT "It looks like this."
210 CALL Pager
220 CALL Picture1
230 CALL Pager
240 PRINT "There are TWO cases of interest: E>Vo and E<Vo."
250 PRINT "We will look at the FORMER case here. The later case is"
252 PRINT "covered in program BAND. Let's enter some "
260 PRINT "numerical data. We will use the MKS system of units."
270 PRINT
280 PRINT "HINT: Good trail values are A = 1nm, B = .5 nm  Vo = 40 ev and M=0.0"
290 PRINT
300 PRINT "Input the barrier WIDTH a in nanometers."
310 INPUT A
320 A = A * 1E-09
330 PRINT "Input the well WIDTH b in nanometers."
340 INPUT B
350 B = B * 1E-09
360 PRINT "What about the barrier height? We are in a unit system "
370 PRINT "such that energy is measured in ELECTRON VOLTS. "
380 PRINT "Input V0";
390 INPUT V0
400 V0 = ABS(V0)
410 C1 = 1.6021892D-19
420 V0P = V0 * C1
430 PRINT " Mass of an electron in AMUs ="; .5110041 / 931.5016
440 PRINT "What about the mass of the particle? Input a mass in AMU";
450 INPUT M
460 M = M * 1.6605655D-27
470 CALL Pager
480 H = 6.626176E-34
490 HBAR = H / (2! * PI)
500 PRINT "Schroedinger equation is: u'' +[2m/hbar^2](E-V(x))u(x) =0."
510 PRINT
520 PRINT "where u is the wave function. The potential is periodic,"
530 PRINT " V(x) = V(x+d)  where d + a+b."
540 PRINT " Thus the wavefunctions will exhibit the SAME periodicity."
550 PRINT "  u(x+d) = C u(x)   and"
560 PRINT "            n"
570 PRINT " u(x+nd) = C u(x)"
580 PRINT " If ABS(C) not equal to 1, u ----> oo or -oo."
590 PRINT " Thus C is a PHASE FACTOR. We let it be represented by"
600 PRINT "  C = exp(i phi)."
610 PRINT " Then u(x+d) = exp(i phi)u(x), etc."
620 PRINT
630 CALL Pager
640 PRINT "Let's look at the diagram again."
650 CALL Picture1
660 CALL Pager
670 PRINT "There are two REGIONS of x to worry about. Region I is in"
680 PRINT "hole,0<x<a for 1st cell, and region II, -b<x<0, the barrier."
690 PRINT "These repeat along the x axis an infinite number of times "
700 PRINT "for this one-dimensional crystal. The wavefunctions in the "
710 PRINT "regions are:"
720 PRINT " I      u(x) = A exp(i beta x) + B exp( -i beta x)  -b<x<0"
730 PRINT " II     u(x) = D exp(i alpha x) +F exp(-i alpha x)  0<x<a"
740 PRINT
750 PRINT "where alpha = SQR[ 2mE/hbar^2] and"
760 PRINT "      beta  = SQR[ 2m(E-V0)/hbar^2]."
770 PRINT " In the next period, u(x) is given by"
780 PRINT
785 PRINT "               A exp(i beta(x-d))+B exp( -i beta(x-d))"
790 PRINT "u(x) =exp(i phi)                                                    "
792 PRINT "               Dexp(i alpha(x-d))+Fexp(-i alpha(x-d))"
800 PRINT " where the top option is a<x<d , and bottom d,x,a+d."
810 CALL Pager
820 PRINT "At x=0 and at x=a we apply the boundary conditions, such"
830 PRINT "that u(I)=u(II) and u'(I)=u'(II). This yields 4 conditions."
840 PRINT "       A   +    B   = D   + F                        1"
850 PRINT "       beta(A-B)    = alpha(D-f)                     2"
860 PRINT "exp(i phi)[ A exp(-i beta b) +B exp( i beta b)] =     "
870 PRINT " D exp(i alpha a ) + F exp( -i alpha a)              3"
890 PRINT "                                                                    "
900 PRINT "beta exp(i phi)[ A exp(-i beta b) - B exp( i beta b)] ="
910 PRINT "     alpha[D exp(i alpha a) -F exp(-i alpha a)]      4"
920 PRINT
930 PRINT "Solution is obtained by setting the determinate of the"
932 PRINT "coefficients of A, B, D, and F to zero. Reduction is NOT"
934 PRINT " simple!!!"
935 CALL Pager
936 PRINT " For the case that E>V0, we obtain after much algebra"
940 PRINT
950 PRINT " cos(phi) = cos( alpha a) cos (beta b)  -"
960 PRINT "                2       2"
970 PRINT "           alpha + beta "
980 PRINT "           -------------  sin(alpha a) sin( beta b)"
990 PRINT "           2 alpha beta"
1000 PRINT "                                                                    "
1020 PRINT " alpha and beta BOTH contain energy and the function"
1022 PRINT "on the RHS of cos(phi) is a mixture of them. Note that the"
1023 PRINT "cosine function is bounded to be -1<= cos <= +1. This "
1024 PRINT "forces energy QUANTIZATION. But here, we do NOT get DISCRETE"
1025 PRINT "ENERGY LEVELS. Instead we obtain separate regions where "
1026 PRINT " energy is continuous and not allowed. Thus we have BAND"
1027 PRINT "structure. This is best seen if one plots the RHS vs energy"
1028 PRINT "and also the maximum and minimum values of cosine. This we"
1029 PRINT "will do below."
1030 CALL Pager
1040 PRINT " The function that we have called the RHS will now be "
1050 PRINT "evaluated for many values of (alpha* a). If we call that"
1060 PRINT "function F(E), we will plot it vs alpha a. Values between"
1062 PRINT "-1< F(E) <+1 yield possible energies."
1072 NPTS = 0
1080 FOR ENG = V0 * 1.001 TO 4 * V0 STEP V0 / 89
1082 NPTS = NPTS + 1
1084 EN = ENG * C1
1086 ENERG(NPTS) = ENG
1090 BETA = SQR(2 * M * (EN - V0P)) / HBAR
1100 ALPHA = SQR(2 * M * EN) / HBAR
1110 ANG1 = ALPHA * A
1120 ANG2 = BETA * B
1130 COEF = (ALPHA ^ 2 + BETA ^ 2) / (2 * ALPHA * BETA)
1140 Y(NPTS) = COS(ANG1) * COS(ANG2) - COEF * SIN(ANG1) * SIN(ANG2)
1150 x(NPTS) = ANG1
1160 NEXT ENG
1170 REM FIND THE BANDS                                                   
1180 PRINT "The function F(E) has been evaluated for"; NPTS; " energies."
1199 NO = 0
1200 FOR i = 1 TO NPTS
1202 IF Y(i) <= 1 AND Y(i) >= -1 THEN NO = NO + 1
1204 NEXT i
1210 IF NO <> 0 THEN PRINT "Of these energies, we find"; NO; " NOT";
1212 PRINT " permitted."
1370 CALL Pager
1380 REM  do a graph                                                      
1390 PRINT "Do you wish to graph the RHS function  vs the incident"
1400 PRINT "energy E?<Yes or No>"
1410 INPUT PLOT$
1420 IF PLOT$ = "YES" OR PLOT$ = "yes" OR PLOT$ = "Yes" THEN 1830
1430 IF PLOT$ = "NO" OR PLOT$ = "no" OR PLOT$ = "No" THEN STOP
1440 REM                                                                  

1770 STOP
1830 REM TEKTRONIX OUTPUT                                               
1840 OPEN "PLTTEK.DAT" FOR OUTPUT AS #2
1850 PRINT #2, 1
1860 PRINT #2, 3
1870 PRINT #2, NPTS
1880 Y$ = "    VALUE OF RHS"
1890 x$ = " ALPHA * A   "
1900 PRINT #2, Y$
1910 PRINT #2, x$
1920 PRINT #2, "RHS Vs E from BAND A="; A * 1E+09; " NM B="; B * 1E+09; " NM"
1930 PRINT #2, "Mass = "; M / 1.6605655D-27; " AMU. E>V0 & valid Energy inside +/-1"
1940 FOR i = 1 TO NPTS
1950 PRINT #2, x(i), Y(i)
1960 NEXT i
FOR i = 1 TO NPTS
PRINT #2, x(i), 1
NEXT i
FOR i = 1 TO NPTS
PRINT #2, x(i), -1
NEXT i

1970 CLOSE #2
1980 PRINT "Now return to windows and run WPLOT."
CHAIN "Tek2wplt.exe"
1990 STOP

2310 END

SUB Pager
 PRINT , "Push RETURN to continue.";
 INPUT T$
 PRINT CHR$(12);
END SUB

SUB Picture1
2130 PRINT , , " ^"; "V(x)"
2140 PRINT , , "I"
2150 PRINT , "-/\/\/\/\/\/\/\/\/\/\/\/\/\/-> E"
2160 PRINT , , "I"
2170 PRINT "-       ------       ------       ------       ------";
2180 PRINT "        ------       ------"
2190 FOR JJ = 1 TO 5
2200 PRINT "|       |    |       |    |       |    |       |    |";
2210 PRINT "        |    |       |    |"
2220 NEXT JJ
2230 FOR i = 1 TO 75
2240 PRINT "-";
2250 NEXT i
2260 PRINT "> x"
2270 PRINT "|  a    | b  |   a   | b  |   a   | b  |   a   | b  |";
2280 PRINT "   a    | b  |   a   | b  |"
2290 PRINT
END SUB