' E-TELSA - A Program to Calculate the Resonant Frequency of Tesla Coils
' Date 1-29-1999  V3.0  Final
' This program is released to the public domain

      pi = 3.141593
      Cap1 = 1E-13
      COLOR 14, 1: CLS : PRINT : PRINT
      PRINT " E-TESLA - Tesla Coil Resonant Frequency Calculator"
      PRINT
      PRINT " Terry Fritz    January 29, 1999"
      PRINT
      PRINT " Version 3.0"
      SLEEP 3

1000  CLS : PRINT : PRINT
2000  PRINT " Input an ACCURACY factor 1 (min) to 5 (max)."
      PRINT
      PRINT " The higher the number, the longer it takes to calculate."
      PRINT
      INPUT Acc
      Acc = INT(Acc)
      IF Acc < 1 OR Acc > 5 THEN PRINT " Must be between 1 and 5!": GOTO 2000
      Acc = Acc / 5

      CLS : PRINT
      PRINT
      PRINT " Input the measurement UNITs you are using."
      PRINT
      Mult = 0
      PRINT " 1 = mm"
      PRINT " 2 = cm"
      PRINT " 3 = inches"
      PRINT " 4 = feet"
      PRINT " 5 = meters"
      PRINT
      PRINT " All dimension inputs will be in these units."
      PRINT
2100  INPUT " Units "; Mult
      IF Mult < 1 OR Mult > 5 OR Mult <> INT(Mult) THEN PRINT " Try again!": GOTO 2100
      IF Mult = 1 THEN Mult = .03937
      IF Mult = 2 THEN Mult = .3937
      IF Mult = 3 THEN Mult = 1
      IF Mult = 4 THEN Mult = 12
      IF Mult = 5 THEN Mult = 39.37

      CLS : PRINT
2200  PRINT " How HIGH is the CEILING?"
      PRINT
      INPUT Ceil
      IF Ceil <= 0 THEN PRINT " Must be greater than zero!": GOTO 2200
      AcF = 120 / (Ceil) * Acc

      CLS : PRINT
2300  PRINT " How FAR is the coil from the WALLs?"
      PRINT
      INPUT Wall
      IF Wall <= 0 THEN PRINT " Must be greater than zero!": GOTO 2300
      IF 2 * Wall < Ceil THEN PRINT " The walls are too close to the coil!": GOTO 2300

      CLS : PRINT
2400  PRINT " What is the HEIGHT, from the floor, to the INNER TURN of the PRIMARY coil? "
      PRINT
      INPUT PriH
      IF PriH < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 2400

      CLS : PRINT
2500  PRINT " What is the INNER DIAMETER of your PRIMARY coil?"
      PRINT
      INPUT PriD0
      IF PriD0 < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 2500
 
      CLS : PRINT
2600  PRINT " What is the OUTER DIAMETER of your PRIMARY coil?"
      PRINT
      INPUT PriD
      IF PriD < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 2600

      CLS : PRINT
2700  PRINT " What is the HEIGHT, from the floor, to the OUTER TURN of the PRIMARY?"
      PRINT
      INPUT PriH1
      IF PriH1 < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 2700

      CLS : PRINT
2800  PRINT " What is the HEIGHT, from the floor, to the FIRST WINDING of the SECONDARY coil?"
      PRINT
      INPUT SecB
      IF SecB < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 2800
     
      CLS : PRINT
2900  PRINT " What is the DIAMETER of your SECONDARY coil?"
      PRINT
      INPUT SecD
      IF SecD < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 2900
     
      CLS : PRINT
3000  PRINT " What is the LENGTH of your SECONDARY coil?"
      PRINT
      INPUT SecL
      IF SecL < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 3000

      CLS : PRINT
3100  PRINT " What is the INDUCTANCE in MILLI-HENRIES of the SECONDARY?"
      PRINT
      INPUT SecInd
      SecInd = SecInd / 1000
      IF SecInd <= 0 THEN PRINT " Must be greater than zero!": GOTO 3100

      CLS : PRINT
3200  PRINT " What is the HEIGHT, from the floor, to the CENTER of the FIRST TOP TERMINAL?"
      PRINT
      INPUT TermH
      IF TermH < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 3200
      IF TermH = 0 THEN TermH = SecB + SecL

      CLS : PRINT
3300  PRINT " What is the RADIUS of the FIRST TERMINAL (sphere or torroid section)?"
      PRINT
      INPUT TermR
      IF TermR < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 3300

      CLS : PRINT
3400  PRINT " What is the DIAMETER of the TERMINAL (center to center if it is a TORROID)?"
      PRINT " Enter 0 if it is a sphere."
      PRINT
      INPUT TermD
      IF TermD < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 3400

      CLS : PRINT
