DECLARE SUB Pager ()
DIM f(800), V(800)
PRINT CHR$(12);
REM GE Copeland
REM Department of Physics
REM Old Dominion University
REM Norfolk VA
REM 3/2001
REM Boltzmann Distrbution for P-P collisions
REM T = 2.32e9
REM vo=speed of proton rem for t=150e+6 K



PRINT "Fusion Process for p-p chain"
PRINT " We calculate the Maxwell-Boltzman Speed Distribution Function"
PRINT "for protons at a given temperature."
REM MKS Sytem of Units
qp = 1.60217733D-19
mp = 1.6726231D-27
k = 1.380658E-23
h = 6.6260755D-34
ke = 8987551787#
pi = 4! * ATN(1!)
PRINT "The temperature at center of the sun is about 15e6 K."
PRINT "Keep the value of T near that value."
PRINT "Input the temperature in Kelvin"
INPUT t
REM  t = 1.5E+07
PRINT " We deduce what classical fraction of the protons will get close enough to each"
PRINT "other for the fusion process  p+p = d + e(+) + ve to happen."


con = (mp / (2! * pi * k * t)) ^ 1.5
   c1 = 4! * pi
PRINT " We select to focus on a relative approach speed of ";
V0 = 7580000!
PRINT V0; "m/s"

W = (.5 * mp * V0 ^ 2)
Teq = W * 2! / (3! * k)
PRINT " This correspondes to a classical temperature of "; Teq; " K"
rclose = 8.999999E+09 * qp ^ 2 / W
PRINT " The classical turning point is at "; rclose * 1E+15; " fermi."
 CALL Pager
 PRINT "We now calculate the speed distribution and determine what fraction"
 PRINT "has speed larger than "; V0

vend = 5! * V0
vstep = V0 * .005
sum = 0
      iprint = 8
        j = 0
        ip = 0
        PRINT "V/V0             f(v)"
        FOR vi = vstep TO vend STEP vstep
     
        f1 = c1 * con * vi ^ 2 * EXP(-(mp * vi ^ 2) / (2 * k * t))
       
       IF f1 = 0! THEN GOTO 4
           ip = ip + 1
           j = j + 1
           V(j) = vi
           f(j) = f1
       IF ip = iprint THEN
       PRINT vi / V0, f(j)
       ip = 0
       ELSE
       END IF
        sum = sum + f(j) * vstep
    NEXT vi
4  PRINT "Sum of the distribution function over all v ="; sum
    CALL Pager
Peak = -.1
FOR i = 1 TO j
IF f(i) > Peak THEN
Peak = f(i)
isave = i
REM PRINT "Isave "; isave
ELSE
END IF
NEXT i

REM now find how much has speed greater than V0
sumv = 0
FOR i = 1 TO j
REM PRINT i;
IF V(i) > V0 THEN
sumv = sumv + f(i) * vstep
REM PRINT "Inside IF in loop"
ELSE
END IF
      NEXT i
OPEN "plttek.dat" FOR OUTPUT AS #2
PRINT "Peak value of f(v) happens when V/v0 = "; V(isave) / V0
PRINT "Fraction of distribution function > v0 = "; sumv
PRINT "Does that surprise you?"
CALL Pager
tq = (4! * (mp / 2!))
t1 = (qp ^ 2) * ke
tq1 = (qp ^ 2 / (3! * k)) * ke
tq = tq / h
t1 = t1 / h
tempq = tq * t1 * tq1
PRINT "The temperature -adjusted for tunneling is just"; tempq; " Kelvin"

        PRINT #2, 1
        PRINT #2, 1
        PRINT #2, j
        PRINT #2, "f(v)"
        PRINT #2, "V in v0 units"
        PRINT #2, "Maxwell v dist  v0="; V0; " m/s p+p fusion"
        PRINT #2, "Temp classical = "; t; " K & Peak Value @"; V(isave) / V0
        FOR i = 1 TO j
        PRINT #2, V(i) / V0, f(i)
        NEXT i
        CLOSE #2
     PRINT "Data file will be written and we chain to TEK2WPLT so you run in windows"
     PRINT "the WPLOT program."
     CALL Pager
     CHAIN "TEK2WPLT.EXE"
        END

 SUB Pager
PRINT " Press Enter to continue";
INPUT Dun$
PRINT CHR$(12);
END SUB