program field_sc ! space charge fields for gaussian distrib.
        implicit double precision (a-h,o-z)
c
        common/position/r,s
        common/beam_para/sx,s3
        external f1,f2,f3,f4
        dimension coord_r(500),coord_s(500),
     +  field_s(500,500),field_r(500,500)
c
        character*5 fileunits(1),fileunitr(1)
        data xmass/105.65/,clight/2.9979d8/,
     +  epsi_0/8.8542d-12/,q/1.6022d-19/
        data fileunitr(1)/'rad01'/        
        data fileunits(1)/'lon01'/
c   OPEN FILE 
        open(unit=28,file='longitudinal.dat',status='unknown',
     +  form='formatted')
        open(unit=29,file='radial.dat',status='unknown',
     +  form='formatted')

        pi=2*asin(1.0d0) ! number pi
c       xmass     ! mass of muons MeV
c       clight    ! speed of light m/s
        write(6,*)'enter muon momentum MeV/c'
        read*,xmom
        energ=sqrt(xmom*xmom+xmass*xmass)
        gamma= energ/xmass
        gbeta=sqrt((gamma-1.)*(gamma+1.))  ! gamma * beta
        beta=gbeta/gamma
        cbeta=clight*beta 
c
c        write(6,*)'enter: Number of particles in the bunch'
c        read*,np
        xnp=4.d12
        write(6,*)'enter: rms radius of the beam (m)'
        read*,sx
        write(6,*)'enter: rms length of the bunch (m)'
        read*,s3
        write(6,*)'enter: No. of longitudinal points in cavity (<500)'
        read*,n_long_points
        write(6,*)'enter:number of radial points in cavity (<500)'
        read*,n_radial_points
        write(28,*)n_radial_points
        write(29,*)n_long_points
        xl_c=4*s3
        r_c=4*sx
        ds=xl_c/real(n_long_points)
        dr=r_c/real(n_radial_points)
        dl=sx/s3 !ratio
        du=1.
        qe=q/epsi_0  ! C/....
        coeff1= qe/sqrt(pi)*4*xnp  
        coeff1=coeff1/(2*(s3-sx)*(s3+sx))**(3/2)
           coeff2=coeff1/(2*(s3-sx)*(s3+sx))
        icount=1
        open(unit=16,file=fileunitr(icount),status='unknown',
     +  form='unformatted')
        open(unit=17,file=fileunits(icount),status='unknown',
     +  form='unformatted')        
        dsxds=0.001
        ds3ds=0.001
        s=-real(xl_c/2)
        do 35 is=1,n_long_points+1 ! loop long. variable
        r=0
c        dsxds=0.001
c        ds3ds=0.0325
        do 40 ir=1,n_radial_points ! loop over radial variable

           call qromb(f1,dl,du,result1)
           field_r(ir,is)=result1*r*coeff1*2
           call qromb(f2,dl,du,result2)
           call qromb(f3,dl,du,result3)
           call qromb(f4,dl,du,result4)
         e1=-(result2*dsxds-result3*ds3ds)*coeff2
         e2=1/gamma**2*result4*coeff1*s*2
         field_s(ir,is)=e2
           coord_r(ir)=r
        r=r+dr
40      end do
        coord_s(is)=s
        s=s+ds
35      end do

        write(16)field_r
        write(17)field_s

        do 70 i=1,n_radial_points
       write(28,*)coord_r(i),
     +  (coord_s(j),field_s(i,j),j=1,n_long_points+1)
70      end do
        do 71 i=1,n_long_points+1
       write(29,*)coord_s(i),
     +  (coord_r(j),field_r(j,i),j=1,n_radial_points)
 71    end do
       write(6,*)'longitudinal.dat contains: r(i),(s(j),es(i,j) )'
       write(6,*)'radial.dat contains:s(i),(r(j),er(j,i))'
        stop
        end               
c   FUNCTIONS
      real*8 function f1(u) ! radial field
      implicit double precision (a-h,o-z)
      common/position/r,s
      common/beam_para/sx,s3
      aux1=sqrt((1-u)*(1+u))/u/u/u
      aux2=1.
      aux3=0.5*(1-u)*(1+u)*((r/u)**2+s**2)/(s3-sx)/(s3+sx)
      f1=aux1*exp(-aux3)
      end
c
      real*8 function f2(u) ! long field sigma_x component
      implicit double precision (a-h,o-z)
      common/position/r,s
      common/beam_para/sx,s3
      aux1=sqrt((1-u)*(1+u))/u/u/u
      aux2=(r/u)**2-r**2-2*(s3-sx)*(s3+sx)
      aux3=0.5*(1-u)*(1+u)*((r/u)**2+s**2)/(s3-sx)/(s3+sx)
      f2=aux1*aux2*exp(-aux3)
      end
c
      real*8 function f3(u) ! long field sigma_3 component
      implicit double precision (a-h,o-z)
      common/position/r,s
      common/beam_para/sx,s3
      aux1=sqrt((1-u)*(1+u))/u
      aux2=2*(s**2*(1-u)*(1+u)-(s3-sx)*(s3+sx))
      aux3=0.5*(1-u)*(1+u)*((r/u)**2+s**2)/(s3-sx)/(s3+sx)
      f3=aux1*aux2*exp(-aux3)
      end
c
      real*8 function f4(u) ! longitudinal field
      implicit double precision (a-h,o-z)
      common/position/r,s
      common/beam_para/sx,s3
      aux1=sqrt((1-u)*(1+u))/u
      aux2=1.
      aux3=0.5*(1-u)*(1+u)*((r/u)**2+s**2)/(s3-sx)/(s3+sx)
      f4=aux1*exp(-aux3)
      end