DEFSNG A-H, O-Z
DECLARE SUB EVENFUN (N2!, B!(), A!, U!(), K1!, A3!(), X!(), J!)

DECLARE SUB AskforInput (A, V0, M)
DECLARE SUB pict1 ()


DECLARE SUB PAGER ()
DECLARE SUB EVENFUN (N2, B(), A, U(), K1, A3(), X(), J)
DECLARE SUB ODDFUN (N2, B(), A, U(), K2, A3(), X(), J)
DECLARE SUB ODDONE (PI, R, R2, E(), A3(), B(), N1, C2, C3, V0)
DECLARE SUB EVENONE (PI, R, R2, E(), A3(), B(), N1, C2, C3, V0)


20 REM WELL PROBLEM                                                
30 REM G E COPELAND                                                
40 REM DEPT OF PHYSICS                                             
50 REM OLD DOMINION UNIVERSITY                                     
60 REM NORFOLK VA 23508                                            
70 REM                                                             
80 REM                                                           
REM updated 3/6/2001

90 DIM SHARED E(50), A3(50), B(50)
100 DIM SHARED X(900), U(900)
110 DIM P1(100), P2(100)
120 REM                                                           
PI = 4! * ATN(1!)
200 REM NORMALIZATION CONSTANT =1                                   
210  REMH = H / (2! * PI)
220 H = 1.0546E-34
230 REM      2=H*H                                                
    C1 = 1.6021917D-19
240 PRINT CHR$(12);
250 PRINT "THE POTENTIAL WELL PROBLEM"
260 PRINT
270 PRINT "We will examine the QUANTUM MECHANICAL character of a potential"
280 PRINT "well that is of WIDTH 2*A and of depth V0."
290 PRINT "It looks like this:"
   CALL PAGER
   CALL pict1
   CALL PAGER
   CALL AskforInput(A, V0, M)

730 CALL PAGER
740 REM                                                             
750 PRINT " As you certainly know by now, there are 4 BOUNDARY CONDITIONS"
760 PRINT "that must be applied to the SCHROEDINGER EQUATION for this"
770 PRINT "POTENTIAL. We do such and get a 4 simultaneous equations.          "
780 PRINT "Using Kramers' rule, we set the determinate of the coefficents"
782 PRINT "equal to zero. This forces the quantization of energy."
790 PRINT "It turns out that there are TWO ways that DET=0 can"
800 PRINT "be obtained and these correspond to two TYPES OF ALLOWED"
810 PRINT "WAVEFUNCTIONS. The two types correspond to ODD and EVEN PARITY."
820 PRINT "Which TYPE are you interested in studying?"
825 INPUT A$
840 IF A$ = "ODD" OR A$ = "odd" THEN GOTO 920
850 IF A$ = "EVEN" OR A$ = "even" THEN GOTO 870
860 GOTO 820
870 REM EVEN PARITY                                                 
880 REM PRINT H,A,M,V0                                              
890 IF A$ = "ODD" OR A$ = "odd" THEN GOTO 920
900 PRINT "We must solve the equation   tan(BETA*A)=ALPHA/BETA"
910 GOTO 930
920 PRINT "We must solve the equation   cot(BETA*A)=-ALPHA/BETA"
930 PRINT "where ALPHA = SQR( 2*M*E/HBAR^2)"
940 PRINT " and BETA = SQR( 2*M*(V0-E)/HBAR^2) "
950 PRINT
960 PRINT "Let X= BETA *A   and Y = ALPHA *A , then a solution"
970 PRINT "will be done NUMERICALY by SEARCHING for "
980 IF A$ = "ODD" OR A$ = "odd" THEN GOTO 1010
990 PRINT "cases where X*TAN(X) =Y1 and X^2+Y^2 =R^2"
1000 GOTO 1020
1010 PRINT "cases where X*COT(X) =-Y1 and X^2+Y^2 = R^2"
1020 PRINT "have intersections, i.e. common values when plotted as a FUNCTION"
1030 PRINT "of X."
1040 R = SQR(2 * M) * SQR(V0) * (A / H)
1050 PRINT "The RADIUS of CIRCLE <X^2+Y^2 =R^2> is ="; R; " RADIANS."
1052 R2 = R * R
1062 CALL PAGER
1064 NROOT = INT(2 * R / PI)
1065 PRINT "MAXIMUM NUMBER of ODD and EVEN solutions ="; NROOT
1066 PRINT "IF this number > about 40, best to start over since solution"
1067 PRINT "time will be VERY LONG."
1070 PRINT A$; "  S O L U T I O N S"
1080 PRINT "NUMBER      Y1       Y2        Y2/Y1          X"
1110 REM WE SEARCH AGAIN.                                            
1130 N1 = 0
1140 REM N1 = QUANTUM NUMBER =1,2,3,.... FOR EVEN PARITY CASE.       
1160 C2 = SQR(2! * M) / H
1170 C3 = (H / (2! * M * A)) * (H / A)
1175 IF A$ = "ODD" OR A$ = "odd" THEN
         CALL ODDONE(PI, R, R2, E(), A3(), B(), N1, C2, C3, V0)
         ELSE
                IF A$ = "EVEN" OR A$ = "even" THEN
                CALL EVENONE(PI, R, R2, E(), A3(), B(), N1, C2, C3, V0)
         ELSE
         END IF
         END IF

