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

       real function Bz(x,z,polom,znord,zsud)
       Bn=polom/((x)**2+(z-znord)**2)
       Bzn=Bn*(z-znord)/sqrt((x)**2+(z-znord)**2)
       Bs= - polom/((x)**2+(z-zsud)**2)
       Bzs=Bs*(z-zsud)/sqrt((x)**2+(z-zsud)**2)
       Bz=Bzn+Bzs
       return
       end

       real function Bx(x,z,polom,znord,zsud)
       Bn=polom/((x)**2+(z-znord)**2)
       Bxn=Bn*x/sqrt((x)**2+(z-znord)**2)
       Bs= - polom/((x)**2+(z-zsud)**2)
       Bxs=Bs*x/sqrt((x)**2+(z-zsud)**2)
       Bx=Bxn+Bxs
       return
       end 

       real function flx(x,z,polom,znord,zsud,omega)
       vy= - omega*z
       flx =  - vy*Bz(x,z,polom,znord,zsud)
*       print *, "flx = ", flx
       return
       end

       real function flz(x,z,polom,znord,zsud,omega)
       vy= - omega*z
       flz =  vy*Bx(x,z,polom,znord,zsud)
*       print *, "flz = ", flz
       return
       end

       real function flxint(xmin,xmax,dx,z,omega,polom,znord,zsud)
       f = 0.
       x=xmin
       do while (x.le.xmax) 
        f=f+dx*flx((x+dx/2.),z,polom,znord,zsud,omega)
        x=x+dx
       end do
       flxint=f
       return
       end

       real function flzint(zmin,zmax,dz,x,omega,polom,znord,zsud)
       f = 0.
       z=zmin
       do while (z.le.zmax)
        f=f+dz*flz(x,(z+dz/2.),polom,znord,zsud,omega)
        z=z+dz
       end do
       flzint=f
       return
       end

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

       real function Ez(x,z,polom,znord,zsud,omega)
*       print*, 'Ez input: ',x,z,polom,znord,zsud,omega
       En = - omega*znord *polom/((x)**2+(z-znord)**2)
       Ezn=En* x  /sqrt((x)**2+(z-znord)**2)
       Es= - (-omega*zsud) *polom/((x)**2+(z-zsud)**2)
       Ezs=Es* x /sqrt((x)**2+(z-zsud)**2)
       Ez=Ezn+Ezs
*       print *, "Ez = ", Ez
       return
       end

       real function Ex(x,z,polom,znord,zsud,omega)
*       print*, 'Ex input: ',x,z,polom,znord,zsud,omega
       En= omega*znord *polom/((x)**2+(z-znord)**2)
       Exn=   En*(z-znord)/sqrt((x)**2+(z-znord)**2)
       Es= - omega*zsud *polom/((x)**2+(z-zsud)**2)
       Exs=   Es*(z-zsud)/sqrt((x)**2+(z-zsud)**2)
       Ex=Exn+Exs
*       print *, "Ex = ", Ex
       return
       end 

       real function Exint(xmin,xmax,dx,z,omega,polom,znord,zsud)
       f = 0.
       x=xmin
       do while (x.le.xmax)
        f=f+dx*Ex((x+dx/2.),z,polom,znord,zsud,omega)
        x=x+dx
       end do
       Exint=f
       return
       end

       real function Ezint(zmin,zmax,dz,x,omega,polom,znord,zsud)
       f = 0.
       z=zmin
       do while (z.le.zmax)
        f=f+dz*Ez(x,(z+dz/2.),polom,znord,zsud,omega)
        z=z+dz
       end do
       Ezint=f
       return
       end

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

       subroutine compute(dx,dz,omega,zs,xs,znord,zsud,polom)

       print *, "******************************************"
       print *, "Semialtezza spira (asse x)   ", xs
       print *, "Semilarghezza spira (asse z) ", zs
       print *, "z del polo puntiforme Nord   ", znord
       print *, "z del polo puntiforme Sud    ", zsud
       print *, "velocita'' di rotazione omega", omega
       print *, "carica magnetica poli N e S  ", polom
       print *, "dx per integrazione          ", dx
       print *, "dz per integrazione          ", dz


*     forze di Lorentz:
       print *, "   caso magnete fermo e spira in rotazione"

       fem1 = flxint(-xs,xs,dx,-zs,omega,polom,znord,zsud)
       fem2 = flxint(-xs,xs,dx,zs,omega,polom,znord,zsud)

       fem3 = flzint(-zs,zs,dz,xs,omega,polom,znord,zsud)
       fem4 = flzint(-zs,zs,dz,-xs,omega,polom,znord,zsud)

       print '(2x,5(f10.6,A))', fem1," - ", fem2, 
     * " + ", fem3," - ", fem4, " = ", 
     *          fem1 - fem2 + fem3 - fem4

*     campi elettrici indotti (con principio di dualita')
       print *, "   caso spira ferma e magnete in rotazione"

       fem1 = Exint(-xs,xs,dx,-zs,omega,polom,znord,zsud)
       fem2 = Exint(-xs,xs,dx,zs,omega,polom,znord,zsud)

       fem3 = Ezint(-zs,zs,dz,xs,omega,polom,znord,zsud)
       fem4 = Ezint(-zs,zs,dz,-xs,omega,polom,znord,zsud)


       print '(2x,5(f10.6,A))', fem1," - ", fem2, 
     * " + ", fem3," - ", fem4, " = ", 
     *          fem1 - fem2 + fem3 - fem4

       return
       end

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

       program fems
       dz=.0001
       dx=.0001
       omega = -1.
       polom=10.

       zs=1.
       xs=1.
       znord=-2.
       zsud=-3.

       call  compute(dx,dz,omega,zs,xs,znord,zsud,polom)

       znord=-1.1
       zsud=-2.1
       call  compute(dx,dz,omega,zs,xs,znord,zsud,polom)

       znord=-2.
       zsud=-3.
       zs=1.9
       xs=1.
       call  compute(dx,dz,omega,zs,xs,znord,zsud,polom)


       znord=-2.
       zsud=-3.
       zs= .1
       xs= 1.
       call  compute(dx,dz,omega,zs,xs,znord,zsud,polom)

       znord=-2.
       zsud=-3.
       zs= 1.
       xs= .1
       call  compute(dx,dz,omega,zs,xs,znord,zsud,polom)

       znord=-2.
       zsud=-30000.
       zs= 1.
       xs= 1.
       call  compute(dx,dz,omega,zs,xs,znord,zsud,polom)  

       stop
       end

