        program errorfldscoil ! vec pot, bx, by and bs for single coil.
        implicit double precision (a-h,o-z)
c
        dimension vector(501,201),brfld(501,201),bsfld(501,201),
     +  b1(501,201,360), b2(501,201,360), b3(501,201,360)
        dimension rad(201),xlong(501),bx(501,201),by(501,201),
     +   bs(501,201)
      common/parameters/XL,R1   ! XL and R1 length of the field region
c   eta is position of center and rho radius
        external a1,a2,a3,a4,a5,a6,a7,a8
       open(unit=29,file='pot_sol_coil.dat',status='unknown',
     +  form='unformatted')
       open(unit=30,file='bx_sol_coil.dat',status='unknown',
     +  form='unformatted')
       open(unit=31,file='by_sol_coil.dat',status='unknown',
     +  form='unformatted')
       open(unit=32,file='bs_sol_coil.dat',status='unknown',
     +  form='unformatted')
       piover2=asin(1.0d0)
       perm_mu0=8*piover2*1.d-7  ! permeability
       write(6,*)'enter: current in [A]'
       read(5,*)curr
       write(6,*)'enter: radius of coil in [m]'
       read(5,*)R1
c        R1=0.7
        XL=2*R1  ! total length is 2*XL; XL is 2 times the radius
        n_radial_points=200
        n_long_points=500
        ds=2*XL/real(n_long_points)
        dr=0.8*R1/real(n_radial_points)   ! 80 % inside solenoid
        delta=0  ! radial position of coil ( i the paper is R)
        eta=0    ! long position of the coil (in the paper is S)
c potential as a power series expansion in r
c        b00=mu_o*I/2
        b00=0.5*perm_mu0*curr
c order of the expansion
c        type*,'up to which order in r? (1,3,5,7...(le.7))'
c        read*,iorder
        write(6,*)'ENTER:tilt angle(psi/rad) and dir.of rotation(degree)'
        read(5,*)psi,roteta
        roteta=roteta*piover2/90.d0
        theta=0.
c        write(6,*)'do you want to ENTER number of angles? Y(1);n(0)'
c        read(5,*)iloop
c        if(iloop.eq.1)then
c           write(6,*)'ENTER number of angles'
c           read(5,*)iangle
           write(6,*)'ENTER angle number for output (< 200)'
c           read(5,*)iout
c           else
c        iangle=360
c       iout=1
c        endif
        iangle=200
        dtheta=1.d0/real(iangle)  ! arbitrarily chosen
        do m=1,iangle+1
           write(6,*)'angle value=',m, theta
          do k=1,n_radial_points+1
           do i=1,n_long_points+1
              vector(i,k)=0
              brfld(i,k)=0
              bsfld(i,k)=0
              end do
              end do
c  loop over positions
            r=0.
              do k=1,n_radial_points+1
           rad(k)=r
           s=-XL
           do i=1,n_long_points+1
          xlong(i)=s
c Now the rotated quantities
           x=r*cos(theta)
           y=r*sin(theta)
           aux1=x*cos(roteta)+y*sin(roteta)
           aux2=x*sin(roteta)-y*cos(roteta)
           pris=s*cos(psi)+aux2*sin(psi)
c           prir=sqrt(x**2+(y*cos(psi))**2+(s*sin(psi))**2-
c     +     2d0*s*y*cos(psi)*sin(psi))
           prir=sqrt((r*cos(psi))**2-(aux1*(1d0-cos(psi)))**2+
     +     (s*sin(psi))**2+2.*cos(psi)*(1d0-cos(psi))*aux1**2-
     +     2.*cos(psi)*sin(psi)*s*aux2)
c
           vec1=0.5*a1(s,eta)
           vec3=1d0/16.*a3(s,eta)
           vec5=1d0/384.*a5(s,eta)
           vec7=1d0/18432.*a7(s,eta)
c
           br1=0.5*a2(s,eta)
           br2=1d0/16.*a4(s,eta)
           br3=1d0/384.*a6(s,eta)
           br4=1d0/18432.*a8(s,eta)
c
           pribr1=0.5*a2(pris,eta)
           pribr2=1d0/16.*a4(pris,eta)
           pribr3=1d0/384.*a6(pris,eta)
           pribr4=1d0/18432.*a8(pris,eta)
c
           bs1=a1(s,eta)
           bs2=1d0/8.*a3(s,eta)
           bs3=1d0/192.*a5(s,eta)
           bs4=1d0/9216.*a7(s,eta)