1410 IF N1 = 0 THEN
         PRINT "This well has NO DISCRETE LEVELS! At least with the CURRENT"
         PRINT "SEARCH PARAMETER ."
 STOP
ELSE
END IF
 REM PRINT OUT OF VALUES                                        

 CALL PAGER
1480 PRINT "QUANTUM NO.        BETA        ENERGY(EV)       ALPHA"
1490 FOR I = 1 TO N1
1580 PRINT I, B(I), E(I) / C1, A3(I)
1590 NEXT I
1600 REM GO DO THE WAVEFUNCTIONS                                     
1610 IF A$ = "ODD" OR A$ = "odd" THEN GOTO 1640
1620 PRINT " There are "; N1; " EVEN PARITY wave functions."
1630 GOTO 1650
1640 PRINT " There are "; N1; " ODD PARITY wave functions."
1650 CALL PAGER
1660 PRINT "These wave functions will be obtained for -3*A <x< 3*A."
1670 PRINT "Output can be in one of three options:"
1700 PRINT "1 graphical OUTPUT OF THE PROB DENSITY FUNCTIONS on the"
1710 PRINT "for use in WPLOT; "
1711 PRINT "2 graphical OUTPUT OF THE WAVEFUNCTION on the"
1712 PRINT "for use in WPLOT; "
1713 PRINT "3 NUMERICAL DATA FOR WAVEFUNCTION.                                  "
1730 PRINT "Input the NUMBER of the one you wish to have?<1,2,3>"
1740 INPUT N8
1750 IF N8 <> INT(N8) THEN GOTO 1730
1760 IF N8 > 3 THEN GOTO 1730
1770 IF N8 < 1 THEN GOTO 1730
1780 REM MUST BE 1 2 3                                          
1800 PRINT
1860 PRINT "Input the NUMBER of the wavefunction to be examined?<1 to"; N1; ">"
1870 INPUT N2
1930 J = 0
1940 IF A$ = "ODD" OR A$ = "odd" THEN
        CALL ODDFUN(N2, B(), A, U(), K2, A3(), X(), J)
        ELSE
         IF A$ = "EVEN" OR A$ = "even" THEN
         CALL EVENFUN(N2, B(), A, U(), K1, A3(), X(), J)
         END IF
         
         END IF


1948 REM ------------------------NORMALIZE THE FUNCTIONS             
1952 REM NORMALIZE SOULUTIONS                                        
1954  SUMIT = 0!
1956  FOR I = 1 TO J
1958  SUMIT = SUMIT + U(I) * U(I)
1960  NEXT I
1961  SUMIT = SQR(SUMIT)
1962  FOR I = 1 TO J
1964  U(I) = U(I) / SUMIT
1966  NEXT I
2200 REM J= # OF DATA POINTS                                         