3500  PRINT " What is the HEIGHT, from the floor, to the CENTER of the SECOND TOP TERMINAL?"
      PRINT
      INPUT TermH2
      IF TermH2 < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 3500
      IF TermH2 < TermH THEN TermH2 = TermH

      CLS : PRINT
3600  PRINT " What is the RADIUS of the SECOND TERMINAL (sphere or torroid section)?"
      PRINT
      INPUT TermR2
      IF TermR2 < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 3600

      CLS : PRINT
3700  PRINT " What is the DIAMETER (center to center if it is a TORROID)?"
      PRINT " Enter 0 if it is a sphere."
      PRINT
      INPUT TermD2
      IF TermD2 < 0 THEN PRINT " Must be greater than or equal to zero!": GOTO 3700

      CLS : PRINT
3800  PRINT " What is the peak TERMINAL VOLTAGE?"
      PRINT " Enter 100 if you don't know."
      PRINT
      INPUT TermV
      IF TermV <= 0 THEN PRINT " Must be greater than zero!": GOTO 3800

      CLS : PRINT
      PRINT " Grid Size = "; Acc * 120
      PRINT " Scale Multiplier = "; Mult
      PRINT " Ceiling Height = "; Ceil
      PRINT " Wall Distance = "; Wall
      PRINT " Primary Height (Inner Turn) = "; PriH
      PRINT " Primary Inner Diameter = "; PriD0
      PRINT " Primary Outer Diameter = "; PriD
      PRINT " Primary Height (Outer Turn) = "; PriH1
      PRINT " Secondary Base Height ="; SecB
      PRINT " Secondary Diameter = "; SecD
      PRINT " Secondary Length = "; SecL
      PRINT " Secondary Inductance = "; SecInd
      PRINT " Terminal-1 Height (center) = "; TermH
      PRINT " Terminal-1 Radius = "; TermR
      PRINT " Terminal-1 Diameter (torroid only)= "; TermD
      PRINT " Terminal-2 Height (center) = "; TermH2
      PRINT " Terminal-2 Radius = "; TermR2
      PRINT " Terminal-2 Diameter (torroid only)= "; TermD2
      PRINT " Terminal Voltage = "; TermV
      PRINT
      PRINT " Some numbers may be changed to default values."
      PRINT
      PRINT " Is this correct?"
      INPUT " y/n "; q$
      IF q$ = "N" OR q$ = "n" GOTO 1000

' Factor Acf so coil diameter falls directly on grid points
    
      AcF = AcF * INT(AcF * SecD / 2) / (AcF * SecD / 2)
     
' Store Inputs to File
      OPEN "output.txt" FOR OUTPUT AS #2
      PRINT #2, " E-TESLA  Version 3.0   Terry Fritz   1-29-1999  "
      PRINT #2, "  ": PRINT #2, TIME$, DATE$
      PRINT #2, "  "
      PRINT #2, " Grid Size = "; Acc * 120
      PRINT #2, " Scale Multiplier = "; Mult
      PRINT #2, " Ceiling Distance = "; Ceil: Ceil = INT(Ceil * AcF + .5)
      PRINT #2, " Wall Distance = "; Wall: Wall = INT(Wall * AcF + .5)
      PRINT #2, " Primary Height (Inner Turn) = "; PriH: PriH = INT(PriH * AcF + .5)
      PRINT #2, " Primary Inner Diameter = "; PriD0: PriD0 = INT(PriD0 * AcF + .5)
      PRINT #2, " Primary  Outer Diameter = "; PriD: PriD = INT(PriD * AcF + .5)
      PRINT #2, " Primary Height (Outer Turn) = "; PriH1: PriH1 = INT(PriH1 * AcF + .5)
      PRINT #2, " Secondary Base Height = "; SecB: SecB = INT(SecB * AcF + .5)
      PRINT #2, " Secondary Diameter = "; SecD: SecD = SecD * AcF
      PRINT #2, " Secondary Length = "; SecL: SecL = INT(SecL * AcF + .5)
      PRINT #2, " Secondary Inductance = "; SecInd
      PRINT #2, " Terminal-1 Height (center) = "; TermH: TermH = INT(TermH * AcF + .5)
      PRINT #2, " Terminal-1 Radius = "; TermR: TermR = TermR * AcF
      PRINT #2, " Terminal-1 Diameter (torroid only)= "; TermD: TermD = TermD * AcF
      PRINT #2, " Terminal-2 Height (center) = "; TermH2: TermH2 = INT(TermH2 * AcF + .5)
      PRINT #2, " Terminal-2 Radius = "; TermR2: TermR2 = TermR2 * AcF
      PRINT #2, " Terminal-2 Diameter (torroid only)= "; TermD2: TermD2 = TermD2 * AcF
      PRINT #2, " Terminal Voltage = "; TermV
      PRINT #2, "  "
      PRINT #2, "  "
      PRINT #2, " Pass"; TAB(10); "Ccalc (pF)"; TAB(30); "% Change"; TAB(50); "Fcalc (Hz)"
      PRINT #2, "  "
      CLOSE #2

      ' Htot is the total effective height of the coil
      ' g0 is the radius of the Gaussian sphere
      Htot = TermR + TermD / 2 + TermH
      IF Htot < TermH2 + TermR2 + TermD2 / 2 THEN Htot = TermH2 + TermR2 + TermD2 / 2
      g0 = INT(((Htot - SecB) + (Ceil - SecB)) / 4 + .5)
                                    
      CLS : PRINT
      PRINT " Calculations are now under way!"
      PRINT
      PRINT " Press Ctrl+Pause to return to BASIC. (F5 to resume)"
      PRINT
      PRINT
      PRINT " Pass"; TAB(10); "Ccalc (pF)"; TAB(30); "% Change"; TAB(50); "Fcalc (Hz)"
      PRINT
     
      DIM a(Wall + 1, Ceil + 1), p(Wall + 1, Ceil + 1)
      ' BASIC automatically sets these arrays to zero

