from __future__ import print_function, division
from __future__ import division
from visual import *

# Faraday's Law: current in coil increases, then stays constant, then decreases, etc.
# Click to advance.

# Ruth Chabay 2007-08-07
# Revised Bruce Sherwood 2007-11-10

# Choose whether to show a ring or a disk through which passes magnetic flux
# Choose whether to show B at many places, not just in flux region
# Magenta arrow represents the vector dB_vec/dt
# Black arrow represents direction and magnitude of conventional current in coil

showring = True # show ring if True, disk if False
showallB = True # if True, show B at many places, not just within circle

scene.height = 800
scene.width = 800
scene.x = scene.y = 0
scene.background = color.white
scene.lights = [0.7*norm(vector(1,0,.5)), 0.7*norm(vector(-1,0.5,0.5))]
Bcolor = (0,.5,.5)
Bcolorother = (.8,1,1)
sw = 0.01
swother = 0.005

def Bfield(source,obsloc):
    r = source-obsloc
    return kmag*(3*dot(mu,norm(r))*norm(r) - mu)/mag(r)**3

def showB():
    for arr in Bother:
        arr.axis=Bscale*Bfield(coil.pos, arr.pos)

xhat = vector(1,0,0)

Rdisk = 0.3
f = cos(pi/4.)
rcoil = 0.12
Lcoil = 0.1
dpole = 0.03

coil = frame(pos=(0,0,0))
loops = helix(frame=coil, pos=(-Lcoil/2,0,0), axis=(Lcoil,0,0), color=color.orange,
             radius=rcoil, coils=9.6)
ampbox = box(frame=coil, pos=(0,-2.5*rcoil,0),
                 size=(0.2,0.1,0.1), color=(0.7,0.7,0.7))
plus = cylinder(frame=coil, pos=ampbox.pos+vector(0.05,0,0.05),
                axis=(0,0,0.01), radius=0.01, color=color.red)
minus = cylinder(frame=coil, pos=ampbox.pos+vector(-0.05,0,0.05),
                axis=(0,0,0.01), radius=0.01, color=color.black)
ammeter = label(frame=coil, pos=(0.03,-.13,0.05), text="I = %0.1f A" % 0., color=color.black)
current = arrow(frame=coil, pos=(0,0,rcoil), axis=(0,0,0), fixedwidth=1,
                shaftwidth=0.01, color=color.black)
pluslead = curve(frame=coil,
        pos=[(0.05,-0.07,-0.1), (.13,-.1,.05), (.13,-.2,.05), plus.pos],
        color=color.red, radius=0.005)
minuslead = curve(frame=coil,
        pos=[(loops.pos.x,0,rcoil), (loops.pos.x-.06,-.1,.08), (loops.pos.x-.06,-.2,.08), minus.pos],
        color=color.white, radius=0.005)

if showring:
    surface = ring(pos=vector(0.4,0,0), radius=Rdisk, thickness=0.001)
    surface.center = surface.pos
else:
    surface = cylinder(pos=vector(0.4,0,0), radius=Rdisk, axis=(0.002,0,0), color=(0.9,0.9,0.9))
    surface.center = surface.pos+surface.axis/2
    
deltax = Lcoil/5
kmag = 1e-7
mu = vector(0,0,0) # magnetic dipole moment of coil

#B = Bfield(coil.pos, surface.pos) # typical 2e-6 tesla

Bscale = 0.15*Rdisk/2e-6
Escale = 0.3*Rdisk/2e-7
xmax = 0.4*Rdisk

#E on perimeter of surface
Earr=[]
for theta in arange(0,2*pi,pi/6):
    a=arrow(pos=(surface.center.x, surface.radius*cos(theta),surface.radius*sin(theta)),
                           axis=(0,0,0), color=color.orange, shaftwidth=0.01)
    a.vv = norm(a.pos - surface.pos)
    Earr.append(a)