2240 A1$ = " Psi<"
2250 A2$ = VALUE$(N2)
2260 A3$ = "> Normalized for Vo ="
2270 A4$ = STR$(V0 / C1)
2280 A5$ = " & 2*A ="
2290 A6$ = STR$(2 * A * 1E+09)
2300 C$ = A$ + A1$ + A2$ + A3$ + A4$ + A5$ + A6$
2310 IF N8 = 1 OR N8 = 2 THEN GOTO 2430
2320 IF N8 = 3 THEN 2610


2430 REM Graphics file output                                       
2440 OPEN "PLTTEK.DAT" FOR OUTPUT AS #2
2450 PRINT #2, 1
2460 PRINT #2, 1
2470 PRINT #2, J
2480 IF N8 = 1 THEN Y$ = "    Probability"
2482 IF N8 = 2 THEN Y$ = "    Value"
2490 X$ = " Distance <nm>------>"
2500 PRINT #2, Y$
2510 PRINT #2, X$
2520 PRINT #2, C$
2530 E$ = " Parity Wavefunctions"
2540 N$ = STR$(N2)
2550 N1$ = " N ="
2560 PRINT #2, A$ + E$ + N1$ + N$
2562 IF N8 = 2 THEN 2580
2570 FOR I = 1 TO J
2572 PRINT #2, 1E+09 * X(I), U(I) * U(I)
2574 NEXT I
2575 GOTO 2592
2580 FOR I = 1 TO J
2582 PRINT #2, 1E+09 * X(I), U(I)
2584 NEXT I

2592 CLOSE #2
CHAIN "Tek2wplt.bas"
2594 PRINT "Now return to windows and run WPLOT."
2600 STOP
REM
REM

2610 REM NUMERICAL OUTPUT                                        
2620 PRINT
2630 PRINT "   NUMERICAL OUTPUT OF PROBABILITY FUNCTION"
2640 PRINT
2650 PRINT "STATE ="; N2; " FOR "; A$; " PARITY WAVEFUNCTIONS NORMALIZED"
2660 PRINT
2670 PRINT "DISTANCE   WAVEFUNCTION       PROBABILITY"
2680 FOR I = 1 TO J STEP INT(J / 20)
2690 PRINT USING "##.########  "; X(I) * 1E+09; U(I); U(I) ^ 2
2700 NEXT I
2710 STOP
4600 END

   SUB AskforInput (A, V0, M)
 PRINT "There are two cases of interest: E>V0 and E<V0."
 PRINT "We will look at the later case here. Let's enter some"
 PRINT "numerical data. We will use the MKS system of units."
 PRINT
 PRINT "HINT: GOOD TRIAL VALUES ARE A = 1 NM,V0=40 EV AND M=0.001 AMU."
 PRINT
 PRINT "Input the well HALF WIDTH A in nanometers. I will figure out 2*A";
 INPUT A
 A9 = A
 A = A * 1E-09
 PRINT "What about the WELL DEPTH? We are in a unit system "
650 PRINT "such that energy is measured in ELECTRON VOLTS. "
660 PRINT "INPUT V0";
670 INPUT V0
671 V0 = ABS(V0)
680 V0 = V0 * 1.6021917D-19
700 PRINT "What about the MASS of the particle? Input a mass in AMU";
710 INPUT M
720 M = M * 1.660531E-27

END SUB

 SUB EVENFUN (N2, B(), A, U(), K1, A3(), X(), J)
4305 K1 = COS(B(N2) * A)
4310 FOR X1 = -3 * A TO 3 * A STEP A / 40
4315 J = J + 1
4320 IF X1 > -A THEN GOTO 4345
4325 REM -3A TO -A                                                 
4330 U(J) = K1 * EXP(A3(N2) * (X1 + A))
4335 GOTO 4380
4340 REM X1>-A                                                     
4345 IF X1 > A THEN GOTO 4370
4350 U(J) = COS(B(N2) * X1)
4365 GOTO 4380
4370 REM X1>A                                                      
4375 U(J) = K1 * EXP(-A3(N2) * (X1 - A))
4380 X(J) = X1
4390 NEXT X1

