Toroidal Coil Field

@Karan Chatrath has posted a detailed solution. I will post my code, because it shows how to derive the path displacement vector by taking something like numerical difference quotients, thus avoiding messy analytical work. In the middle of the toroid coil, the field matches the standard formula. At the origin, the field largely cancels, but not entirely. I ran the simulation with a wide range of spatial resolutions to ensure that this lack of cancellation was physically real, and not just a numerical artifact.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118

import math

# Scan resolution

N = 10**7

######################################################

# Constants

R = 4.0
r = 1.0

u0 = 1.0
I = 1.0
Bmult = u0*I/(4.0*math.pi)

thetaf = 2.0*math.pi

dtheta = thetaf/N

######################################################

# Test point

Nx = 0.0
Ny = 0.0
Nz = 0.0

######################################################
######################################################

# Determine B field using Biot Savart

Bx = 0.0
By = 0.0
Bz = 0.0

theta = 0.0

while theta <= thetaf - dtheta:         

    # coordinates

    theta1 = theta
    theta2 = theta + dtheta

    x1 = (R + r*math.cos(1000.0*theta1))*math.cos(theta1)
    y1 = (R + r*math.cos(1000.0*theta1))*math.sin(theta1)             
    z1 = r*math.sin(1000.0*theta1)

    x2 = (R + r*math.cos(1000.0*theta2))*math.cos(theta2)
    y2 = (R + r*math.cos(1000.0*theta2))*math.sin(theta2)             
    z2 = r*math.sin(1000.0*theta2)

    dx = x2 - x1
    dy = y2 - y1        # path displacement (components of dL)
    dz = z2 - z1

    Dx = x2 - Nx                  # vector to test point
    Dy = y2 - Ny
    Dz = z2 - Nz

    D = math.sqrt(Dx**2.0 + Dy**2.0 + Dz**2.0)  # rmag 

    crossx = dy*Dz - dz*Dy      # dL cross r from Biot Savart
    crossy = -(dx*Dz - dz*Dx)
    crossz = dx*Dy - dy*Dx

    dBx = Bmult*crossx/(D**3.0)  # infinitesimal B field contributions
    dBy = Bmult*crossy/(D**3.0)
    dBz = Bmult*crossz/(D**3.0)

    Bx = Bx + dBx
    By = By + dBy
    Bz = Bz + dBz

    theta = theta + dtheta

######################################################
######################################################

# Print results

print N
print ""
print Bx
print By
print Bz
print ""
B = math.sqrt(Bx**2.0 + By**2.0 + Bz**2.0)

print B

######################################################
######################################################

# Result for (R,0,0)

#10000000

#1.17915938536e-05
#39.7887330646
#-0.0587505108876

#39.7887764389

######################################################

# Result for (0,0,0)

#10000000

#5.13950703891e-14
#-3.6167629907e-09
#-0.12302337116

#0.12302337116

Similar solution outline as the solenoid exercise. Specific details are left out.

Consider a point on the solenoid as such:

$\vec{r}_p = x(\theta)\hat{i} + y(\theta)\hat{j} + z(\theta)\hat{k}=\left(R+r\cos(1000\theta)\right)\cos{\theta} \hat{i} + \left(R+r\cos(1000\theta)\right)\sin{\theta} \hat{j} + r\sin(1000\theta)\hat{k}$

An arc length element tangential to the point is:

$d\vec{r}_p = \left(\frac{dx}{d\theta}\hat{i} + \frac{dy}{d\theta} \hat{j} + \frac{dz}{d\theta}\hat{k}\right) \ d\theta$

Consider the point of interest at which field is to be computed to be:

$\vec{r}_{c1} = R \hat{i} + 0 \hat{j} +0\hat{k}$ $\vec{r}_{c2} = 0 \hat{i} + 0 \hat{j} +0\hat{k}$

In general, the point of interest will be denoted by $\vec{r}_{c}$ . The vector joining a point on the solenoid and the point of interest, and directed towards the point of interest, is:

$\vec{r} = \vec{r}_c - \vec{r}_p$

Now, Biot-Savart Law is applied:

$d\vec{B} = \frac{\mu_o I}{4 \pi}\left(\frac{d\vec{r}_p \times \vec{r}}{\lvert \vec{r} \rvert^3}\right)= f_x(\theta) \hat{i} + f_y(\theta) \hat{j} +f_z(\theta) \hat{k}$

$d\vec{B}$ has three components each of which is a function of $\theta$ only. This gives three integrals:

$B_x = \int_{0}^{2000 \pi} f_x(\theta) \ d\theta$
$B_y = \int_{0}^{2000 \pi} f_y(\theta) \ d\theta$ $B_z = \int_{0}^{2000 \pi} f_z(\theta) \ d\theta$

$B_{1,2} = \lvert \vec{B} \rvert = \sqrt{B_x^2 + B_y^2 + B_z^2}$

I have written a script of code to compute the answer for each case which turns out to be:

$\boxed{\frac{B_1}{B_2} \approx 323.4}$

To numerically integrate to obtain the fields, it was required to set the numerical resolution quite high. In fact, visualising the toroid in a 3D plot itself required a high resolution. Step sizes such as even $\pi/1000$ were not adequate. Also, I did not get this right in my initial attempts. I made algebraic errors while computing derivatives and in other intermediate steps.