
**********************************************

       function scalar(x1,y1,z1, x2,y2,z2)
       scalar = x1*x2+y1*y2+z1*z2
       return
       end

*      vecprod: a vec b = [a2b3 - a3b2, a3b1 - a1b3, a1b2 - a2b1].
       function vecprodx(a1,a2,a3, b1,b2,b3)
       vecprodx = a2*b3 - a3*b2
       return
       end
       function vecprody(a1,a2,a3, b1,b2,b3)
       vecprody = a3*b1 - a1*b3
       return
       end
       function vecprodz(a1,a2,a3, b1,b2,b3)
       vecprodz = a1*b2 - a2*b1
       return
       end

**********************************************

       subroutine computeB (x,y,z, Bx,By,Bz)
       common /poles/poleN,xpn,ypn,zpn,poleS,xps,yps,zps
       parameter (numgrid=2*10+1)
       Bx=0.
       By=0.
       Bz=0.
       
       pole=poleN/(numgrid**2)
       dx = xpn/numgrid
       dy = ypn/numgrid 
       do ix = -numgrid/2, numgrid/2
        do iy = -numgrid/2, numgrid/2       
         xp = ix*dx
         yp = iy*dy
         call computeBsinglepole(x,y,z,xp,yp,zpn,pole, 
     *                           partBx,partBy,partBz)
         Bx=Bx+partBx
         By=By+partBy
         Bz=Bz+partBz
        end do
       end do

       pole=poleS/(numgrid**2)
       dx = xps/numgrid
       dy = yps/numgrid 
       do ix = -numgrid/2, numgrid/2
        do iy = -numgrid/2, numgrid/2       
         xp = ix*dx
         yp = iy*dy
         call computeBsinglepole(x,y,z,xp,yp,zps,pole, 
     *                           partBx,partBy,partBz)
         Bx=Bx+partBx
         By=By+partBy
         Bz=Bz+partBz
        end do
       end do

       return
       end

       subroutine computeBsinglepole(x,y,z,xp,yp,zp,pole, Bx,By,Bz)
       squaredist = scalar(x-xp,y-yp,z-zp, x-xp,y-yp,z-zp)
       B =  pole/squaredist
       dist = sqrt(squaredist)
       Bx =  B*(x-xp)/dist
       By =  B*(y-yp)/dist                 
       Bz =  B*(z-zp)/dist
       return
       end

       subroutine computeLorentz (x,y,z,vx,vy,vz, flx,fly,flz)
       call computeB(x,y,z, Bx,By,Bz)
       flx = vecprodx(vx,vy,vz, Bx,By,Bz)
       fly = vecprody(vx,vy,vz, Bx,By,Bz)
       flz = vecprodz(vx,vy,vz, Bx,By,Bz)
       return
       end

       subroutine computev(x,y,z,omega, vx,vy,vz)
*      coil rotating around x axis
       vx = 0.
       vy = omega*z
       vz = -omega*y
       return
       end   

       function fem(x1,y1,z1, x2,y2,z2, ds, omega)
       dist = sqrt(scalar(x2-x1,y2-y1,z2-z1, x2-x1,y2-y1,z2-z1))
       nsteps = dist/ds
       dx = ds*(x2-x1)/dist
       dy = ds*(y2-y1)/dist
       dz = ds*(z2-z1)/dist
       fem = 0.
       do n=0,nsteps-1
        x = x1+(n+.5)*dx
        y = y1+(n+.5)*dy
        z = z1+(n+.5)*dz          
        call computev(x,y,z,omega, vx,vy,vz)
        call computeLorentz (x,y,z,vx,vy,vz, flx,fly,flz)
        fem = fem + scalar(flx,fly,flz, dx,dy,dz)
       end do
       return
       end

       function totalfem(t, FiloLunghezza,distanzaDaAsse,ds, omega)
* dopo rotazione omega*t a partire dal piano xz - x non varia
       theta = omega*t
* 1mo lato di connessione
       x1 = -FiloLunghezza/2.
       y1 = 0.
       z1 = 0.
       x2 = -FiloLunghezza/2.
       y2 = 0.
       z2 = -distanzaDaAsse
       y1=z1*sin(theta)
       z1=z1*cos(theta)
       y2=z2*sin(theta)
       z2=z2*cos(theta)
* l' ordine e' importante
       totalfem = 0.
       totalfem = totalfem + fem(x1,y1,z1, x2,y2,z2, ds, omega) 
* lato parallelo all'asse
       x1 = -FiloLunghezza/2.
       y1 = 0.
       z1 = -distanzaDaAsse
       x2 = + FiloLunghezza/2.
       y2 = 0.
       z2 = -distanzaDaAsse
       y1=z1*sin(theta)
       z1=z1*cos(theta)
       y2=z2*sin(theta)
       z2=z2*cos(theta)
* l' ordine e' importante
       totalfem = totalfem + fem(x1,y1,z1, x2,y2,z2, ds, omega)
* 2ndo lato di connessione
       x1 = FiloLunghezza/2.
       y1 = 0.
       z1 = -distanzaDaAsse
       x2 = FiloLunghezza/2.
       y2 = 0.
       z2 = 0.
       y1=z1*sin(theta)
       z1=z1*cos(theta)
       y2=z2*sin(theta)
       z2=z2*cos(theta)
* l' ordine e' importante
       totalfem = totalfem + fem(x1,y1,z1, x2,y2,z2, ds, omega) 
       return
       end

**********************************************

       program computeFEM
       common /poles/poleN,xpn,ypn,zpn,poleS,xps,yps,zps
       parameter (nsamples = 100, npositions=5)
       dimension out(0:nsamples, npositions)

       poleN = 4.
       xpn = 50.
       ypn = 25.
       zpn = -50.

       poleS = -4.
       xps = 50.
       yps = 25.
       zps = -25.

       ds=.1
       omega = 1.
       pi = 4.*atan(1.)

       filoLunghezza = 230.

       TT = 2*pi/omega
       Dt=TT/nsamples

       do npos=1,npositions 
        distanzaDaAsse = 20. - 10.*float(npos-1)
        print *, 'distanza da asse ', distanzaDaAsse
        do nsample=0,nsamples
         t = 0. + nsample*Dt
         out(nsample,npos) =  totalfem
     *           (t, FiloLunghezza,distanzaDaAsse, ds, omega)   
        enddo
       enddo


       do nsample=0,nsamples
        write (10,*) nsample, (out(nsample,npos),npos=1,npositions)
       end do

       stop 
       end