END SUB

 SUB EVENONE (PI, R, R2, E(), A3(), B(), N1, C2, C3, V0)
3151 XINC = .0001
3155 X1 = .4 * PI
3160 IF X1 > R THEN EXIT SUB
3250 Y2 = X1 * TAN(X1)
3260 IF Y2 < 0 THEN GOTO 3392
3262 Y1 = SQR(R2 - X1 * X1)
3310 REM IF Y1>Y2+.5  THEN  GOTO 3392                              
3314 REM  IF Y2>Y1+.5 THEN GOTO 3392                               
3315 IF Y1 > Y2 + .1 THEN GOTO 3392
3316 IF Y2 > Y1 + .1 THEN GOTO 3392
3320 REM FOUND A ROOT                                              
3340 N1 = N1 + 1
3350 E(N1) = -Y1 * Y1 * C3
3360 A3(N1) = C2 * SQR(ABS(E(N1)))
3370 B(N1) = C2 * SQR(V0 - ABS(E(N1)))
3380 PRINT USING " ###.#######  "; N1; Y1; Y2; Y1 / Y2; X1
3382 X1 = X1 + .8 * PI
3384 GOTO 3160
3392 X1 = X1 + XINC
3395 GOTO 3160

END SUB

 SUB ODDFUN (N2, B(), A, U(), K2, A3(), X(), J)
4110 K2 = SIN(B(N2) * A)
4120 FOR X1 = -3 * A TO 3 * A STEP A / 40
4130 J = J + 1
4140 IF X1 > -A THEN GOTO 4180
4150 REM -3A TO -A                                                 
4160 U(J) = -K2 * EXP(A3(N2) * (X1 + A))
4170 GOTO 4230
4180 REM X1>-A                                                     
4190 IF X1 > A THEN GOTO 4210
4200 U(J) = SIN(B(N2) * X1)
4205 GOTO 4230
4210 REM X1>A                                                      
4220 U(J) = K2 * EXP(-A3(N2) * (X1 - A))
4230 X(J) = X1
4240 NEXT X1

END SUB

3500 SUB ODDONE (PI, R, R2, E(), A3(), B(), N1, C2, C3, V0)
3580 X1 = PI * .8
3582 XINC = .001
3630 IF X1 >= R THEN EXIT SUB
3680 REM
REM  y2=-x1*cot(x1)
     Y2 = -X1 / TAN(X1)
3690  IF Y2 < 0 THEN GOTO 3900
3692 Y1 = SQR(R2 - X1 * X1)
3693 IF Y1 > Y2 + .1 THEN GOTO 3900
3700 IF Y2 > Y1 + .1 THEN GOTO 3900
3820 REM  FOUND A ROOT                                             
3840 N1 = N1 + 1
3850 E(N1) = -Y1 * Y1 * C3
3860 A3(N1) = C2 * SQR(ABS(E(N1)))
3870 B(N1) = C2 * SQR(V0 - ABS(E(N1)))
3880 PRINT USING " #####.##### "; N1; Y1; Y2; Y1 / Y2; X1
3890 X1 = X1 + PI * .8
3892 GOTO 3630
3900 X1 = X1 + XINC
3902 GOTO 3630

END SUB

 SUB PAGER
 PRINT , "Push  RETURN to continue."
 INPUT dum$
 PRINT CHR$(12);

END SUB

     SUB pict1
310 PRINT , , "^"; "V(X)"
320 FOR I = 1 TO 5
330 PRINT , , "I"
340 NEXT I
341 PRINT "0";
350 FOR I = 1 TO 14
360 PRINT "-";
370 NEXT I
380 PRINT , "I",
385                                                               
390 FOR I = 1 TO 14
400 PRINT "-";
410 NEXT I
411 PRINT " X->"
420                                                               
430 FOR I = 1 TO 5
440 PRINT , "|", "I", "|"
450 NEXT I
451 PRINT "===    -V0   ";
460 FOR I = 1 TO 59
470 PRINT "=";
480 NEXT I
490 PRINT
500 PRINT , "-A", "0", "+A";
520 PRINT

END SUB