DECLARE SUB SHOCK (S!)
DECLARE SUB CIRCLE1 (xj!, yj!, cxj!, cyj!, R!, color1!)
DECLARE SUB CIRCLE2 (xj!, yj!, rj!, anglj!, color1!)
DECLARE FUNCTION gT! (T!)
DECLARE FUNCTION fT! (T!)
DIM SHARED x(20), y(20), R(20), color1(20), cx(20), cy(20), angl(20)
DIM SHARED pi, xdmax, ydmax, speed, nT, scale, c, scalec

'                     RADIATE.BAS
'      Revision 1.0   Original                               14 Apr 1999
'      Revision 1.1   Added angle for movement 3             15 Apr 1999
'      Revision 1.2   Added shock wave example               13 May 1999

'      For program description select:
'             <A HREF="progrm15.htm">Program 15: Field Radiation</A>
'      For program implementation guide lines select:
'             <A HREF="implemnt.htm">Implementation Details</A>

SCREEN 12
xdmax = 639: ydmax = 479       ' display maximum

PSET (xdmax, ydmax)
pi = 4 * ATN(1)
basespeed = 1.7
rrbase = 100         ' radius in movement = 3
deltax = 4.2         ' delta x particle
deltar = speed * deltax  ' delta r field

DO:
  COLOR 7
  LOCATE 1, 1
  PRINT "  0 = End "
  PRINT "  1 = MTW page 111 "
  PRINT "  2 = Origin centered"
  PRINT "  3 = 4 shock wave display "
  COLOR 15: INPUT "Enter theory "; Theory
  IF Theory = 0 THEN END
  IF Theory = 3 THEN GOTO prog1
  COLOR 7
  PRINT "  0 = End "
  PRINT "  1 = Straight Line"
  PRINT "  2 = 180 Degrees Reflection"
  PRINT "  3 = Slowly Bended"
  COLOR 15: INPUT "Enter particle movement "; movement
  IF movement = 0 THEN END
  IF movement = 3 THEN
    INPUT "Enter radius (Standard = 100)"; rr1
    IF rr1 = 0 THEN rr = rrbase ELSE rr = rr1
  END IF
    COLOR 7: PRINT "  Minimum 1.1, > 50 is infinite "
  COLOR 15: INPUT "Enter speed of field (Standard 1.7) "; speed1
  IF speed1 = 0 THEN speed = basespeed ELSE speed = speed1
  IF speed < 1.1 THEN speed = 1.1
  IF speed > 50 THEN speed = 50              ' Infinite

  deltar = speed * deltax      ' Delta radius r
  x = INT(xdmax / 2): y = INT(ydmax / 2)
  xmax = x
  ' Set color code for field
  FOR i = 0 TO 20
    color1(i) = 6
    IF i MOD 4 = 0 THEN color1(i) = 15
  NEXT i
  Time1 = TIMER
  FOR i = 0 TO 20
    CLS
    COLOR 7
    LOCATE 1, 45
      PRINT "Theory   = ";
      IF Theory = 1 THEN PRINT "MTW page 111 "
      IF Theory = 2 THEN PRINT "Origin centered"
    LOCATE 2, 45
      PRINT "Movement = ";
      IF movement = 1 THEN PRINT "Straight Line"
      IF movement = 2 THEN PRINT "180 Degrees Reflection"
      IF movement = 3 THEN PRINT "Slowly Bended"
    LOCATE 3, 55
      PRINT "Speed = "; speed
    LOCATE 4, 55
      IF movement = 3 THEN PRINT "Radius = "; rr
    FOR j = i TO 0 STEP -1
      IF j = i THEN           ' create new circle
        y(j) = y
        angl(j) = 0
        SELECT CASE movement
        CASE 1
          x(j) = x - deltax * j
        CASE 2
          IF j < 10 THEN
            xmax = xmax - deltax
            x(j) = xmax
          ELSE
            xmax = xmax + deltax
            x(j) = xmax
          END IF
        CASE 3
          x0 = xdmax / 2: y0 = ydmax / 2 - rr
          dalp = deltax * 360 / (2 * pi * rr)
          alint = 90 - 10 * dalp
          'PRINT dalp: END
          alp = alint + j * dalp
          x(j) = x0 + rr * COS(alp * pi / 180)
          y(j) = y0 + rr * SIN(alp * pi / 180)
          angl(j) = alp
          ' PRINT j; alp; x(j); y(j)
        END SELECT
        R(j) = 0
        cx(j) = x(j)
        cy(j) = y(j)
      ELSE
        R(j) = R(j) + deltar
        SELECT CASE Theory
        CASE 1
          SELECT CASE movement
          CASE 1
            cx(j) = cx(j) - deltax
          CASE 2
            IF j < 10 THEN
              cx(j) = cx(j) - deltax
            ELSE
              cx(j) = cx(j) + deltax
            END IF
          CASE 3
            cx(j) = cx(j) - deltax
          END SELECT
        CASE 2
        END SELECT
      END IF
      FOR k = 0 TO j
        PSET (x(j), ydmax - y(j)), 2
      NEXT k
      SELECT CASE Theory
      CASE 1
        IF speed <> 50 OR (speed = 50 AND j = i - 1) THEN
          CIRCLE1 x(j), y(j), cx(j), cy(j), R(j), color1(j)
        END IF
      CASE 2
        IF speed <> 50 OR (speed = 50 AND j = i - 1) THEN
          CIRCLE2 x(j), y(j), R(j), angl(j), color1(j)
        END IF
      END SELECT
    NEXT j
    DO:  LOOP UNTIL INT(TIMER) <> INT(Time1)
    Time1 = TIMER
  NEXT i