' Set initial conditions and lock bondaries
      GOSUB 51000

30000 'Main Loop
      FOR zz = 1 TO 100
      FOR y = 1 TO Ceil - 1: FOR x = 1 TO Wall - 1
      IF p(x, y) GOTO 31000
      a(x, y) = .25 * (a(x + 1, y) + a(x - 1, y) + a(x, y + 1) + a(x, y - 1))
      a(x, y) = a(x, y) + (a(x + 1, y) - a(x - 1, y)) / (8 * x)
31000 NEXT x, y
      FOR z = 1 TO Ceil - 1
      IF p(0, z) GOTO 32000
      a(0, z) = (4 * a(1, z) + a(0, z - 1) + a(0, z + 1)) / 6
32000 NEXT z
      NEXT zz
      count = count + 100

' Calculate E-field by summing all flux through a sphere

      Et = 0
      y0 = g0 + PriH + .5
      FOR ts = -1.57 TO 1.57 STEP .001
      'CircFac is the area of the section of the sphere
      CircFac = (SIN(ts + .0005) - SIN(ts - .0005))
      xr = COS(ts) * (g0 - .5)
      yr = SIN(ts) * (g0 - .5) + y0
      GOSUB 50000
      vs1 = vr
      xr = COS(ts) * (g0 + .5)
      yr = SIN(ts) * (g0 + .5) + y0
      GOSUB 50000
      vs2 = vr
      vs = (vs1 - vs2)
      Et = Et + vs * CircFac
      NEXT ts

      Et = Et * 5.4 * g0 * g0 / AcF / 39.37
     
' Calculate the capacitance
     
      Qt = 8.8542E-12 * Et * Mult
      Cap = Qt / TermV
      Erp = (Cap - Cap1) / Cap1 * 100
      IF Cap < 0 THEN Fo = 0: GOTO 45000
      Fo = 1 / (2 * pi * SQR(SecInd * Cap))

45000 PRINT count; TAB(10);
      PRINT "C="; Cap * 1E+12; TAB(30); Erp; TAB(50); "F ="; Fo
      OPEN "output.txt" FOR APPEND AS #3
      PRINT #3, count; TAB(10);
      PRINT #3, "C="; Cap * 1E+12; TAB(30); Erp; TAB(50); "F ="; Fo
      CLOSE #3
      Cap1 = Cap

' Output the whole voltage grid space to a file for Excel
      OPEN "Voutput.xls" FOR OUTPUT AS #1
      FOR y = Ceil TO 0 STEP -1: FOR x = Wall TO 0 STEP -1
      PRINT #1, a(x, y);
      NEXT x
      ' Mirror left side
      FOR x = 1 TO Wall
      PRINT #1, a(x, y);
      NEXT x
      PRINT #1, "  "
      NEXT y
      CLOSE #1