Bsurface = []  ## arrows on surface at which to calculate field, flux
dR = 0.2*Rdisk
for y in arange(-0.8*Rdisk,0.9*Rdisk,dR):
    a = arrow(pos=(surface.center.x, y, 0),axis=(0,0,0),
              color=Bcolor, fixedwidth=1, shaftwidth=sw)
    Bsurface.append(a)

## locations at which to display magnetic field around coil
Bother = []
dtheta = pi/6
phi = pi/4
for theta in arange(dtheta, pi-dtheta/2, dtheta):
    x = Rdisk*cos(theta)
    y = Rdisk*sin(theta)
    z = 0
    Bother.append( arrow(pos=(x,y,z), axis=(0,0,0)) )
    Bother.append( arrow(pos=(x,-y,z), axis=(0,0,0)) )
    Bother.append( arrow(pos=(x,z,y), axis=(0,0,0)) )
    Bother.append( arrow(pos=(x,z-y), axis=(0,0,0)) )
    Bother.append( arrow(pos=(x,y,z), axis=(0,0,0)) )
    Bother.append( arrow(pos=(x,y,z), axis=(0,0,0)) )
    Bother.append( arrow(pos=(x,y,z), axis=(0,0,0)) )
    Bother.append( arrow(pos=(x,y,z), axis=(0,0,0)) )
    Bother.append( arrow(pos=(x,y,z), axis=(0,0,0)) )

    a = vector(x,y,z)
    b = rotate(a,angle=phi, axis=(1,0,0))
    Bother.append( arrow(pos=(b.x,b.y,b.z), axis=(0,0,0)) )
    Bother.append( arrow(pos=(b.x,-b.y,b.z), axis=(0,0,0)) )
    b = rotate(a, angle=3*phi, axis=(1,0,0))
    Bother.append( arrow(pos=(b.x,b.y,b.z), axis=(0,0,0)) )
    Bother.append( arrow(pos=(b.x,-b.y,b.z), axis=(0,0,0)) )

for arr in Bother:
    arr.color = Bcolorother
    arr.shaftwidth = swother
    arr.fixedwidth = 1
    
scene.center = surface.pos/2
scene.forward = -vector(1,0,2)
scene.autoscale = 0

flux = 0
dt = 0.01
count = 0
limit = 50
changemu = 1/limit
# 0 to 7: rising, flat, falling, flat, falling, flat, rising, flat, repeat
changes = [1,0,-1,0,-1,0,1,0]
stage = 0
scene.mouse.getclick()
dBdtarr = arrow(pos=surface.pos+vector(0,-0.1*Rdisk,0.2*Rdisk), axis=(0,0,0), color=color.magenta,
                fixedwidth=1, shaftwidth=sw)
while 1:
    rate(limit/5)
    ammeter.text = "I = %0.1f A" % (abs(int(100*mu.x)/10))
    if showallB: showB()
    oldflux = flux
    flux = 0
    for arr in Bsurface:
        B = Bfield(coil.pos, arr.pos)
        arr.axis = B*Bscale
        flux += dot(B,xhat)*pi*abs(arr.y)*dR # 2 B's per y
    dflux = flux - oldflux
    dBdtarr.axis = (0.06*dflux/1e-8,0,0)
    E = (dflux/dt)/(2*pi*surface.radius)
    for a in Earr:
        a.axis = -E*Escale*cross(xhat,a.vv)
    count += 1
    if count > limit:
        count = 0
        stage += 1
        dBdtarr.visible = 0
        if stage >= len(changes):
            stage = 0
        for a in Earr:
            a.axis = vector(0,0,0)
        if changes[stage] == 0: # constant current
            if stage == 3 or stage == 7:
                for arr in Bsurface+Bother:
                    arr.visible = 0
            scene.mouse.getclick()
            for arr in Bsurface+Bother:
                arr.visible = 1
            
            stage += 1 # skip a stage with constant current in the coil
            if stage >= len(changes):
                stage = 0
        dBdtarr.visible = 1
    else:
        mu.x = mu.x+changes[stage]*changemu
        current.axis = vector(0,-mu.x/5,0)