Private Sub Sect6

Private Sub Sect6_Click()

'*******************************************************************

'This program computes the densities of particles in frame K',

'at the single instant t'=0. It does this both for the case where the

'particles are at rest in K and the case where they move in K.

'*******************************************************************

Const pi As Double = 3.141592654

Const R As Double = 1 'Radius of circular contour in frame K

Const N As Double = 36000 'Number of particles

Const c As Double = 300000000# 'Speed of light

Const v As Double = 0.8 * c 'Speed of frame K1 relative to K

Const u As Double = v 'Particle speed in K

Const omega As Double = u / R 'Angular rate in K

Const dphi As Double = 2 * pi / 360 'in radians

Dim deltat As Double

Dim t, oldt, tp, oldtp As Double

Dim x, oldxp As Double

Dim xp(N), xpResting(N) As Double

Dim y, oldyp As Double

Dim yp(N), ypResting(N) As Double

Dim theta As Double

Dim gamma As Double

Dim index, index1 As Long

Dim phi(360), xmin(360), xmax(360), ymin(360), ymax(360), rho(360), rhoResting(360) As Double

gamma = 1 / Sqr(1 - v ^ 2 / c ^ 2) 'Used in transformations

For index = 0 To N - 1

theta = 2 * pi * index / N

x = R * Cos(theta)

y = R * Sin(theta)

'K'-coordinates are analytic when particles are at rest in K.

xpResting(index) = x / gamma

ypResting(index) = y

Next index

'In order to save compute time, the creeping up strategy in this case is to

'start out with a coarse deltat, and to keep decreasing it until t' is

'acceptably close to zero.

For index = 0 To N - 1

deltat = 4 * pi * R / (N * u) 'Start out with large deltat

theta = 2 * pi * index / N 'Angle with positive x-axis

x = R * Cos(theta) 'x-coordinate

xp(index) = gamma * x 'x1 coordinate (length contracted x)

y = R * Sin(theta) 'y-coordinate

yp(index) = y 'y' coordinate (no length contraction)

t = 0 'Reading of clock coinciding with particle in K

tp = -gamma * v * x / c ^ 2 'Reading of K' clock

deltat = deltat / 2

oldxp = xp(index) 'Used to retrench when t1 changes sign

oldyp = yp(index)

oldtp = tp

oldt = t

If tp <= 0 Then

t = t + deltat

Else

t = t - deltat

End If

theta = 2 * pi * index / N + omega * t

x = R * Cos(theta) 'Adjusted coordinates in K at new time t

y = R * Sin(theta)

xp(index) = gamma * (x - v * t) 'Adjusted coordinates in K'

yp(index) = y

tp = gamma * (t - v * x / c ^ 2)

Loop Until Sgn(tp) <> Sgn(oldtp) 'When sign of t' changes, scale back or quit

If Abs(tp) < 1E-23 Then Exit Do

xp(index) = oldxp

yp(index) = oldyp

tp = oldtp

t = oldt

Loop

Next index

'phi is the angle (in whole degrees) with the positive x'-axis. xmin, xmax,

'ymin and ymax are boundaries between successive angles phi, for

'tallying number of particles and hence density.

For index1 = 0 To 359

phi(index1) = index1 * dphi 'phi initially in radians, for computing trig functions

Next index1

xmax(0) = 1 / gamma * R

xmin(0) = Sin(Atn(1 / gamma * Tan(phi(1)))) * R / Tan(phi(1))

For index1 = 1 To 88

xmax(index1) = xmin(index1 - 1)

xmin(index1) = Sin(Atn(1 / gamma * Tan(phi(index1 + 1)))) * R / Tan(phi(index1 + 1))

Next index1

xmax(89) = xmin(88)

xmin(89) = 0

xmax(90) = 0

xmin(90) = -Sin(Atn(1 / gamma * Tan(phi(91)))) * R / Tan(phi(91))

For index1 = 91 To 178

xmax(index1) = xmin(index1 - 1)

xmin(index1) = -Sin(Atn(1 / gamma * Tan(phi(index1 + 1)))) * R / Tan(phi(index1 + 1))

Next index1

xmax(179) = xmin(178)

xmin(179) = -1 / gamma * R

xmin(180) = -1 / gamma * R

xmax(180) = -Sin(Atn(1 / gamma * Tan(phi(181)))) * R / Tan(phi(181))

For index1 = 181 To 268

xmin(index1) = xmax(index1 - 1)

xmax(index1) = -Sin(Atn(1 / gamma * Tan(phi(index1 + 1)))) * R / Tan(phi(index1 + 1))

Next index1

xmin(269) = xmax(268)

xmax(269) = 0

xmin(270) = 0

xmax(270) = Sin(Atn(1 / gamma * Tan(phi(271)))) * R / Tan(phi(271))

For index1 = 271 To 358

xmin(index1) = xmax(index1 - 1)

xmax(index1) = Sin(Atn(1 / gamma * Tan(phi(index1 + 1)))) * R / Tan(phi(index1 + 1))

Next index1

xmin(359) = xmax(358)

xmax(359) = 0

For index1 = 0 To 88

ymin(index1) = Sin(Atn(1 / gamma * Tan(phi(index1)))) * R

ymax(index1) = Sin(Atn(1 / gamma * Tan(phi(index1 + 1)))) * R

Next index1

ymin(89) = ymax(88)

ymax(89) = R

ymin(90) = -Sin(Atn(1 / gamma * Tan(phi(91)))) * R

ymax(90) = R

For index1 = 91 To 178

ymax(index1) = -Sin(Atn(1 / gamma * Tan(phi(index1)))) * R

ymin(index1) = -Sin(Atn(1 / gamma * Tan(phi(index1 + 1)))) * R

Next index1

ymax(179) = ymin(178)

ymin(179) = 0

ymax(180) = 0

ymin(180) = -Sin(Atn(1 / gamma * Tan(phi(180)))) * R

For index1 = 181 To 268

ymax(index1) = -Sin(Atn(1 / gamma * Tan(phi(index1)))) * R

ymin(index1) = -Sin(Atn(1 / gamma * Tan(phi(index1 + 1)))) * R

Next index1

ymax(269) = ymin(268)

ymin(269) = -R

ymax(270) = ymax(269)

ymin(270) = -R

For index1 = 271 To 358

ymin(index1) = Sin(Atn(1 / gamma * Tan(phi(index1)))) * R

ymax(index1) = Sin(Atn(1 / gamma * Tan(phi(index1 + 1)))) * R

Next index1

ymin(359) = ymax(358)

ymax(359) = 0

For index1 = 0 To 359

phi(index1) = index1 'Express phi now in whole degrees

rho(index1) = 0 'Initialize density to zero

rhoResting(index1) = 0

For index = 0 To N - 1

If xp(index) >= xmin(index1) And xp(index) < xmax(index1) And yp(index) >= ymin(index1) And yp(index) < ymax(index1) Then

rho(index1) = rho(index1) + 1

End If

If xpResting(index) >= xmin(index1) And xpResting(index) < xmax(index1) And ypResting(index) >= ymin(index1) And ypResting(index) < ymax(index1) Then

rhoResting(index1) = rhoResting(index1) + 1

End If

Next index

Next index1

Open "c:\WINMCAD\Physics\PolarizationRho.prn" For Output As #1

For index1 = 0 To 359

Write #1, phi(index1), rho(index1), rhoResting(index1)

Next index1

MsgBox ("Ready for plotting")

End Sub