c
           pribs1=a1(pris,eta)
           pribs2=1d0/8.*a3(pris,eta)
           pribs3=1d0/192.*a5(pris,eta)
           pribs4=1d0/9216.*a7(pris,eta)
c
c           vector(i,k)=b00*r*(vec1+r**2*(vec3+r**2*(vec5+r**2*vec7)))
c           brfld(i,k)=-b00*r*(br1+r**2*(br2+r**2*(br3+r**2*br4)))
c           bsfld(i,k)=b00*(bs1+r**2*(bs2+r**2*(bs3+r**2*bs4)))
c
           vector(i,k)=r*(vec1+r**2*(vec3+r**2*(vec5+r**2*vec7)))
           brfld(i,k)=-r*(br1+r**2*(br2+r**2*(br3+r**2*br4)))
           bsfld(i,k)=(bs1+r**2*(bs2+r**2*(bs3+r**2*bs4)))
c
c           privector(i,k)=prir*(privec1+prir**2*(privec3+prir**2*
c     +     (privec5+prir**2*privec7)))

           pribrfld=-prir*(pribr1+prir**2*(pribr2+prir**2*
     +     (pribr3+prir**2*pribr4)))
           pribsfld=(pribs1+prir**2*(pribs2+prir**2*
     +     (pribs3+prir**2*pribs4)))

c rotated fields - the unrotated fields === ERROR FIELDS
           b1(i,k,m)=(pribrfld-brfld(i,k))*cos(roteta)+
     +     pribrfld*sin(psi)*sin(roteta)-pribsfld*sin(psi)
           b2(i,k,m)=(pribrfld*cos(psi)-brfld(i,k))*cos(roteta)-
     +     pribrfld*sin(roteta)+sin(psi)*pribsfld
           b3(i,k,m)=-sin(psi)*pribrfld+ pribsfld*cos(psi)-bsfld(i,k)
c           b1(i,k,m)=pribrfld*sin(theta)
c           b2(i,k,m)=pribrfld*cos(theta)+psi*pribsfld
c           b3(i,k,m)=-psi*pribrfld*cos(theta)+ pribsfld
c
           if(i.eq.300.and.m.eq.iout)then
           write(20,*)rad(k),b1(i,k,m)
           write(21,*)rad(k),b2(i,k,m)
           write(22,*)rad(k),b3(i,k,m)
           endif
          s=s+ds
        end do  ! end s do loop
        if(k.eq.1.and.m.eq.iout)then
           do i=1,n_long_points+1
           write(10,*)xlong(i),b1(i,k,m)
           write(11,*)xlong(i),b2(i,k,m)
           write(12,*)xlong(i),b3(i,k,m)
           end do
           endif
        r=r+dr
        end do  ! end r do loop
        theta=theta+dtheta*4.*piover2
        end do   ! end theta do loop
        write(29)n_radial_points+1,n_long_points+1
       write(29)vector
        write(30)n_radial_points+1,n_long_points+1
        do m=1,iangle
       if(m.eq.iout)then
          do k=1,n_radial_points+1
          do i=1,n_long_points+1 
          bx(i,k)=b1(i,k,m)
       end do
       end do
       write(30)bx
       endif
       end do
        write(31)n_radial_points+1,n_long_points+1
        do m=1,iangle
       if(m.eq.iout)then 
          do k=1,n_radial_points+1 
          do i=1,n_long_points+1
          by(i,k)=b2(i,k,m)
       end do
       end do
       write(31)by
       endif
       end do
         write(32)n_radial_points+1,n_long_points+1
        do m=1,iangle
       if(m.eq.iout)then
          do k=1,n_radial_points+1 
             do i=1,n_long_points+1 
         bs(i,k)=b3(i,k,m)
       end do
       end do
       write(32)bs
       endif
       end do
       write(6,*)'pot_sol_coil.dat contains:POT; NORMALIZE',
     +      n_radial_points+1,n_long_points+1
       write(6,*)'bx_sol_coil.dat contains:Bx; NORMALIZE',
     +      n_radial_points+1,n_long_points+1
       write(6,*)'by_sol_coil.dat contains:By; NORMALIZE',
     +      n_radial_points+1,n_long_points+1
       write(6,*)'bs_sol_coil.dat contains:Bs; NORMALIZE',
     +      n_radial_points+1,n_long_points+1
       write(6,*)'fort.10 contains: bx(1,s) vs s'
       write(6,*)'fort.11 contains: by(1,s) vs s'
       write(6,*)'fort.12 contains: bs(1,s) vs s'
       write(6,*)'fort.20 contains: bx(r,250) vs r'
       write(6,*)'fort.21 contains: by(r,250) vs r'
       write(6,*)'fort.22 contains: bs(r,250) vs r'
