program error_sheet ! errors in cylindrical current sheets.
        implicit double precision (a-h,o-z)
c
      common/parameters/XL,R1,R2
c   u is tild angle; eta is position of center and rho radius
        external f1,f2,f3,f4,f5
       open(unit=29,file='data_solenoid.dat',status='unknown',
     +  form='formatted')
        XL=0.2  ! total length is 2*L=0.4 m
        R1=0.7
        R2=R1+0.15
        n_long_points=500
        ds=60*XL/n_long_points
        s=-30*XL
        psi=0.01 ! tilde angle (rad)
        d=0.1  ! translation along y-axis (cm)
        delta=0.01 ! rms power fluct.( 1%)
c field at center of solenoid
c        b00=mu_o*I/(R2-R1)*Ln(R2+sqrt(R2**2+XL**2))-Ln(R1+sqrt(R1**2+XL**2))
        b00=log((R2+sqrt(R2**2+XL**2))/(R1+sqrt(R1**2+XL**2)))/(2*XL)
        do i=1,n_long_points
           by=f1(psi,s,XL,R2)-f1(psi,s,XL,R1)-(f1(psi,s,-XL,R2)-
     +  f1(psi,s,-XL,R1))
           bs=f2(psi,s,XL,R2)-f2(psi,s,XL,R1)-(f2(psi,s,-XL,R2)-
     +  f2(psi,s,-XL,R1))  
           bty=f3(d,s,XL,R2)-f3(d,s,XL,R1)-(f3(d,s,-XL,R2)-
     +  f3(d,s,-XL,R1))   
           bts=f4(d,s,XL,R2)-f4(d,s,XL,R1)-(f4(d,s,-XL,R2)-
     +  f4(d,s,-XL,R1))   
           bp=f5(delta,s,XL,R2)-f5(delta,s,XL,R1)-(f5(delta,s,-XL,R2)-
     +  f5(delta,s,-XL,R1))      
       write(29,50)s, by/b00,bs/b00,bty/b00,bts/b00,bp/b00
        s=s+ds
        end do
       write(6,*)'data_solenoid.dat contains:s, by bs,bty,bts,bp'
 50    format(1x,6(g12.5,1x))
        stop
        end               
c   FUNCTIONS
      real*8 function f1(u,ss,eta,rho) ! tilded solenoid
      implicit double precision (a-h,o-z)
c   u is tild angle; eta is position of center and rho radius
      aux0=sqrt((ss-eta)**2+rho**2)
      aux1=u*((2*eta-ss)*log(rho+sqrt((ss-eta)**2+rho**2))-
     +     ss*rho/aux0)
      aux2=4
      f1=aux1/aux2
      end
c

      real*8 function f2(u,ss,eta,rho) ! tilded solenoid
      implicit double precision (a-h,o-z)
c   u is tild angle; eta is position of center and rho radius
      aux0=sqrt((ss-eta)**2+rho**2)
      aux1=-u*u*(2*ss*log(rho+aux0)+rho*(ss-eta)*(eta**2+(ss-eta)**2)
     +     /aux0**3+rho/aux0*(3*(ss-eta)-eta*(2*ss-eta)/(ss-eta)))
      aux2=8
      f2=aux1/aux2
      end
c
      real*8 function f3(u,ss,eta,rho) ! translated solenoid
      implicit double precision (a-h,o-z)
c   u is translation distance; eta is position of center and rho radius
      aux0=sqrt((ss-eta)**2+rho**2)
      aux1=u*(log(rho+aux0)-rho/aux0)
      aux2=4
      f3=aux1/aux2
      end
c
      real*8 function f4(u,ss,eta,rho) ! translated solenoid
      implicit double precision (a-h,o-z)
c   u is translation distance; eta is position of center and rho radius
      aux0=sqrt((ss-eta)**2+rho**2)
      aux1=u*u*rho*(ss-eta)/aux0*(1/(ss-eta)**2-1/aux0)
      aux2=8
      f4=aux1/aux2
      end
c
      real*8 function f5(u,ss,eta,rho) ! power supply fluctuation
      implicit double precision (a-h,o-z)
c   u is rms power fluct.; eta is position of center and rho radius
      aux0=sqrt((ss-eta)**2+rho**2)
      aux1=u*(ss-eta)*log(rho+aux0)
      aux2=2
      f5=aux1/aux2
      end