' Output the whole stress grid space to a file for Excel
      OPEN "Soutput.xls" FOR OUTPUT AS #1
      FOR y = Ceil - 1 TO 1 STEP -1: FOR x = Wall - 1 TO 1 STEP -1
      Str = 0
      IF ABS(a(x, y) - a(x - 1, y)) > Str THEN Str = ABS(a(x, y) - a(x - 1, y))
      IF ABS(a(x, y) - a(x + 1, y)) > Str THEN Str = ABS(a(x, y) - a(x + 1, y))
      IF ABS(a(x, y) - a(x, y + 1)) > Str THEN Str = ABS(a(x, y) - a(x, y + 1))
      IF ABS(a(x, y) - a(x, y - 1)) > Str THEN Str = ABS(a(x, y) - a(x, y - 1))
      PRINT #1, Str;
      NEXT x
      ' Mirror left side
      FOR x = 1 TO Wall - 1
      Str = 0
      IF ABS(a(x, y) - a(x - 1, y)) > Str THEN Str = ABS(a(x, y) - a(x - 1, y))
      IF ABS(a(x, y) - a(x + 1, y)) > Str THEN Str = ABS(a(x, y) - a(x + 1, y))
      IF ABS(a(x, y) - a(x, y + 1)) > Str THEN Str = ABS(a(x, y) - a(x, y + 1))
      IF ABS(a(x, y) - a(x, y - 1)) > Str THEN Str = ABS(a(x, y) - a(x, y - 1))
     
      PRINT #1, Str;
      NEXT x
      PRINT #1, "  "
      NEXT y
      CLOSE #1
     
49000 GOTO 30000

49999 END
     
50000 REM RESOLVER (interpolates any fractional x,y point)
      Vr1 = a(INT(xr), INT(yr))
      Vr2 = a(INT(xr + 1), INT(yr))
      Vr3 = a(INT(xr), INT(yr + 1))
      Vr4 = a(INT(xr + 1), INT(yr + 1))
      Dr1 = SQR(INT(xr) ^ 2 + (INT(yr) - y0) ^ 2)
      Dr2 = SQR(INT(xr + 1) ^ 2 + (INT(yr) - y0) ^ 2)
      Dr3 = SQR(INT(xr) ^ 2 + (INT(yr + 1) - y0) ^ 2)
      Dr4 = SQR(INT(xr + 1) ^ 2 + (INT(yr + 1) - y0) ^ 2)

      Qtt = (Vr1 / Dr1 + Vr2 / Dr2 + Vr3 / Dr3 + Vr4 / Dr4) / 4
    
      vr = Qtt * SQR(xr ^ 2 + (yr - y0) ^ 2)

RETURN

51000 ' Set Initial Conditions and Lock Boundaries

' Primary Coil
      FOR z = PriD0 / 2 TO PriD / 2
      q = PriD / 2 - SecD / 2 - 1
      a(z, PriH + (PriH1 - PriH) * (z - SecD / 2 - 1) / q) = 0
      p(z, PriH + (PriH1 - PriH) * (z - SecD / 2 - 1) / q) = 1
      NEXT z

' Voltage Along Secondary
      FOR z = SecB TO SecB + SecL
      'boosted voltage for radius not on the grid
      a(SecD / 2, z) = TermV * SIN(pi / 2 * (z - SecB) / SecL)
      p(SecD / 2, z) = 1
      NEXT z
 
' Top of Secondary to Top Terminal-1
      FOR z = SecB + SecL TO TermH
      a(0, z) = TermV
      p(0, z) = 1
      NEXT z

' Out to Radius-1
      FOR z = 0 TO TermD / 2
      a(z, TermH) = TermV
      p(z, TermH) = 1
      NEXT z

' Terminal-1 Radius
      FOR z = -pi / 2 TO pi / 2 STEP .0785
      x = TermD / 2 + TermR * COS(z): y = TermH + TermR * SIN(z)
      a(x, y) = TermV
      p(x, y) = 1
      IF TermD / 2 + 1 < TermR THEN 51100
      x = TermD / 2 - TermR * COS(z)
      a(x, y) = TermV
      p(x, y) = 1
51100 NEXT z
      a(TermD / 2, TermH + TermR) = TermV
      a(TermD / 2, TermH - TermR) = TermV
      p(TermD / 2, TermH + TermR) = 1
      p(TermD / 2, TermH - TermR) = 1
   
' Top of Terminal-1 to Top of Terminal-2
      FOR z = TermH TO TermH2
      a(0, z) = TermV
      p(0, z) = 1
      NEXT z

' Out to Radius-1
      FOR z = 0 TO TermD2 / 2
      a(z, TermH2) = TermV
      p(z, TermH2) = 1
      NEXT z

' Terminal-2 Radius
      FOR z = -pi / 2 TO pi / 2 STEP .0785
      x = TermD2 / 2 + TermR2 * COS(z): y = TermH2 + TermR2 * SIN(z)
      a(x, y) = TermV
      p(x, y) = 1
      IF TermD2 / 2 + 1 < TermR2 THEN 51200
      x = TermD2 / 2 - TermR2 * COS(z)
      a(x, y) = TermV
      p(x, y) = 1
51200 NEXT z
      a(TermD2 / 2, TermH2 + TermR2) = TermV
      a(TermD2 / 2, TermH2 - TermR2) = TermV
      p(TermD2 / 2, TermH2 + TermR2) = 1
      p(TermD2 / 2, TermH2 - TermR2) = 1
 
RETURN