Option Explicit

Private Sub cmdComputeF_Click()

'******

'This procedure computes the magnetic force on Loop B when its center lies

'to the right of loop A and both loop planes coincide with the xz-plane.

'TA is short for theta Loop A, TB is short for theta Loop B.

'******

Const pi As Double = 3.131492654

Const Radius As Double = 1 'common loop radius

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

'Select the appropriate values for dmin and dmax, and comment out the

'other values. dmin and dmax refer to the minimum and maximum values

'of d, the separation of the current loop centers.

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

'Select the following values for Figs. 2_1a and 2_2a.

'Const dmin As Double = 2.01 * Radius

'Const dmax As Double = 2.1 * Radius

'Select the following values for Figs. 2_3a and b.

'Const dmin As Double = -1.75 * Radius

'Const dmax As Double = 1.75 * Radius

'Select the following values for Fig. 2_3c.

Const dmin As Double = 0.25 * Radius

Const dmax As Double = 3.75 * Radius

Const nd As Long = 100 'number of divisions of d (loop center separation)

Const dd As Double = (dmax - dmin) / nd 'd increment

Const nTA As Double = 180 'number of TA values

Const dTA As Double = 2 * pi / nTA 'TA increment

Const qB As Double = 1 'circulating charge of Loop B

Const nTB As Double = 500 'number of qB increments

Const dTB As Double = 2 * pi / nTB 'TB increment

Const dqB As Double = qB / nTB 'size of qB increment

Const c As Double = 300000000# 'speed of light

Const eps0 As Double = 0.00000000000885 'permittivity constant

Const amps As Double = 1 'loop current

Const lambda As Double = qB / (2 * pi * Radius) 'line density, circ chge, B

Const v As Double = amps / lambda 'speed of circ chge, B

Const C1 As Double = amps / (4 * pi * eps0 * c ^ 2) 'saves cpu time

Dim xB, yB, zB As Double 'point on Loop B

Dim xA, yA, zA As Double 'point on LoopA

Dim TB, TA As Double 'values of TB and TA

Dim dlxA, dlyA, dlzA As Double 'dl increment components on LoopA

Dim rx, ry, rz, r As Double 'displacement from point on LoopA to pt on LoopB

Dim Bx, By, Bz As Double 'magnetic fld of LoopA at point on LoopB

Dim indexTB, indexTA As Long 'loop counters for TB and TA

Dim indexd As Long 'loop counter for d

Dim dlxrx, dlxry, dlxrz As Double 'dl(LoopA) cross r components

Dim d(nd) As Double 'values of d, the separation of Loop A and Loop B centers

Dim vx, vy, vz As Double 'loop B circulating chge increment velocity

Dim dFx(nTB), dFy(nTB), dFz(nTB) As Double 'force increments on parts of Loop B

Dim Fx(nTB), Fy(nTB), Fz(nTB) As Double 'total force on Loop B increment

Dim NetFx(nd) As Double 'net force on Loop B output for plotting vs. d

'For each value of d (the separation of the loop centers)...

For indexd = 0 To nd - 1

d(indexd) = dmin + indexd * dd

Debug.Print indexd

'...and for each charge increment in Loop B...

For indexTB = 0 To nTB / 2

TB = indexTB * dTB

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

'Select the first value for when the dipole moments are parallel.

'Select the second value for when the dipole moments are antiparallel.

'vz = -v * Cos(TB)

vz = v * Cos(TB)

'...compute the B field of Loop A at the B increment position;

xB = d(indexd) + Radius * Cos(TB)

yB = 0

zB = -Radius * Sin(TB)

Bx = 0

By = 0

Bz = 0

dFx(indexTB) = 0

dFy(indexTB) = 0

dFz(indexTB) = 0

For indexTA = 0 To nTA - 1

TA = indexTA * dTA

xA = Radius * Cos(TA)

yA = 0

zA = -Radius * Sin(TA)

'At non-overlap points...

If Sqr((xB - xA) ^ 2 + (zB - zA) ^ 2) > dd Then

dlxA = -Radius * dTA * Sin(TA)

dlyA = 0

dlzA = -Radius * dTA * Cos(TA)

rx = xB - xA

ry = 0

rz = zB - zA

r = Sqr(rx ^ 2 + ry ^ 2 + rz ^ 2)

dlxrx = dlyA * rz - dlzA * ry

dlxry = dlzA * rx - dlxA * rz

dlxrz = dlxA * ry - dlyA * rx

Bx = 0

