#To run this program install Python and Vpython
#see: http://vpython.org/
#for Windows download and install Python and Vpython from: 
# http://vpython.org/contents/download_windows.html
#Copy this program and paste it in IDLE (for example), save it as
# Travellinwavelinac.py and run it.
#
#Simple travelling wave Linear accelerator model
#Author:  Manuel Jose Paez
#Instituto de Fisica
#Universidad de Antioquia
#Medellin, Colombia
#mpaez@fisica.udea.edu.co
'''
The tube is sectioned to  see through.
The electric field is represented by blue arrows
The accelerated particle is an electron
Is accelerated in the cavities, the electric field is opposite to its direction
Notice that the particle finds the same phase of the electric field (wave)
in each cavity. This is the principle of the accelerator
The white curves represent the electric field lines. Magnetic field lines not
drawn, they are circles tanget to the longitudinal direction
The cavities are fed with klystrons that produce electromagnetic oscillations
in each cavity. There is a phase difference between a cavity and the following
is pi/6 in this model.
The whole wavelength contains 12 cavities.
'''
from visual.graph import * #graphics and math classes
scene=display(x=0,y=0,width=800,height=300,range=800,
           title='Linear Accelerator. E: electric field white arrows')
ball=sphere(pos=(-660,0,0),radius=10,color=color.green)  #the electron
bola=sphere(radius=2,color=color.white)
bola.trail=curve(color=color.white,radius=2)
bola2=sphere(radius=2,color=color.white)
bola2.trail=curve(color=color.white,radius=2)
bola3=sphere(radius=2,color=color.white)
bola3.trail=curve(color=color.white,radius=2) 
arrow1=arrow(color=color.blue)          #arrows for electric field
arrow2=arrow(color=color.blue)          #In this way, arrows will be erased
arrow3=arrow(color=color.blue)
arrow4=arrow(color=color.blue)
arrow5=arrow(color=color.blue)
arrow6=arrow(color=color.blue)           #arrows for electric field
arrow7=arrow(color=color.blue)           #In this way, arrows will be erased
arrow8=arrow(color=color.blue)           #and repainted as the phase changes
arrow9=arrow(color=color.blue)
arrow10=arrow(color=color.blue)
arrow11=arrow(color=color.blue)
arrow12=arrow(color=color.blue)
def iris():                              #circular apertures to increase
    xpos=-720                            #the impedance of the cavities
    length=2
    for i in range(0,13):
        iris1=helix(pos=(xpos,0,0),radius=50,axis=(length,0,0),coils=1.5,
                thickness=12,color=color.orange)
        iris2=helix(pos=(xpos,0,0),radius=45,axis=(length,0,0),coils=1.5,
                thickness=12,color=color.orange)
        iris3=helix(pos=(xpos,0,0),radius=55,axis=(length,0,0),coils=1.5,
                thickness=12,color=color.orange)
        xpos +=120
def tube(xpos,length):              #the tube is plotted with lines
    thincr=math.pi/40.0             #forming a cylinder. No lines in front
    th=0.0                          #to see the motion of the electron
    for i in range(1,40):
        th=th+thincr
        y=55*math.sin(th-math.pi/2)
        z=55*math.cos(th-math.pi/2)
        line=curve(pos=[(xpos,y,-z),(xpos+length,y,-z)],radius=5,color=color.orange)
def Efield(xp):                     #plots the field lines as white curves
    th=0.0                          #only few are drawn.
    thincr=math.pi/10.
   
    x=xp
    xincr=12.0
    zincr=thincr
    for i in range(0,11):
        y=25.*math.cos(math.pi/4.)*(1.0-cos(th-math.pi/2.0))
        z =25.0*sin(th)-55.*cos(math.pi/4.)
        bola.pos=vector(x,y+10,z)
        bola.trail.append(pos=(x,y+10,z))
        bola2.pos=vector(x,-y-10,z)
        bola2.trail.append(pos=(x,-y-10,z))
        bola3.pos=vector(x,0,z)
        bola3.trail.append(pos=(x,0,z))
        
        th +=thincr
        x +=xincr
       

def Ecavity(phase):                #shows the arrows that represnet the E field
    ax=zeros((13),Float)           #gives at each time the size of the arrows
    posarrow=zeros((13),Float)
    xx=-720
    Eth=math.pi/12.
    for i in range(1,13):
        ax[i]=-90.*cos(Eth-phase)
        if abs(ax[i])<50:
            ax[i]=50.0
        if ax[i]<=0:    #to center in cavity the arrow of the Electric field
            posarrow[i]=xx+120 - (120-abs(ax[i]))/2.0
        else:
            posarrow[i]=xx+(120-abs(ax[i]))/2.0
        arrow1.pos=vector(posarrow[1],0,0)  #arrows  represent the electric field
        arrow1.axis=(ax[1],0,0)             #at the center of the cavity
        arrow2.pos=vector(posarrow[2],0,0)
        arrow2.axis=(ax[2],0,0)
        arrow3.pos=vector(posarrow[3],0,0)
        arrow3.axis=(ax[3],0,0)
        arrow4.pos=vector(posarrow[4],0,0)
        arrow4.axis=(ax[4],0,0)
        arrow5.pos=vector(posarrow[5],0,0)
        arrow5.axis=(ax[5],0,0)
        arrow6.pos=vector(posarrow[6],0,0)
        arrow6.axis=(ax[6],0,0)
        arrow7.pos=vector(posarrow[7],0,0)
        arrow7.axis=(ax[7],0,0)
        arrow8.pos=vector(posarrow[8],0,0)
        arrow8.axis=(ax[8],0,0)
        arrow9.pos=vector(posarrow[9],0,0)
        arrow9.axis=(ax[9],0,0)
        arrow10.pos=vector(posarrow[10],0,0)
        arrow10.axis=(ax[10],0,0)
        arrow11.pos=vector(posarrow[11],0,0)
        arrow11.axis=(ax[11],0,0)
        arrow12.pos=vector(posarrow[12],0,0)
        arrow12.axis=(ax[12],0,0)
        #print xx
        if phase==0:
             Efield(xx)
        Eth +=math.pi/6.0       #phase difference between consecutive cavities
        #print xx
        xx +=120               #longitudinal size of each cavity
        
iris()
tube(-720,1440)
xb=-660
imax=12
#xf=-720
for i in range(1,imax):
    #rate(2)
    phase=(i-1)*math.pi/6.0
    #Efield(xf)
    Ecavity(phase)
   
    dt=1.0/(12.0-i)
    time=0.0
    #xf +=120
    if i==0:
        xb=-660
    if i==1:  
        for j in range(0,12):
          rate(10)
          xbincr=120.0*(time)**2
          xb =-660+120+xbincr
          ball.pos=vector(xb,0,0)
          
          time +=dt
    else:
        time=dt
        for j in range(1,imax+1-i):
            rate(10)
            xbincr=120.0*(time)**2
            xb=-660+i*120.0+xbincr
            ball.pos=vector(xb,0,0)               
            time +=dt
        if j==imax:                
            ball.pos=vector(800,0,0)    #out of the accelerator