DECLARE FUNCTION sinh! (x!)
DECLARE SUB picture1 ()
DECLARE SUB pager ()

60 REM                                                                    
80 REM BARRIER PROBLEM                                                    
100 REM G E COPELAND                                                       
120 REM DEPT OF PHYSICS                                                    
140 REM OLD DOMINION UNIVERSITY                                            
160 REM NORFOLK VA 23508                                                   
REM   Updated Fall 2003                                                                
200 REM                                                                    
240 DIM x(300), Y(300)
280 REM                                                                    
300 PRINT CHR$(12);
320 PRINT "THE POTENTIAL BARRIER PROBLEM"
340 PRINT
360 PRINT "We will examine the QUANTUM MECHANICAL character of a POTENTIAL"
380 PRINT "BARRIER that is of width A and of height  V0."
400 PRINT "It looks like this."
420 CALL pager
440 CALL picture1
460 CALL pager
480 PRINT "There are TWO cases of interest: E>Vo and E<Vo."
500 PRINT "We will look at the later case first. Let's enter some"
520 PRINT "numerical data. We will use the MKS system of units."
540 PRINT
560 PRINT "HINT: Good trail values are A = 1nm, Vo = 40 ev and M=0.001 AMU."
580 PRINT
600 PRINT "Input the barrier WIDTH A in nanometers."
620 INPUT A
660 A = A * 1E-09
680 PRINT "What about the barrier height? We are in a unit system "
700 PRINT "such that energy is measured in ELECTRON VOLTS. "
720 PRINT "Input V0";
740 INPUT V0
760 V0 = ABS(V0)
770 C1 = 1.6021892D-19
780 V0P = V0 * C1
782 PRINT " Mass of an electron in AMUs ="; .5110041 / 931.5016
820 PRINT "What about the mass of the particle? Input a mass in AMU";
840 INPUT M
860 M = M * 1.6605655D-27
880 CALL pager
900 H = 6.626176E-34
pi = 4! * ATN(1!)
920 HBAR = H / (2! * pi)
940 PRINT " We will solve this problem in reduced variables. That is,"
960 PRINT "energy of the incoming particle and the barrier height will"
980 PRINT "be rescaled.  Let E = 2m E'/hbar**2 and Vo = 2m Vo'/hbar**2"
1000 PRINT "where hbar = h/(2pi), E' and Vo' are the energies in ev."
1020 CALL pager
1040 PRINT "Let's look at the diagram again."
1060 CALL picture1
1064 CALL pager
1080 PRINT " The solutions in the THREE regions must be matched at the"
1100 PRINT "boundries. The Schrodinger wave equations and its solutions"
1120 PRINT "in these regions are:"
1140 PRINT "Region1:  u'' = -Eu     &   u = A exp(ik1 x) + B exp(-ik1 x)"
1160 PRINT "Region2:  u'' = -(E-V)u &   u = C exp( k2 x) + D exp(- k2 x)"
1180 PRINT "Region3:  u'' = -Eu     &   u = G exp(ik1 x) + F exp(-ik1 x)"
1200 PRINT  'where:'                                                         
1220 PRINT " the symbol ' = d/dx and k1 = SQRT(E) and k2 = SQRT(Vo-E)."
1222 PRINT " The constant F =0, since it represents particles moving in the"
1224 PRINT "the -x direction. Next we apply the boundry conditions at x=0 and"
1226 PRINT "x=A. u1(0)=u2(0), u1'(0)=u2'(0), u2(A)=u3(A) & u2'(A)=u3'(A)."
1228 PRINT
1230 PRINT "We get 4 equations with 5 constants."
1232 PRINT "Solving for ratios of the constants A/G,B/G, etc. can be done."
1234 PRINT " The probability of transmission of the incident particle is:"
1236 PRINT "            transmitted current (x>a)   v3 G^* G   G^* G"
1238 PRINT " T   = ------------------------------ = -------- = -----"
1240 PRINT "           incident current (x<0 )      v1 A^* A   A^* A "
1242 PRINT " ^* means complex conjugate."
1244 PRINT "The velocities v3 and v1 are equal, and after much algebra we"
1246 PRINT "find a value of T."
1250 CALL pager
1251 PRINT "                  2"
1252 PRINT "               Vo'            2                           ^-1"
1254 PRINT " T = {    --------------  sinh [ 2m(Vo' -E')]a/hbar } +1 }"
1256 PRINT "        4 E'(Vo'-E')"
1258 PRINT
1260 PRINT "The problem can be worked again for the case of E'>Vo' !!!"
1262 PRINT "To do such set u2 = C exp(i K x) + D exp(-i K x)"
1264 PRINT "where K = SQRT( E-Vo), i.e. using reduced values. This is "
1266 PRINT "very messy, but luckily sinh(ix)=i sin(x), so we can get T."
1271 PRINT "                  2"
1272 PRINT "               Vo'            2                           ^-1"
1274 PRINT " T = {    --------------  sin  [ 2m(E' -Vo')]a/hbar } +1 }"
1276 PRINT "        4 E'(E'-Vo')"
1278 PRINT
1280 CALL pager
1300 PRINT "You have already input the barrier height Vo', the width A,"
1310 PRINT " and the mass of the incident particle. Now select an incident"
1312 PRINT "energy?(in ev)"
1314 INPUT EP
1316 EP = EP * C1
1320 IF EP > V0P THEN PRINT " E > V0 CASE."
1322 IF EP <= V0P THEN PRINT " E < V0 CASE"
1324 IF EP > V0P THEN 1400
1330 REM E < V0 CASE                                                        
1332 INSIDE = SQR(2 * M * (V0P - EP)) * A / HBAR
1334 COEF = V0P ^ 2 / (4! * EP * (V0P - EP))
1340 T = 1! + COEF * sinh(INSIDE) ^ 2
1342 T = 1! / T
1344 PRINT "The transmission coefficient T ="; T
1350 PRINT " Do you want to input another energy?(Yes or No)"
1352 INPUT ENG$
1354 IF ENG$ = "YES" OR ENG$ = "yes" OR ENG$ = "Yes" THEN GOTO 1300
1356 IF ENG$ = "NO" OR ENG$ = "no" OR ENG$ = "No" THEN GOTO 1500
1358  PRINT "I don't understand your answer."
1359 GOTO 1350
1400 REM E > V0 CASE                                                        
1432 INSIDE = SQR(2 * M * (EP - V0P)) * A / HBAR
1434 COEF = V0P ^ 2 / (4! * EP * (EP - V0P))
1440 T = COEF * SIN(INSIDE) ^ 2 + 1!
1442 T = 1! / T
1444 PRINT "The transmission coefficient T ="; T
1446 PRINT
1448 PRINT " NOTE: when   SQRT( 2m(E'-Vo')) A/hbar = n Pi   n = 1,2,3,4,... "
1449 PRINT "  when n is integer, T =1. Between cases where n is interger"
1450 PRINT "T will fall to some value less than one. This could NEVER"
1452 PRINT "happen in a classical system."
1455 PRINT
1499 GOTO 1350
1500 CALL pager
1510 REM  do a graph                                                        
1520 PRINT "Do you wish to graph the transmission function T vs the incident"
1530 PRINT "energy E?<Yes or No>"
1540 INPUT PLOT$
1542 IF PLOT$ = "YES" OR PLOT$ = "yes" OR PLOT$ = "Yes" THEN 1550
1544 IF PLOT$ = "NO" OR PLOT$ = "no" OR PLOT$ = "No" THEN STOP
1550 NPTS = 0
1551 V0P2 = V0P ^ 2
1552 FOR E = .8 * EP TO 1.4 * EP STEP EP / 200
1554 NPTS = NPTS + 1
1555 x(NPTS) = E / C1
1556 IF E > V0P THEN 1580
1558 REM E < V0 CASE                                                        
1560 INSIDE = SQR(2 * M * (V0P - E)) * A / HBAR
1562 COEF = V0P2 / (4! * E * (V0P - E))
1564 T = 1! + COEF * sinh(INSIDE) ^ 2
1566 Y(NPTS) = 1! / T
1568 GOTO 1595
1580 REM E > V0 CASE                                                        
1582 INSIDE = SQR(2 * M * (E - V0P)) * A / HBAR
1584 COEF = V0P2 / (4! * E * (E - V0P))
1586 T = COEF * SIN(INSIDE) ^ 2 + 1!
1588 Y(NPTS) = 1! / T
1595 NEXT E
4560 REM TEKTRONIX OUTPUT                                                 
PRINT "Opening PLTTEK.DAT"
4580 OPEN "PLTTEK.dat" FOR OUTPUT AS #2
4600 PRINT #2, 1
4620 PRINT #2, 1
4640 PRINT #2, NPTS
4660 Y$ = "    TRANSMISSION"
4680 x$ = " ENERGY IN EV"
4700 PRINT #2, Y$
4720 PRINT #2, x$
4740 PRINT #2, "T VS E FROM TUNNEL A="; A * 1E+09; " NM"
4820 PRINT #2, " MASS = "; M / 1.6605655D-27; " AMU."
4840 FOR I = 1 TO NPTS
4860 PRINT #2, x(I), Y(I)
4880 NEXT I
4900 CLOSE #2
4920 PRINT "Now return to chain to TEK2WPLT for graphics"
CHAIN "tek2wplt.bas"
4940 STOP
4960 REM NUMERICAL OUTPUT                                                   
4980 PRINT
5000 PRINT "   NUMERICAL OUTPUT OF PROBABILITY FUNCTION"
5020 PRINT
5040 PRINT "STATE ="; N2; " FOR "; A$; " PARITY WAVEFUNCTIONS UNNORMALIZED"
5060 PRINT
5080 PRINT "DISTANCE           PROBABILITY"
5100 FOR I = 1 TO J
5120 PRINT x(I) * 1E+09, U(I) ^ 2
5140 NEXT I
5160 STOP
5180 REM ODD PARITY CALCULATIONS                                            

5540 END

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

SUB picture1
5220 PRINT , , "^"; "V(x)"
5240 FOR I = 1 TO 4
5260 PRINT , , "I"
5280 NEXT I
5300 PRINT , , "---------------------"
5320 FOR JJ = 1 TO 10
5340 PRINT , , "|                   |"
5360 NEXT JJ
5380 PRINT "0";
5400 FOR I = 1 TO 73
5420 PRINT "-";
5440 NEXT I
5460 PRINT "> X"
5480 PRINT , , "0                   A"
5500 PRINT "region 1", , "      region 2 ", "         region 3"

END SUB

FUNCTION sinh (x)
ans = -.5 * (EXP(x) - EXP(-x))
sinh = ans
END FUNCTION