By = By + C1 * dlxry / r ^ 3

Bz = 0

dFx(indexTB) = dFx(indexTB) + dqB * (vy * Bz - vz * By)

dFy(indexTB) = dFy(indexTB) + dqB * (vz * Bx - vx * Bz)

dFz(indexTB) = dFz(indexTB) + dqB * (vx * By - vy * Bx)

Else

'At overlapping points...

dFx(indexTB) = 0

End If

Next indexTA

Fx(indexTB) = dFx(indexTB)

Fy(indexTB) = dFy(indexTB)

Fz(indexTB) = dFz(indexTB)

Next indexTB

'Set the total force for the symmetric Loop B chge increments.

For indexTB = 1 To nTB / 2 - 1

Fx(nTB - indexTB) = Fx(indexTB)

Next indexTB

'Now add up all the forces on the individual Loop B increments

'to obtain the net force on Loop B at separation d.

NetFx(indexd) = 0

For indexTB = 0 To nTB - 1

NetFx(indexd) = NetFx(indexd) + Fx(indexTB)

Next indexTB

Next indexd

'Output the values of Fx(d) for plotting purposes.

Open "C:\WINMCAD\Physics\Nucleons.PRN" For Output As #1

For indexd = 0 To nd - 1

Write #1, d(indexd), NetFx(indexd)

Next indexd

Close

MsgBox ("Ready for plotting Fx")

Stop

End Sub

Private Sub ComputeF1_Click()

'******

'This procedure computes the force on Loop B when its center lies above that

'of Loop A, and when their dipole moments are parallel or antiparallel. The

'planes of both loops are parallel to the xz-plane.

'******

Const pi As Double = 3.131492654

Const Radius As Double = 1

Const dmin As Double = 0.01 * Radius

Const dmax As Double = 0.1 * Radius

Const nd As Long = 100

Const dd As Double = (dmax - dmin) / nd

Const nTA As Double = 360

Const dTA As Double = 2 * pi / nTA

Const qB As Double = 1

Const nTB As Double = 500

Const dqB As Double = qB / nTB

Const c As Double = 300000000#

Const eps0 As Double = 0.00000000000885

Const amps As Double = 1

Const lambda As Double = qB / (2 * pi * Radius)

Const v As Double = amps / lambda

Const C1 As Double = amps / (4 * pi * eps0 * c ^ 2)

Dim xB, yB, zB As Double

Dim xA, yA, zA As Double

Dim TA As Double

Dim dlxA, dlyA, dlzA As Double

Dim rx, ry, rz, r As Double

Dim Bx, By, Bz As Double

Dim indexTA As Long

Dim indexd As Long

Dim dlxrx, dlxry, dlxrz As Double

Dim d(nd) As Double

Dim vx, vy, vz As Double

Dim dFx, dFy, dFz As Double

Dim NetFy(nd) As Double

For indexd = 0 To nd - 1

d(indexd) = dmin + indexd * dd

Debug.Print indexd

vx = 0

vy = 0

'vz = -v

vz = v

xB = Radius

yB = d(indexd)

zB = 0

Bx = 0

By = 0

Bz = 0

dFx = 0

dFy = 0

dFz = 0

For indexTA = 0 To nTA - 1

TA = indexTA * dTA

xA = Radius * Cos(TA)

yA = 0

zA = -Radius * Sin(TA)

dlxA = -Radius * dTA * Sin(TA)

dlyA = 0

dlzA = -Radius * dTA * Cos(TA)

rx = xB - xA

ry = yB

rz = zB - zA

r = Sqr(rx ^ 2 + ry ^ 2 + rz ^ 2)

dlxrx = dlyA * rz - dlzA * ry

dlxry = dlzA * rx - dlxA * rz

dlxrz = dlxA * ry - dlyA * rx

Bx = Bx + C1 * dlxrx / r ^ 3

By = By + C1 * dlxry / r ^ 3

Bz = Bz + C1 * dlxrz / r ^ 3

dFx = dFx + dqB * (vy * Bz - vz * By)

dFy = dFy + dqB * (vz * Bx - vx * Bz)

dFz = dFz + dqB * (vx * By - vy * Bx)

Next indexTA

NetFy(indexd) = nTA * dFy

Next indexd

Open "C:\WINMCAD\Physics\Nucleons.PRN" For Output As #1

For indexd = 0 To nd - 1

Write #1, d(indexd), NetFy(indexd)

Next indexd

Close

MsgBox ("Ready for plotting NetFy(d)")

Stop

End Sub