program random_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=mfile,file=fname,status='unknown',iostat=ioc)
          character*10 fname
c        open(unit=29,file='data_solenoid.dat',status='unknown',
c     +  iostat=ioc)
c  set seed random number generator
          irandom=-777
        XL=0.2  ! total length is 2*L=0.4 m
        R1=0.7
        R2=R1+0.15
        n_long_points=500
c
        mx=1
        my=1
        mz=n_long_points
c
        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))
c
        mfile=29
        fname='data_solenoid.dat'
        write(mfile,*)
        write(mfile,*)mx,my,mz
c  find out about these vectors
c	    write(mfile,*) ordtrkbkg
c	    write(mfile,*) (htrkbkg(i),i=1,mz)
c	    write(mfile,*) (hptrkbkg(i),i=1,mz)
          nrecords = mx*my*mz
!
          close(mfile)
!
        b00=log((R2+sqrt(R2**2+XL**2))/(R1+sqrt(R1**2+XL**2)))/(2*XL)

        do i=1,n_long_points
           bx=0
           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))      
c   call to gaussian random number generator
           rx1=gasdev(irandom)
           by=by*b00*rx1
           bs=bs*b00*rx1
           rx2=gasdev(irandom)
           bty=bty*b00*rx1
           bts=bts*b00*rx1
           rx3=gasdev(irandom)
           bp=bp*b00*rx3
           write(mfile,*)mx,my,1,bx,by+bty,bs+bts+bp
        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

      FUNCTION gasdev(idum)
      INTEGER idum
      DOUBLE PRECISION gasdev
CU    USES ran1
      INTEGER iset
      DOUBLE PRECISION fac,gset,rsq,v1,v2,ran1
      SAVE iset,gset
      DATA iset/0/
      if (iset.eq.0) then
1       v1=2.d0*ran1(idum)-1.d0
        v2=2.d0*ran1(idum)-1.d0
        rsq=v1**2+v2**2
        if(rsq.ge.1.d0.or.rsq.eq.0.d0)goto 1
        fac=sqrt(-2.d0*log(rsq)/rsq)
        gset=v1*fac
        gasdev=v2*fac
        iset=1
      else
        gasdev=gset
        iset=0
      endif
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software *#'".,.
      FUNCTION ran1(idum)
      INTEGER idum,IA,IM,IQ,IR,NTAB,NDIV
      DOUBLE PRECISION ran1,AM,EPS,RNMX
      PARAMETER (IA=16807,IM=2147483647,AM=1.d0/IM,IQ=127773,IR=2836,
     *NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2d-7,RNMX=1.d0-EPS)
      INTEGER j,k,iv(NTAB),iy
      SAVE iv,iy
      DATA iv /NTAB*0/, iy /0/
      if (idum.le.0.or.iy.eq.0) then
        idum=max(-idum,1)
        do 11 j=NTAB+8,1,-1
          k=idum/IQ
          idum=IA*(idum-k*IQ)-IR*k
          if (idum.lt.0) idum=idum+IM
          if (j.le.NTAB) iv(j)=idum
11      continue
        iy=iv(1)
      endif
      k=idum/IQ
      idum=IA*(idum-k*IQ)-IR*k
      if (idum.lt.0) idum=idum+IM
      j=1+iy/NDIV
      iy=iv(j)
      iv(j)=idum
      ran1=min(AM*iy,RNMX)
      return
      END
C  (C) Copr. 1986-92 Numerical Recipes Software *#'".,.