LOOP
END

prog1:

' Program by Chris Hillman

'% Define trajectory
'f(T) = .5 * EXP(-(T - 1) ^ 2) * COS(3 * T)
'g(T) = .5 * EXP(-(T - 3) ^ 2) * SIN(3 * T)

'% Set time interval between wavefronts
COLOR 7
nT = .075
scale = 60
cbase = 1        'speed of propagation
   
INPUT "Enter speed of propagation (Standard = 1)"; c1
IF c1 = 0 THEN c = cbase ELSE c = c1
scalec = scale / c

'% Plot the speed of trajectory
  'showspeed1 = SQR(fT(T) ^ 2 + g(T) ^ 2)  ' {t,0,4}]
  ceiling = 5 / nT
    PRINT TAB(3); "nT"; TAB(12); "fT(n * T)"; TAB(27); "gT(n * T)"; TAB(45); "R"; TAB(60); "vx"
  FOR n = 0 TO ceiling
    R = SQR(fT(n * nT) ^ 2 + gT(n * nT) ^ 2)
    ' Select distance for speed calculation
    IF n * nT <= 2 THEN
       Dist0 = fT((n - 1) * nT): Dist1 = fT(n * nT)
    ELSE
       Dist0 = gT((n - 1) * nT): Dist1 = gT(n * nT)
    END IF
    PRINT n * nT; TAB(10); fT(n * nT); TAB(25); gT(n * nT); TAB(40); R; TAB(55); (Dist1 - Dist0) / nT
    IF n = 26 OR n = 52 THEN
      INPUT A$
      CLS
      PRINT TAB(3); "nT"; TAB(12); "fT(n * T)"; TAB(27); "gT(n * T)"; TAB(45); "R"; TAB(60); "vy"
    END IF
  NEXT n
  INPUT A$

'% Make the wavefronts for time S
  'wavefronts(S) = Table(Graphics(Thickness(.001))):
  'CIRCLE (f(n * T), g(n * T)), S - n * T
  '{n,0,Ceiling[S/T]}]

'% Draw picture of wavefronts at time S
'shock(S) = Show(wavefronts(S))
                  ' PlotRange -> {{-6,6},{-6,6}},
                  ' AspectRatio-> Automatic];

'% Draw four pictures
'shockpics = Show(GraphicsArray(shock(1), shock(2), shock(4), shock(3)))

FOR n = 1 TO 4
  CLS
  SHOCK n
  INPUT A$
NEXT n

SUB CIRCLE1 (xj, yj, cxj, cyj, R, color1)
COLOR color1

FOR angle = 0 TO 360 STEP 30
  angle1 = angle * pi / 180
  dx = xj - cxj
  IF angle <> 90 AND angle <> 270 THEN
    tga = TAN(angle1)
    A = tga * tga + 1
    b = 2 * dx * tga * tga
    c = tga * tga * dx * dx - R * R
    a1 = (-b + SQR(b * b - 4 * A * c)) / (2 * A)
    a2 = (-b - SQR(b * b - 4 * A * c)) / (2 * A)
  ELSE
    a1 = -dx: a2 = -dx
  END IF
  IF angle <= 90 OR angle >= 270 THEN a3 = a1 ELSE a3 = a2
  b = SQR(R * R - a3 * a3)             ' vertical offset
  IF angle <= 180 THEN y1 = cyj + b:  ELSE y1 = cyj - b
  xk = dx + a3                         ' horizontal offset
  x1 = cxj + xk
  IF speed <> 50 THEN
    PSET (x1, ydmax - y1), color1
  ELSE
    x1 = xj + R * COS(angle1)
    y1 = yj + R * SIN(angle1)
    LINE (xj, ydmax - yj)-(x1, ydmax - y1), color1
  END IF
NEXT angle
'INPUT a$
END SUB

SUB CIRCLE2 (xj, yj, rj, anglj, color1)
COLOR color1
FOR angle = 0 TO 360 STEP 30
  angle1 = (anglj + angle) * pi / 180
  x1 = xj + rj * COS(angle1)
  y1 = yj + rj * SIN(angle1)
  IF speed <> 50 THEN
    PSET (x1, ydmax - y1), color1
  ELSE
    LINE (xj, ydmax - yj)-(x1, ydmax - y1), color1
  END IF
NEXT angle
END SUB

FUNCTION fT (T)
   fT = .5 * EXP(-(T - 1) ^ 2) * COS(3 * T)
END FUNCTION

FUNCTION gT (T)
gT = .5 * EXP(-(T - 3) ^ 2) * SIN(3 * T)
END FUNCTION

SUB SHOCK (S)
  ceiling = INT((S - .001) / nT)
  PRINT "Shock "; S; TAB(15); "Speed "; c
  FOR n = 0 TO ceiling
    IF n = ceiling THEN COLOR 15 ELSE COLOR 7
    CIRCLE (xdmax / 2 + scalec * fT(n * nT), ydmax / 2 - scalec * gT(n * nT)), scale * (S - n * nT)
  NEXT n
END SUB