c 50    format(1x,6(g12.5,1x))
        stop
        end
               
c   FUNCTIONS
      real*8 function a1(ss,eta) ! 1st order
      implicit double precision (a-h,o-z)
      common/parameters/XL,R1
c   eta is position of center and R1 radius
      aux0=sqrt(R1**2+(ss-eta)**2)
      aux1=1/aux0
      aux2=-(ss-eta)**2/aux0**3
      a1=aux1+aux2
      end
c
      real*8 function a2(ss,eta) ! 2nd order
      implicit double precision (a-h,o-z)
      common/parameters/XL,R1
c   eta is position of center and R1 radius
      aux0=sqrt(R1**2+(ss-eta)**2)
      aux1=-3d0*(ss-eta)/aux0**3
      aux2=3d0*(ss-eta)**3/aux0**5
      a2=aux1+aux2
      end
c
      real*8 function a3(ss,eta) ! 3rd order
      implicit double precision (a-h,o-z)
c   eta is position of center and R1 radius
      common/parameters/XL,R1
      aux0=sqrt(R1**2+(ss-eta)**2)
      aux1=-3d0/aux0**3
      aux2=18d0*(ss-eta)**2/aux0**5
      aux3=-15d0*(ss-eta)**4/aux0**7
      a3=aux1+aux2+aux3
      end
c
      real*8 function a4(ss,eta) ! 4th order
      implicit double precision (a-h,o-z)
      common/parameters/XL,R1
c   eta is position of center and R1 radius
      aux0=sqrt(R1**2+(ss-eta)**2)
      aux1=45d0*(ss-eta)/aux0**5
      aux2=-150d0*(ss-eta)**3/aux0**7
      aux3=105d0*(ss-eta)**5/aux0**9
      a4=aux1+aux2+aux3
      end
c
      real*8 function a5(ss,eta) ! 5th order
      implicit double precision (a-h,o-z)
c   eta is position of center and R1 radius
      common/parameters/XL,R1
      aux0=sqrt(R1**2+(ss-eta)**2)
      aux1=-945d0*(ss-eta)**6/aux0**11
      aux2=1575d0*(ss-eta)**4/aux0**9
      aux3=-675d0*(ss-eta)**2/aux0**7
      aux4=45d0/aux0**5
      a5=aux1+aux2+aux3+aux4
      end
c
      real*8 function a6(ss,eta) ! 6th order
      implicit double precision (a-h,o-z)
      common/parameters/XL,R1
c   eta is position of center and R1 radius
      aux0=sqrt(R1**2+(ss-eta)**2)
      aux1=-1575d0*(ss-eta)/aux0**7
      aux2=11025d0*(ss-eta)**3/aux0**9
      aux3=-19845d0*(ss-eta)**5/aux0**11
      aux4=10395d0*(ss-eta)**7/aux0**13
      a6=aux1+aux2+aux3+aux4
      end
c
      real*8 function a7(ss,eta) ! 7th order
      implicit double precision (a-h,o-z)
c  eta is position of center and R1 radius
      common/parameters/XL,R1
      aux0=sqrt(R1**2+(ss-eta)**2)
      aux1=-135135*(ss-eta)**8/aux0**15
      aux2=291060*(ss-eta)**6/aux0**13
      aux3=-198450*(ss-eta)**4/aux0**11
      aux4=44100*(ss-eta)**2/aux0**9
      aux5=-1575/aux0**7
      a7=aux1+aux2+aux3+aux4+aux5
      end
c
      real*8 function a8(ss,eta) ! 8th order
      implicit double precision (a-h,o-z)
      common/parameters/XL,R1
c   eta is position of center and R1 radius
      aux0=sqrt(R1**2+(ss-eta)**2)
      aux1=2027025d0*(ss-eta)**9/aux0**17
      aux2=-4864860d0*(ss-eta)**7/aux0**15
      aux3=3929310d0*(ss-eta)**5/aux0**13
      aux4=-1190700d0*(ss-eta)**3/aux0**11
      aux5=99225d0*(ss-eta)/aux0**9
      a8=aux1+aux2+aux3+aux4+aux5
      end
c 
